Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

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 
00287 
00288 void EMAN1Ctf::compute_2d_complex(EMData * image, CtfType type, XYData * sf)
00289 {
00290         if (!image) {
00291                 LOGERR("image is null. cannot computer 2D complex CTF");
00292                 return;
00293         }
00294 
00295         if (image->is_complex() == false) {
00296                 LOGERR("compute_2d_complex can only work on complex images");
00297                 return;
00298         }
00299 
00300         int nx = image->get_xsize();
00301         int ny = image->get_ysize();
00302 
00303         if (nx != ny + 2) {
00304                 LOGERR("compute_2d_complex only works on (nx, nx-2) images");
00305                 return;
00306         }
00307 
00308         float ds = 1.0f / (apix * ny);
00309         image->to_one();
00310 
00311         float *d = image->get_data();
00312         float g1 = calc_g1();
00313         float g2 = calc_g2();
00314 
00315         if (type == CTF_BACKGROUND) {
00316                 for (int y = 0; y < ny; y++) {
00317                         int ynx = y * nx;
00318 
00319                         for (int x = 0; x < nx / 2; x++) {
00320 #ifdef  _WIN32
00321                                 float s = (float) _hypot(x, y - ny / 2.0f) * ds;
00322 #else
00323                                 float s = (float) hypot(x, y - ny / 2.0f) * ds;
00324 #endif
00325                                 d[x * 2 + ynx] = calc_noise(s);
00326                                 d[x * 2 + ynx + 1] = 0;                 // The phase is somewhat arbitrary
00327                         }
00328                 }
00329         }
00330         else if (type == CTF_AMP) {
00331                 for (int y = 0; y < ny; y++) {
00332                         int ynx = y * nx;
00333 
00334                         for (int x = 0; x < nx / 2; x++) {
00335 #ifdef  _WIN32
00336                                 float s = (float)_hypot((float) x, (float) y - ny / 2) * ds;
00337 #else
00338                                 float s = (float)hypot((float) x, (float) y - ny / 2) * ds;
00339 #endif  //_WIN32
00340                                 float gamma = calc_gamma(g1, g2, s);
00341                                 float v = fabs(calc_amplitude(gamma));
00342                                 d[x * 2 + ynx] = v;
00343                                 d[x * 2 + ynx + 1] = 0;
00344                         }
00345                 }
00346         }
00347         else if (type == CTF_SIGN) {
00348                 for (int y = 0; y < ny; y++) {
00349                         int ynx = y * nx;
00350                         for (int x = 0; x < nx / 2; x++) {
00351 #ifdef  _WIN32
00352                                 float s = (float)_hypot(x, y - ny / 2.0f) * ds;
00353 #else
00354                                 float s = (float)hypot(x, y - ny / 2.0f) * ds;
00355 #endif
00356                                 float gamma = calc_gamma(g1, g2, s);
00357                                 float v = calc_amplitude(gamma);
00358                                 d[x * 2 + ynx] = v > 0 ? 1.0f : -1.0f;
00359                                 d[x * 2 + ynx + 1] = 0;
00360                         }
00361                 }
00362 
00363         }
00364         else if (type == CTF_SNR || type == CTF_SNR_SMOOTH) {
00365                 float amp1 = calc_amp1();
00366 
00367                 for (int y = 0; y < ny; y++) {
00368                         int ynx = y * nx;
00369 
00370                         for (int x = 0; x < nx / 2; x++) {
00371 
00372 #ifdef  _WIN32
00373                                 float s = (float)_hypot(x, y - ny / 2.0f) * ds;
00374 #else
00375                                 float s = (float)hypot(x, y - ny / 2.0f) * ds;
00376 #endif
00377                                 float gamma = calc_gamma(g1, g2, s);
00378                                 float f = calc_ctf1(amp1, gamma, s);
00379                                 float noise = calc_noise(s);
00380                                 f = f * f / noise;
00381 
00382                                 if (s && sf) {
00383                                         f *= sf->get_yatx(s);
00384                                 }
00385                                 d[x * 2 + ynx] *= f;
00386                                 d[x * 2 + ynx + 1] = 0;
00387                         }
00388                 }
00389         }
00390         else if (type == CTF_WIENER_FILTER) {
00391                 float amp1 = calc_amp1();
00392 
00393                 for (int y = 0; y < ny; y++) {
00394                         int ynx = y * nx;
00395 
00396                         for (int x = 0; x < nx / 2; x++) {
00397 
00398 #ifdef  _WIN32
00399                                 float s = (float)_hypot(x, y - ny / 2.0f) * ds;
00400 #else
00401                                 float s = (float)hypot(x, y - ny / 2.0f) * ds;
00402 #endif
00403                                 float gamma = calc_gamma(g1, g2, s);
00404                                 float f = calc_ctf1(amp1, gamma, s);
00405                                 float noise = calc_noise(s);
00406                                 f = f * f / noise;
00407 
00408                                 if (s) {
00409                                         f *= sf->get_yatx(s);
00410                                 }
00411                                 f = 1.0f / (1.0f + 1.0f / f);
00412                                 d[x * 2 + ynx] *= f;
00413                                 d[x * 2 + ynx + 1] = 0;
00414                         }
00415                 }
00416         }
00417         else if (type == CTF_TOTAL) {
00418                 float amp1 = calc_amp1();
00419 
00420                 for (int y = 0; y < ny; y++) {
00421                         int ynx = y * nx;
00422 
00423                         for (int x = 0; x < nx / 2; x++) {
00424 
00425 #ifdef  _WIN32
00426                                 float s = (float)_hypot(x, y - ny / 2.0f) * ds;
00427 #else
00428                                 float s = (float)hypot(x, y - ny / 2.0f) * ds;
00429 #endif
00430                                 float gamma = calc_gamma(g1, g2, s);
00431                                 float f = calc_ctf1(amp1, gamma, s);
00432                                 float noise = calc_noise(s);
00433                                 f = f * f;
00434 
00435                                 if (sf && s) {
00436                                         f *= sf->get_yatx(s);
00437                                 }
00438                                 f+=noise;
00439 
00440                                 d[x * 2 + ynx] *= f;
00441                                 d[x * 2 + ynx + 1] = 0;
00442                         }
00443                 }
00444         }
00445 
00446         image->update();
00447 }
00448 
00449 
00450 
00451 bool EMAN1Ctf::equal(const Ctf * ctf1) const
00452 {
00453         if (ctf1) {
00454                 const EMAN1Ctf *c = static_cast<const EMAN1Ctf *>(ctf1);
00455                 if (defocus == c->defocus &&
00456                         bfactor == c->bfactor &&
00457                         amplitude == c->amplitude &&
00458                         ampcont == c->ampcont &&
00459                         noise1 == c->noise1 &&
00460                         noise2 == c->noise2 &&
00461                         noise3 == c->noise3 &&
00462                         noise4 == c->noise4 && voltage == c->voltage && cs == c->cs && apix == c->apix) {
00463                         return true;
00464                 }
00465         }
00466         return false;
00467 }
00468 
00469 /*************************************
00470 EMAN2Ctf
00471 *************************************/
00472 
00473 EMAN2Ctf::EMAN2Ctf()
00474 {
00475         defocus = 0;
00476         dfdiff = 0;
00477         dfang = 0;
00478         bfactor = 0;
00479         ampcont = 0;
00480         voltage = 0;
00481         cs = 0;
00482         apix = 1.0;
00483         dsbg=-1;
00484         background.clear();
00485         snr.clear();
00486 }
00487 
00488 
00489 EMAN2Ctf::~EMAN2Ctf()
00490 {
00491 }
00492 
00493 
00494 int EMAN2Ctf::from_string(const string & ctf)
00495 {
00496         Assert(ctf != "");
00497         char type=' ';
00498         int pos=-1,i,j;
00499         int bglen=0,snrlen=0;
00500         float v;
00501         const char *s=ctf.c_str();
00502 
00503         sscanf(s, "%c%f %f %f %f %f %f %f %f %f %d%n",
00504                                    &type,&defocus, &dfdiff,&dfang,&bfactor,&ampcont,&voltage, &cs, &apix,&dsbg,&bglen,&pos);
00505         if (type!='E') throw InvalidValueException(type,"Trying to initialize Ctf object with bad string");
00506         if (pos==-1) throw InvalidValueException(s," Invalid CTF string");
00507         
00508 
00509         background.resize(bglen);
00510         for (i=0; i<bglen; i++) {
00511                 if (sscanf(s+pos,",%f%n",&v,&j)<1) return(1);
00512                 background[i]=v;
00513                 pos+=j;
00514         }
00515 
00516         sscanf(s+pos," %d%n",&snrlen,&j);
00517         pos+=j;
00518         snr.resize(snrlen);
00519         for (i=0; i<snrlen; i++) {
00520                 if (sscanf(s+pos,",%f%n",&v,&j)<1) return(1);
00521                 snr[i]=v;
00522                 pos+=j;
00523         }
00524 
00525         return 0;
00526 
00527 }
00528 
00529 string EMAN2Ctf::to_string() const
00530 {
00531         char ctf[256];
00532         sprintf(ctf, "E%1.4g %1.4g %1.4g %1.4g %1.4g %1.4g %1.4g %1.4g %1.4g %d",
00533                         defocus, dfdiff, dfang, bfactor, ampcont, voltage, cs, apix, dsbg,(int)background.size());
00534 
00535         string ret=ctf;
00536         for (int i=0; i<(int)background.size(); i++) {
00537                 sprintf(ctf,",%1.3g",background[i]);
00538                 ret+=ctf;
00539         }
00540 
00541         sprintf(ctf, " %d",(int)snr.size());
00542         ret+=ctf;
00543         for (int i=0; i<(int)snr.size(); i++) {
00544                 sprintf(ctf,",%1.3g",snr[i]);
00545                 ret+=ctf;
00546         }
00547 
00548 
00549         return ret;
00550 }
00551 
00552 void EMAN2Ctf::from_dict(const Dict & dict)
00553 {
00554         defocus = (float)dict["defocus"];
00555         dfdiff = (float)dict["dfdiff"];
00556         dfang = (float)dict["dfang"];
00557         bfactor = (float)dict["bfactor"];
00558         ampcont = (float)dict["ampcont"];
00559         voltage = (float)dict["voltage"];
00560         cs = (float)dict["cs"];
00561         apix = (float)dict["apix"];
00562         dsbg = (float)dict["dsbg"];
00563         background = dict["background"];
00564         snr = dict["snr"];
00565 }
00566 
00567 Dict EMAN2Ctf::to_dict() const
00568 {
00569         Dict dict;
00570         dict["defocus"] = defocus;
00571         dict["dfdiff"] = dfdiff;
00572         dict["dfang"] = dfang;
00573         dict["bfactor"] = bfactor;
00574         dict["ampcont"] = ampcont;
00575         dict["voltage"] = voltage;
00576         dict["cs"] = cs;
00577         dict["apix"] = apix;
00578         dict["dsbg"] = dsbg;
00579         dict["background"] = background;
00580         dict["snr"] = snr;
00581 
00582         return dict;
00583 }
00584 
00585 void EMAN2Ctf::from_vector(const vector<float>& vctf)
00586 {
00587         int i;
00588         defocus = vctf[0];
00589         dfdiff = vctf[1];
00590         dfang = vctf[2];
00591         bfactor = vctf[3];
00592         ampcont = vctf[4];
00593         voltage = vctf[5];
00594         cs = vctf[6];
00595         apix = vctf[7];
00596         dsbg = vctf[8];
00597         background.resize((int)vctf[9]);
00598         for (i=0; i<(int)vctf[9]; i++) background[i]=vctf[i+10];
00599         snr.resize((int)vctf[i+10]);
00600         for (int j=0; j<(int)vctf[i+10]; j++) snr[j]=vctf[i+j+11];
00601 }
00602 
00603 vector<float> EMAN2Ctf::to_vector() const
00604 {
00605         vector<float> vctf;
00606 
00607         vctf.push_back(defocus);
00608         vctf.push_back(dfdiff);
00609         vctf.push_back(dfang);
00610         vctf.push_back(bfactor);
00611         vctf.push_back(ampcont);
00612         vctf.push_back(voltage);
00613         vctf.push_back(cs);
00614         vctf.push_back(apix);
00615         vctf.push_back(dsbg);
00616         vctf.push_back((float)background.size());
00617         for (unsigned int i=0; i<background.size(); i++) vctf.push_back(background[i]);
00618         vctf.push_back((float)snr.size());
00619         for (unsigned int j=0; j<snr.size(); j++) vctf.push_back(snr[j]);
00620 
00621         return vctf;
00622 }
00623 
00624 
00625 
00626 void EMAN2Ctf::copy_from(const Ctf * new_ctf)
00627 {
00628         if (new_ctf) {
00629                 const EMAN2Ctf *c = static_cast<const EMAN2Ctf *>(new_ctf);
00630                 defocus = c->defocus;
00631                 dfdiff = c->dfdiff;
00632                 dfang = c->dfang;
00633                 bfactor = c->bfactor;
00634                 ampcont = c->ampcont;
00635                 voltage = c->voltage;
00636                 cs = c->cs;
00637                 apix = c->apix;
00638                 dsbg = c->dsbg;
00639                 background = c->background;
00640                 snr = c->snr;
00641         }
00642 }
00643 
00644 inline int max_int(int a,int b) { return a>b?a:b; }
00645 inline int min_int(int a,int b) { return a<b?a:b; }
00646 
00647 vector < float >EMAN2Ctf::compute_1d(int size,float ds, CtfType type, XYData * sf)
00648 {
00649         Assert(size > 0);
00650 
00651 //      float tmp_f1 =  sqrt((float) 2) * size / 2;
00652 //      int np = (int) ceil(tmp_f1) + 2;
00653         int np=size/2;
00654         vector < float >r;
00655 
00656         r.resize(np);
00657 
00658 //      float ds = 1 / (apix * size);
00659         float s = 0;
00660         float g1 = calc_g1();
00661         float g2 = calc_g2();
00662         float amp1 = calc_amp1();
00663 
00664         switch (type) {
00665         case CTF_AMP:
00666                 for (int i = 0; i < np; i++) {
00667                         float gamma = calc_gamma(g1, g2, s);
00668                         r[i] = calc_ctf1(amp1, gamma, s);
00669                         s += ds;
00670                 }
00671                 break;
00672 
00673         case CTF_SIGN:
00674                 for (int i = 0; i < np; i++) {
00675                         float gamma = calc_gamma(g1, g2, s);
00676                         r[i] = calc_ctf1(amp1, gamma, s)>=0?1.0f:-1.0f;
00677                         s += ds;
00678                 }
00679                 break;
00680 
00681         case CTF_BACKGROUND:
00682                 for (int i = 0; i < np; i++) {
00683                         float f = s/dsbg;
00684                         int j = (int)floor(f);
00685                         f-=j;
00686                         if (j>(int)background.size()-2) r[i]=background.back();
00687                         else r[i]=background[j]*(1.0f-f)+background[j+1]*f;
00688                         s+=ds;
00689                 }
00690                 break;
00691 
00692         case CTF_SNR:
00693                 if (snr.size()<1) {
00694                         throw InvalidValueException("CTF_SNR"," ERROR: No SNR info in CTF header.");
00695                         break;
00696                 }
00697                 for (int i = 0; i < np; i++) {
00698                         float f = s/dsbg;
00699                         int j = (int)floor(f);
00700                         f-=j;
00701                         if (j>(int)snr.size()-2) r[i]=snr.back();
00702                         else r[i]=snr[j]*(1.0f-f)+snr[j+1]*f;
00703 //                      printf("%d\t%f\n",j,snr[j]);
00704                         s+=ds;
00705                 }
00706                 break;
00707         case CTF_SNR_SMOOTH:
00708                 // This apparently complicated routine tries to make a nice smooth and accurate SNR curve. It does this
00709                 // by fitting local regions of the SNR vs the theoretical SNR (theoretical CTF^2/measured background),
00710                 // then taking the slope of the result times the theoretical SNR to produce a local SNR estimate
00711 
00712                 { // <- is to permit new temporary value allocation
00713                         vector < float >tsnr;   // theoretical SNR
00714                         tsnr.resize(np);
00715                         vector < float >dsnr;   // data SNR
00716                         dsnr.resize(np);
00717                         
00718                         float s0=s;
00719                         
00720                         for (int i = 0; i < np; i++) {
00721                                 float gamma = calc_gamma(g1, g2, s);
00722                                 tsnr[i] = calc_ctf1(amp1, gamma, s);            // ctf amp
00723 
00724                                 // background value
00725                                 float f = s/dsbg;
00726                                 int j = (int)floor(f);
00727                                 f-=j;
00728                                 float bg;
00729                                 if (j>(int)background.size()-2) bg=background.back();
00730                                 else bg=background[j]*(1.0f-f)+background[j+1]*f;
00731                                 if (bg <=0) bg=.001f;
00732 
00733                                 tsnr[i] = tsnr[i]*tsnr[i]/bg;           // This is now a SNR curve
00734                                 if (sf && s) {
00735                                         tsnr[i] *= sf->get_yatx(s);
00736                                 }
00737 
00738                                 
00739                                 // This is the SNR computed from the data without fitting
00740                                 if (j>(int)snr.size()-2) dsnr[i]=snr.back();
00741                                 else dsnr[i]=snr[j]*(1.0f-f)+snr[j+1]*f;
00742                                 
00743                                 s+=ds;
00744                         }
00745 
00746                         int npsm=np/25;                 // 1/2 number of points to smooth over, 25 is arbitrary
00747                         if (npsm<2) npsm=2;
00748                         
00749                         s=s0;
00750                         for (int i = 1; i < np; i++) {
00751                                 // simple linear regression embedded here
00752                                 double sum = 0;
00753                                 double sum_x = 0;
00754                                 double sum_y = 0;
00755                                 double sum_xx = 0;
00756                                 double sum_xy = 0;
00757 
00758                                 for (int k=max_int(i-npsm,1); k<=min_int(i+npsm,np-1); k++) {
00759                                         double y = dsnr[k];
00760                                         double x = tsnr[k];
00761 
00762                                         sum_x += x;
00763                                         sum_y += y;
00764                                         sum_xx += x * x;
00765                                         sum_xy += x * y;
00766                                         sum++;
00767                                 }
00768 
00769                                 double div = sum * sum_xx - sum_x * sum_x;
00770 //                              if (div == 0) {
00771 //                                      div = 0.0000001f;
00772 //                              }
00773 
00774         //                      *intercept = (float) ((sum_xx * sum_y - sum_x * sum_xy) / div);
00775         //                      *slope = (float) ((sum * sum_xy - sum_x * sum_y) / div);
00776 
00777                                 if (div!=0.0) r[i]=(float) ((sum * sum_xy - sum_x * sum_y) / div)*tsnr[i];
00778                                 else r[i]=0.0;
00779                                 if (r[i]<0) r[i]=0;
00780                                 
00781                                 s+=ds;
00782                         }
00783                         r[0]=0;
00784                 }
00785                 break;
00786 
00787         case CTF_WIENER_FILTER:
00788 //              if (!sf) {
00789 //                      LOGERR("CTF computation error, no SF found\n");
00790 //                      return r;
00791 //              }
00792 
00793                 for (int i = 0; i < np; i++) {
00794                         float f = s/dsbg;
00795                         int j = (int)floor(f);
00796                         float bg;
00797                         f-=j;
00798                         if (j>(int)snr.size()-2) {
00799 /*                              r[i]=snr.back();
00800                                 bg=background.back();*/
00801                                 r[i]=0;
00802                         }
00803                         else {
00804                                 r[i]=snr[j]*(1.0f-f)+snr[j+1]*f;
00805                                 bg=background[j]*(1.0f-f)+background[j+1]*f;
00806                         }
00807                         if (r[i]<0) r[i]=0;
00808                         r[i]=r[i]/(r[i]+1.0f);          // Full Wiener filter assuming perfect image with noise
00809 //                      r[i]=sqrt(r[i]/bg)/(r[i]+1.0);  // Wiener filter with 1/CTF term (sort of) to restore image then filter
00810                         s+=ds;
00811                 }
00812                 r[0]=0;
00813                 break;
00814 
00815         case CTF_TOTAL:
00816 
00817                 for (int i = 0; i < np; i++) {
00818                         float gamma = calc_gamma(g1, g2, s);
00819                         if (sf) {
00820                                 r[i] = calc_ctf1(amp1, gamma, s);
00821                                 r[i] = r[i] * r[i] * sf->get_yatx(s) + calc_noise(s);
00822                         }
00823                         else {
00824                                 r[i] = calc_ctf1(amp1, gamma, s);
00825                                 r[i] = r[i] * r[i] + calc_noise(s);
00826                         }
00827                         s += ds;
00828                 }
00829                 break;
00830         default:
00831                 break;
00832         }
00833 
00834         return r;
00835 }
00836 
00837 
00838 void EMAN2Ctf::compute_2d_real(EMData *, CtfType, XYData *)
00839 {
00840 
00841 
00842 }
00843 
00844 
00845 
00846 void EMAN2Ctf::compute_2d_complex(EMData * image, CtfType type, XYData * sf)
00847 {
00848         if (!image) {
00849                 LOGERR("image is null. cannot computer 2D complex CTF");
00850                 return;
00851         }
00852 
00853         if (image->is_complex() == false) {
00854                 LOGERR("compute_2d_complex can only work on complex images");
00855                 return;
00856         }
00857 
00858         int nx = image->get_xsize();
00859         int ny = image->get_ysize();
00860 
00861         if ((ny%2==1 && nx!=ny+1) || (ny%2==0 && nx != ny + 2)) {
00862                 printf("nx=%d  ny=%d\n",nx,ny);
00863                 LOGERR("compute_2d_complex only works on (nx, nx-2) images");
00864                 return;
00865         }
00866 
00867         float ds = 1.0f / (apix * ny);
00868         image->to_one();
00869 
00870         float *d = image->get_data();
00871         float g1 = calc_g1();
00872         float g2 = calc_g2();
00873         float amp1 = calc_amp1();
00874 
00875         if (type == CTF_BACKGROUND) {
00876                 for (int y = -ny/2; y < ny/2; y++) {
00877                         int y2=(y+ny)%ny;
00878                         int ynx = y2 * nx;
00879 
00880                         for (int x = 0; x < nx / 2; x++) {
00881                                 float s = (float) Util::hypot_fast(x, y ) * ds;
00882                                 d[x * 2 + ynx] = calc_noise(s);
00883                                 d[x * 2 + ynx + 1] = 0;                 // The phase is somewhat arbitrary
00884                         }
00885                 }
00886         }
00887         else if (type == CTF_AMP) {
00888                 for (int y = -ny/2; y < ny/2; y++) {
00889                         int y2=(y+ny)%ny;
00890                         int ynx = y2 * nx;
00891 
00892                         for (int x = 0; x < nx / 2; x++) {
00893                                 float s = (float)Util::hypot_fast(x,y ) * ds;
00894                                 float gamma = calc_gamma(g1, g2, s);
00895                                 float v = calc_ctf1(amp1, gamma, s);
00896 //                              float v = calc_amplitude(gamma);
00897                                 d[x * 2 + ynx] = v;
00898                                 d[x * 2 + ynx + 1] = 0;
00899                         }
00900                 }
00901         }
00902         else if (type == CTF_SIGN) {
00903                 for (int y = -ny/2; y < ny/2; y++) {
00904                         int y2=(y+ny)%ny;
00905                         int ynx = y2 * nx;
00906                         for (int x = 0; x < nx / 2; x++) {
00907                                 float s = (float)Util::hypot_fast(x,y ) * ds;
00908                                 float gamma = calc_gamma(g1, g2, s);
00909                                 float v = calc_amplitude(gamma);
00910                                 d[x * 2 + ynx] = v >= 0 ? 1.0f : -1.0f;
00911                                 d[x * 2 + ynx + 1] = 0;
00912                         }
00913                 }
00914         }
00915         else if (type == CTF_SNR) {
00916 
00917                 for (int y = -ny/2; y < ny/2; y++) {
00918                         int y2=(y+ny)%ny;
00919                         int ynx = y2 * nx;
00920 
00921                         for (int x = 0; x < nx / 2; x++) {
00922 
00923                                 float s = (float)Util::hypot_fast(x,y ) * ds;
00924                                 float f = s/dsbg;
00925                                 int j = (int)floor(f);
00926                                 f-=j;
00927                                 if (j>(int)snr.size()-2) d[x*2+ynx]=snr.back();
00928                                 else d[x*2+ynx]=snr[j]*(1.0f-f)+snr[j+1]*f;
00929                                 d[x * 2 + ynx + 1] = 0;
00930                         }
00931                 }
00932                 d[0]=0;
00933         }
00934         else if (type == CTF_SNR_SMOOTH) {
00935                 for (int y = -ny/2; y < ny/2; y++) {
00936                         int y2=(y+ny)%ny;
00937                         int ynx = y2 * nx;
00938 
00939                         for (int x = 0; x < nx / 2; x++) {
00940 
00941                                 float s = (float)Util::hypot_fast(x,y ) * ds;
00942                                 float f = s/dsbg;
00943                                 int j = (int)floor(f);
00944                                 f-=j;
00945                                 if (j>(int)snr.size()-2) d[x*2+ynx]=snr.back();
00946                                 else d[x*2+ynx]=snr[j]*(1.0f-f)+snr[j+1]*f;
00947                                 d[x * 2 + ynx + 1] = 0;
00948                         }
00949                 }
00950                 d[0]=0;
00951         }
00952         else if (type == CTF_WIENER_FILTER) {
00953                 if (dsbg==0) printf("Warning, DSBG set to 0\n");
00954                 for (int y = -ny/2; y < ny/2; y++) {
00955                         int y2=(y+ny)%ny;
00956                         int ynx = y2 * nx;
00957 
00958                         for (int x = 0; x < nx / 2; x++) {
00959 
00960                                 float s = (float)Util::hypot_fast(x,y ) * ds;
00961                                 float f = s/dsbg;
00962                                 int j = (int)floor(f);
00963                                 float bg,snrf;
00964                                 f-=j;
00965                                 if (j>(int)snr.size()-2) {
00966 /*                                      bg=background.back();
00967                                         d[x*2+ynx]=snr.back()/(snr.back()+1.0);*/
00968                                         d[x*2+ynx]=0;
00969                                 }
00970                                 else {
00971                                         bg=background[j]*(1.0f-f)+background[j+1]*f;
00972                                         snrf=snr[j]*(1.0f-f)+snr[j+1]*f;
00973 
00974 //                                      printf("%d\t%f\n",x,sqrt(snrf/bg)/(snrf+1.0));
00975                                         if (snrf<0) snrf=0.0;
00976 //                                      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
00977                                         d[x*2+ynx]=snrf/(snrf+1);       // This is just the simple Wiener filter
00978 
00979                                 }
00980                                 d[x * 2 + ynx + 1] = 0;
00981                         }
00982                 }
00983                 d[0]=0;
00984         }
00985         else if (type == CTF_TOTAL) {
00986                 float amp1 = calc_amp1();
00987 
00988                 for (int y = -ny/2; y < ny/2; y++) {
00989                         int y2=(y+ny)%ny;
00990                         int ynx = y2 * nx;
00991 
00992                         for (int x = 0; x < nx / 2; x++) {
00993 
00994                                 float s = (float)Util::hypot_fast(x,y ) * ds;
00995                                 float gamma = calc_gamma(g1, g2, s);
00996                                 float f = calc_ctf1(amp1, gamma, s);
00997                                 float noise = 0;
00998                                 f = f * f;
00999 
01000                                 if (sf && s) {
01001                                         f *= sf->get_yatx(s);
01002                                 }
01003                                 f+=noise;
01004 
01005                                 d[x * 2 + ynx] *= f;
01006                                 d[x * 2 + ynx + 1] = 0;
01007                         }
01008                 }
01009         }
01010         else printf("Unknown CTF image mode\n");
01011 
01012         image->update();
01013 }
01014 
01015 
01016 
01017 bool EMAN2Ctf::equal(const Ctf * ctf1) const
01018 {
01019         if (ctf1) {
01020                 const EMAN2Ctf *c = static_cast<const EMAN2Ctf *>(ctf1);
01021                 if (defocus != c->defocus ||
01022                         dfdiff != c->dfdiff ||
01023                         dfang != c->dfang ||
01024                         bfactor != c->bfactor ||
01025                         ampcont != c->ampcont ||
01026                         voltage != c->voltage ||
01027                         cs != c->cs ||
01028                         apix != c->apix ||
01029                         dsbg != c->dsbg ||
01030                         background.size() != c->background.size() ||
01031                         snr.size() != c->snr.size()
01032                         ) return false;
01033 
01034                 for (unsigned int i=0; i<background.size(); i++) {
01035                         if (background[i]!=c->background[i]) return false;
01036                 }
01037                 for (unsigned int i=0; i<snr.size(); i++) {
01038                         if (snr[i]!=c->snr[i]) return false;
01039                 }
01040                 return true;
01041         }
01042         return false;
01043 }

Generated on Tue Jun 11 13:40:37 2013 for EMAN2 by  doxygen 1.3.9.1