Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

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 IterationAverager::NAME = "iteration";
00049 const string WeightingAverager::NAME = "snr_weight";
00050 const string CtfCAverager::NAME = "ctfc";
00051 const string CtfCWAverager::NAME = "ctfcw";
00052 const string CtfCWautoAverager::NAME = "ctfw.auto";
00053 const string CtfCAutoAverager::NAME = "ctf.auto";
00054 
00055 template <> Factory < Averager >::Factory()
00056 {
00057         force_add<ImageAverager>();
00058         force_add<MinMaxAverager>();
00059         force_add<IterationAverager>();
00060 //      force_add<WeightingAverager>();
00061         // These commented out until we're happy they're working. (d.woolford, Feb 3rd 2009)
00062 //      force_add(&CtfCAverager::NEW);
00063 //      force_add(&CtfCWAverager::NEW);
00064         force_add<CtfCWautoAverager>();
00065         force_add<CtfCAutoAverager>();
00066         force_add<TomoAverager>();
00067 //      force_add<XYZAverager>();
00068 }
00069 
00070 void Averager::mult(const float& s)
00071 {
00072         if ( result != 0 )
00073         {
00074                 result->mult(s);
00075         }
00076         else throw NullPointerException("Error, attempted to multiply the result image, which is NULL");
00077 }
00078 
00079 
00080 void Averager::add_image_list(const vector<EMData*> & image_list)
00081 {
00082         for (size_t i = 0; i < image_list.size(); i++) {
00083                 add_image(image_list[i]);
00084         }
00085 }
00086 
00087 TomoAverager::TomoAverager()
00088         : norm_image(0)
00089 {
00090 
00091 }
00092 
00093 void TomoAverager::add_image(EMData * image)
00094 {
00095         if (!image) {
00096                 return;
00097         }
00098 
00099         if (!image->is_complex()) {
00100                 image=image->do_fft();
00101                 image->set_attr("free_me",(int)1);
00102         }
00103                 
00104         if (result!=0 && !EMUtil::is_same_size(image, result)) {
00105                 LOGERR("%s Averager can only process same-size Images",
00106                            get_name().c_str());
00107                 return;
00108         }
00109 
00110         int nx = image->get_xsize();
00111         int ny = image->get_ysize();
00112         int nz = image->get_zsize();
00113         size_t image_size = (size_t)nx * ny * nz;
00114 
00115         if (norm_image == 0) {
00116                 printf("init average %d %d %d",nx,ny,nz);
00117                 result = image->copy_head();
00118                 result->to_zero();
00119                 
00120                 norm_image = image->copy_head();
00121                 norm_image->to_zero();
00122                 
00123                 thresh_sigma = (float)params.set_default("thresh_sigma", 0.01);
00124         }
00125 
00126         float *result_data = result->get_data();
00127         float *norm_data = norm_image->get_data();
00128         float *data = image->get_data();
00129         float thresh=image->get_attr("sigma");
00130         thresh=2.0f*thresh*thresh*thresh_sigma;
00131         
00132         // Add any values above threshold to the result image, and add 1 to the corresponding pixels in the norm image
00133         for (size_t j = 0; j < image_size; j+=2) {
00134                 float f=data[j];        // real
00135                 float g=data[j+1];      // imag
00136                 float inten=f*f+g*g;
00137                 
00138                 if (inten<thresh) continue;
00139                 
00140                 result_data[j]  +=f;
00141                 result_data[j+1]+=g;
00142                 
00143                 norm_data[j]  +=1.0;
00144                 norm_data[j+1]+=1.0;
00145         }
00146         
00147         if (image->has_attr("free_me")) delete image;
00148 }
00149 
00150 EMData * TomoAverager::finish()
00151 {
00152         if (norm_image==0 || result==0) return NULL;
00153         
00154         int nx = result->get_xsize();
00155         int ny = result->get_ysize();
00156         int nz = result->get_zsize();
00157         size_t image_size = (size_t)nx * ny * nz;
00158         
00159         float *result_data = result->get_data();
00160         float *norm_data = norm_image->get_data();
00161         
00162         printf("finish average %d %d %d",nx,ny,nz);
00163         // normalize the average
00164         for (size_t j = 0; j < image_size; j++) {
00165                 if (norm_data[j]==0.0) result_data[j]=0.0;
00166                 else result_data[j]/=norm_data[j];
00167         }
00168         
00169         norm_image->update();
00170         result->update();
00171         
00172         EMData *ret = result->do_ift();
00173         ret->set_attr("ptcl_repr",norm_image->get_attr("maximum"));
00174         if ((int)params.set_default("save_norm", 0)) 
00175                 norm_image->write_image("norm.hdf");
00176         
00177         delete result;
00178         delete norm_image;
00179         result=0;
00180         norm_image=0;
00181         
00182         return ret;
00183 }
00184 
00185 
00186 ImageAverager::ImageAverager()
00187         : sigma_image(0), nimg_n0(0), ignore0(0), nimg(0)
00188 {
00189 
00190 }
00191 
00192 void ImageAverager::add_image(EMData * image)
00193 {
00194         if (!image) {
00195                 return;
00196         }
00197 
00198         if (nimg >= 1 && !EMUtil::is_same_size(image, result)) {
00199                 LOGERR("%sAverager can only process same-size Image",
00200                            get_name().c_str());
00201                 return;
00202         }
00203 
00204         nimg++;
00205 
00206         int nx = image->get_xsize();
00207         int ny = image->get_ysize();
00208         int nz = image->get_zsize();
00209         size_t image_size = (size_t)nx * ny * nz;
00210 
00211         if (nimg == 1) {
00212                 result = new EMData();
00213                 result->set_size(nx, ny, nz);
00214                 sigma_image = params.set_default("sigma", (EMData*)0);
00215                 ignore0 = params["ignore0"];
00216 
00217                 nimg_n0 = new int[image_size];
00218                 for (size_t i = 0; i < image_size; ++i) {
00219                         nimg_n0[i] = 0;
00220                 }
00221         }
00222 
00223         float *result_data = result->get_data();
00224         float *sigma_image_data = 0;
00225         if (sigma_image) {
00226                 sigma_image->set_size(nx, ny, nz);
00227                 sigma_image_data = sigma_image->get_data();
00228         }
00229 
00230         float * image_data = image->get_data();
00231 
00232         if (!ignore0) {
00233                 for (size_t j = 0; j < image_size; ++j) {
00234                         float f = image_data[j];
00235                         result_data[j] += f;
00236                         if (sigma_image_data) {
00237                                 sigma_image_data[j] += f * f;
00238                         }
00239                 }
00240         }
00241         else {
00242                 for (size_t j = 0; j < image_size; ++j) {
00243                         float f = image_data[j];
00244                         if (f) {
00245                                 result_data[j] += f;
00246                                 if (sigma_image_data) {
00247                                         sigma_image_data[j] += f * f;
00248                                 }
00249                                 nimg_n0[j]++;
00250                         }
00251                 }
00252         }
00253 }
00254 
00255 EMData * ImageAverager::finish()
00256 {
00257         if (result && nimg > 1) {
00258                 size_t image_size = (size_t)result->get_xsize() * result->get_ysize() * result->get_zsize();
00259                 float * result_data = result->get_data();
00260 
00261                 if (!ignore0) {
00262                         for (size_t j = 0; j < image_size; ++j) {
00263                                 result_data[j] /= nimg;
00264                         }
00265 
00266                         if (sigma_image) {
00267                                 float * sigma_image_data = sigma_image->get_data();
00268 
00269                                 for (size_t j = 0; j < image_size; ++j) {
00270                                         float f1 = sigma_image_data[j] / nimg;
00271                                         float f2 = result_data[j];
00272                                         sigma_image_data[j] = sqrt(f1 - f2 * f2);
00273                                 }
00274 
00275                                 sigma_image->update();
00276                         }
00277                 }
00278                 else {
00279                         for (size_t j = 0; j < image_size; ++j) {
00280                                 result_data[j] /= nimg_n0[j];
00281                         }
00282                         if (sigma_image) {
00283                                 float * sigma_image_data = sigma_image->get_data();
00284 
00285                                 for (size_t j = 0; j < image_size; ++j) {
00286                                         float f1 = sigma_image_data[j] / nimg_n0[j];
00287                                         float f2 = result_data[j];
00288                                         sigma_image_data[j] = sqrt(f1 - f2 * f2);
00289                                 }
00290 
00291                                 sigma_image->update();
00292                         }
00293                 }
00294 
00295                 result->update();
00296 
00297         }
00298         result->set_attr("ptcl_repr",nimg);
00299 
00300         if( nimg_n0 )
00301         {
00302                 delete [] nimg_n0;
00303                 nimg_n0 = 0;
00304         }
00305 
00306         return result;
00307 }
00308 
00309 #if 0
00310 EMData *ImageAverager::average(const vector < EMData * >&image_list) const
00311 {
00312         if (image_list.size() == 0) {
00313                 return 0;
00314         }
00315 
00316         EMData *sigma_image = params["sigma"];
00317         int ignore0 = params["ignore0"];
00318 
00319         EMData *image0 = image_list[0];
00320 
00321         int nx = image0->get_xsize();
00322         int ny = image0->get_ysize();
00323         int nz = image0->get_zsize();
00324         size_t image_size = (size_t)nx * ny * nz;
00325 
00326         EMData *result = new EMData();
00327         result->set_size(nx, ny, nz);
00328 
00329         float *result_data = result->get_data();
00330         float *sigma_image_data = 0;
00331 
00332         if (sigma_image) {
00333                 sigma_image->set_size(nx, ny, nz);
00334                 sigma_image_data = sigma_image->get_data();
00335         }
00336 
00337         int c = 1;
00338         if (ignore0) {
00339                 for (size_t j = 0; j < image_size; ++j) {
00340                         int g = 0;
00341                         for (size_t i = 0; i < image_list.size(); i++) {
00342                                 float f = image_list[i]->get_value_at(j);
00343                                 if (f) {
00344                                         g++;
00345                                         result_data[j] += f;
00346                                         if (sigma_image_data) {
00347                                                 sigma_image_data[j] += f * f;
00348                                         }
00349                                 }
00350                         }
00351                         if (g) {
00352                                 result_data[j] /= g;
00353                         }
00354                 }
00355         }
00356         else {
00357                 float *image0_data = image0->get_data();
00358                 if (sigma_image_data) {
00359                         memcpy(sigma_image_data, image0_data, image_size * sizeof(float));
00360 
00361                         for (size_t j = 0; j < image_size; ++j) {
00362                                 sigma_image_data[j] *= sigma_image_data[j];
00363                         }
00364                 }
00365 
00366                 image0->update();
00367                 memcpy(result_data, image0_data, image_size * sizeof(float));
00368 
00369                 for (size_t i = 1; i < image_list.size(); i++) {
00370                         EMData *image = image_list[i];
00371 
00372                         if (EMUtil::is_same_size(image, result)) {
00373                                 float *image_data = image->get_data();
00374 
00375                                 for (size_t j = 0; j < image_size; ++j) {
00376                                         result_data[j] += image_data[j];
00377                                 }
00378 
00379                                 if (sigma_image_data) {
00380                                         for (size_t j = 0; j < image_size; ++j) {
00381                                                 sigma_image_data[j] += image_data[j] * image_data[j];
00382                                         }
00383                                 }
00384 
00385                                 image->update();
00386                                 c++;
00387                         }
00388                 }
00389 
00390                 for (size_t j = 0; j < image_size; ++j) {
00391                         result_data[j] /= static_cast < float >(c);
00392                 }
00393         }
00394 
00395         if (sigma_image_data) {
00396                 for (size_t j = 0; j < image_size; ++j) {
00397                         float f1 = sigma_image_data[j] / c;
00398                         float f2 = result_data[j];
00399                         sigma_image_data[j] = sqrt(f1 - f2 * f2);
00400                 }
00401         }
00402 
00403         if (sigma_image_data) {
00404                 sigma_image->update();
00405         }
00406 
00407         result->update();
00408         return result;
00409 }
00410 #endif
00411 
00412 MinMaxAverager::MinMaxAverager()
00413         : nimg(0)
00414 {
00415         /*move max out of initializer list, since this max(0) is considered as a macro
00416          * in Visual Studio, which we defined somewhere else*/
00417         max = 0;
00418 }
00419 
00420 void MinMaxAverager::add_image(EMData * image)
00421 {
00422         if (!image) {
00423                 return;
00424         }
00425 
00426         if (nimg >= 1 && !EMUtil::is_same_size(image, result)) {
00427                 LOGERR("%sAverager can only process same-size Image",
00428                            get_name().c_str());
00429                 return;
00430         }
00431 
00432         nimg++;
00433 
00434         int nx = image->get_xsize();
00435         int ny = image->get_ysize();
00436         int nz = image->get_zsize();
00437 
00438         if (nimg == 1) {
00439                 result = image->copy();
00440                 max = params["max"];
00441                 return;
00442         }
00443 
00444         for (int z=0; z<nz; z++) {
00445                 for (int y=0; y<ny; y++) {
00446                         for (int x=0; x<nx; x++) {
00447                                 if (result->get_value_at(x,y,z)>image->get_value_at(x,y,z))
00448                                         { if (!max) result->set_value_at(x,y,z,image->get_value_at(x,y,z)); }
00449                                 else { if (max) result->set_value_at(x,y,z,image->get_value_at(x,y,z)); }
00450                         }
00451                 }
00452         }
00453 
00454 }
00455 
00456 EMData *MinMaxAverager::finish()
00457 {
00458         result->update();
00459         result->set_attr("ptcl_repr",nimg);
00460         
00461         if (result && nimg > 1) return result;
00462 
00463         return NULL;
00464 }
00465 
00466 
00467 IterationAverager::IterationAverager() : nimg(0)
00468 {
00469 
00470 }
00471 
00472 void IterationAverager::add_image( EMData * image)
00473 {
00474         if (!image) {
00475                 return;
00476         }
00477 
00478         if (nimg >= 1 && !EMUtil::is_same_size(image, result)) {
00479                 LOGERR("%sAverager can only process same-size Image",
00480                                                          get_name().c_str());
00481                 return;
00482         }
00483 
00484         nimg++;
00485 
00486         int nx = image->get_xsize();
00487         int ny = image->get_ysize();
00488         int nz = image->get_zsize();
00489         size_t image_size = (size_t)nx * ny * nz;
00490 
00491         if (nimg == 1) {
00492                 result = new EMData();
00493                 result->set_size(nx, ny, nz);
00494                 sigma_image = new EMData();
00495                 sigma_image->set_size(nx, ny, nz);
00496         }
00497 
00498         float *image_data = image->get_data();
00499         float *result_data = result->get_data();
00500         float *sigma_image_data = sigma_image->get_data();
00501 
00502         for (size_t j = 0; j < image_size; ++j) {
00503                 float f = image_data[j];
00504                 result_data[j] += f;
00505                 sigma_image_data[j] += f * f;
00506         }
00507 
00508 
00509 }
00510 
00511 EMData * IterationAverager::finish()
00512 {
00513         if (nimg < 1) {
00514                 return result;
00515         }
00516 
00517         int nx = result->get_xsize();
00518         int ny = result->get_ysize();
00519         int nz = result->get_zsize();
00520         size_t image_size = (size_t)nx * ny * nz;
00521 
00522         float *result_data = result->get_data();
00523         float *sigma_image_data = sigma_image->get_data();
00524 
00525         for (size_t j = 0; j < image_size; ++j) {
00526                 result_data[j] /= nimg;
00527                 float f1 = sigma_image_data[j] / nimg;
00528                 float f2 = result_data[j];
00529                 sigma_image_data[j] = sqrt(f1 - f2 * f2) / sqrt((float)nimg);
00530         }
00531 
00532         result->update();
00533         sigma_image->update();
00534 
00535         result->append_image("iter.hed");
00536         float sigma = sigma_image->get_attr("sigma");
00537         float *sigma_image_data2 = sigma_image->get_data();
00538         float *result_data2 = result->get_data();
00539         float *d2 = new float[nx * ny];
00540         size_t sec_size = nx * ny * sizeof(float);
00541 
00542         memcpy(d2, result_data2, sec_size);
00543         memcpy(sigma_image_data2, result_data2, sec_size);
00544 
00545         printf("Iter sigma=%f\n", sigma);
00546 
00547         for (int k = 0; k < 1000; k++) {
00548                 for (int i = 1; i < nx - 1; i++) {
00549                         for (int j = 1; j < ny - 1; j++) {
00550                                 int l = i + j * nx;
00551                                 float c1 = (d2[l - 1] + d2[l + 1] + d2[l - nx] + d2[l + nx]) / 4.0f - d2[l];
00552                                 float c2 = fabs(result_data2[l] - sigma_image_data2[l]) / sigma;
00553                                 result_data2[l] += c1 * Util::eman_erfc(c2) / 100.0f;
00554                         }
00555                 }
00556 
00557                 memcpy(d2, result_data2, sec_size);
00558         }
00559 
00560         if( d2 )
00561         {
00562                 delete[]d2;
00563                 d2 = 0;
00564         }
00565 
00566         sigma_image->update();
00567         if( sigma_image )
00568         {
00569                 delete sigma_image;
00570                 sigma_image = 0;
00571         }
00572 
00573         result->update();
00574         result->append_image("iter.hed");
00575 
00576 
00577         return result;
00578 }
00579 
00580 CtfCWautoAverager::CtfCWautoAverager()
00581         : nimg(0)
00582 {
00583 
00584 }
00585 
00586 
00587 void CtfCWautoAverager::add_image(EMData * image)
00588 {
00589         if (!image) {
00590                 return;
00591         }
00592 
00593 
00594 
00595         EMData *fft=image->do_fft();
00596 
00597         if (nimg >= 1 && !EMUtil::is_same_size(fft, result)) {
00598                 LOGERR("%s Averager can only process images of the same size", get_name().c_str());
00599                 return;
00600         }
00601 
00602         nimg++;
00603         if (nimg == 1) {
00604                 result = fft->copy_head();
00605                 result->to_zero();
00606         }
00607 
00608         Ctf *ctf = (Ctf *)image->get_attr("ctf");
00609 //string cc=ctf->to_string();
00610 //FILE *out=fopen("ctf.txt","a");
00611 //fprintf(out,"%s\n",cc.c_str());
00612 //fclose(out);
00613         float b=ctf->bfactor;
00614         ctf->bfactor=100.0;                     // FIXME - this is a temporary fixed B-factor which does a (very) little sharpening
00615 
00616 //      if (nimg==1) unlink("snr.hdf");
00617 
00618         EMData *snr = result -> copy();
00619         ctf->compute_2d_complex(snr,Ctf::CTF_SNR);
00620 //      snr->write_image("snr.hdf",-1);
00621         EMData *ctfi = result-> copy();
00622         ctf->compute_2d_complex(ctfi,Ctf::CTF_AMP);
00623 
00624         ctf->bfactor=b; // return to its original value
00625 
00626         float *outd = result->get_data();
00627         float *ind = fft->get_data();
00628         float *snrd = snr->get_data();
00629         float *ctfd = ctfi->get_data();
00630 
00631         size_t sz=snr->get_xsize()*snr->get_ysize();
00632         for (size_t i = 0; i < sz; i+=2) {
00633                 if (snrd[i]<0) snrd[i]=0;
00634                 ctfd[i]=fabs(ctfd[i]);
00635                 if (ctfd[i]<.05) ctfd[i]=0.05f;
00636 //              {
00637 //                      if (snrd[i]<=0) ctfd[i]=.05f;
00638 //                      else ctfd[i]=snrd[i]*10.0f;
00639 //              }
00640                 outd[i]+=ind[i]*snrd[i]/ctfd[i];
00641                 outd[i+1]+=ind[i+1]*snrd[i]/ctfd[i];
00642         }
00643 
00644         if (nimg==1) {
00645                 snrsum=snr->copy_head();
00646                 float *ssnrd=snrsum->get_data();
00647                 // we're only using the real component, and we need to start with 1.0
00648                 for (size_t i = 0; i < sz; i+=2) { ssnrd[i]=1.0; ssnrd[i+1]=0.0; }
00649         }
00650 //      snr->write_image("snr.hdf",-1);
00651         snr->process_inplace("math.absvalue");
00652         snrsum->add(*snr);
00653 
00654         delete ctf;
00655         delete fft;
00656         delete snr;
00657         delete ctfi;
00658 }
00659 
00660 EMData * CtfCWautoAverager::finish()
00661 {
00662 /*      EMData *tmp=result->do_ift();
00663         tmp->write_image("ctfcw.hdf",0);
00664         delete tmp;
00665 
00666         tmp=snrsum->do_ift();
00667         tmp->write_image("ctfcw.hdf",1);
00668         delete tmp;*/
00669 
00670 //snrsum->write_image("snrsum.hdf",-1);
00671         //size_t sz=result->get_xsize()*result->get_ysize();
00672         int nx=result->get_xsize();
00673         int ny=result->get_ysize();     
00674         float *snrsd=snrsum->get_data();
00675         float *outd=result->get_data();
00676 
00677         int rm=(ny-2)*(ny-2)/4;
00678         for (int j=0; j<ny; j++) {
00679                 for (int i=0; i<nx; i+=2) {
00680                         size_t ii=i+j*nx;
00681                         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; }
00682                         outd[ii]/=snrsd[ii];            // snrsd contains total SNR
00683                         outd[ii+1]/=snrsd[ii];
00684                 }
00685         }
00686 
00687         result->update();
00688         result->set_attr("ptcl_repr",nimg);
00689         result->set_attr("ctf_snr_total",snrsum->calc_radial_dist(snrsum->get_ysize()/2,0,1,false));
00690         result->set_attr("ctf_wiener_filtered",1);
00691         
00692         delete snrsum;
00693         EMData *ret=result->do_ift();
00694         delete result;
00695         result=NULL;
00696         return ret;
00697 }
00698 
00699 CtfCAutoAverager::CtfCAutoAverager()
00700         : nimg(0)
00701 {
00702 
00703 }
00704 
00705 
00706 void CtfCAutoAverager::add_image(EMData * image)
00707 {
00708         if (!image) {
00709                 return;
00710         }
00711 
00712 
00713 
00714         EMData *fft=image->do_fft();
00715 
00716         if (nimg >= 1 && !EMUtil::is_same_size(fft, result)) {
00717                 LOGERR("%s Averager can only process images of the same size", get_name().c_str());
00718                 return;
00719         }
00720 
00721         nimg++;
00722         if (nimg == 1) {
00723                 result = fft->copy_head();
00724                 result->to_zero();
00725         }
00726 
00727         Ctf *ctf = (Ctf *)image->get_attr("ctf");
00728         float b=ctf->bfactor;
00729         ctf->bfactor=0;                 // NO B-FACTOR CORRECTION !
00730 
00731         EMData *snr = result -> copy();
00732         ctf->compute_2d_complex(snr,Ctf::CTF_SNR);
00733         EMData *ctfi = result-> copy();
00734         ctf->compute_2d_complex(ctfi,Ctf::CTF_AMP);
00735 
00736         ctf->bfactor=b; // return to its original value
00737 
00738         float *outd = result->get_data();
00739         float *ind = fft->get_data();
00740         float *snrd = snr->get_data();
00741         float *ctfd = ctfi->get_data();
00742 
00743         size_t sz=snr->get_xsize()*snr->get_ysize();
00744         for (size_t i = 0; i < sz; i+=2) {
00745                 if (snrd[i]<0) snrd[i]=0;
00746                 ctfd[i]=fabs(ctfd[i]);
00747                 
00748                 // This limits the maximum possible amplification in CTF correction to 10x
00749                 if (ctfd[i]<.05)  ctfd[i]=0.05f;
00750 //              {
00751 //                      if (snrd[i]<=0) ctfd[i]=.05f;
00752 //                      else ctfd[i]=snrd[i]*10.0f;
00753 //              }
00754                 
00755                 // SNR weight and CTF correction
00756                 outd[i]+=ind[i]*snrd[i]/ctfd[i];
00757                 outd[i+1]+=ind[i+1]*snrd[i]/ctfd[i];
00758         }
00759 
00760         if (nimg==1) {
00761                 snrsum=snr->copy_head();
00762                 float *ssnrd=snrsum->get_data();
00763                 // we're only using the real component, for Wiener filter we put 1.0 in R, but for just SNR weight we use 0
00764                 for (size_t i = 0; i < sz; i+=2) { ssnrd[i]=0.0; ssnrd[i+1]=0.0; }
00765         }
00766         snr->process_inplace("math.absvalue");
00767         snrsum->add(*snr);
00768 
00769         delete ctf;
00770         delete fft;
00771         delete snr;
00772         delete ctfi;
00773 }
00774 
00775 EMData * CtfCAutoAverager::finish()
00776 {
00777 /*      EMData *tmp=result->do_ift();
00778         tmp->write_image("ctfcw.hdf",0);
00779         delete tmp;
00780 
00781         tmp=snrsum->do_ift();
00782         tmp->write_image("ctfcw.hdf",1);
00783         delete tmp;*/
00784 
00785 //      snrsum->write_image("snrsum.hdf",-1);
00786         //size_t sz=result->get_xsize()*result->get_ysize();
00787         int nx=result->get_xsize();
00788         int ny=result->get_ysize();     
00789         float *snrsd=snrsum->get_data();
00790         float *outd=result->get_data();
00791 
00792         int rm=(ny-2)*(ny-2)/4;
00793         for (int j=0; j<ny; j++) {
00794                 for (int i=0; i<nx; i+=2) {
00795                         size_t ii=i+j*nx;
00796                         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; }
00797                         // we aren't wiener filtering, but if the total SNR is too low, we don't want TOO much exaggeration of noise
00798                         if (snrsd[ii]<.05) {            
00799                                 outd[ii]*=20.0;         // 1/0.05
00800                                 outd[ii+1]*=20.0;
00801                         }
00802                         else {
00803                                 outd[ii]/=snrsd[ii];            // snrsd contains total SNR
00804                                 outd[ii+1]/=snrsd[ii];
00805                         }
00806                 }
00807         }
00808         result->update();
00809         result->set_attr("ptcl_repr",nimg);
00810         result->set_attr("ctf_snr_total",snrsum->calc_radial_dist(snrsum->get_ysize()/2,0,1,false));
00811         result->set_attr("ctf_wiener_filtered",0);
00812         
00813 /*      snrsum->write_image("snr.hdf",-1);
00814         result->write_image("avg.hdf",-1);*/
00815         
00816         delete snrsum;
00817         EMData *ret=result->do_ift();
00818         delete result;
00819         result=NULL;
00820         return ret;
00821 }
00822 
00823 #if 0
00824 EMData *IterationAverager::average(const vector < EMData * >&image_list) const
00825 {
00826         if (image_list.size() == 0) {
00827                 return 0;
00828         }
00829 
00830         EMData *image0 = image_list[0];
00831 
00832         int nx = image0->get_xsize();
00833         int ny = image0->get_ysize();
00834         int nz = image0->get_zsize();
00835         size_t image_size = (size_t)nx * ny * nz;
00836 
00837         EMData *result = new EMData();
00838         result->set_size(nx, ny, nz);
00839 
00840         EMData *sigma_image = new EMData();
00841         sigma_image->set_size(nx, ny, nz);
00842 
00843         float *image0_data = image0->get_data();
00844         float *result_data = result->get_data();
00845         float *sigma_image_data = sigma_image->get_data();
00846 
00847         memcpy(result_data, image0_data, image_size * sizeof(float));
00848         memcpy(sigma_image_data, image0_data, image_size * sizeof(float));
00849 
00850         for (size_t j = 0; j < image_size; ++j) {
00851                 sigma_image_data[j] *= sigma_image_data[j];
00852         }
00853 
00854         image0->update();
00855 
00856         int nc = 1;
00857         for (size_t i = 1; i < image_list.size(); i++) {
00858                 EMData *image = image_list[i];
00859 
00860                 if (EMUtil::is_same_size(image, result)) {
00861                         float *image_data = image->get_data();
00862 
00863                         for (size_t j = 0; j < image_size; j++) {
00864                                 result_data[j] += image_data[j];
00865                                 sigma_image_data[j] += image_data[j] * image_data[j];
00866                         }
00867 
00868                         image->update();
00869                         nc++;
00870                 }
00871         }
00872 
00873         float c = static_cast < float >(nc);
00874         for (size_t j = 0; j < image_size; ++j) {
00875                 float f1 = sigma_image_data[j] / c;
00876                 float f2 = result_data[j] / c;
00877                 sigma_image_data[j] = sqrt(f1 - f2 * f2) / sqrt(c);
00878         }
00879 
00880 
00881         for (size_t j = 0; j < image_size; ++j) {
00882                 result_data[j] /= c;
00883         }
00884 
00885         result->update();
00886         sigma_image->update();
00887 
00888         result->append_image("iter.hed");
00889 
00890         float sigma = sigma_image->get_attr("sigma");
00891         float *sigma_image_data2 = sigma_image->get_data();
00892         float *result_data2 = result->get_data();
00893         float *d2 = new float[nx * ny];
00894         size_t sec_size = nx * ny * sizeof(float);
00895 
00896         memcpy(d2, result_data2, sec_size);
00897         memcpy(sigma_image_data2, result_data2, sec_size);
00898 
00899         printf("Iter sigma=%f\n", sigma);
00900 
00901         for (int k = 0; k < 1000; k++) {
00902                 for (int i = 1; i < nx - 1; i++) {
00903                         for (int j = 1; j < ny - 1; j++) {
00904                                 int l = i + j * nx;
00905                                 float c1 = (d2[l - 1] + d2[l + 1] + d2[l - nx] + d2[l + nx]) / 4.0f - d2[l];
00906                                 float c2 = fabs(result_data2[l] - sigma_image_data2[l]) / sigma;
00907                                 result_data2[l] += c1 * Util::eman_erfc(c2) / 100.0f;
00908                         }
00909                 }
00910 
00911                 memcpy(d2, result_data2, sec_size);
00912         }
00913 
00914         if( d2 )
00915         {
00916                 delete[]d2;
00917                 d2 = 0;
00918         }
00919 
00920         sigma_image->update();
00921         if( sigma_image )
00922         {
00923                 delete sigma_image;
00924                 sigma_image = 0;
00925         }
00926 
00927         result->update();
00928         result->append_image("iter.hed");
00929 
00930         return result;
00931 }
00932 #endif
00933 
00934 
00935 CtfAverager::CtfAverager() :
00936         sf(0), curves(0), need_snr(false), outfile(0),
00937         image0_fft(0), image0_copy(0), snri(0), snrn(0),
00938         tdr(0), tdi(0), tn(0),
00939         filter(0), nimg(0), nx(0), ny(0), nz(0)
00940 {
00941 
00942 }
00943 
00944 void CtfAverager::add_image(EMData * image)
00945 {
00946         if (!image) {
00947                 return;
00948         }
00949 
00950         if (nimg >= 1 && !EMUtil::is_same_size(image, result)) {
00951                 LOGERR("%sAverager can only process same-size Image",
00952                                                          get_name().c_str());
00953                 return;
00954         }
00955 
00956         if (image->get_zsize() != 1) {
00957                 LOGERR("%sAverager: Only 2D images are currently supported",
00958                                                          get_name().c_str());
00959         }
00960 
00961         string alg_name = get_name();
00962 
00963         if (alg_name == "CtfCW" || alg_name == "CtfCWauto") {
00964                 if (image->get_ctf() != 0 && !image->has_ctff()) {
00965                         LOGERR("%sAverager: Attempted CTF Correction with no ctf parameters",
00966                                                                  get_name().c_str());
00967                 }
00968         }
00969         else {
00970                 if (image->get_ctf() != 0) {
00971                         LOGERR("%sAverager: Attempted CTF Correction with no ctf parameters",
00972                                                                  get_name().c_str());
00973                 }
00974         }
00975 
00976         nimg++;
00977 
00978 
00979         if (nimg == 1) {
00980                 image0_fft = image->do_fft();
00981 
00982                 nx = image0_fft->get_xsize();
00983                 ny = image0_fft->get_ysize();
00984                 nz = image0_fft->get_zsize();
00985 
00986                 result = new EMData();
00987                 result->set_size(nx - 2, ny, nz);
00988 
00989 
00990                 if (alg_name == "Weighting" && curves) {
00991                         if (!sf) {
00992                                 LOGWARN("CTF curve in file will contain relative, not absolute SNR!");
00993                         }
00994                         curves->set_size(Ctf::CTFOS * ny / 2, 3, 1);
00995                         curves->to_zero();
00996                 }
00997 
00998 
00999                 if (alg_name == "CtfC") {
01000                         filter = params["filter"];
01001                         if (filter == 0) {
01002                                 filter = 22.0f;
01003                         }
01004                         float apix_y = image->get_attr_dict().get("apix_y");
01005                         float ds = 1.0f / (apix_y * ny * Ctf::CTFOS);
01006                         filter = 1.0f / (filter * ds);
01007                 }
01008 
01009                 if (alg_name == "CtfCWauto") {
01010                         int nxy2 = nx * ny/2;
01011 
01012                         snri = new float[ny / 2];
01013                         snrn = new float[ny / 2];
01014                         tdr = new float[nxy2];
01015                         tdi = new float[nxy2];
01016                         tn = new float[nxy2];
01017 
01018                         for (int i = 0; i < ny / 2; i++) {
01019                                 snri[i] = 0;
01020                                 snrn[i] = 0;
01021                         }
01022 
01023                         for (int i = 0; i < nxy2; i++) {
01024                                 tdr[i] = 1;
01025                                 tdi[i] = 1;
01026                                 tn[i] = 1;
01027                         }
01028                 }
01029 
01030                 image0_copy = image0_fft->copy_head();
01031                 image0_copy->ap2ri();
01032                 image0_copy->to_zero();
01033         }
01034 
01035         Ctf::CtfType curve_type = Ctf::CTF_AMP;
01036         if (alg_name == "CtfCWauto") {
01037                 curve_type = Ctf::CTF_AMP;
01038         }
01039 
01040         float *src = image->get_data();
01041         image->ap2ri();
01042         Ctf *image_ctf = image->get_ctf();
01043         int ny2 = image->get_ysize();
01044 
01045         vector<float> ctf1 = image_ctf->compute_1d(ny2,1.0f/(image_ctf->apix*image->get_ysize()), curve_type);
01046 
01047         if (ctf1.size() == 0) {
01048                 LOGERR("Unexpected CTF correction problem");
01049         }
01050 
01051         ctf.push_back(ctf1);
01052 
01053         vector<float> ctfn1;
01054         if (sf) {
01055                 ctfn1 = image_ctf->compute_1d(ny2,1.0f/(image_ctf->apix*image->get_ysize()), Ctf::CTF_SNR, sf);
01056         }
01057         else {
01058                 ctfn1 = image_ctf->compute_1d(ny2,1.0f/(image_ctf->apix*image->get_ysize()), Ctf::CTF_SNR);
01059         }
01060 
01061         ctfn.push_back(ctfn1);
01062 
01063         if (alg_name == "CtfCWauto") {
01064                 int j = 0;
01065                 for (int y = 0; y < ny; y++) {
01066                         for (int x = 0; x < nx / 2; x++, j += 2) {
01067 #ifdef  _WIN32
01068                                 float r = (float)_hypot((float)x, (float)(y - ny / 2.0f));
01069 #else
01070                                 float r = (float)hypot((float)x, (float)(y - ny / 2.0f));
01071 #endif  //_WIN32
01072                                 int l = static_cast < int >(Util::fast_floor(r));
01073 
01074                                 if (l >= 0 && l < ny / 2) {
01075                                         int k = y*nx/2 + x;
01076                                         tdr[k] *= src[j];
01077                                         tdi[k] *= src[j + 1];
01078 #ifdef  _WIN32
01079                                         tn[k] *= (float)_hypot(src[j], src[j + 1]);
01080 #else
01081                                         tn[k] *= (float)hypot(src[j], src[j + 1]);
01082 #endif  //_WIN32
01083                                 }
01084                         }
01085                 }
01086         }
01087 
01088 
01089         float *tmp_data = image0_copy->get_data();
01090 
01091         int j = 0;
01092         for (int y = 0; y < ny; y++) {
01093                 for (int x = 0; x < nx / 2; x++, j += 2) {
01094                         float r = Ctf::CTFOS * sqrt(x * x + (y - ny / 2.0f) * (y - ny / 2.0f));
01095                         int l = static_cast < int >(Util::fast_floor(r));
01096                         r -= l;
01097 
01098                         float f = 0;
01099                         if (l <= Ctf::CTFOS * ny / 2 - 2) {
01100                                 f = (ctf1[l] * (1 - r) + ctf1[l + 1] * r);
01101                         }
01102                         tmp_data[j] += src[j] * f;
01103                         tmp_data[j + 1] += src[j + 1] * f;
01104                 }
01105         }
01106 
01107         EMData *image_fft = image->do_fft();
01108         image_fft->update();
01109         if(image_ctf) {delete image_ctf; image_ctf=0;}
01110 }
01111 
01112 EMData * CtfAverager::finish()
01113 {
01114         int j = 0;
01115         for (int y = 0; y < ny; y++) {
01116                 for (int x = 0; x < nx / 2; x++, j += 2) {
01117 #ifdef  _WIN32
01118                         float r = (float) _hypot(x, y - ny / 2.0f);
01119 #else
01120                         float r = (float) hypot(x, y - ny / 2.0f);
01121 #endif
01122                         int l = static_cast < int >(Util::fast_floor(r));
01123                         if (l >= 0 && l < ny / 2) {
01124                                 int k = y*nx/2 + x;
01125                                 snri[l] += (tdr[k] + tdi[k]/tn[k]);
01126                                 snrn[l] += 1;
01127                         }
01128                 }
01129         }
01130 
01131         for (int i = 0; i < ny / 2; i++) {
01132                 snri[i] *= nimg / snrn[i];
01133         }
01134 
01135         if(strcmp(outfile, "") != 0) {
01136                 Util::save_data(0, 1, snri, ny / 2, outfile);
01137         }
01138 
01139 
01140         float *cd = 0;
01141         if (curves) {
01142                 cd = curves->get_data();
01143         }
01144 
01145         for (int i = 0; i < Ctf::CTFOS * ny / 2; i++) {
01146                 float ctf0 = 0;
01147                 for (int j = 0; j < nimg; j++) {
01148                         ctf0 += ctfn[j][i];
01149                         if (ctf[j][i] == 0) {
01150                                 ctf[j][i] = 1.0e-12f;
01151                         }
01152 
01153                         if (curves) {
01154                                 cd[i] += ctf[j][i] * ctfn[j][i];
01155                                 cd[i + Ctf::CTFOS * ny / 2] += ctfn[j][i];
01156                                 cd[i + 2 * Ctf::CTFOS * ny / 2] += ctfn[j][i];
01157                         }
01158                 }
01159 
01160                 string alg_name = get_name();
01161 
01162                 if (alg_name == "CtfCW" && need_snr) {
01163                         snr[i] = ctf0;
01164                 }
01165 
01166                 float ctf1 = ctf0;
01167                 if (alg_name == "CtfCWauto") {
01168                         ctf1 = snri[i / Ctf::CTFOS];
01169                 }
01170 
01171                 if (ctf1 <= 0.0001f) {
01172                         ctf1 = 0.1f;
01173                 }
01174 
01175                 if (alg_name == "CtfC") {
01176                         for (int j = 0; j < nimg; j++) {
01177                                 ctf[j][i] = exp(-i * i / (filter * filter)) * ctfn[j][i] / (fabs(ctf[j][i]) * ctf1);
01178                         }
01179                 }
01180                 else if (alg_name == "Weighting") {
01181                         for (int j = 0; j < nimg; j++) {
01182                                 ctf[j][i] = ctfn[j][i] / ctf1;
01183                         }
01184                 }
01185                 else if (alg_name == "CtfCW") {
01186                         for (int j = 0; j < nimg; j++) {
01187                                 ctf[j][i] = (ctf1 / (ctf1 + 1)) * ctfn[j][i] / (ctf[j][i] * ctf1);
01188                         }
01189                 }
01190                 else if (alg_name == "CtfCWauto") {
01191                         for (int j = 0; j < nimg; j++) {
01192                                 ctf[j][i] = ctf1 * ctfn[j][i] / (fabs(ctf[j][i]) * ctf0);
01193                         }
01194                 }
01195         }
01196 
01197 
01198         if (curves) {
01199                 for (int i = 0; i < Ctf::CTFOS * ny / 2; i++) {
01200                         cd[i] /= cd[i + Ctf::CTFOS * ny / 2];
01201                 }
01202                 curves->update();
01203         }
01204 
01205         image0_copy->update();
01206 
01207         float *result_data = result->get_data();
01208         EMData *tmp_ift = image0_copy->do_ift();
01209         float *tmp_ift_data = tmp_ift->get_data();
01210         memcpy(result_data, tmp_ift_data, (nx - 2) * ny * sizeof(float));
01211 
01212         tmp_ift->update();
01213         result->update();
01214         result->set_attr("ptcl_repr",nimg);
01215 
01216         if( image0_copy )
01217         {
01218                 delete image0_copy;
01219                 image0_copy = 0;
01220         }
01221 
01222         if (snri) {
01223                 delete[]snri;
01224                 snri = 0;
01225         }
01226 
01227         if (snrn) {
01228                 delete[]snrn;
01229                 snrn = 0;
01230         }
01231 
01232         if( snri )
01233         {
01234                 delete [] snri;
01235                 snri = 0;
01236         }
01237         if( snrn )
01238         {
01239                 delete [] snrn;
01240                 snrn = 0;
01241         }
01242         if( tdr )
01243         {
01244                 delete [] tdr;
01245                 tdr = 0;
01246         }
01247         if( tdi )
01248         {
01249                 delete [] tdi;
01250                 tdi = 0;
01251         }
01252         if( tn )
01253         {
01254                 delete [] tn;
01255                 tn = 0;
01256         }
01257 
01258         return result;
01259 }
01260 
01261 #if 0
01262 EMData *CtfAverager::average(const vector < EMData * >&image_list) const
01263 {
01264         if (image_list.size() == 0) {
01265                 return 0;
01266         }
01267 
01268         EMData *image0 = image_list[0];
01269         if (image0->get_zsize() != 1) {
01270                 LOGERR("Only 2D images are currently supported");
01271                 return 0;
01272         }
01273 
01274         string alg_name = get_name();
01275 
01276         if (alg_name == "CtfCW" || alg_name == "CtfCWauto") {
01277                 if (image0->get_ctf() != 0 && !image0->has_ctff()) {
01278                         LOGERR("Attempted CTF Correction with no ctf parameters");
01279                         return 0;
01280                 }
01281         }
01282         else {
01283                 if (image0->get_ctf() != 0) {
01284                         LOGERR("Attempted CTF Correction with no ctf parameters");
01285                         return 0;
01286                 }
01287         }
01288 
01289         size_t num_images = image_list.size();
01290         vector < float >*ctf = new vector < float >[num_images];
01291         vector < float >*ctfn = new vector < float >[num_images];
01292         float **src = new float *[num_images];
01293 
01294         Ctf::CtfType curve_type = Ctf::CTF_ABS_AMP;
01295         if (alg_name == "CtfCWauto") {
01296                 curve_type = Ctf::CTF_AMP;
01297         }
01298 
01299         for (size_t i = 0; i < num_images; i++) {
01300                 EMData *image = image_list[i]->do_fft();
01301                 image->ap2ri();
01302                 src[i] = image->get_data();
01303                 Ctf *image_ctf = image->get_ctf();
01304                 int ny = image->get_ysize();
01305                 ctf[i] = image_ctf->compute_1d(ny, curve_type);
01306 
01307                 if (ctf[i].size() == 0) {
01308                         LOGERR("Unexpected CTF correction problem");
01309                         return 0;
01310                 }
01311 
01312                 if (sf) {
01313                         ctfn[i] = image_ctf->compute_1d(ny, Ctf::CTF_ABS_SNR, sf);
01314                 }
01315                 else {          result->set_attr("ptcl_repr",nimg);
01316                         ctfn[i] = image_ctf->compute_1d(ny, Ctf::CTF_RELATIVE_SNR);
01317                 }
01318 
01319                 if(image_ctf) {delete image_ctf; image_ctf=0;}
01320         }
01321 
01322         EMData *image0_fft = image0->do_fft();
01323 
01324         int nx = image0_fft->get_xsize();
01325         int ny = image0_fft->get_ysize();
01326         int nz = image0_fft->get_zsize();
01327 
01328         EMData *result = new EMData();
01329         result->set_size(nx - 2, ny, nz);
01330 
01331         float *cd = 0;
01332         if (alg_name == "Weighting" && curves) {
01333                 if (!sf) {
01334                         LOGWARN("CTF curve in file will contain relative, not absolute SNR!");
01335                 }
01336                 curves->set_size(Ctf::CTFOS * ny / 2, 3, 1);
01337                 curves->to_zero();
01338                 cd = curves->get_data();
01339         }
01340 
01341         float filter = 0;
01342         if (alg_name == "CtfC") {
01343                 filter = params["filter"];
01344                 if (filter == 0) {
01345                         filter = 22.0f;
01346                 }
01347                 float apix_y = image0->get_attr_dict().get("apix_y");
01348                 float ds = 1.0f / (apix_y * ny * Ctf::CTFOS);
01349                 filter = 1.0f / (filter * ds);
01350         }
01351 
01352         float *snri = 0;
01353         float *snrn = 0;
01354 
01355         if (alg_name == "CtfCWauto") {
01356                 snri = new float[ny / 2];
01357                 snrn = new float[ny / 2];
01358 
01359                 for (int i = 0; i < ny / 2; i++) {
01360                         snri[i] = 0;
01361                         snrn[i] = 0;
01362                 }
01363 
01364                 int j = 0;
01365                 for (int y = 0; y < ny; y++) {
01366                         for (int x = 0; x < nx / 2; x++, j += 2) {
01367                                 float r = hypot(x, y - ny / 2.0f);
01368                                 int l = static_cast < int >(Util::fast_floor(r));
01369 
01370                                 if (l >= 0 && l < ny / 2) {
01371                                         float tdr = 1;
01372                                         float tdi = 1;
01373                                         float tn = 1;
01374 
01375                                         for (size_t i = 0; i < num_images; i++) {
01376                                                 tdr *= src[i][j];
01377                                                 tdi *= src[i][j + 1];
01378                                                 tn *= hypot(src[i][j], src[i][j + 1]);
01379                                         }
01380 
01381                                         tdr += tdi / tn;
01382                                         snri[l] += tdr;
01383                                         snrn[l] += 1;
01384                                 }
01385                         }
01386                 }
01387 
01388                 for (int i = 0; i < ny / 2; i++) {
01389                         snri[i] *= num_images / snrn[i];
01390                 }
01391                 if (outfile != "") {
01392                         Util::save_data(0, 1, snri, ny / 2, outfile);
01393                 }
01394         }
01395 
01396         for (int i = 0; i < Ctf::CTFOS * ny / 2; i++) {
01397                 float ctf0 = 0;
01398                 for (size_t j = 0; j < num_images; j++) {
01399                         ctf0 += ctfn[j][i];
01400                         if (ctf[j][i] == 0) {
01401                                 ctf[j][i] = 1.0e-12;
01402                         }
01403 
01404                         if (curves) {
01405                                 cd[i] += ctf[j][i] * ctfn[j][i];
01406                                 cd[i + Ctf::CTFOS * ny / 2] += ctfn[j][i];
01407                                 cd[i + 2 * Ctf::CTFOS * ny / 2] += ctfn[j][i];
01408                         }
01409                 }
01410 
01411                 if (alg_name == "CtfCW" && need_snr) {
01412                         snr[i] = ctf0;
01413                 }
01414 
01415                 float ctf1 = ctf0;
01416                 if (alg_name == "CtfCWauto") {
01417                         ctf1 = snri[i / Ctf::CTFOS];
01418                 }
01419 
01420                 if (ctf1 <= 0.0001f) {
01421                         ctf1 = 0.1f;
01422                 }
01423 
01424                 if (alg_name == "CtfC") {
01425                         for (size_t j = 0; j < num_images; j++) {
01426                                 ctf[j][i] = exp(-i * i / (filter * filter)) * ctfn[j][i] / (fabs(ctf[j][i]) * ctf1);
01427                         }
01428                 }
01429                 else if (alg_name == "Weighting") {
01430                         for (size_t j = 0; j < num_images; j++) {
01431                                 ctf[j][i] = ctfn[j][i] / ctf1;
01432                         }
01433                 }
01434                 else if (alg_name == "CtfCW") {
01435                         for (size_t j = 0; j < num_images; j++) {
01436                                 ctf[j][i] = (ctf1 / (ctf1 + 1)) * ctfn[j][i] / (ctf[j][i] * ctf1);
01437                         }
01438                 }
01439                 else if (alg_name == "CtfCWauto") {
01440                         for (size_t j = 0; j < num_images; j++) {
01441                                 ctf[j][i] = ctf1 * ctfn[j][i] / (fabs(ctf[j][i]) * ctf0);
01442                         }
01443                 }
01444         }
01445 
01446 
01447         if (curves) {
01448                 for (int i = 0; i < Ctf::CTFOS * ny / 2; i++) {
01449                         cd[i] /= cd[i + Ctf::CTFOS * ny / 2];
01450                 }
01451                 curves->update();
01452         }
01453 
01454         EMData *image0_copy = image0_fft->copy_head();
01455         image0_copy->ap2ri();
01456 
01457         float *tmp_data = image0_copy->get_data();
01458 
01459         int j = 0;
01460         for (int y = 0; y < ny; y++) {
01461                 for (int x = 0; x < nx / 2; x++, j += 2) {
01462                         float r = Ctf::CTFOS * sqrt(x * x + (y - ny / 2.0f) * (y - ny / 2.0f));
01463                         int l = static_cast < int >(Util::fast_floor(r));
01464                         r -= l;
01465 
01466                         tmp_data[j] = 0;
01467                         tmp_data[j + 1] = 0;
01468 
01469                         for (size_t i = 0; i < num_images; i++) {
01470                                 float f = 0;
01471                                 if (l <= Ctf::CTFOS * ny / 2 - 2) {
01472                                         f = (ctf[i][l] * (1 - r) + ctf[i][l + 1] * r);
01473                                 }
01474                                 tmp_data[j] += src[i][j] * f;
01475                                 tmp_data[j + 1] += src[i][j + 1] * f;
01476                         }
01477                 }
01478         }
01479 
01480         image0_copy->update();
01481 
01482         float *result_data = result->get_data();
01483         EMData *tmp_ift = image0_copy->do_ift();
01484         float *tmp_ift_data = tmp_ift->get_data();
01485         memcpy(result_data, tmp_ift_data, (nx - 2) * ny * sizeof(float));
01486 
01487         tmp_ift->update();
01488 
01489         if( image0_copy )
01490         {
01491                 delete image0_copy;
01492                 image0_copy = 0;
01493         }
01494 
01495         for (size_t i = 0; i < num_images; i++) {
01496                 EMData *img = image_list[i]->do_fft();
01497                 img->update();
01498         }
01499 
01500         if( src )
01501         {
01502                 delete[]src;
01503                 src = 0;
01504         }
01505 
01506         if( ctf )
01507         {
01508                 delete[]ctf;
01509                 ctf = 0;
01510         }
01511 
01512         if( ctfn )
01513         {
01514                 delete[]ctfn;
01515                 ctfn = 0;
01516         }
01517 
01518         if (snri) {
01519                 delete[]snri;
01520                 snri = 0;
01521         }
01522 
01523         if (snrn) {
01524                 delete[]snrn;
01525                 snrn = 0;
01526         }
01527 
01528         result->update();
01529         return result;
01530 }
01531 #endif
01532 
01533 
01534 void EMAN::dump_averagers()
01535 {
01536         dump_factory < Averager > ();
01537 }
01538 
01539 map<string, vector<string> > EMAN::dump_averagers_list()
01540 {
01541         return dump_factory_list < Averager > ();
01542 }

Generated on Thu Nov 17 12:43:44 2011 for EMAN2 by  doxygen 1.3.9.1