EMAN2
ctf.cpp
Go to the documentation of this file.
00001 
00005 /*
00006  * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
00007  * Copyright (c) 2000-2006 Baylor College of Medicine
00008  *
00009  * This software is issued under a joint BSD/GNU license. You may use the
00010  * source code in this file under either license. However, note that the
00011  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00012  * so you are responsible for compliance with the licenses of these packages
00013  * if you opt to use BSD licensing. The warranty disclaimer below holds
00014  * in either instance.
00015  *
00016  * This complete copyright notice must be included in any revised version of the
00017  * source code. Additional authorship citations may be added, but existing
00018  * author citations must be preserved.
00019  *
00020  * This program is free software; you can redistribute it and/or modify
00021  * it under the terms of the GNU General Public License as published by
00022  * the Free Software Foundation; either version 2 of the License, or
00023  * (at your option) any later version.
00024  *
00025  * This program is distributed in the hope that it will be useful,
00026  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00027  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00028  * GNU General Public License for more details.
00029  *
00030  * You should have received a copy of the GNU General Public License
00031  * along with this program; if not, write to the Free Software
00032  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
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, &amplitude, &ampcont, &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 //      float ds = 1 / (apix * size * CTFOS);
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 //              if (!sf) {
00220 //                      LOGERR("CTF computation error, no SF found\n");
00221 //                      return r;
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 vector <float> EMAN1Ctf::compute_1d_fromimage(int size, float ds, EMData *image)
00287 {
00288         vector < float >r;
00289 
00290         printf("Sorry, this function hasn't been implemented for EMAN1Ctf\n");
00291         
00292         int np=size/2;
00293         r.resize(np);
00294 
00295         return r;
00296 }
00297 
00298 
00299 void EMAN1Ctf::compute_2d_complex(EMData * image, CtfType type, XYData * sf)
00300 {
00301         if (!image) {
00302                 LOGERR("image is null. cannot computer 2D complex CTF");
00303                 return;
00304         }
00305 
00306         if (image->is_complex() == false) {
00307                 LOGERR("compute_2d_complex can only work on complex images");
00308                 return;
00309         }
00310 
00311         int nx = image->get_xsize();
00312         int ny = image->get_ysize();
00313 
00314         if (nx != ny + 2) {
00315                 LOGERR("compute_2d_complex only works on (nx, nx-2) images");
00316                 return;
00317         }
00318 
00319         float ds = 1.0f / (apix * ny);
00320         image->to_one();
00321 
00322         float *d = image->get_data();
00323         float g1 = calc_g1();
00324         float g2 = calc_g2();
00325 
00326         if (type == CTF_BACKGROUND) {
00327                 for (int y = 0; y < ny; y++) {
00328                         int ynx = y * nx;
00329 
00330                         for (int x = 0; x < nx / 2; x++) {
00331 #ifdef  _WIN32
00332                                 float s = (float) _hypot(x, y - ny / 2.0f) * ds;
00333 #else
00334                                 float s = (float) hypot(x, y - ny / 2.0f) * ds;
00335 #endif
00336                                 d[x * 2 + ynx] = calc_noise(s);
00337                                 d[x * 2 + ynx + 1] = 0;                 // The phase is somewhat arbitrary
00338                         }
00339                 }
00340         }
00341         else if (type == CTF_AMP) {
00342                 for (int y = 0; y < ny; y++) {
00343                         int ynx = y * nx;
00344 
00345                         for (int x = 0; x < nx / 2; x++) {
00346 #ifdef  _WIN32
00347                                 float s = (float)_hypot((float) x, (float) y - ny / 2) * ds;
00348 #else
00349                                 float s = (float)hypot((float) x, (float) y - ny / 2) * ds;
00350 #endif  //_WIN32
00351                                 float gamma = calc_gamma(g1, g2, s);
00352                                 float v = fabs(calc_amplitude(gamma));
00353                                 d[x * 2 + ynx] = v;
00354                                 d[x * 2 + ynx + 1] = 0;
00355                         }
00356                 }
00357         }
00358         else if (type == CTF_SIGN) {
00359                 for (int y = 0; y < ny; y++) {
00360                         int ynx = y * nx;
00361                         for (int x = 0; x < nx / 2; x++) {
00362 #ifdef  _WIN32
00363                                 float s = (float)_hypot(x, y - ny / 2.0f) * ds;
00364 #else
00365                                 float s = (float)hypot(x, y - ny / 2.0f) * ds;
00366 #endif
00367                                 float gamma = calc_gamma(g1, g2, s);
00368                                 float v = calc_amplitude(gamma);
00369                                 d[x * 2 + ynx] = v > 0 ? 1.0f : -1.0f;
00370                                 d[x * 2 + ynx + 1] = 0;
00371                         }
00372                 }
00373 
00374         }
00375         else if (type == CTF_SNR || type == CTF_SNR_SMOOTH) {
00376                 float amp1 = calc_amp1();
00377 
00378                 for (int y = 0; y < ny; y++) {
00379                         int ynx = y * nx;
00380 
00381                         for (int x = 0; x < nx / 2; x++) {
00382 
00383 #ifdef  _WIN32
00384                                 float s = (float)_hypot(x, y - ny / 2.0f) * ds;
00385 #else
00386                                 float s = (float)hypot(x, y - ny / 2.0f) * ds;
00387 #endif
00388                                 float gamma = calc_gamma(g1, g2, s);
00389                                 float f = calc_ctf1(amp1, gamma, s);
00390                                 float noise = calc_noise(s);
00391                                 f = f * f / noise;
00392 
00393                                 if (s && sf) {
00394                                         f *= sf->get_yatx(s);
00395                                 }
00396                                 d[x * 2 + ynx] *= f;
00397                                 d[x * 2 + ynx + 1] = 0;
00398                         }
00399                 }
00400         }
00401         else if (type == CTF_WIENER_FILTER) {
00402                 float amp1 = calc_amp1();
00403 
00404                 for (int y = 0; y < ny; y++) {
00405                         int ynx = y * nx;
00406 
00407                         for (int x = 0; x < nx / 2; x++) {
00408 
00409 #ifdef  _WIN32
00410                                 float s = (float)_hypot(x, y - ny / 2.0f) * ds;
00411 #else
00412                                 float s = (float)hypot(x, y - ny / 2.0f) * ds;
00413 #endif
00414                                 float gamma = calc_gamma(g1, g2, s);
00415                                 float f = calc_ctf1(amp1, gamma, s);
00416                                 float noise = calc_noise(s);
00417                                 f = f * f / noise;
00418 
00419                                 if (s) {
00420                                         f *= sf->get_yatx(s);
00421                                 }
00422                                 f = 1.0f / (1.0f + 1.0f / f);
00423                                 d[x * 2 + ynx] *= f;
00424                                 d[x * 2 + ynx + 1] = 0;
00425                         }
00426                 }
00427         }
00428         else if (type == CTF_TOTAL) {
00429                 float amp1 = calc_amp1();
00430 
00431                 for (int y = 0; y < ny; y++) {
00432                         int ynx = y * nx;
00433 
00434                         for (int x = 0; x < nx / 2; x++) {
00435 
00436 #ifdef  _WIN32
00437                                 float s = (float)_hypot(x, y - ny / 2.0f) * ds;
00438 #else
00439                                 float s = (float)hypot(x, y - ny / 2.0f) * ds;
00440 #endif
00441                                 float gamma = calc_gamma(g1, g2, s);
00442                                 float f = calc_ctf1(amp1, gamma, s);
00443                                 float noise = calc_noise(s);
00444                                 f = f * f;
00445 
00446                                 if (sf && s) {
00447                                         f *= sf->get_yatx(s);
00448                                 }
00449                                 f+=noise;
00450 
00451                                 d[x * 2 + ynx] *= f;
00452                                 d[x * 2 + ynx + 1] = 0;
00453                         }
00454                 }
00455         }
00456 
00457         image->update();
00458 }
00459 
00460 
00461 
00462 bool EMAN1Ctf::equal(const Ctf * ctf1) const
00463 {
00464         if (ctf1) {
00465                 const EMAN1Ctf *c = static_cast<const EMAN1Ctf *>(ctf1);
00466                 if (defocus == c->defocus &&
00467                         bfactor == c->bfactor &&
00468                         amplitude == c->amplitude &&
00469                         ampcont == c->ampcont &&
00470                         noise1 == c->noise1 &&
00471                         noise2 == c->noise2 &&
00472                         noise3 == c->noise3 &&
00473                         noise4 == c->noise4 && voltage == c->voltage && cs == c->cs && apix == c->apix) {
00474                         return true;
00475                 }
00476         }
00477         return false;
00478 }
00479 
00480 float EMAN1Ctf::zero(int n) const { return 0; }
00481 
00482 /*************************************
00483 EMAN2Ctf
00484 *************************************/
00485 
00486 EMAN2Ctf::EMAN2Ctf()
00487 {
00488         defocus = 0;
00489         dfdiff = 0;
00490         dfang = 0;
00491         bfactor = 0;
00492         ampcont = 0;
00493         voltage = 0;
00494         cs = 0;
00495         apix = 1.0;
00496         dsbg=-1;
00497         background.clear();
00498         snr.clear();
00499 }
00500 
00501 
00502 EMAN2Ctf::~EMAN2Ctf()
00503 {
00504 }
00505 
00506 
00507 int EMAN2Ctf::from_string(const string & ctf)
00508 {
00509         Assert(ctf != "");
00510         char type=' ';
00511         int pos=-1,i,j;
00512         int bglen=0,snrlen=0;
00513         float v;
00514         const char *s=ctf.c_str();
00515 
00516         sscanf(s, "%c%f %f %f %f %f %f %f %f %f %d%n",
00517                                    &type,&defocus, &dfdiff,&dfang,&bfactor,&ampcont,&voltage, &cs, &apix,&dsbg,&bglen,&pos);
00518         if (type!='E') throw InvalidValueException(type,"Trying to initialize Ctf object with bad string");
00519         if (pos==-1) throw InvalidValueException(s," Invalid CTF string");
00520         
00521 
00522         background.resize(bglen);
00523         for (i=0; i<bglen; i++) {
00524                 if (sscanf(s+pos,",%f%n",&v,&j)<1) return(1);
00525                 background[i]=v;
00526                 pos+=j;
00527         }
00528 
00529         sscanf(s+pos," %d%n",&snrlen,&j);
00530         pos+=j;
00531         snr.resize(snrlen);
00532         for (i=0; i<snrlen; i++) {
00533                 if (sscanf(s+pos,",%f%n",&v,&j)<1) return(1);
00534                 snr[i]=v;
00535                 pos+=j;
00536         }
00537 
00538         return 0;
00539 
00540 }
00541 
00542 string EMAN2Ctf::to_string() const
00543 {
00544         char ctf[256];
00545         sprintf(ctf, "E%1.4g %1.4g %1.4g %1.4g %1.4g %1.4g %1.4g %1.4g %1.4g %d",
00546                         defocus, dfdiff, dfang, bfactor, ampcont, voltage, cs, apix, dsbg,(int)background.size());
00547 
00548         string ret=ctf;
00549         for (int i=0; i<(int)background.size(); i++) {
00550                 sprintf(ctf,",%1.3g",background[i]);
00551                 ret+=ctf;
00552         }
00553 
00554         sprintf(ctf, " %d",(int)snr.size());
00555         ret+=ctf;
00556         for (int i=0; i<(int)snr.size(); i++) {
00557                 sprintf(ctf,",%1.3g",snr[i]);
00558                 ret+=ctf;
00559         }
00560 
00561 
00562         return ret;
00563 }
00564 
00565 void EMAN2Ctf::from_dict(const Dict & dict)
00566 {
00567 //      printf("st\n"); 
00568         defocus = (float)dict["defocus"];
00569 //      printf("df\n");
00570         dfdiff = (float)dict["dfdiff"];
00571 //      printf("dfd\n");
00572         dfang = (float)dict["dfang"];
00573 //      printf("dfa\n");
00574         bfactor = (float)dict["bfactor"];
00575 //      printf("bf\n");
00576         ampcont = (float)dict["ampcont"];
00577 //      printf("ac\n");
00578         voltage = (float)dict["voltage"];
00579 //      printf("vo\n");
00580         cs = (float)dict["cs"];
00581 //      printf("cs\n");
00582         apix = (float)dict["apix"];
00583 //      printf("ap\n");
00584         dsbg = (float)dict["dsbg"];
00585 //      printf("ds\n");
00586         background = dict["background"];
00587 //      printf("bg\n");
00588         snr = dict["snr"];
00589 //      printf("done\n");
00590 }
00591 
00592 Dict EMAN2Ctf::to_dict() const
00593 {
00594         Dict dict;
00595         dict["defocus"] = defocus;
00596         dict["dfdiff"] = dfdiff;
00597         dict["dfang"] = dfang;
00598         dict["bfactor"] = bfactor;
00599         dict["ampcont"] = ampcont;
00600         dict["voltage"] = voltage;
00601         dict["cs"] = cs;
00602         dict["apix"] = apix;
00603         dict["dsbg"] = dsbg;
00604         dict["background"] = background;
00605         dict["snr"] = snr;
00606 
00607         return dict;
00608 }
00609 
00610 void EMAN2Ctf::from_vector(const vector<float>& vctf)
00611 {
00612         int i;
00613         defocus = vctf[0];
00614         dfdiff = vctf[1];
00615         dfang = vctf[2];
00616         bfactor = vctf[3];
00617         ampcont = vctf[4];
00618         voltage = vctf[5];
00619         cs = vctf[6];
00620         apix = vctf[7];
00621         dsbg = vctf[8];
00622         background.resize((int)vctf[9]);
00623         for (i=0; i<(int)vctf[9]; i++) background[i]=vctf[i+10];
00624         snr.resize((int)vctf[i+10]);
00625         for (int j=0; j<(int)vctf[i+10]; j++) snr[j]=vctf[i+j+11];
00626 }
00627 
00628 vector<float> EMAN2Ctf::to_vector() const
00629 {
00630         vector<float> vctf;
00631 
00632         vctf.push_back(defocus);
00633         vctf.push_back(dfdiff);
00634         vctf.push_back(dfang);
00635         vctf.push_back(bfactor);
00636         vctf.push_back(ampcont);
00637         vctf.push_back(voltage);
00638         vctf.push_back(cs);
00639         vctf.push_back(apix);
00640         vctf.push_back(dsbg);
00641         vctf.push_back((float)background.size());
00642         for (unsigned int i=0; i<background.size(); i++) vctf.push_back(background[i]);
00643         vctf.push_back((float)snr.size());
00644         for (unsigned int j=0; j<snr.size(); j++) vctf.push_back(snr[j]);
00645 
00646         return vctf;
00647 }
00648 
00649 
00650 
00651 void EMAN2Ctf::copy_from(const Ctf * new_ctf)
00652 {
00653         if (new_ctf) {
00654                 const EMAN2Ctf *c = static_cast<const EMAN2Ctf *>(new_ctf);
00655                 defocus = c->defocus;
00656                 dfdiff = c->dfdiff;
00657                 dfang = c->dfang;
00658                 bfactor = c->bfactor;
00659                 ampcont = c->ampcont;
00660                 voltage = c->voltage;
00661                 cs = c->cs;
00662                 apix = c->apix;
00663                 dsbg = c->dsbg;
00664                 background = c->background;
00665                 snr = c->snr;
00666         }
00667 }
00668 
00669 
00670 vector <float> EMAN2Ctf::compute_1d_fromimage(int size, float ds, EMData *image)
00671 {
00672         size/=2;        // to be consistent with the next routine
00673         vector <float> ret(size);
00674         vector <float> norm(size);
00675         
00676         if (!image->is_complex() || !image->is_ri()) ImageFormatException("EMAN2Ctf::compute_1d_fromimage requires a complex, ri image");
00677         int isinten=image->get_attr_default("is_intensity",0);
00678         if (isinten) ImageFormatException("EMAN2Ctf::compute_1d_fromimage does not take intensity images");
00679         
00680         int nx=image->get_attr("nx");
00681         int ny=image->get_attr("ny");
00682         float *data=image->get_data();
00683         
00684         int x,y,i;
00685         for (y=i=0; y<ny; y++) {
00686                 for (x=0; x<nx; x+=2,i+=2) {
00687                         if (x==0 && y>ny/2) continue;
00688                         float r=(float)(Util::hypot_fast(x/2,y<ny/2?y:ny-y));           // origin at 0,0; periodic
00689                         float a=(float)atan2((float)(y<ny/2?y:ny-y),(float)(x/2));              // angle to point (radians)
00690                         r=stos2(r*ds,-dfdiff/2*cos(2.0*a-2.0*M_PI/180.0*dfang))/ds;             // angle-based defocus -> convert s to remove astigmatism, dfang in degrees
00691                         int f=int(r);   // safe truncation, so floor isn't needed
00692                         r-=float(f);    // r is now the fractional spacing between bins
00693                         
00694                         float v=Util::square_sum(data[i],data[i+1]);                            //intensity
00695                         if (f>=0 && f<size) {
00696                                 ret[f]+=v*(1.0f-r);
00697                                 norm[f]+=(1.0f-r);
00698                                 if (f<size-1) {
00699                                         ret[f+1]+=v*r;
00700                                         norm[f+1]+=r;
00701                                 }
00702                         }
00703                 }
00704         }
00705 
00706         for (i=0; i<size; i++) ret[i]/=norm[i]?norm[i]:1.0f;    // Normalize
00707 
00708         
00709         return ret;
00710 }
00711 
00712 
00713 inline int max_int(int a,int b) { return a>b?a:b; }
00714 inline int min_int(int a,int b) { return a<b?a:b; }
00715 
00716 vector <float> EMAN2Ctf::compute_1d(int size,float ds, CtfType type, XYData * sf)
00717 {
00718         Assert(size > 0);
00719 
00720 //      float tmp_f1 =  sqrt((float) 2) * size / 2;
00721 //      int np = (int) ceil(tmp_f1) + 2;
00722         int np=size/2;
00723         vector < float >r;
00724 
00725         r.resize(np);
00726 
00727         float s = 0;
00728         float g1=M_PI/2.0*cs*1.0e7*pow(lambda(),3.0f);  // s^4 coefficient for gamma, cached in a variable for simplicity (maybe speed? depends on the compiler)
00729         float g2=M_PI*lambda()*defocus*10000.0;                 // s^2 coefficient for gamma 
00730         float acac=acos(ampcont/100.0);                                 // instead of ac*cos(g)+sqrt(1-ac^2)*sin(g), we can use cos(g-acos(ac)) and save a trig op
00731         switch (type) {
00732         case CTF_AMP:
00733                 for (int i = 0; i < np; i++) {
00734                         float gam=-g1*s*s*s*s+g2*s*s;
00735                         r[i] = cos(gam-acac)*exp(-(bfactor/4.0f * s * s));
00736                         s += ds;
00737                 }
00738                 break;
00739 
00740         case CTF_INTEN:
00741                 for (int i = 0; i < np; i++) {
00742                         float gam=-g1*s*s*s*s+g2*s*s;
00743                         r[i] = cos(gam-acac);
00744                         r[i]*=r[i];
00745                         s += ds;
00746                 }
00747                 break;
00748                 
00749         case CTF_SIGN:
00750                 for (int i = 0; i < np; i++) {
00751                         float gam=-g1*s*s*s*s+g2*s*s;
00752                         r[i] = cos(gam-acac)<0?-1.0:1.0;
00753                         s += ds;
00754                 }
00755                 break;
00756 
00757         case CTF_BACKGROUND:
00758                 for (int i = 0; i < np; i++) {
00759                         float f = s/dsbg;
00760                         int j = (int)floor(f);
00761                         f-=j;
00762                         if (j>(int)background.size()-2) r[i]=background.back();
00763                         else r[i]=background[j]*(1.0f-f)+background[j+1]*f;
00764                         s+=ds;
00765                 }
00766                 break;
00767 
00768         case CTF_SNR:
00769                 if (snr.size()<1) {
00770                         throw InvalidValueException("CTF_SNR"," ERROR: No SNR info in CTF header.");
00771                         break;
00772                 }
00773                 for (int i = 0; i < np; i++) {
00774                         float f = s/dsbg;
00775                         int j = (int)floor(f);
00776                         f-=j;
00777                         if (j>(int)snr.size()-2) r[i]=snr.back();
00778                         else r[i]=snr[j]*(1.0f-f)+snr[j+1]*f;
00779 //                      printf("%d\t%f\n",j,snr[j]);
00780                         s+=ds;
00781                 }
00782                 break;
00783         case CTF_SNR_SMOOTH:
00784                 // This apparently complicated routine tries to make a nice smooth and accurate SNR curve. It does this
00785                 // by fitting local regions of the SNR vs the theoretical SNR (theoretical CTF^2/measured background),
00786                 // then taking the slope of the result times the theoretical SNR to produce a local SNR estimate
00787 
00788                 { // <- is to permit new temporary value allocation
00789                         vector < float >tsnr;   // theoretical SNR
00790                         tsnr.resize(np);
00791                         vector < float >dsnr;   // data SNR
00792                         dsnr.resize(np);
00793                         
00794                         float s0=s;
00795                         
00796                         for (int i = 0; i < np; i++) {
00797                                 float gam=-g1*s*s*s*s+g2*s*s;
00798                                 tsnr[i] = cos(gam-acac)*exp(-(bfactor/4.0f * s * s));           // ctf amp
00799 
00800                                 // background value
00801                                 float f = s/dsbg;
00802                                 int j = (int)floor(f);
00803                                 f-=j;
00804                                 float bg;
00805                                 if (j>(int)background.size()-2) bg=background.back();
00806                                 else bg=background[j]*(1.0f-f)+background[j+1]*f;
00807                                 if (bg <=0) bg=.001f;
00808 
00809                                 tsnr[i] = tsnr[i]*tsnr[i]/bg;           // This is now a SNR curve
00810                                 if (sf && s) {
00811                                         tsnr[i] *= sf->get_yatx(s);
00812                                 }
00813 
00814                                 
00815                                 // This is the SNR computed from the data without fitting
00816                                 if (j>(int)snr.size()-2) dsnr[i]=snr.back();
00817                                 else dsnr[i]=snr[j]*(1.0f-f)+snr[j+1]*f;
00818                                 
00819                                 s+=ds;
00820                         }
00821 
00822                         int npsm=np/25;                 // 1/2 number of points to smooth over, 25 is arbitrary
00823                         if (npsm<2) npsm=2;
00824                         
00825                         s=s0;
00826                         for (int i = 1; i < np; i++) {
00827                                 // simple linear regression embedded here
00828                                 double sum = 0;
00829                                 double sum_x = 0;
00830                                 double sum_y = 0;
00831                                 double sum_xx = 0;
00832                                 double sum_xy = 0;
00833 
00834                                 for (int k=max_int(i-npsm,1); k<=min_int(i+npsm,np-1); k++) {
00835                                         double y = dsnr[k];
00836                                         double x = tsnr[k];
00837 
00838                                         sum_x += x;
00839                                         sum_y += y;
00840                                         sum_xx += x * x;
00841                                         sum_xy += x * y;
00842                                         sum++;
00843                                 }
00844 
00845                                 double div = sum * sum_xx - sum_x * sum_x;
00846 //                              if (div == 0) {
00847 //                                      div = 0.0000001f;
00848 //                              }
00849 
00850         //                      *intercept = (float) ((sum_xx * sum_y - sum_x * sum_xy) / div);
00851         //                      *slope = (float) ((sum * sum_xy - sum_x * sum_y) / div);
00852 
00853                                 if (div!=0.0) r[i]=(float) ((sum * sum_xy - sum_x * sum_y) / div)*tsnr[i];
00854                                 else r[i]=0.0;
00855                                 if (r[i]<0) r[i]=0;
00856                                 
00857                                 s+=ds;
00858                         }
00859                         r[0]=0;
00860                 }
00861                 break;
00862 
00863         case CTF_WIENER_FILTER:
00864 //              if (!sf) {
00865 //                      LOGERR("CTF computation error, no SF found\n");
00866 //                      return r;
00867 //              }
00868 
00869                 for (int i = 0; i < np; i++) {
00870                         float f = s/dsbg;
00871                         int j = (int)floor(f);
00872                         float bg;
00873                         f-=j;
00874                         if (j>(int)snr.size()-2) {
00875 /*                              r[i]=snr.back();
00876                                 bg=background.back();*/
00877                                 r[i]=0;
00878                         }
00879                         else {
00880                                 r[i]=snr[j]*(1.0f-f)+snr[j+1]*f;
00881                                 bg=background[j]*(1.0f-f)+background[j+1]*f;
00882                         }
00883                         if (r[i]<0) r[i]=0;
00884                         r[i]=r[i]/(r[i]+1.0f);          // Full Wiener filter assuming perfect image with noise
00885 //                      r[i]=sqrt(r[i]/bg)/(r[i]+1.0);  // Wiener filter with 1/CTF term (sort of) to restore image then filter
00886                         s+=ds;
00887                 }
00888                 r[0]=0;
00889                 break;
00890 
00891         case CTF_TOTAL:
00892 
00893                 for (int i = 0; i < np; i++) {
00894                         float gam=-g1*s*s*s*s+g2*s*s;
00895                         if (sf) {
00896                                 r[i] = cos(gam-acac)*exp(-(bfactor/4.0f * s*s));
00897                                 r[i] = r[i] * r[i] * sf->get_yatx(s) + calc_noise(s);
00898                         }
00899                         else {
00900                                 r[i] = cos(gam-acac)*exp(-(bfactor/4.0f * s*s));
00901                                 r[i] = r[i] * r[i] + calc_noise(s);
00902                         }
00903                         s += ds;
00904                 }
00905                 break;
00906         case CTF_NOISERATIO:
00907                 for (int i = 0; i < np; i++) {
00908                         float gam=-g1*s*s*s*s+g2*s*s;
00909                         if (sf) {
00910                                 r[i] = cos(gam-acac)*exp(-(bfactor/4.0f * s*s));
00911                                 r[i] = r[i] * r[i] * sf->get_yatx(s);
00912                                 r[i] = r[i]/(r[i]+calc_noise(s));
00913                         }
00914                         else {
00915                                 r[i] = cos(gam-acac)*exp(-(bfactor/4.0f * s*s));
00916                                 r[i] = r[i] * r[i];
00917                                 r[i] = r[i]/(r[i]+calc_noise(s));
00918                         }
00919                         
00920                         s+=ds;
00921                 }
00922                         
00923                 break;
00924         default:
00925                 break;
00926         }
00927 
00928         return r;
00929 }
00930 
00931 
00932 void EMAN2Ctf::compute_2d_real(EMData *, CtfType, XYData *)
00933 {
00934 
00935 
00936 }
00937 
00938 
00939 
00940 void EMAN2Ctf::compute_2d_complex(EMData * image, CtfType type, XYData * sf)
00941 {
00942         if (!image) {
00943                 LOGERR("image is null. cannot computer 2D complex CTF");
00944                 return;
00945         }
00946 
00947         if (image->is_complex() == false) {
00948                 LOGERR("compute_2d_complex can only work on complex images");
00949                 return;
00950         }
00951 
00952         int nx = image->get_xsize();
00953         int ny = image->get_ysize();
00954 
00955         if ((ny%2==1 && nx!=ny+1) || (ny%2==0 && nx != ny + 2)) {
00956                 printf("nx=%d  ny=%d\n",nx,ny);
00957                 LOGERR("compute_2d_complex only works on (nx, nx-2) images");
00958                 return;
00959         }
00960 
00961         float ds = 1.0f / (apix * ny);
00962         image->to_one();
00963 
00964         float *d = image->get_data();
00965         float g1=M_PI/2.0*cs*1.0e7*pow(lambda(),3.0f);  // s^4 coefficient for gamma, cached in a variable for simplicity (maybe speed? depends on the compiler)
00966         float g2=M_PI*lambda()*10000.0;                                 // s^2 coefficient for gamma 
00967         float acac=acos(ampcont/100.0);                                 // instead of ac*cos(g)+sqrt(1-ac^2)*sin(g), we can use cos(g-acos(ac)) and save a trig op
00968 
00969         if (type == CTF_BACKGROUND) {
00970                 for (int y = -ny/2; y < ny/2; y++) {
00971                         int y2=(y+ny)%ny;
00972                         int ynx = y2 * nx;
00973 
00974                         for (int x = 0; x < nx / 2; x++) {
00975                                 float s = (float) Util::hypot_fast(x, y ) * ds;
00976                                 d[x * 2 + ynx] = calc_noise(s);
00977                                 d[x * 2 + ynx + 1] = 0;                 // The phase is somewhat arbitrary
00978                         }
00979                 }
00980         }
00981         else if (type == CTF_AMP) {
00982                 for (int y = -ny/2; y < ny/2; y++) {
00983                         int y2=(y+ny)%ny;
00984                         int ynx = y2 * nx;
00985 
00986                         for (int x = 0; x < nx / 2; x++) {
00987                                 float s = (float)Util::hypot_fast(x,y ) * ds;
00988                                 float gam;
00989                                 if (dfdiff==0) gam=-g1*s*s*s*s+g2*defocus*s*s;
00990                                 else gam=-g1*s*s*s*s+g2*df(atan2((float)y,(float)x))*s*s;
00991                                 float v = cos(gam-acac)*exp(-(bfactor/4.0f * s*s));
00992                                 d[x * 2 + ynx] = v;
00993                                 d[x * 2 + ynx + 1] = 0;
00994                         }
00995                 }
00996         }
00997         else if (type == CTF_INTEN) {
00998                 for (int y = -ny/2; y < ny/2; y++) {
00999                         int y2=(y+ny)%ny;
01000                         int ynx = y2 * nx;
01001 
01002                         for (int x = 0; x < nx / 2; x++) {
01003                                 float s = (float)Util::hypot_fast(x,y ) * ds;
01004                                 float gam;
01005                                 if (dfdiff==0) gam=-g1*s*s*s*s+g2*defocus*s*s;
01006                                 else gam=-g1*s*s*s*s+g2*df(atan2((float)y,(float)x))*s*s;
01007                                 float v = cos(gam-acac);
01008                                 d[x * 2 + ynx] = v*v;
01009                                 d[x * 2 + ynx + 1] = 0;
01010                         }
01011                 }
01012         }
01013         else if (type == CTF_SIGN) {
01014                 for (int y = -ny/2; y < ny/2; y++) {
01015                         int y2=(y+ny)%ny;
01016                         int ynx = y2 * nx;
01017                         for (int x = 0; x < nx / 2; x++) {
01018                                 float s = (float)Util::hypot_fast(x,y ) * ds;
01019                                 float gam;
01020                                 if (dfdiff==0) gam=-g1*s*s*s*s+g2*defocus*s*s;
01021                                 else gam=-g1*s*s*s*s+g2*df(atan2((float)y,(float)x))*s*s;
01022                                 float v = cos(gam-acac);
01023                                 d[x * 2 + ynx] = v<0?-1.0:1.0;
01024                                 d[x * 2 + ynx + 1] = 0;
01025                         }
01026                 }
01027         }
01028         else if (type == CTF_FITREF) {
01029                 for (int y = -ny/2; y < ny/2; y++) {
01030                         int y2=(y+ny)%ny;
01031                         int ynx = y2 * nx;
01032 
01033                         for (int x = 0; x < nx / 2; x++) {
01034                                 float s = (float)Util::hypot_fast(x,y ) * ds;
01035                                 // We exclude very low frequencies from consideration due to the strong structure factor
01036                                 if (s<.04) {
01037                                         d[x * 2 + ynx] = 0;
01038                                         d[x * 2 + ynx + 1] = 0;
01039                                         continue;
01040                                 }
01041                                 // Rather than suddenly "turning on" the CTF, we do it gradually
01042                                 float mf=1.0;
01043                                 if (s<.05) {
01044                                         mf=1.0-exp(-pow((s-.04f)*300.0f,2.0f));
01045                                 }
01046                                 float gam;
01047                                 if (dfdiff==0) gam=-g1*s*s*s*s+g2*defocus*s*s;
01048                                 else gam=-g1*s*s*s*s+g2*df(atan2((float)y,(float)x))*s*s;
01049                                 float v = cos(gam-acac)*exp(-(bfactor/4.0f * s*s));
01050                                 d[x * 2 + ynx] = mf*v*v;
01051                                 d[x * 2 + ynx + 1] = 0;
01052                         }
01053                 }
01054         }
01055         else if (type == CTF_SNR) {
01056 
01057                 for (int y = -ny/2; y < ny/2; y++) {
01058                         int y2=(y+ny)%ny;
01059                         int ynx = y2 * nx;
01060 
01061                         for (int x = 0; x < nx / 2; x++) {
01062 
01063                                 float s = (float)Util::hypot_fast(x,y ) * ds;
01064                                 float f = s/dsbg;
01065                                 int j = (int)floor(f);
01066                                 f-=j;
01067                                 if (j>(int)snr.size()-2) d[x*2+ynx]=snr.back();
01068                                 else d[x*2+ynx]=snr[j]*(1.0f-f)+snr[j+1]*f;
01069                                 d[x * 2 + ynx + 1] = 0;
01070                         }
01071                 }
01072                 d[0]=0;
01073         }
01074         else if (type == CTF_SNR_SMOOTH) {
01075                 for (int y = -ny/2; y < ny/2; y++) {
01076                         int y2=(y+ny)%ny;
01077                         int ynx = y2 * nx;
01078 
01079                         for (int x = 0; x < nx / 2; x++) {
01080 
01081                                 float s = (float)Util::hypot_fast(x,y ) * ds;
01082                                 float f = s/dsbg;
01083                                 int j = (int)floor(f);
01084                                 f-=j;
01085                                 if (j>(int)snr.size()-2) d[x*2+ynx]=snr.back();
01086                                 else d[x*2+ynx]=snr[j]*(1.0f-f)+snr[j+1]*f;
01087                                 d[x * 2 + ynx + 1] = 0;
01088                         }
01089                 }
01090                 d[0]=0;
01091         }
01092         else if (type == CTF_WIENER_FILTER) {
01093                 if (dsbg==0) printf("Warning, DSBG set to 0\n");
01094                 for (int y = -ny/2; y < ny/2; y++) {
01095                         int y2=(y+ny)%ny;
01096                         int ynx = y2 * nx;
01097 
01098                         for (int x = 0; x < nx / 2; x++) {
01099 
01100                                 float s = (float)Util::hypot_fast(x,y ) * ds;
01101                                 float f = s/dsbg;
01102                                 int j = (int)floor(f);
01103                                 float bg,snrf;
01104                                 f-=j;
01105                                 if (j>(int)snr.size()-2) {
01106 /*                                      bg=background.back();
01107                                         d[x*2+ynx]=snr.back()/(snr.back()+1.0);*/
01108                                         d[x*2+ynx]=0;
01109                                 }
01110                                 else {
01111                                         bg=background[j]*(1.0f-f)+background[j+1]*f;
01112                                         snrf=snr[j]*(1.0f-f)+snr[j+1]*f;
01113 
01114 //                                      printf("%d\t%f\n",x,sqrt(snrf/bg)/(snrf+1.0));
01115                                         if (snrf<0) snrf=0.0;
01116 //                                      d[x*2+ynx]=sqrt(snrf/bg)/(snrf+1.0);    // Note that this is a Wiener filter with a 1/CTF term to compensate for the filtration already applied, but needs to be multiplied by the structure factor
01117                                         d[x*2+ynx]=snrf/(snrf+1);       // This is just the simple Wiener filter
01118 
01119                                 }
01120                                 d[x * 2 + ynx + 1] = 0;
01121                         }
01122                 }
01123                 d[0]=0;
01124         }
01125         else if (type == CTF_TOTAL) {
01126 
01127                 for (int y = -ny/2; y < ny/2; y++) {
01128                         int y2=(y+ny)%ny;
01129                         int ynx = y2 * nx;
01130 
01131                         for (int x = 0; x < nx / 2; x++) {
01132                                 float s = (float)Util::hypot_fast(x,y ) * ds;
01133                                 float gam;
01134                                 if (dfdiff==0) gam=-g1*s*s*s*s+g2*defocus*s*s;
01135                                 else gam=-g1*s*s*s*s+g2*df(atan2((float)y,(float)x))*s*s;
01136                                 float v = cos(gam-acac)*exp(-(bfactor/4.0f * s*s));
01137                                 d[x * 2 + ynx] = v*v+calc_noise(s);
01138                                 d[x * 2 + ynx + 1] = 0;
01139 //                              if ((x==101||x==100) && y2==79) printf("%f %f %f %f\n",s,gam,v,d[x * 2 + ynx]);
01140                         }
01141                 }
01142         }
01143         else printf("Unknown CTF image mode\n");
01144 
01145         image->update();
01146 }
01147 
01148 
01149 
01150 bool EMAN2Ctf::equal(const Ctf * ctf1) const
01151 {
01152         if (ctf1) {
01153                 const EMAN2Ctf *c = static_cast<const EMAN2Ctf *>(ctf1);
01154                 if (defocus != c->defocus ||
01155                         dfdiff != c->dfdiff ||
01156                         dfang != c->dfang ||
01157                         bfactor != c->bfactor ||
01158                         ampcont != c->ampcont ||
01159                         voltage != c->voltage ||
01160                         cs != c->cs ||
01161                         apix != c->apix ||
01162                         dsbg != c->dsbg ||
01163                         background.size() != c->background.size() ||
01164                         snr.size() != c->snr.size()
01165                         ) return false;
01166 
01167                 for (unsigned int i=0; i<background.size(); i++) {
01168                         if (background[i]!=c->background[i]) return false;
01169                 }
01170                 for (unsigned int i=0; i<snr.size(); i++) {
01171                         if (snr[i]!=c->snr[i]) return false;
01172                 }
01173                 return true;
01174         }
01175         return false;
01176 }
01177 
01178 float EMAN2Ctf::zero(int n) const {
01179 vector <float>zeroes;
01180 float lam=lambda();
01181 
01182 int m=n*2;
01183 if (m<15) m=15;         // A bit arbitrary. I believe this will always get us the zeroes we need in any practical situation
01184 for (int i=-m; i<m; i++) {
01185         float r1=defocus*defocus*1.0e8 + cs*lam*1.0e7 - 2.0*cs*i*lam*1.0e7 - 2*cs*1.0e7*lam*acos(ampcont/100.0)/M_PI;
01186 //      printf("%f\n",r1);
01187         if (r1<0) continue;
01188         float r2=defocus*1.0e4+sqrt(r1);
01189         if (r2>0) zeroes.push_back(sqrt(r2/(cs*lam*lam*1.0e7)));
01190         float r2a=defocus*1.0e4-sqrt(r1);
01191         if (r2a>0) zeroes.push_back(sqrt(r2a/(cs*lam*lam*1.0e7)));
01192 }
01193 if (zeroes.size()==0) return 0.0f;
01194 std::sort(zeroes.begin(),zeroes.end());
01195 
01196 return zeroes[n];
01197 }
01198 
01199