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 "ctf.h"
00037 #include "emdata.h"
00038 #include "xydata.h"
00039 #include "emassert.h"
00040
00041 using namespace EMAN;
00042
00043 EMAN1Ctf::EMAN1Ctf()
00044 {
00045 defocus = 0;
00046 bfactor = 0;
00047 amplitude = 0;
00048 ampcont = 0;
00049 noise1 = 0;
00050 noise2 = 0;
00051 noise3 = 0;
00052 noise4 = 0;
00053 voltage = 0;
00054 cs = 0;
00055 apix = 0;
00056 }
00057
00058
00059 EMAN1Ctf::~EMAN1Ctf()
00060 {
00061 }
00062
00063
00064 int EMAN1Ctf::from_string(const string & ctf)
00065 {
00066 Assert(ctf != "");
00067 char type;
00068 int pos;
00069 int i = sscanf(ctf.c_str(), "%c%f %f %f %f %f %f %f %f %f %f %f%n",
00070 &type,&defocus, &bfactor, &litude, &cont, &noise1,
00071 &noise2, &noise3, &noise4, &voltage, &cs, &apix, &pos);
00072 if (pos==-1) {
00073 const char *s=ctf.c_str();
00074 throw InvalidValueException(s," Invalid CTF string");
00075 }
00076
00077 return 0;
00078 }
00079
00080 void EMAN1Ctf::from_dict(const Dict & dict)
00081 {
00082 defocus = dict["defocus"];
00083 bfactor = dict["bfactor"];
00084 amplitude = dict["amplitude"];
00085 ampcont = dict["ampcont"];
00086 noise1 = dict["noise1"];
00087 noise2 = dict["noise2"];
00088 noise3 = dict["noise3"];
00089 noise4 = dict["noise4"];
00090 voltage = dict["voltage"];
00091 cs = dict["cs"];
00092 apix = dict["apix"];
00093 }
00094
00095 Dict EMAN1Ctf::to_dict() const
00096 {
00097 Dict dict;
00098 dict["defocus"] = defocus;
00099 dict["bfactor"] = bfactor;
00100 dict["amplitude"] = amplitude;
00101 dict["ampcont"] = ampcont;
00102 dict["noise1"] = noise1;
00103 dict["noise2"] = noise2;
00104 dict["noise3"] = noise3;
00105 dict["noise4"] = noise4;
00106 dict["voltage"] = voltage;
00107 dict["cs"] = cs;
00108 dict["apix"] = apix;
00109
00110 return dict;
00111 }
00112
00113 void EMAN1Ctf::from_vector(const vector<float>& vctf)
00114 {
00115 defocus = vctf[0];
00116 bfactor = vctf[1];
00117 amplitude = vctf[2];
00118 ampcont = vctf[3];
00119 noise1 = vctf[4];
00120 noise2 = vctf[5];
00121 noise3 = vctf[6];
00122 noise4 = vctf[7];
00123 voltage = vctf[8];
00124 cs = vctf[9];
00125 apix = vctf[10];
00126 }
00127
00128 vector<float> EMAN1Ctf::to_vector() const
00129 {
00130 vector<float> vctf;
00131
00132 vctf.push_back(defocus);
00133 vctf.push_back(bfactor);
00134 vctf.push_back(amplitude);
00135 vctf.push_back(ampcont);
00136 vctf.push_back(noise1);
00137 vctf.push_back(noise2);
00138 vctf.push_back(noise3);
00139 vctf.push_back(noise4);
00140 vctf.push_back(voltage);
00141 vctf.push_back(cs);
00142 vctf.push_back(apix);
00143
00144 return vctf;
00145 }
00146
00147
00148 string EMAN1Ctf::to_string() const
00149 {
00150 char ctf[1024];
00151 sprintf(ctf, "O%1.3g %1.3g %1.3g %1.3g %1.3g %1.3g %1.3g %1.3g %1.3g %1.3g %1.3g",
00152 defocus, bfactor, amplitude, ampcont, noise1, noise2, noise3, noise4, voltage, cs,
00153 apix);
00154
00155 return string(ctf);
00156 }
00157
00158 void EMAN1Ctf::copy_from(const Ctf * new_ctf)
00159 {
00160 if (new_ctf) {
00161 const EMAN1Ctf *c = static_cast<const EMAN1Ctf *>(new_ctf);
00162 defocus = c->defocus;
00163 bfactor = c->bfactor;
00164 amplitude = c->amplitude;
00165 ampcont = c->ampcont;
00166 noise1 = c->noise1;
00167 noise2 = c->noise2;
00168 noise3 = c->noise3;
00169 noise4 = c->noise4;
00170 voltage = c->voltage;
00171 cs = c->cs;
00172 apix = c->apix;
00173 }
00174 }
00175
00176
00177 vector < float >EMAN1Ctf::compute_1d(int size, float ds, CtfType type, XYData * sf)
00178 {
00179 Assert(size > 0);
00180
00181 float tmp_f1 = CTFOS * sqrt((float) 2) * size / 2;
00182 int np = (int) ceil(tmp_f1) + 2;
00183 vector < float >r;
00184
00185 r.resize(np);
00186
00187
00188 float s = 0;
00189 float g1 = calc_g1();
00190 float g2 = calc_g2();
00191 float amp1 = calc_amp1();
00192
00193 switch (type) {
00194 case CTF_AMP:
00195 for (int i = 0; i < np; i++) {
00196 float gamma = calc_gamma(g1, g2, s);
00197 r[i] = calc_ctf1(amp1, gamma, s);
00198 s += ds;
00199 }
00200 break;
00201
00202 case CTF_SIGN:
00203 for (int i = 0; i < np; i++) {
00204 float gamma = calc_gamma(g1, g2, s);
00205 r[i] = calc_ctf1(amp1, gamma, s)>0?1.0f:-1.0f;
00206 s += ds;
00207 }
00208 break;
00209
00210 case CTF_BACKGROUND:
00211 for (int i = 0; i < np; i++) {
00212 r[i] = calc_noise(s);
00213 s += ds;
00214 }
00215 break;
00216
00217 case CTF_SNR:
00218 case CTF_SNR_SMOOTH:
00219
00220
00221
00222
00223
00224 for (int i = 0; i < np; i++) {
00225 float gamma = calc_gamma(g1, g2, s);
00226 r[i] = calc_snr(amp1, gamma, s);
00227 if (s && sf) {
00228 r[i] *= sf->get_yatx(s);
00229 }
00230 s += ds;
00231 }
00232
00233 break;
00234
00235 case CTF_WIENER_FILTER:
00236 if (!sf) {
00237 LOGERR("CTF computation error, no SF found\n");
00238 return r;
00239 }
00240
00241 for (int i = 0; i < np; i++) {
00242 float gamma = calc_gamma(g1, g2, s);
00243 r[i] = calc_snr(amp1, gamma, s);
00244 if (s && sf) {
00245 r[i] *= sf->get_yatx(s);
00246 }
00247
00248 r[i] = 1.0f / (1.0f + 1.0f / r[i]);
00249 s += ds;
00250 }
00251 break;
00252
00253 case CTF_TOTAL:
00254 if (!sf) {
00255 LOGERR("CTF computation error, no SF found\n");
00256 return r;
00257 }
00258
00259 for (int i = 0; i < np; i++) {
00260 float gamma = calc_gamma(g1, g2, s);
00261 if (sf) {
00262 r[i] = calc_ctf1(amp1, gamma, s);
00263 r[i] = r[i] * r[i] * sf->get_yatx(s) + calc_noise(s);
00264 }
00265 else {
00266 r[i] = calc_ctf1(amp1, gamma, s);
00267 r[i] = r[i] * r[i] + calc_noise(s);
00268 }
00269 s += ds;
00270 }
00271 break;
00272 default:
00273 break;
00274 }
00275
00276 return r;
00277 }
00278
00279
00280 void EMAN1Ctf::compute_2d_real(EMData *, CtfType, XYData *)
00281 {
00282
00283
00284 }
00285
00286
00287
00288 void EMAN1Ctf::compute_2d_complex(EMData * image, CtfType type, XYData * sf)
00289 {
00290 if (!image) {
00291 LOGERR("image is null. cannot computer 2D complex CTF");
00292 return;
00293 }
00294
00295 if (image->is_complex() == false) {
00296 LOGERR("compute_2d_complex can only work on complex images");
00297 return;
00298 }
00299
00300 int nx = image->get_xsize();
00301 int ny = image->get_ysize();
00302
00303 if (nx != ny + 2) {
00304 LOGERR("compute_2d_complex only works on (nx, nx-2) images");
00305 return;
00306 }
00307
00308 float ds = 1.0f / (apix * ny);
00309 image->to_one();
00310
00311 float *d = image->get_data();
00312 float g1 = calc_g1();
00313 float g2 = calc_g2();
00314
00315 if (type == CTF_BACKGROUND) {
00316 for (int y = 0; y < ny; y++) {
00317 int ynx = y * nx;
00318
00319 for (int x = 0; x < nx / 2; x++) {
00320 #ifdef _WIN32
00321 float s = (float) _hypot(x, y - ny / 2.0f) * ds;
00322 #else
00323 float s = (float) hypot(x, y - ny / 2.0f) * ds;
00324 #endif
00325 d[x * 2 + ynx] = calc_noise(s);
00326 d[x * 2 + ynx + 1] = 0;
00327 }
00328 }
00329 }
00330 else if (type == CTF_AMP) {
00331 for (int y = 0; y < ny; y++) {
00332 int ynx = y * nx;
00333
00334 for (int x = 0; x < nx / 2; x++) {
00335 #ifdef _WIN32
00336 float s = (float)_hypot((float) x, (float) y - ny / 2) * ds;
00337 #else
00338 float s = (float)hypot((float) x, (float) y - ny / 2) * ds;
00339 #endif //_WIN32
00340 float gamma = calc_gamma(g1, g2, s);
00341 float v = fabs(calc_amplitude(gamma));
00342 d[x * 2 + ynx] = v;
00343 d[x * 2 + ynx + 1] = 0;
00344 }
00345 }
00346 }
00347 else if (type == CTF_SIGN) {
00348 for (int y = 0; y < ny; y++) {
00349 int ynx = y * nx;
00350 for (int x = 0; x < nx / 2; x++) {
00351 #ifdef _WIN32
00352 float s = (float)_hypot(x, y - ny / 2.0f) * ds;
00353 #else
00354 float s = (float)hypot(x, y - ny / 2.0f) * ds;
00355 #endif
00356 float gamma = calc_gamma(g1, g2, s);
00357 float v = calc_amplitude(gamma);
00358 d[x * 2 + ynx] = v > 0 ? 1.0f : -1.0f;
00359 d[x * 2 + ynx + 1] = 0;
00360 }
00361 }
00362
00363 }
00364 else if (type == CTF_SNR || type == CTF_SNR_SMOOTH) {
00365 float amp1 = calc_amp1();
00366
00367 for (int y = 0; y < ny; y++) {
00368 int ynx = y * nx;
00369
00370 for (int x = 0; x < nx / 2; x++) {
00371
00372 #ifdef _WIN32
00373 float s = (float)_hypot(x, y - ny / 2.0f) * ds;
00374 #else
00375 float s = (float)hypot(x, y - ny / 2.0f) * ds;
00376 #endif
00377 float gamma = calc_gamma(g1, g2, s);
00378 float f = calc_ctf1(amp1, gamma, s);
00379 float noise = calc_noise(s);
00380 f = f * f / noise;
00381
00382 if (s && sf) {
00383 f *= sf->get_yatx(s);
00384 }
00385 d[x * 2 + ynx] *= f;
00386 d[x * 2 + ynx + 1] = 0;
00387 }
00388 }
00389 }
00390 else if (type == CTF_WIENER_FILTER) {
00391 float amp1 = calc_amp1();
00392
00393 for (int y = 0; y < ny; y++) {
00394 int ynx = y * nx;
00395
00396 for (int x = 0; x < nx / 2; x++) {
00397
00398 #ifdef _WIN32
00399 float s = (float)_hypot(x, y - ny / 2.0f) * ds;
00400 #else
00401 float s = (float)hypot(x, y - ny / 2.0f) * ds;
00402 #endif
00403 float gamma = calc_gamma(g1, g2, s);
00404 float f = calc_ctf1(amp1, gamma, s);
00405 float noise = calc_noise(s);
00406 f = f * f / noise;
00407
00408 if (s) {
00409 f *= sf->get_yatx(s);
00410 }
00411 f = 1.0f / (1.0f + 1.0f / f);
00412 d[x * 2 + ynx] *= f;
00413 d[x * 2 + ynx + 1] = 0;
00414 }
00415 }
00416 }
00417 else if (type == CTF_TOTAL) {
00418 float amp1 = calc_amp1();
00419
00420 for (int y = 0; y < ny; y++) {
00421 int ynx = y * nx;
00422
00423 for (int x = 0; x < nx / 2; x++) {
00424
00425 #ifdef _WIN32
00426 float s = (float)_hypot(x, y - ny / 2.0f) * ds;
00427 #else
00428 float s = (float)hypot(x, y - ny / 2.0f) * ds;
00429 #endif
00430 float gamma = calc_gamma(g1, g2, s);
00431 float f = calc_ctf1(amp1, gamma, s);
00432 float noise = calc_noise(s);
00433 f = f * f;
00434
00435 if (sf && s) {
00436 f *= sf->get_yatx(s);
00437 }
00438 f+=noise;
00439
00440 d[x * 2 + ynx] *= f;
00441 d[x * 2 + ynx + 1] = 0;
00442 }
00443 }
00444 }
00445
00446 image->update();
00447 }
00448
00449
00450
00451 bool EMAN1Ctf::equal(const Ctf * ctf1) const
00452 {
00453 if (ctf1) {
00454 const EMAN1Ctf *c = static_cast<const EMAN1Ctf *>(ctf1);
00455 if (defocus == c->defocus &&
00456 bfactor == c->bfactor &&
00457 amplitude == c->amplitude &&
00458 ampcont == c->ampcont &&
00459 noise1 == c->noise1 &&
00460 noise2 == c->noise2 &&
00461 noise3 == c->noise3 &&
00462 noise4 == c->noise4 && voltage == c->voltage && cs == c->cs && apix == c->apix) {
00463 return true;
00464 }
00465 }
00466 return false;
00467 }
00468
00469
00470
00471
00472
00473 EMAN2Ctf::EMAN2Ctf()
00474 {
00475 defocus = 0;
00476 dfdiff = 0;
00477 dfang = 0;
00478 bfactor = 0;
00479 ampcont = 0;
00480 voltage = 0;
00481 cs = 0;
00482 apix = 1.0;
00483 dsbg=-1;
00484 background.clear();
00485 snr.clear();
00486 }
00487
00488
00489 EMAN2Ctf::~EMAN2Ctf()
00490 {
00491 }
00492
00493
00494 int EMAN2Ctf::from_string(const string & ctf)
00495 {
00496 Assert(ctf != "");
00497 char type=' ';
00498 int pos=-1,i,j;
00499 int bglen=0,snrlen=0;
00500 float v;
00501 const char *s=ctf.c_str();
00502
00503 sscanf(s, "%c%f %f %f %f %f %f %f %f %f %d%n",
00504 &type,&defocus, &dfdiff,&dfang,&bfactor,&cont,&voltage, &cs, &apix,&dsbg,&bglen,&pos);
00505 if (type!='E') throw InvalidValueException(type,"Trying to initialize Ctf object with bad string");
00506 if (pos==-1) throw InvalidValueException(s," Invalid CTF string");
00507
00508
00509 background.resize(bglen);
00510 for (i=0; i<bglen; i++) {
00511 if (sscanf(s+pos,",%f%n",&v,&j)<1) return(1);
00512 background[i]=v;
00513 pos+=j;
00514 }
00515
00516 sscanf(s+pos," %d%n",&snrlen,&j);
00517 pos+=j;
00518 snr.resize(snrlen);
00519 for (i=0; i<snrlen; i++) {
00520 if (sscanf(s+pos,",%f%n",&v,&j)<1) return(1);
00521 snr[i]=v;
00522 pos+=j;
00523 }
00524
00525 return 0;
00526
00527 }
00528
00529 string EMAN2Ctf::to_string() const
00530 {
00531 char ctf[256];
00532 sprintf(ctf, "E%1.4g %1.4g %1.4g %1.4g %1.4g %1.4g %1.4g %1.4g %1.4g %d",
00533 defocus, dfdiff, dfang, bfactor, ampcont, voltage, cs, apix, dsbg,(int)background.size());
00534
00535 string ret=ctf;
00536 for (int i=0; i<(int)background.size(); i++) {
00537 sprintf(ctf,",%1.3g",background[i]);
00538 ret+=ctf;
00539 }
00540
00541 sprintf(ctf, " %d",(int)snr.size());
00542 ret+=ctf;
00543 for (int i=0; i<(int)snr.size(); i++) {
00544 sprintf(ctf,",%1.3g",snr[i]);
00545 ret+=ctf;
00546 }
00547
00548
00549 return ret;
00550 }
00551
00552 void EMAN2Ctf::from_dict(const Dict & dict)
00553 {
00554 defocus = (float)dict["defocus"];
00555 dfdiff = (float)dict["dfdiff"];
00556 dfang = (float)dict["dfang"];
00557 bfactor = (float)dict["bfactor"];
00558 ampcont = (float)dict["ampcont"];
00559 voltage = (float)dict["voltage"];
00560 cs = (float)dict["cs"];
00561 apix = (float)dict["apix"];
00562 dsbg = (float)dict["dsbg"];
00563 background = dict["background"];
00564 snr = dict["snr"];
00565 }
00566
00567 Dict EMAN2Ctf::to_dict() const
00568 {
00569 Dict dict;
00570 dict["defocus"] = defocus;
00571 dict["dfdiff"] = dfdiff;
00572 dict["dfang"] = dfang;
00573 dict["bfactor"] = bfactor;
00574 dict["ampcont"] = ampcont;
00575 dict["voltage"] = voltage;
00576 dict["cs"] = cs;
00577 dict["apix"] = apix;
00578 dict["dsbg"] = dsbg;
00579 dict["background"] = background;
00580 dict["snr"] = snr;
00581
00582 return dict;
00583 }
00584
00585 void EMAN2Ctf::from_vector(const vector<float>& vctf)
00586 {
00587 int i;
00588 defocus = vctf[0];
00589 dfdiff = vctf[1];
00590 dfang = vctf[2];
00591 bfactor = vctf[3];
00592 ampcont = vctf[4];
00593 voltage = vctf[5];
00594 cs = vctf[6];
00595 apix = vctf[7];
00596 dsbg = vctf[8];
00597 background.resize((int)vctf[9]);
00598 for (i=0; i<(int)vctf[9]; i++) background[i]=vctf[i+10];
00599 snr.resize((int)vctf[i+10]);
00600 for (int j=0; j<(int)vctf[i+10]; j++) snr[j]=vctf[i+j+11];
00601 }
00602
00603 vector<float> EMAN2Ctf::to_vector() const
00604 {
00605 vector<float> vctf;
00606
00607 vctf.push_back(defocus);
00608 vctf.push_back(dfdiff);
00609 vctf.push_back(dfang);
00610 vctf.push_back(bfactor);
00611 vctf.push_back(ampcont);
00612 vctf.push_back(voltage);
00613 vctf.push_back(cs);
00614 vctf.push_back(apix);
00615 vctf.push_back(dsbg);
00616 vctf.push_back((float)background.size());
00617 for (unsigned int i=0; i<background.size(); i++) vctf.push_back(background[i]);
00618 vctf.push_back((float)snr.size());
00619 for (unsigned int j=0; j<snr.size(); j++) vctf.push_back(snr[j]);
00620
00621 return vctf;
00622 }
00623
00624
00625
00626 void EMAN2Ctf::copy_from(const Ctf * new_ctf)
00627 {
00628 if (new_ctf) {
00629 const EMAN2Ctf *c = static_cast<const EMAN2Ctf *>(new_ctf);
00630 defocus = c->defocus;
00631 dfdiff = c->dfdiff;
00632 dfang = c->dfang;
00633 bfactor = c->bfactor;
00634 ampcont = c->ampcont;
00635 voltage = c->voltage;
00636 cs = c->cs;
00637 apix = c->apix;
00638 dsbg = c->dsbg;
00639 background = c->background;
00640 snr = c->snr;
00641 }
00642 }
00643
00644 inline int max_int(int a,int b) { return a>b?a:b; }
00645 inline int min_int(int a,int b) { return a<b?a:b; }
00646
00647 vector < float >EMAN2Ctf::compute_1d(int size,float ds, CtfType type, XYData * sf)
00648 {
00649 Assert(size > 0);
00650
00651
00652
00653 int np=size/2;
00654 vector < float >r;
00655
00656 r.resize(np);
00657
00658
00659 float s = 0;
00660 float g1 = calc_g1();
00661 float g2 = calc_g2();
00662 float amp1 = calc_amp1();
00663
00664 switch (type) {
00665 case CTF_AMP:
00666 for (int i = 0; i < np; i++) {
00667 float gamma = calc_gamma(g1, g2, s);
00668 r[i] = calc_ctf1(amp1, gamma, s);
00669 s += ds;
00670 }
00671 break;
00672
00673 case CTF_SIGN:
00674 for (int i = 0; i < np; i++) {
00675 float gamma = calc_gamma(g1, g2, s);
00676 r[i] = calc_ctf1(amp1, gamma, s)>=0?1.0f:-1.0f;
00677 s += ds;
00678 }
00679 break;
00680
00681 case CTF_BACKGROUND:
00682 for (int i = 0; i < np; i++) {
00683 float f = s/dsbg;
00684 int j = (int)floor(f);
00685 f-=j;
00686 if (j>(int)background.size()-2) r[i]=background.back();
00687 else r[i]=background[j]*(1.0f-f)+background[j+1]*f;
00688 s+=ds;
00689 }
00690 break;
00691
00692 case CTF_SNR:
00693 if (snr.size()<1) {
00694 throw InvalidValueException("CTF_SNR"," ERROR: No SNR info in CTF header.");
00695 break;
00696 }
00697 for (int i = 0; i < np; i++) {
00698 float f = s/dsbg;
00699 int j = (int)floor(f);
00700 f-=j;
00701 if (j>(int)snr.size()-2) r[i]=snr.back();
00702 else r[i]=snr[j]*(1.0f-f)+snr[j+1]*f;
00703
00704 s+=ds;
00705 }
00706 break;
00707 case CTF_SNR_SMOOTH:
00708
00709
00710
00711
00712 {
00713 vector < float >tsnr;
00714 tsnr.resize(np);
00715 vector < float >dsnr;
00716 dsnr.resize(np);
00717
00718 float s0=s;
00719
00720 for (int i = 0; i < np; i++) {
00721 float gamma = calc_gamma(g1, g2, s);
00722 tsnr[i] = calc_ctf1(amp1, gamma, s);
00723
00724
00725 float f = s/dsbg;
00726 int j = (int)floor(f);
00727 f-=j;
00728 float bg;
00729 if (j>(int)background.size()-2) bg=background.back();
00730 else bg=background[j]*(1.0f-f)+background[j+1]*f;
00731 if (bg <=0) bg=.001f;
00732
00733 tsnr[i] = tsnr[i]*tsnr[i]/bg;
00734 if (sf && s) {
00735 tsnr[i] *= sf->get_yatx(s);
00736 }
00737
00738
00739
00740 if (j>(int)snr.size()-2) dsnr[i]=snr.back();
00741 else dsnr[i]=snr[j]*(1.0f-f)+snr[j+1]*f;
00742
00743 s+=ds;
00744 }
00745
00746 int npsm=np/25;
00747 if (npsm<2) npsm=2;
00748
00749 s=s0;
00750 for (int i = 1; i < np; i++) {
00751
00752 double sum = 0;
00753 double sum_x = 0;
00754 double sum_y = 0;
00755 double sum_xx = 0;
00756 double sum_xy = 0;
00757
00758 for (int k=max_int(i-npsm,1); k<=min_int(i+npsm,np-1); k++) {
00759 double y = dsnr[k];
00760 double x = tsnr[k];
00761
00762 sum_x += x;
00763 sum_y += y;
00764 sum_xx += x * x;
00765 sum_xy += x * y;
00766 sum++;
00767 }
00768
00769 double div = sum * sum_xx - sum_x * sum_x;
00770
00771
00772
00773
00774
00775
00776
00777 if (div!=0.0) r[i]=(float) ((sum * sum_xy - sum_x * sum_y) / div)*tsnr[i];
00778 else r[i]=0.0;
00779 if (r[i]<0) r[i]=0;
00780
00781 s+=ds;
00782 }
00783 r[0]=0;
00784 }
00785 break;
00786
00787 case CTF_WIENER_FILTER:
00788
00789
00790
00791
00792
00793 for (int i = 0; i < np; i++) {
00794 float f = s/dsbg;
00795 int j = (int)floor(f);
00796 float bg;
00797 f-=j;
00798 if (j>(int)snr.size()-2) {
00799
00800
00801 r[i]=0;
00802 }
00803 else {
00804 r[i]=snr[j]*(1.0f-f)+snr[j+1]*f;
00805 bg=background[j]*(1.0f-f)+background[j+1]*f;
00806 }
00807 if (r[i]<0) r[i]=0;
00808 r[i]=r[i]/(r[i]+1.0f);
00809
00810 s+=ds;
00811 }
00812 r[0]=0;
00813 break;
00814
00815 case CTF_TOTAL:
00816
00817 for (int i = 0; i < np; i++) {
00818 float gamma = calc_gamma(g1, g2, s);
00819 if (sf) {
00820 r[i] = calc_ctf1(amp1, gamma, s);
00821 r[i] = r[i] * r[i] * sf->get_yatx(s) + calc_noise(s);
00822 }
00823 else {
00824 r[i] = calc_ctf1(amp1, gamma, s);
00825 r[i] = r[i] * r[i] + calc_noise(s);
00826 }
00827 s += ds;
00828 }
00829 break;
00830 default:
00831 break;
00832 }
00833
00834 return r;
00835 }
00836
00837
00838 void EMAN2Ctf::compute_2d_real(EMData *, CtfType, XYData *)
00839 {
00840
00841
00842 }
00843
00844
00845
00846 void EMAN2Ctf::compute_2d_complex(EMData * image, CtfType type, XYData * sf)
00847 {
00848 if (!image) {
00849 LOGERR("image is null. cannot computer 2D complex CTF");
00850 return;
00851 }
00852
00853 if (image->is_complex() == false) {
00854 LOGERR("compute_2d_complex can only work on complex images");
00855 return;
00856 }
00857
00858 int nx = image->get_xsize();
00859 int ny = image->get_ysize();
00860
00861 if ((ny%2==1 && nx!=ny+1) || (ny%2==0 && nx != ny + 2)) {
00862 printf("nx=%d ny=%d\n",nx,ny);
00863 LOGERR("compute_2d_complex only works on (nx, nx-2) images");
00864 return;
00865 }
00866
00867 float ds = 1.0f / (apix * ny);
00868 image->to_one();
00869
00870 float *d = image->get_data();
00871 float g1 = calc_g1();
00872 float g2 = calc_g2();
00873 float amp1 = calc_amp1();
00874
00875 if (type == CTF_BACKGROUND) {
00876 for (int y = -ny/2; y < ny/2; y++) {
00877 int y2=(y+ny)%ny;
00878 int ynx = y2 * nx;
00879
00880 for (int x = 0; x < nx / 2; x++) {
00881 float s = (float) Util::hypot_fast(x, y ) * ds;
00882 d[x * 2 + ynx] = calc_noise(s);
00883 d[x * 2 + ynx + 1] = 0;
00884 }
00885 }
00886 }
00887 else if (type == CTF_AMP) {
00888 for (int y = -ny/2; y < ny/2; y++) {
00889 int y2=(y+ny)%ny;
00890 int ynx = y2 * nx;
00891
00892 for (int x = 0; x < nx / 2; x++) {
00893 float s = (float)Util::hypot_fast(x,y ) * ds;
00894 float gamma = calc_gamma(g1, g2, s);
00895 float v = calc_ctf1(amp1, gamma, s);
00896
00897 d[x * 2 + ynx] = v;
00898 d[x * 2 + ynx + 1] = 0;
00899 }
00900 }
00901 }
00902 else if (type == CTF_SIGN) {
00903 for (int y = -ny/2; y < ny/2; y++) {
00904 int y2=(y+ny)%ny;
00905 int ynx = y2 * nx;
00906 for (int x = 0; x < nx / 2; x++) {
00907 float s = (float)Util::hypot_fast(x,y ) * ds;
00908 float gamma = calc_gamma(g1, g2, s);
00909 float v = calc_amplitude(gamma);
00910 d[x * 2 + ynx] = v >= 0 ? 1.0f : -1.0f;
00911 d[x * 2 + ynx + 1] = 0;
00912 }
00913 }
00914 }
00915 else if (type == CTF_SNR) {
00916
00917 for (int y = -ny/2; y < ny/2; y++) {
00918 int y2=(y+ny)%ny;
00919 int ynx = y2 * nx;
00920
00921 for (int x = 0; x < nx / 2; x++) {
00922
00923 float s = (float)Util::hypot_fast(x,y ) * ds;
00924 float f = s/dsbg;
00925 int j = (int)floor(f);
00926 f-=j;
00927 if (j>(int)snr.size()-2) d[x*2+ynx]=snr.back();
00928 else d[x*2+ynx]=snr[j]*(1.0f-f)+snr[j+1]*f;
00929 d[x * 2 + ynx + 1] = 0;
00930 }
00931 }
00932 d[0]=0;
00933 }
00934 else if (type == CTF_SNR_SMOOTH) {
00935 for (int y = -ny/2; y < ny/2; y++) {
00936 int y2=(y+ny)%ny;
00937 int ynx = y2 * nx;
00938
00939 for (int x = 0; x < nx / 2; x++) {
00940
00941 float s = (float)Util::hypot_fast(x,y ) * ds;
00942 float f = s/dsbg;
00943 int j = (int)floor(f);
00944 f-=j;
00945 if (j>(int)snr.size()-2) d[x*2+ynx]=snr.back();
00946 else d[x*2+ynx]=snr[j]*(1.0f-f)+snr[j+1]*f;
00947 d[x * 2 + ynx + 1] = 0;
00948 }
00949 }
00950 d[0]=0;
00951 }
00952 else if (type == CTF_WIENER_FILTER) {
00953 if (dsbg==0) printf("Warning, DSBG set to 0\n");
00954 for (int y = -ny/2; y < ny/2; y++) {
00955 int y2=(y+ny)%ny;
00956 int ynx = y2 * nx;
00957
00958 for (int x = 0; x < nx / 2; x++) {
00959
00960 float s = (float)Util::hypot_fast(x,y ) * ds;
00961 float f = s/dsbg;
00962 int j = (int)floor(f);
00963 float bg,snrf;
00964 f-=j;
00965 if (j>(int)snr.size()-2) {
00966
00967
00968 d[x*2+ynx]=0;
00969 }
00970 else {
00971 bg=background[j]*(1.0f-f)+background[j+1]*f;
00972 snrf=snr[j]*(1.0f-f)+snr[j+1]*f;
00973
00974
00975 if (snrf<0) snrf=0.0;
00976
00977 d[x*2+ynx]=snrf/(snrf+1);
00978
00979 }
00980 d[x * 2 + ynx + 1] = 0;
00981 }
00982 }
00983 d[0]=0;
00984 }
00985 else if (type == CTF_TOTAL) {
00986 float amp1 = calc_amp1();
00987
00988 for (int y = -ny/2; y < ny/2; y++) {
00989 int y2=(y+ny)%ny;
00990 int ynx = y2 * nx;
00991
00992 for (int x = 0; x < nx / 2; x++) {
00993
00994 float s = (float)Util::hypot_fast(x,y ) * ds;
00995 float gamma = calc_gamma(g1, g2, s);
00996 float f = calc_ctf1(amp1, gamma, s);
00997 float noise = 0;
00998 f = f * f;
00999
01000 if (sf && s) {
01001 f *= sf->get_yatx(s);
01002 }
01003 f+=noise;
01004
01005 d[x * 2 + ynx] *= f;
01006 d[x * 2 + ynx + 1] = 0;
01007 }
01008 }
01009 }
01010 else printf("Unknown CTF image mode\n");
01011
01012 image->update();
01013 }
01014
01015
01016
01017 bool EMAN2Ctf::equal(const Ctf * ctf1) const
01018 {
01019 if (ctf1) {
01020 const EMAN2Ctf *c = static_cast<const EMAN2Ctf *>(ctf1);
01021 if (defocus != c->defocus ||
01022 dfdiff != c->dfdiff ||
01023 dfang != c->dfang ||
01024 bfactor != c->bfactor ||
01025 ampcont != c->ampcont ||
01026 voltage != c->voltage ||
01027 cs != c->cs ||
01028 apix != c->apix ||
01029 dsbg != c->dsbg ||
01030 background.size() != c->background.size() ||
01031 snr.size() != c->snr.size()
01032 ) return false;
01033
01034 for (unsigned int i=0; i<background.size(); i++) {
01035 if (background[i]!=c->background[i]) return false;
01036 }
01037 for (unsigned int i=0; i<snr.size(); i++) {
01038 if (snr[i]!=c->snr[i]) return false;
01039 }
01040 return true;
01041 }
01042 return false;
01043 }