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 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
00060
00061
00062 force_add<CtfCWautoAverager>();
00063
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
00311
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;
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;
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
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
00546
00547
00548
00549
00550
00551
00552
00553
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];
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 }