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