EMAN2
averager.cpp
Go to the documentation of this file.
00001 
00005 /*
00006  * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
00007  * Copyright (c) 2000-2006 Baylor College of Medicine
00008  *
00009  * This software is issued under a joint BSD/GNU license. You may use the
00010  * source code in this file under either license. However, note that the
00011  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00012  * so you are responsible for compliance with the licenses of these packages
00013  * if you opt to use BSD licensing. The warranty disclaimer below holds
00014  * in either instance.
00015  *
00016  * This complete copyright notice must be included in any revised version of the
00017  * source code. Additional authorship citations may be added, but existing
00018  * author citations must be preserved.
00019  *
00020  * This program is free software; you can redistribute it and/or modify
00021  * it under the terms of the GNU General Public License as published by
00022  * the Free Software Foundation; either version 2 of the License, or
00023  * (at your option) any later version.
00024  *
00025  * This program is distributed in the hope that it will be useful,
00026  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00027  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00028  * GNU General Public License for more details.
00029  *
00030  * You should have received a copy of the GNU General Public License
00031  * along with this program; if not, write to the Free Software
00032  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00033  *
00034  * */
00035 
00036 #include "averager.h"
00037 #include "emdata.h"
00038 #include "xydata.h"
00039 #include "ctf.h"
00040 #include <cstring>
00041 #include "plugins/averager_template.h"
00042 
00043 using namespace EMAN;
00044 
00045 const string ImageAverager::NAME = "mean";
00046 const string TomoAverager::NAME = "mean.tomo";
00047 const string MinMaxAverager::NAME = "minmax";
00048 const string AbsMaxMinAverager::NAME = "absmaxmin";
00049 const string IterationAverager::NAME = "iteration";
00050 const string CtfCWautoAverager::NAME = "ctfw.auto";
00051 const string CtfCAutoAverager::NAME = "ctf.auto";
00052 const string FourierWeightAverager::NAME = "weightedfourier";
00053 
00054 template <> Factory < Averager >::Factory()
00055 {
00056         force_add<ImageAverager>();
00057         force_add<MinMaxAverager>();
00058         force_add<AbsMaxMinAverager>();
00059         force_add<IterationAverager>();
00060         force_add<CtfCWautoAverager>();
00061         force_add<CtfCAutoAverager>();
00062         force_add<TomoAverager>();
00063         force_add<FourierWeightAverager>();
00064 //      force_add<XYZAverager>();
00065 }
00066 
00067 void Averager::mult(const float& s)
00068 {
00069         if ( result != 0 )
00070         {
00071                 result->mult(s);
00072         }
00073         else throw NullPointerException("Error, attempted to multiply the result image, which is NULL");
00074 }
00075 
00076 
00077 void Averager::add_image_list(const vector<EMData*> & image_list)
00078 {
00079         for (size_t i = 0; i < image_list.size(); i++) {
00080                 add_image(image_list[i]);
00081         }
00082 }
00083 
00084 TomoAverager::TomoAverager()
00085         : norm_image(0)
00086 {
00087 
00088 }
00089 
00090 void TomoAverager::add_image(EMData * image)
00091 {
00092         if (!image) {
00093                 return;
00094         }
00095 
00096         if (!image->is_complex()) {
00097                 image=image->do_fft();
00098                 image->set_attr("free_me",(int)1);
00099         }
00100                 
00101         if (result!=0 && !EMUtil::is_same_size(image, result)) {
00102                 LOGERR("%s Averager can only process same-size Images",
00103                            get_name().c_str());
00104                 return;
00105         }
00106 
00107         int nx = image->get_xsize();
00108         int ny = image->get_ysize();
00109         int nz = image->get_zsize();
00110         size_t image_size = (size_t)nx * ny * nz;
00111 
00112         if (norm_image == 0) {
00113                 printf("init average %d %d %d",nx,ny,nz);
00114                 result = image->copy_head();
00115                 result->to_zero();
00116                 
00117                 norm_image = image->copy_head();
00118                 norm_image->to_zero();
00119                 
00120                 thresh_sigma = (float)params.set_default("thresh_sigma", 0.01);
00121         }
00122 
00123         float *result_data = result->get_data();
00124         float *norm_data = norm_image->get_data();
00125         float *data = image->get_data();
00126         float thresh=image->get_attr("sigma");
00127         thresh=2.0f*thresh*thresh*thresh_sigma;
00128         
00129         // Add any values above threshold to the result image, and add 1 to the corresponding pixels in the norm image
00130         for (size_t j = 0; j < image_size; j+=2) {
00131                 float f=data[j];        // real
00132                 float g=data[j+1];      // imag
00133                 float inten=f*f+g*g;
00134                 
00135                 if (inten<thresh) continue;
00136                 
00137                 result_data[j]  +=f;
00138                 result_data[j+1]+=g;
00139                 
00140                 norm_data[j]  +=1.0;
00141                 norm_data[j+1]+=1.0;
00142         }
00143         
00144         if (image->has_attr("free_me")) delete image;
00145 }
00146 
00147 EMData * TomoAverager::finish()
00148 {
00149         if (norm_image==0 || result==0) return NULL;
00150         
00151         int nx = result->get_xsize();
00152         int ny = result->get_ysize();
00153         int nz = result->get_zsize();
00154         size_t image_size = (size_t)nx * ny * nz;
00155         
00156         float *result_data = result->get_data();
00157         float *norm_data = norm_image->get_data();
00158         
00159         printf("finish average %d %d %d",nx,ny,nz);
00160         // normalize the average
00161         for (size_t j = 0; j < image_size; j++) {
00162                 if (norm_data[j]==0.0) result_data[j]=0.0;
00163                 else result_data[j]/=norm_data[j];
00164         }
00165         
00166         norm_image->update();
00167         result->update();
00168         
00169         EMData *ret = result->do_ift();
00170         ret->set_attr("ptcl_repr",norm_image->get_attr("maximum"));
00171         if ((int)params.set_default("save_norm", 0)) 
00172                 norm_image->write_image("norm.hdf");
00173         
00174         delete result;
00175         delete norm_image;
00176         result=0;
00177         norm_image=0;
00178         
00179         return ret;
00180 }
00181 
00182 
00183 ImageAverager::ImageAverager()
00184         : sigma_image(0), ignore0(0), normimage(0), freenorm(0), nimg(0)
00185 {
00186 
00187 }
00188 
00189 void ImageAverager::add_image(EMData * image)
00190 {
00191         if (!image) {
00192                 return;
00193         }
00194 
00195         if (nimg >= 1 && !EMUtil::is_same_size(image, result)) {
00196                 LOGERR("%sAverager can only process same-size Image",
00197                            get_name().c_str());
00198                 return;
00199         }
00200 
00201         nimg++;
00202 
00203         int nx = image->get_xsize();
00204         int ny = image->get_ysize();
00205         int nz = image->get_zsize();
00206         size_t image_size = (size_t)nx * ny * nz;
00207 
00208         if (nimg == 1) {
00209                 result = image->copy_head();
00210                 result->set_size(nx, ny, nz);
00211                 sigma_image = params.set_default("sigma", (EMData*)0);
00212                 ignore0 = params["ignore0"];
00213 
00214                 normimage = params.set_default("normimage", (EMData*)0);
00215                 if (ignore0 && normimage==0) { normimage=new EMData(nx,ny,nz); freenorm=1; }
00216                 if (normimage) normimage->to_zero();
00217         }
00218 
00219         float *result_data = result->get_data();
00220         float *sigma_image_data = 0;
00221         if (sigma_image) {
00222                 sigma_image->set_size(nx, ny, nz);
00223                 sigma_image_data = sigma_image->get_data();
00224         }
00225 
00226         float * image_data = image->get_data();
00227 
00228         if (!ignore0) {
00229                 for (size_t j = 0; j < image_size; ++j) {
00230                         float f = image_data[j];
00231                         result_data[j] += f;
00232                         if (sigma_image_data) {
00233                                 sigma_image_data[j] += f * f;
00234                         }
00235                 }
00236         }
00237         else {
00238                 for (size_t j = 0; j < image_size; ++j) {
00239                         float f = image_data[j];
00240                         if (f) {
00241                                 result_data[j] += f;
00242                                 if (sigma_image_data) {
00243                                         sigma_image_data[j] += f * f;
00244                                 }
00245                                 normimage->set_value_at_fast(j,normimage->get_value_at(j)+1.0);
00246                         }
00247                 }
00248         }
00249 }
00250 
00251 EMData * ImageAverager::finish()
00252 {
00253         if (result && nimg > 1) {
00254                 size_t image_size = (size_t)result->get_xsize() * result->get_ysize() * result->get_zsize();
00255                 float * result_data = result->get_data();
00256 
00257                 if (!ignore0) {
00258                         for (size_t j = 0; j < image_size; ++j) {
00259                                 result_data[j] /= nimg;
00260                         }
00261 
00262                         if (sigma_image) {
00263                                 float * sigma_image_data = sigma_image->get_data();
00264 
00265                                 for (size_t j = 0; j < image_size; ++j) {
00266                                         float f1 = sigma_image_data[j] / nimg;
00267                                         float f2 = result_data[j];
00268                                         sigma_image_data[j] = sqrt(f1 - f2 * f2);
00269                                 }
00270 
00271                                 sigma_image->update();
00272                         }
00273                 }
00274                 else {
00275                         for (size_t j = 0; j < image_size; ++j) {
00276                                 if (normimage->get_value_at(j)>0) result_data[j] /= normimage->get_value_at(j);
00277                         }
00278                         if (sigma_image) {
00279                                 float * sigma_image_data = sigma_image->get_data();
00280 
00281                                 for (size_t j = 0; j < image_size; ++j) {
00282                                         float f1 = 0;
00283                                         if (normimage->get_value_at(j)>0) f1=sigma_image_data[j] / normimage->get_value_at(j);
00284                                         float f2 = result_data[j];
00285                                         sigma_image_data[j] = sqrt(f1 - f2 * f2);
00286                                 }
00287 
00288                                 sigma_image->update();
00289                         }
00290                 }
00291 
00292                 result->update();
00293 
00294         }               
00295         result->set_attr("ptcl_repr",nimg);
00296 
00297         if (freenorm) { delete normimage; normimage=(EMData*)0; }
00298 
00299         return result;
00300 }
00301 
00302 FourierWeightAverager::FourierWeightAverager()
00303         : normimage(0), freenorm(0), nimg(0)
00304 {
00305 
00306 }
00307 
00308 void FourierWeightAverager::add_image(EMData * image)
00309 {
00310         if (!image) {
00311                 return;
00312         }
00313 
00314 
00315 
00316         EMData *img=image->do_fft();
00317         if (nimg >= 1 && !EMUtil::is_same_size(img, result)) {
00318                 LOGERR("%sAverager can only process same-size Image",
00319                            get_name().c_str());
00320                 return;
00321         }
00322         
00323         nimg++;
00324 
00325         int nx = img->get_xsize();
00326         int ny = img->get_ysize();
00327         int nz = 1;
00328 //      size_t image_size = (size_t)nx * ny * nz;
00329 
00330         XYData *weight=(XYData *)image->get_attr("avg_weight");
00331         
00332         if (nimg == 1) {
00333                 result = new EMData(nx,ny,nz);
00334                 result->set_complex(true);
00335                 result->to_zero();
00336 
00337                 normimage = params.set_default("normimage", (EMData*)0);
00338                 if (normimage==0) { normimage=new EMData(nx/2,ny,nz); freenorm=1; }
00339                 normimage->to_zero();
00340         }
00341 
00342         // We're using routines that handle complex image wraparound for us, so we iterate over the half-plane
00343         for (int y=-ny/2; y<ny/2; y++) {
00344                 for (int x=0; x<nx/2; x++) {
00345                         std::complex<float> v=img->get_complex_at(x,y);
00346                         float r=Util::hypot2(y/(float)ny,x/(float)nx);
00347                         float wt=weight->get_yatx(r);
00348                         result->set_complex_at(x,y,result->get_complex_at(x,y)+v*wt);
00349                         normimage->set_value_at(x,y+ny/2,normimage->get_value_at(x,y+ny/2)+wt);
00350                 }
00351         }
00352 
00353         delete img;
00354 }
00355 
00356 EMData * FourierWeightAverager::finish()
00357 {
00358         EMData *ret = (EMData *)0;
00359         
00360         if (result && nimg >= 1) {
00361         // We're using routines that handle complex image wraparound for us, so we iterate over the half-plane
00362                 int nx = result->get_xsize();
00363                 int ny = result->get_ysize();
00364                 
00365                 for (int y=-ny/2; y<ny/2; y++) {
00366                         for (int x=0; x<nx/2; x++) {
00367                                 float norm=normimage->get_value_at(x,y+ny/2);
00368                                 if (norm<=0) result->set_complex_at(x,y,0.0f);
00369                                 else result->set_complex_at(x,y,result->get_complex_at(x,y)/norm);
00370                         }
00371                 }
00372 
00373                 result->update();
00374 //              result->mult(1.0f/(float)result->get_attr("sigma"));
00375 //              result->write_image("tmp.hdf",0);
00376 //              printf("%g %g %g\n",(float)result->get_attr("sigma"),(float)result->get_attr("minimum"),(float)result->get_attr("maximum"));
00377                 ret=result->do_ift();
00378                 delete result;
00379                 result=(EMData*) 0;
00380         }
00381         ret->set_attr("ptcl_repr",nimg);
00382 
00383         if (freenorm) { delete normimage; normimage=(EMData*)0; }
00384         nimg=0;
00385 
00386         return ret;
00387 }
00388 
00389 
00390 #if 0
00391 EMData *ImageAverager::average(const vector < EMData * >&image_list) const
00392 {
00393         if (image_list.size() == 0) {
00394                 return 0;
00395         }
00396 
00397         EMData *sigma_image = params["sigma"];
00398         int ignore0 = params["ignore0"];
00399 
00400         EMData *image0 = image_list[0];
00401 
00402         int nx = image0->get_xsize();
00403         int ny = image0->get_ysize();
00404         int nz = image0->get_zsize();
00405         size_t image_size = (size_t)nx * ny * nz;
00406 
00407         EMData *result = new EMData();
00408         result->set_size(nx, ny, nz);
00409 
00410         float *result_data = result->get_data();
00411         float *sigma_image_data = 0;
00412 
00413         if (sigma_image) {
00414                 sigma_image->set_size(nx, ny, nz);
00415                 sigma_image_data = sigma_image->get_data();
00416         }
00417 
00418         int c = 1;
00419         if (ignore0) {
00420                 for (size_t j = 0; j < image_size; ++j) {
00421                         int g = 0;
00422                         for (size_t i = 0; i < image_list.size(); i++) {
00423                                 float f = image_list[i]->get_value_at(j);
00424                                 if (f) {
00425                                         g++;
00426                                         result_data[j] += f;
00427                                         if (sigma_image_data) {
00428                                                 sigma_image_data[j] += f * f;
00429                                         }
00430                                 }
00431                         }
00432                         if (g) {
00433                                 result_data[j] /= g;
00434                         }
00435                 }
00436         }
00437         else {
00438                 float *image0_data = image0->get_data();
00439                 if (sigma_image_data) {
00440                         memcpy(sigma_image_data, image0_data, image_size * sizeof(float));
00441 
00442                         for (size_t j = 0; j < image_size; ++j) {
00443                                 sigma_image_data[j] *= sigma_image_data[j];
00444                         }
00445                 }
00446 
00447                 image0->update();
00448                 memcpy(result_data, image0_data, image_size * sizeof(float));
00449 
00450                 for (size_t i = 1; i < image_list.size(); i++) {
00451                         EMData *image = image_list[i];
00452 
00453                         if (EMUtil::is_same_size(image, result)) {
00454                                 float *image_data = image->get_data();
00455 
00456                                 for (size_t j = 0; j < image_size; ++j) {
00457                                         result_data[j] += image_data[j];
00458                                 }
00459 
00460                                 if (sigma_image_data) {
00461                                         for (size_t j = 0; j < image_size; ++j) {
00462                                                 sigma_image_data[j] += image_data[j] * image_data[j];
00463                                         }
00464                                 }
00465 
00466                                 image->update();
00467                                 c++;
00468                         }
00469                 }
00470 
00471                 for (size_t j = 0; j < image_size; ++j) {
00472                         result_data[j] /= static_cast < float >(c);
00473                 }
00474         }
00475 
00476         if (sigma_image_data) {
00477                 for (size_t j = 0; j < image_size; ++j) {
00478                         float f1 = sigma_image_data[j] / c;
00479                         float f2 = result_data[j];
00480                         sigma_image_data[j] = sqrt(f1 - f2 * f2);
00481                 }
00482         }
00483 
00484         if (sigma_image_data) {
00485                 sigma_image->update();
00486         }
00487 
00488         result->update();
00489         return result;
00490 }
00491 #endif
00492 
00493 MinMaxAverager::MinMaxAverager()
00494         : nimg(0)
00495 {
00496         /*move max out of initializer list, since this max(0) is considered as a macro
00497          * in Visual Studio, which we defined somewhere else*/
00498         max = 0;
00499 }
00500 
00501 void MinMaxAverager::add_image(EMData * image)
00502 {
00503         if (!image) {
00504                 return;
00505         }
00506 
00507         if (nimg >= 1 && !EMUtil::is_same_size(image, result)) {
00508                 LOGERR("%sAverager can only process same-size Image",
00509                            get_name().c_str());
00510                 return;
00511         }
00512 
00513         nimg++;
00514 
00515         int nx = image->get_xsize();
00516         int ny = image->get_ysize();
00517         int nz = image->get_zsize();
00518 
00519         if (nimg == 1) {
00520                 result = image->copy();
00521                 max = params["max"];
00522                 return;
00523         }
00524 
00525         for (int z=0; z<nz; z++) {
00526                 for (int y=0; y<ny; y++) {
00527                         for (int x=0; x<nx; x++) {
00528                                 if (result->get_value_at(x,y,z)>image->get_value_at(x,y,z))
00529                                         { if (!max) result->set_value_at(x,y,z,image->get_value_at(x,y,z)); }
00530                                 else { if (max) result->set_value_at(x,y,z,image->get_value_at(x,y,z)); }
00531                         }
00532                 }
00533         }
00534 
00535 }
00536 
00537 EMData *MinMaxAverager::finish()
00538 {
00539         result->update();
00540         result->set_attr("ptcl_repr",nimg);
00541         
00542         if (result && nimg > 1) return result;
00543 
00544         return NULL;
00545 }
00546 
00547 AbsMaxMinAverager::AbsMaxMinAverager() : nimg(0)
00548 {
00549         /*move max out of initializer list, since this max(0) is considered as a macro
00550          * in Visual Studio, which we defined somewhere else*/
00551         min = 0;
00552 }
00553 
00554 void AbsMaxMinAverager::add_image(EMData * image)
00555 {
00556         if (!image) {
00557                 return;
00558         }
00559 
00560         if (nimg >= 1 && !EMUtil::is_same_size(image, result)) {
00561                 LOGERR("%sAverager can only process same-size Image",
00562                            get_name().c_str());
00563                 return;
00564         }
00565 
00566         nimg++;
00567 
00568         int nx = image->get_xsize();
00569         int ny = image->get_ysize();
00570         int nz = image->get_zsize();
00571 
00572         size_t imgsize = (size_t)nx*ny*nz;
00573 
00574         if (nimg == 1) {
00575                 result = image->copy();
00576                 min = params["min"];
00577                 return;
00578         }
00579 
00580         float * data     = result->get_data();
00581         float * src_data = image->get_data();
00582 
00583         for(size_t i=0; i<imgsize; ++i) {
00584                 if(!min) {      //average to maximum by default
00585                         if (fabs(data[i]) < fabs(src_data[i])) data[i]=src_data[i];
00586                 }
00587                 else {  //average to minimum if set 'min'
00588                         if (fabs(data[i]) > fabs(src_data[i])) data[i]=src_data[i];
00589                 }
00590         }
00591 }
00592 
00593 EMData *AbsMaxMinAverager::finish()
00594 {
00595         result->update();
00596         result->set_attr("ptcl_repr",nimg);
00597 
00598         if (result && nimg > 1) return result;
00599 
00600         return NULL;
00601 }
00602 
00603 IterationAverager::IterationAverager() : nimg(0)
00604 {
00605 
00606 }
00607 
00608 void IterationAverager::add_image( EMData * image)
00609 {
00610         if (!image) {
00611                 return;
00612         }
00613 
00614         if (nimg >= 1 && !EMUtil::is_same_size(image, result)) {
00615                 LOGERR("%sAverager can only process same-size Image",
00616                                                          get_name().c_str());
00617                 return;
00618         }
00619 
00620         nimg++;
00621 
00622         int nx = image->get_xsize();
00623         int ny = image->get_ysize();
00624         int nz = image->get_zsize();
00625         size_t image_size = (size_t)nx * ny * nz;
00626 
00627         if (nimg == 1) {
00628                 result = image->copy_head();
00629                 result->set_size(nx, ny, nz);
00630                 sigma_image = image->copy_head();
00631                 sigma_image->set_size(nx, ny, nz);
00632         }
00633 
00634         float *image_data = image->get_data();
00635         float *result_data = result->get_data();
00636         float *sigma_image_data = sigma_image->get_data();
00637 
00638         for (size_t j = 0; j < image_size; ++j) {
00639                 float f = image_data[j];
00640                 result_data[j] += f;
00641                 sigma_image_data[j] += f * f;
00642         }
00643 
00644 
00645 }
00646 
00647 EMData * IterationAverager::finish()
00648 {
00649         if (nimg < 1) {
00650                 return result;
00651         }
00652 
00653         int nx = result->get_xsize();
00654         int ny = result->get_ysize();
00655         int nz = result->get_zsize();
00656         size_t image_size = (size_t)nx * ny * nz;
00657 
00658         float *result_data = result->get_data();
00659         float *sigma_image_data = sigma_image->get_data();
00660 
00661         for (size_t j = 0; j < image_size; ++j) {
00662                 result_data[j] /= nimg;
00663                 float f1 = sigma_image_data[j] / nimg;
00664                 float f2 = result_data[j];
00665                 sigma_image_data[j] = sqrt(f1 - f2 * f2) / sqrt((float)nimg);
00666         }
00667 
00668         result->update();
00669         sigma_image->update();
00670 
00671         result->append_image("iter.hed");
00672         float sigma = sigma_image->get_attr("sigma");
00673         float *sigma_image_data2 = sigma_image->get_data();
00674         float *result_data2 = result->get_data();
00675         float *d2 = new float[nx * ny];
00676         size_t sec_size = nx * ny * sizeof(float);
00677 
00678         memcpy(d2, result_data2, sec_size);
00679         memcpy(sigma_image_data2, result_data2, sec_size);
00680 
00681         printf("Iter sigma=%f\n", sigma);
00682 
00683         for (int k = 0; k < 1000; k++) {
00684                 for (int i = 1; i < nx - 1; i++) {
00685                         for (int j = 1; j < ny - 1; j++) {
00686                                 int l = i + j * nx;
00687                                 float c1 = (d2[l - 1] + d2[l + 1] + d2[l - nx] + d2[l + nx]) / 4.0f - d2[l];
00688                                 float c2 = fabs(result_data2[l] - sigma_image_data2[l]) / sigma;
00689                                 result_data2[l] += c1 * Util::eman_erfc(c2) / 100.0f;
00690                         }
00691                 }
00692 
00693                 memcpy(d2, result_data2, sec_size);
00694         }
00695 
00696         if( d2 )
00697         {
00698                 delete[]d2;
00699                 d2 = 0;
00700         }
00701 
00702         sigma_image->update();
00703         if( sigma_image )
00704         {
00705                 delete sigma_image;
00706                 sigma_image = 0;
00707         }
00708 
00709         result->update();
00710         result->append_image("iter.hed");
00711 
00712 
00713         return result;
00714 }
00715 
00716 CtfCWautoAverager::CtfCWautoAverager()
00717         : nimg(0)
00718 {
00719 
00720 }
00721 
00722 
00723 void CtfCWautoAverager::add_image(EMData * image)
00724 {
00725         if (!image) {
00726                 return;
00727         }
00728 
00729 
00730 
00731         EMData *fft=image->do_fft();
00732 
00733         if (nimg >= 1 && !EMUtil::is_same_size(fft, result)) {
00734                 LOGERR("%s Averager can only process images of the same size", get_name().c_str());
00735                 return;
00736         }
00737 
00738         nimg++;
00739         if (nimg == 1) {
00740                 result = fft->copy_head();
00741                 result->to_zero();
00742         }
00743 
00744         Ctf *ctf = (Ctf *)image->get_attr("ctf");
00745 //string cc=ctf->to_string();
00746 //FILE *out=fopen("ctf.txt","a");
00747 //fprintf(out,"%s\n",cc.c_str());
00748 //fclose(out);
00749         float b=ctf->bfactor;
00750         ctf->bfactor=100.0;                     // FIXME - this is a temporary fixed B-factor which does a (very) little sharpening
00751 
00752 //      if (nimg==1) unlink("snr.hdf");
00753 
00754         EMData *snr = result -> copy();
00755         ctf->compute_2d_complex(snr,Ctf::CTF_SNR);
00756 //      snr->write_image("snr.hdf",-1);
00757         EMData *ctfi = result-> copy();
00758         ctf->compute_2d_complex(ctfi,Ctf::CTF_AMP);
00759 
00760         ctf->bfactor=b; // return to its original value
00761 
00762         float *outd = result->get_data();
00763         float *ind = fft->get_data();
00764         float *snrd = snr->get_data();
00765         float *ctfd = ctfi->get_data();
00766 
00767         size_t sz=snr->get_xsize()*snr->get_ysize();
00768         for (size_t i = 0; i < sz; i+=2) {
00769                 if (snrd[i]<0) snrd[i]=0.001;   // Used to be 0. See ctfcauto averager
00770                 ctfd[i]=fabs(ctfd[i]);
00771                 if (ctfd[i]<.05) ctfd[i]=0.05f;
00772 //              {
00773 //                      if (snrd[i]<=0) ctfd[i]=.05f;
00774 //                      else ctfd[i]=snrd[i]*10.0f;
00775 //              }
00776                 outd[i]+=ind[i]*snrd[i]/ctfd[i];
00777                 outd[i+1]+=ind[i+1]*snrd[i]/ctfd[i];
00778         }
00779 
00780         if (nimg==1) {
00781                 snrsum=snr->copy_head();
00782                 float *ssnrd=snrsum->get_data();
00783                 // we're only using the real component, and we need to start with 1.0
00784                 for (size_t i = 0; i < sz; i+=2) { ssnrd[i]=1.0; ssnrd[i+1]=0.0; }
00785         }
00786 //      snr->write_image("snr.hdf",-1);
00787         snr->process_inplace("math.absvalue");
00788         snrsum->add(*snr);
00789 
00790         delete ctf;
00791         delete fft;
00792         delete snr;
00793         delete ctfi;
00794 }
00795 
00796 EMData * CtfCWautoAverager::finish()
00797 {
00798 /*      EMData *tmp=result->do_ift();
00799         tmp->write_image("ctfcw.hdf",0);
00800         delete tmp;
00801 
00802         tmp=snrsum->do_ift();
00803         tmp->write_image("ctfcw.hdf",1);
00804         delete tmp;*/
00805 
00806 //snrsum->write_image("snrsum.hdf",-1);
00807         //size_t sz=result->get_xsize()*result->get_ysize();
00808         int nx=result->get_xsize();
00809         int ny=result->get_ysize();     
00810         float *snrsd=snrsum->get_data();
00811         float *outd=result->get_data();
00812 
00813         int rm=(ny-2)*(ny-2)/4;
00814         for (int j=0; j<ny; j++) {
00815                 for (int i=0; i<nx; i+=2) {
00816                         size_t ii=i+j*nx;
00817                         if ((j<ny/2 && i*i/4+j*j>rm) ||(j>=ny/2 && i*i/4+(ny-j)*(ny-j)>rm) || snrsd[ii]==0) { outd[ii]=outd[ii+1]=0; continue; }
00818                         outd[ii]/=snrsd[ii];            // snrsd contains total SNR
00819                         outd[ii+1]/=snrsd[ii];
00820                 }
00821         }
00822 
00823         result->update();
00824         result->set_attr("ptcl_repr",nimg);
00825         result->set_attr("ctf_snr_total",snrsum->calc_radial_dist(snrsum->get_ysize()/2,0,1,false));
00826         result->set_attr("ctf_wiener_filtered",1);
00827         
00828         delete snrsum;
00829         EMData *ret=result->do_ift();
00830         delete result;
00831         result=NULL;
00832         return ret;
00833 }
00834 
00835 CtfCAutoAverager::CtfCAutoAverager()
00836         : nimg(0)
00837 {
00838 
00839 }
00840 
00841 
00842 void CtfCAutoAverager::add_image(EMData * image)
00843 {
00844         if (!image) {
00845                 return;
00846         }
00847 
00848 
00849 
00850         EMData *fft=image->do_fft();
00851 
00852         if (nimg >= 1 && !EMUtil::is_same_size(fft, result)) {
00853                 LOGERR("%s Averager can only process images of the same size", get_name().c_str());
00854                 return;
00855         }
00856 
00857         nimg++;
00858         if (nimg == 1) {
00859                 result = fft->copy_head();
00860                 result->to_zero();
00861         }
00862 
00863         Ctf *ctf = (Ctf *)image->get_attr("ctf");
00864         float b=ctf->bfactor;
00865         ctf->bfactor=0;                 // NO B-FACTOR CORRECTION !
00866 
00867         EMData *snr = result -> copy();
00868         ctf->compute_2d_complex(snr,Ctf::CTF_SNR);
00869         EMData *ctfi = result-> copy();
00870         ctf->compute_2d_complex(ctfi,Ctf::CTF_AMP);
00871 
00872         ctf->bfactor=b; // return to its original value
00873 
00874         float *outd = result->get_data();
00875         float *ind = fft->get_data();
00876         float *snrd = snr->get_data();
00877         float *ctfd = ctfi->get_data();
00878 
00879         size_t sz=snr->get_xsize()*snr->get_ysize();
00880         for (size_t i = 0; i < sz; i+=2) {
00881                 if (snrd[i]<=0) snrd[i]=0.001;          // used to be 0. Trying to insure that there is always at least a little signal used. In cases with dense particles, SNR may be dramatically underestimated
00882                 ctfd[i]=fabs(ctfd[i]);
00883                 
00884                 // This limits the maximum possible amplification in CTF correction to 10x
00885                 if (ctfd[i]<.05)  ctfd[i]=0.05f;
00886 //              {
00887 //                      if (snrd[i]<=0) ctfd[i]=.05f;
00888 //                      else ctfd[i]=snrd[i]*10.0f;
00889 //              }
00890                 
00891                 // SNR weight and CTF correction
00892                 outd[i]+=ind[i]*snrd[i]/ctfd[i];
00893                 outd[i+1]+=ind[i+1]*snrd[i]/ctfd[i];
00894         }
00895 
00896         if (nimg==1) {
00897                 snrsum=snr->copy_head();
00898                 float *ssnrd=snrsum->get_data();
00899                 // we're only using the real component, for Wiener filter we put 1.0 in R, but for just SNR weight we use 0
00900                 for (size_t i = 0; i < sz; i+=2) { ssnrd[i]=0.0; ssnrd[i+1]=0.0; }
00901         }
00902         snr->process_inplace("math.absvalue");
00903         snrsum->add(*snr);
00904 
00905         delete ctf;
00906         delete fft;
00907         delete snr;
00908         delete ctfi;
00909 }
00910 
00911 EMData * CtfCAutoAverager::finish()
00912 {
00913 /*      EMData *tmp=result->do_ift();
00914         tmp->write_image("ctfcw.hdf",0);
00915         delete tmp;
00916 
00917         tmp=snrsum->do_ift();
00918         tmp->write_image("ctfcw.hdf",1);
00919         delete tmp;*/
00920 
00921 //      snrsum->write_image("snrsum.hdf",-1);
00922         //size_t sz=result->get_xsize()*result->get_ysize();
00923         int nx=result->get_xsize();
00924         int ny=result->get_ysize();     
00925         float *snrsd=snrsum->get_data();
00926         float *outd=result->get_data();
00927 
00928         int rm=(ny-2)*(ny-2)/4;
00929         for (int j=0; j<ny; j++) {
00930                 for (int i=0; i<nx; i+=2) {
00931                         size_t ii=i+j*nx;
00932                         if ((j<ny/2 && i*i/4+j*j>rm) ||(j>=ny/2 && i*i/4+(ny-j)*(ny-j)>rm) || snrsd[ii]==0) { outd[ii]=outd[ii+1]=0; continue; }
00933                         // we aren't wiener filtering, but if the total SNR is too low, we don't want TOO much exaggeration of noise
00934                         if (snrsd[ii]<.05) {            
00935                                 outd[ii]*=20.0;         // 1/0.05
00936                                 outd[ii+1]*=20.0;
00937                         }
00938                         else {
00939                                 outd[ii]/=snrsd[ii];            // snrsd contains total SNR
00940                                 outd[ii+1]/=snrsd[ii];
00941                         }
00942                 }
00943         }
00944         result->update();
00945         result->set_attr("ptcl_repr",nimg);
00946         result->set_attr("ctf_snr_total",snrsum->calc_radial_dist(snrsum->get_ysize()/2,0,1,false));
00947         result->set_attr("ctf_wiener_filtered",0);
00948         
00949 /*      snrsum->write_image("snr.hdf",-1);
00950         result->write_image("avg.hdf",-1);*/
00951         
00952         delete snrsum;
00953         EMData *ret=result->do_ift();
00954         delete result;
00955         result=NULL;
00956         return ret;
00957 }
00958 
00959 #if 0
00960 EMData *IterationAverager::average(const vector < EMData * >&image_list) const
00961 {
00962         if (image_list.size() == 0) {
00963                 return 0;
00964         }
00965 
00966         EMData *image0 = image_list[0];
00967 
00968         int nx = image0->get_xsize();
00969         int ny = image0->get_ysize();
00970         int nz = image0->get_zsize();
00971         size_t image_size = (size_t)nx * ny * nz;
00972 
00973         EMData *result = new EMData();
00974         result->set_size(nx, ny, nz);
00975 
00976         EMData *sigma_image = new EMData();
00977         sigma_image->set_size(nx, ny, nz);
00978 
00979         float *image0_data = image0->get_data();
00980         float *result_data = result->get_data();
00981         float *sigma_image_data = sigma_image->get_data();
00982 
00983         memcpy(result_data, image0_data, image_size * sizeof(float));
00984         memcpy(sigma_image_data, image0_data, image_size * sizeof(float));
00985 
00986         for (size_t j = 0; j < image_size; ++j) {
00987                 sigma_image_data[j] *= sigma_image_data[j];
00988         }
00989 
00990         image0->update();
00991 
00992         int nc = 1;
00993         for (size_t i = 1; i < image_list.size(); i++) {
00994                 EMData *image = image_list[i];
00995 
00996                 if (EMUtil::is_same_size(image, result)) {
00997                         float *image_data = image->get_data();
00998 
00999                         for (size_t j = 0; j < image_size; j++) {
01000                                 result_data[j] += image_data[j];
01001                                 sigma_image_data[j] += image_data[j] * image_data[j];
01002                         }
01003 
01004                         image->update();
01005                         nc++;
01006                 }
01007         }
01008 
01009         float c = static_cast < float >(nc);
01010         for (size_t j = 0; j < image_size; ++j) {
01011                 float f1 = sigma_image_data[j] / c;
01012                 float f2 = result_data[j] / c;
01013                 sigma_image_data[j] = sqrt(f1 - f2 * f2) / sqrt(c);
01014         }
01015 
01016 
01017         for (size_t j = 0; j < image_size; ++j) {
01018                 result_data[j] /= c;
01019         }
01020 
01021         result->update();
01022         sigma_image->update();
01023 
01024         result->append_image("iter.hed");
01025 
01026         float sigma = sigma_image->get_attr("sigma");
01027         float *sigma_image_data2 = sigma_image->get_data();
01028         float *result_data2 = result->get_data();
01029         float *d2 = new float[nx * ny];
01030         size_t sec_size = nx * ny * sizeof(float);
01031 
01032         memcpy(d2, result_data2, sec_size);
01033         memcpy(sigma_image_data2, result_data2, sec_size);
01034 
01035         printf("Iter sigma=%f\n", sigma);
01036 
01037         for (int k = 0; k < 1000; k++) {
01038                 for (int i = 1; i < nx - 1; i++) {
01039                         for (int j = 1; j < ny - 1; j++) {
01040                                 int l = i + j * nx;
01041                                 float c1 = (d2[l - 1] + d2[l + 1] + d2[l - nx] + d2[l + nx]) / 4.0f - d2[l];
01042                                 float c2 = fabs(result_data2[l] - sigma_image_data2[l]) / sigma;
01043                                 result_data2[l] += c1 * Util::eman_erfc(c2) / 100.0f;
01044                         }
01045                 }
01046 
01047                 memcpy(d2, result_data2, sec_size);
01048         }
01049 
01050         if( d2 )
01051         {
01052                 delete[]d2;
01053                 d2 = 0;
01054         }
01055 
01056         sigma_image->update();
01057         if( sigma_image )
01058         {
01059                 delete sigma_image;
01060                 sigma_image = 0;
01061         }
01062 
01063         result->update();
01064         result->append_image("iter.hed");
01065 
01066         return result;
01067 }
01068 #endif
01069 
01070 
01071 
01072 
01073 void EMAN::dump_averagers()
01074 {
01075         dump_factory < Averager > ();
01076 }
01077 
01078 map<string, vector<string> > EMAN::dump_averagers_list()
01079 {
01080         return dump_factory_list < Averager > ();
01081 }