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] *= pow(10.0f, 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] *= pow(10.0f, 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] * pow(10.0f, 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 *= pow(10.0f, 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 *= pow(10.0f, 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 *= pow(10.0f, 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 for (int i = 0; i < np; i++) {
00701 float gamma = calc_gamma(g1, g2, s);
00702 float f = s/dsbg;
00703 int j = (int)floor(f);
00704
00705 double sum=0,norm=0;
00706 for (int k=max_int(j-5,1); k<min_int(j+6,np); k++) {
00707 float s2=f-k;
00708 float e=exp(-14.0f*s2*s2/(gamma));
00709 if (k>(int)snr.size()) break;
00710 sum+=e*snr[k];
00711 norm+=e;
00712
00713 }
00714 r[i]=(float)(norm==0?0:sum/norm);
00715 s+=ds;
00716 }
00717 r[0]=0;
00718 break;
00719
00720 case CTF_WIENER_FILTER:
00721
00722
00723
00724
00725
00726 for (int i = 0; i < np; i++) {
00727 float f = s/dsbg;
00728 int j = (int)floor(f);
00729 float bg;
00730 f-=j;
00731 if (j>(int)snr.size()-2) {
00732
00733
00734 r[i]=0;
00735 }
00736 else {
00737 r[i]=snr[j]*(1.0f-f)+snr[j+1]*f;
00738 bg=background[j]*(1.0f-f)+background[j+1]*f;
00739 }
00740 if (r[i]<0) r[i]=0;
00741 r[i]=r[i]/(r[i]+1.0f);
00742
00743 s+=ds;
00744 }
00745 r[0]=0;
00746 break;
00747
00748 case CTF_TOTAL:
00749 if (!sf) {
00750 LOGERR("CTF computation error, no SF found\n");
00751 return r;
00752 }
00753
00754 for (int i = 0; i < np; i++) {
00755 float gamma = calc_gamma(g1, g2, s);
00756 if (sf) {
00757 r[i] = calc_ctf1(amp1, gamma, s);
00758
00759 }
00760 else {
00761 r[i] = calc_ctf1(amp1, gamma, s);
00762
00763 }
00764 s += ds;
00765 }
00766 break;
00767 default:
00768 break;
00769 }
00770
00771 return r;
00772 }
00773
00774
00775 void EMAN2Ctf::compute_2d_real(EMData *, CtfType, XYData *)
00776 {
00777
00778
00779 }
00780
00781
00782
00783 void EMAN2Ctf::compute_2d_complex(EMData * image, CtfType type, XYData * sf)
00784 {
00785 if (!image) {
00786 LOGERR("image is null. cannot computer 2D complex CTF");
00787 return;
00788 }
00789
00790 if (image->is_complex() == false) {
00791 LOGERR("compute_2d_complex can only work on complex images");
00792 return;
00793 }
00794
00795 int nx = image->get_xsize();
00796 int ny = image->get_ysize();
00797
00798 if ((ny%2==1 && nx!=ny+1) || (ny%2==0 && nx != ny + 2)) {
00799 printf("nx=%d ny=%d\n",nx,ny);
00800 LOGERR("compute_2d_complex only works on (nx, nx-2) images");
00801 return;
00802 }
00803
00804 float ds = 1.0f / (apix * ny);
00805 image->to_one();
00806
00807 float *d = image->get_data();
00808 float g1 = calc_g1();
00809 float g2 = calc_g2();
00810 float amp1 = calc_amp1();
00811
00812 if (type == CTF_BACKGROUND) {
00813 for (int y = -ny/2; y < ny/2; y++) {
00814 int y2=(y+ny)%ny;
00815 int ynx = y2 * nx;
00816
00817 for (int x = 0; x < nx / 2; x++) {
00818 float s = (float) Util::hypot_fast(x, y ) * ds;
00819 d[x * 2 + ynx] = calc_noise(s);
00820 d[x * 2 + ynx + 1] = 0;
00821 }
00822 }
00823 }
00824 else if (type == CTF_AMP) {
00825 for (int y = -ny/2; y < ny/2; y++) {
00826 int y2=(y+ny)%ny;
00827 int ynx = y2 * nx;
00828
00829 for (int x = 0; x < nx / 2; x++) {
00830 float s = (float)Util::hypot_fast(x,y ) * ds;
00831 float gamma = calc_gamma(g1, g2, s);
00832 float v = calc_ctf1(amp1, gamma, s);
00833
00834 d[x * 2 + ynx] = v;
00835 d[x * 2 + ynx + 1] = 0;
00836 }
00837 }
00838 }
00839 else if (type == CTF_SIGN) {
00840 for (int y = -ny/2; y < ny/2; y++) {
00841 int y2=(y+ny)%ny;
00842 int ynx = y2 * nx;
00843 for (int x = 0; x < nx / 2; x++) {
00844 float s = (float)Util::hypot_fast(x,y ) * ds;
00845 float gamma = calc_gamma(g1, g2, s);
00846 float v = calc_amplitude(gamma);
00847 d[x * 2 + ynx] = v >= 0 ? 1.0f : -1.0f;
00848 d[x * 2 + ynx + 1] = 0;
00849 }
00850 }
00851 }
00852 else if (type == CTF_SNR) {
00853
00854 for (int y = -ny/2; y < ny/2; y++) {
00855 int y2=(y+ny)%ny;
00856 int ynx = y2 * nx;
00857
00858 for (int x = 0; x < nx / 2; x++) {
00859
00860 float s = (float)Util::hypot_fast(x,y ) * ds;
00861 float f = s/dsbg;
00862 int j = (int)floor(f);
00863 f-=j;
00864 if (j>(int)snr.size()-2) d[x*2+ynx]=snr.back();
00865 else d[x*2+ynx]=snr[j]*(1.0f-f)+snr[j+1]*f;
00866 d[x * 2 + ynx + 1] = 0;
00867 }
00868 }
00869 d[0]=0;
00870 }
00871 else if (type == CTF_SNR_SMOOTH) {
00872 for (int y = -ny/2; y < ny/2; y++) {
00873 int y2=(y+ny)%ny;
00874 int ynx = y2 * nx;
00875
00876 for (int x = 0; x < nx / 2; x++) {
00877
00878 float s = (float)Util::hypot_fast(x,y ) * ds;
00879 float f = s/dsbg;
00880 int j = (int)floor(f);
00881 f-=j;
00882 if (j>(int)snr.size()-2) d[x*2+ynx]=snr.back();
00883 else d[x*2+ynx]=snr[j]*(1.0f-f)+snr[j+1]*f;
00884 d[x * 2 + ynx + 1] = 0;
00885 }
00886 }
00887 d[0]=0;
00888 }
00889 else if (type == CTF_WIENER_FILTER) {
00890 if (dsbg==0) printf("Warning, DSBG set to 0\n");
00891 for (int y = -ny/2; y < ny/2; y++) {
00892 int y2=(y+ny)%ny;
00893 int ynx = y2 * nx;
00894
00895 for (int x = 0; x < nx / 2; x++) {
00896
00897 float s = (float)Util::hypot_fast(x,y ) * ds;
00898 float f = s/dsbg;
00899 int j = (int)floor(f);
00900 float bg,snrf;
00901 f-=j;
00902 if (j>(int)snr.size()-2) {
00903
00904
00905 d[x*2+ynx]=0;
00906 }
00907 else {
00908 bg=background[j]*(1.0f-f)+background[j+1]*f;
00909 snrf=snr[j]*(1.0f-f)+snr[j+1]*f;
00910
00911
00912 if (snrf<0) snrf=0.0;
00913
00914 d[x*2+ynx]=snrf/(snrf+1);
00915
00916 }
00917 d[x * 2 + ynx + 1] = 0;
00918 }
00919 }
00920 d[0]=0;
00921 }
00922 else if (type == CTF_TOTAL) {
00923 float amp1 = calc_amp1();
00924
00925 for (int y = -ny/2; y < ny/2; y++) {
00926 int y2=(y+ny)%ny;
00927 int ynx = y2 * nx;
00928
00929 for (int x = 0; x < nx / 2; x++) {
00930
00931 float s = (float)Util::hypot_fast(x,y ) * ds;
00932 float gamma = calc_gamma(g1, g2, s);
00933 float f = calc_ctf1(amp1, gamma, s);
00934 float noise = 0;
00935 f = f * f;
00936
00937 if (sf && s) {
00938 f *= pow(10.0f, sf->get_yatx(s));
00939 }
00940 f+=noise;
00941
00942 d[x * 2 + ynx] *= f;
00943 d[x * 2 + ynx + 1] = 0;
00944 }
00945 }
00946 }
00947 else printf("Unknown CTF image mode\n");
00948
00949 image->update();
00950 }
00951
00952
00953
00954 bool EMAN2Ctf::equal(const Ctf * ctf1) const
00955 {
00956 if (ctf1) {
00957 const EMAN2Ctf *c = static_cast<const EMAN2Ctf *>(ctf1);
00958 if (defocus != c->defocus ||
00959 dfdiff != c->dfdiff ||
00960 dfang != c->dfang ||
00961 bfactor != c->bfactor ||
00962 ampcont != c->ampcont ||
00963 voltage != c->voltage ||
00964 cs != c->cs ||
00965 apix != c->apix ||
00966 dsbg != c->dsbg ||
00967 background.size() != c->background.size() ||
00968 snr.size() != c->snr.size()
00969 ) return false;
00970
00971 for (unsigned int i=0; i<background.size(); i++) {
00972 if (background[i]!=c->background[i]) return false;
00973 }
00974 for (unsigned int i=0; i<snr.size(); i++) {
00975 if (snr[i]!=c->snr[i]) return false;
00976 }
00977 return true;
00978 }
00979 return false;
00980 }