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

Generated on Fri Aug 10 16:30:37 2012 for EMAN2 by  doxygen 1.4.7