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

Generated on Tue May 25 17:13:31 2010 for EMAN2 by  doxygen 1.4.7