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