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