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 i = sscanf(ctf.c_str(), "%c%f %f %f %f %f %f %f %f %f %f %f",
00069                                    &type,&defocus, &bfactor, &amplitude, &ampcont, &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 //      float ds = 1 / (apix * size * CTFOS);
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 //              if (!sf) {
00217 //                      LOGERR("CTF computation error, no SF found\n");
00218 //                      return r;
00219 //              }
00220 
00221                 for (int i = 0; i < np; i++) {
00222                         float gamma = calc_gamma(g1, g2, s);
00223                         r[i] = calc_snr(amp1, gamma, s);
00224                         if (s && sf) {
00225                                 r[i] *= sf->get_yatx(s);
00226                         }
00227                         s += ds;
00228                 }
00229 
00230                 break;
00231 
00232         case CTF_WIENER_FILTER:
00233                 if (!sf) {
00234                         LOGERR("CTF computation error, no SF found\n");
00235                         return r;
00236                 }
00237 
00238                 for (int i = 0; i < np; i++) {
00239                         float gamma = calc_gamma(g1, g2, s);
00240                         r[i] = calc_snr(amp1, gamma, s);
00241                         if (s && sf) {
00242                                 r[i] *= sf->get_yatx(s);
00243                         }
00244 
00245                         r[i] = 1.0f / (1.0f + 1.0f / r[i]);
00246                         s += ds;
00247                 }
00248                 break;
00249 
00250         case CTF_TOTAL:
00251                 if (!sf) {
00252                         LOGERR("CTF computation error, no SF found\n");
00253                         return r;
00254                 }
00255 
00256                 for (int i = 0; i < np; i++) {
00257                         float gamma = calc_gamma(g1, g2, s);
00258                         if (sf) {
00259                                 r[i] = calc_ctf1(amp1, gamma, s);
00260                                 r[i] = r[i] * r[i] * sf->get_yatx(s) + calc_noise(s);
00261                         }
00262                         else {
00263                                 r[i] = calc_ctf1(amp1, gamma, s);
00264                                 r[i] = r[i] * r[i] + calc_noise(s);
00265                         }
00266                         s += ds;
00267                 }
00268                 break;
00269         default:
00270                 break;
00271         }
00272 
00273         return r;
00274 }
00275 
00276 
00277 void EMAN1Ctf::compute_2d_real(EMData *, CtfType, XYData *)
00278 {
00279 
00280 
00281 }
00282 
00283 
00284 
00285 void EMAN1Ctf::compute_2d_complex(EMData * image, CtfType type, XYData * sf)
00286 {
00287         if (!image) {
00288                 LOGERR("image is null. cannot computer 2D complex CTF");
00289                 return;
00290         }
00291 
00292         if (image->is_complex() == false) {
00293                 LOGERR("compute_2d_complex can only work on complex images");
00294                 return;
00295         }
00296 
00297         int nx = image->get_xsize();
00298         int ny = image->get_ysize();
00299 
00300         if (nx != ny + 2) {
00301                 LOGERR("compute_2d_complex only works on (nx, nx-2) images");
00302                 return;
00303         }
00304 
00305         float ds = 1.0f / (apix * ny);
00306         image->to_one();
00307 
00308         float *d = image->get_data();
00309         float g1 = calc_g1();
00310         float g2 = calc_g2();
00311 
00312         if (type == CTF_BACKGROUND) {
00313                 for (int y = 0; y < ny; y++) {
00314                         int ynx = y * nx;
00315 
00316                         for (int x = 0; x < nx / 2; x++) {
00317 #ifdef  _WIN32
00318                                 float s = (float) _hypot(x, y - ny / 2.0f) * ds;
00319 #else
00320                                 float s = (float) hypot(x, y - ny / 2.0f) * ds;
00321 #endif
00322                                 d[x * 2 + ynx] = calc_noise(s);
00323                                 d[x * 2 + ynx + 1] = 0;                 // The phase is somewhat arbitrary
00324                         }
00325                 }
00326         }
00327         else if (type == CTF_AMP) {
00328                 for (int y = 0; y < ny; y++) {
00329                         int ynx = y * nx;
00330 
00331                         for (int x = 0; x < nx / 2; x++) {
00332 #ifdef  _WIN32
00333                                 float s = (float)_hypot((float) x, (float) y - ny / 2) * ds;
00334 #else
00335                                 float s = (float)hypot((float) x, (float) y - ny / 2) * ds;
00336 #endif  //_WIN32
00337                                 float gamma = calc_gamma(g1, g2, s);
00338                                 float v = fabs(calc_amplitude(gamma));
00339                                 d[x * 2 + ynx] = v;
00340                                 d[x * 2 + ynx + 1] = 0;
00341                         }
00342                 }
00343         }
00344         else if (type == CTF_SIGN) {
00345                 for (int y = 0; y < ny; y++) {
00346                         int ynx = y * nx;
00347                         for (int x = 0; x < nx / 2; x++) {
00348 #ifdef  _WIN32
00349                                 float s = (float)_hypot(x, y - ny / 2.0f) * ds;
00350 #else
00351                                 float s = (float)hypot(x, y - ny / 2.0f) * ds;
00352 #endif
00353                                 float gamma = calc_gamma(g1, g2, s);
00354                                 float v = calc_amplitude(gamma);
00355                                 d[x * 2 + ynx] = v > 0 ? 1.0f : -1.0f;
00356                                 d[x * 2 + ynx + 1] = 0;
00357                         }
00358                 }
00359 
00360         }
00361         else if (type == CTF_SNR || type == CTF_SNR_SMOOTH) {
00362                 float amp1 = calc_amp1();
00363 
00364                 for (int y = 0; y < ny; y++) {
00365                         int ynx = y * nx;
00366 
00367                         for (int x = 0; x < nx / 2; x++) {
00368 
00369 #ifdef  _WIN32
00370                                 float s = (float)_hypot(x, y - ny / 2.0f) * ds;
00371 #else
00372                                 float s = (float)hypot(x, y - ny / 2.0f) * ds;
00373 #endif
00374                                 float gamma = calc_gamma(g1, g2, s);
00375                                 float f = calc_ctf1(amp1, gamma, s);
00376                                 float noise = calc_noise(s);
00377                                 f = f * f / noise;
00378 
00379                                 if (s && sf) {
00380                                         f *= sf->get_yatx(s);
00381                                 }
00382                                 d[x * 2 + ynx] *= f;
00383                                 d[x * 2 + ynx + 1] = 0;
00384                         }
00385                 }
00386         }
00387         else if (type == CTF_WIENER_FILTER) {
00388                 float amp1 = calc_amp1();
00389 
00390                 for (int y = 0; y < ny; y++) {
00391                         int ynx = y * nx;
00392 
00393                         for (int x = 0; x < nx / 2; x++) {
00394 
00395 #ifdef  _WIN32
00396                                 float s = (float)_hypot(x, y - ny / 2.0f) * ds;
00397 #else
00398                                 float s = (float)hypot(x, y - ny / 2.0f) * ds;
00399 #endif
00400                                 float gamma = calc_gamma(g1, g2, s);
00401                                 float f = calc_ctf1(amp1, gamma, s);
00402                                 float noise = calc_noise(s);
00403                                 f = f * f / noise;
00404 
00405                                 if (s) {
00406                                         f *= sf->get_yatx(s);
00407                                 }
00408                                 f = 1.0f / (1.0f + 1.0f / f);
00409                                 d[x * 2 + ynx] *= f;
00410                                 d[x * 2 + ynx + 1] = 0;
00411                         }
00412                 }
00413         }
00414         else if (type == CTF_TOTAL) {
00415                 float amp1 = calc_amp1();
00416 
00417                 for (int y = 0; y < ny; y++) {
00418                         int ynx = y * nx;
00419 
00420                         for (int x = 0; x < nx / 2; x++) {
00421 
00422 #ifdef  _WIN32
00423                                 float s = (float)_hypot(x, y - ny / 2.0f) * ds;
00424 #else
00425                                 float s = (float)hypot(x, y - ny / 2.0f) * ds;
00426 #endif
00427                                 float gamma = calc_gamma(g1, g2, s);
00428                                 float f = calc_ctf1(amp1, gamma, s);
00429                                 float noise = calc_noise(s);
00430                                 f = f * f;
00431 
00432                                 if (sf && s) {
00433                                         f *= sf->get_yatx(s);
00434                                 }
00435                                 f+=noise;
00436 
00437                                 d[x * 2 + ynx] *= f;
00438                                 d[x * 2 + ynx + 1] = 0;
00439                         }
00440                 }
00441         }
00442 
00443         image->update();
00444 }
00445 
00446 
00447 
00448 bool EMAN1Ctf::equal(const Ctf * ctf1) const
00449 {
00450         if (ctf1) {
00451                 const EMAN1Ctf *c = static_cast<const EMAN1Ctf *>(ctf1);
00452                 if (defocus == c->defocus &&
00453                         bfactor == c->bfactor &&
00454                         amplitude == c->amplitude &&
00455                         ampcont == c->ampcont &&
00456                         noise1 == c->noise1 &&
00457                         noise2 == c->noise2 &&
00458                         noise3 == c->noise3 &&
00459                         noise4 == c->noise4 && voltage == c->voltage && cs == c->cs && apix == c->apix) {
00460                         return true;
00461                 }
00462         }
00463         return false;
00464 }
00465 
00466 /*************************************
00467 EMAN2Ctf
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,&ampcont,&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 //      float tmp_f1 =  sqrt((float) 2) * size / 2;
00648 //      int np = (int) ceil(tmp_f1) + 2;
00649         int np=size/2;
00650         vector < float >r;
00651 
00652         r.resize(np);
00653 
00654 //      float ds = 1 / (apix * size);
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 //                      printf("%d\t%f\n",j,snr[j]);
00696                         s+=ds;
00697                 }
00698                 break;
00699         case CTF_SNR_SMOOTH:
00700                 // This apparently complicated routine tries to make a nice smooth and accurate SNR curve. It does this
00701                 // by fitting local regions of the SNR vs the theoretical SNR (theoretical CTF^2/measured background),
00702                 // then taking the slope of the result times the theoretical SNR to produce a local SNR estimate
00703 
00704                 { // <- is to permit new temporary value allocation
00705                         vector < float >tsnr;   // theoretical SNR
00706                         tsnr.resize(np);
00707                         vector < float >dsnr;   // data SNR
00708                         dsnr.resize(np);
00709                         
00710                         float s0=s;
00711                         
00712                         for (int i = 0; i < np; i++) {
00713                                 float gamma = calc_gamma(g1, g2, s);
00714                                 tsnr[i] = calc_ctf1(amp1, gamma, s);            // ctf amp
00715 
00716                                 // background value
00717                                 float f = s/dsbg;
00718                                 int j = (int)floor(f);
00719                                 f-=j;
00720                                 float bg;
00721                                 if (j>(int)background.size()-2) bg=background.back();
00722                                 else bg=background[j]*(1.0f-f)+background[j+1]*f;
00723                                 if (bg <=0) bg=.001f;
00724 
00725                                 tsnr[i] = tsnr[i]*tsnr[i]/bg;           // This is now a SNR curve
00726                                 if (sf && s) {
00727                                         tsnr[i] *= sf->get_yatx(s);
00728                                 }
00729 
00730                                 
00731                                 // This is the SNR computed from the data without fitting
00732                                 if (j>(int)snr.size()-2) dsnr[i]=snr.back();
00733                                 else dsnr[i]=snr[j]*(1.0f-f)+snr[j+1]*f;
00734                                 
00735                                 s+=ds;
00736                         }
00737 
00738                         int npsm=np/25;                 // 1/2 number of points to smooth over, 25 is arbitrary
00739                         if (npsm<2) npsm=2;
00740                         
00741                         s=s0;
00742                         for (int i = 1; i < np; i++) {
00743                                 // simple linear regression embedded here
00744                                 double sum = 0;
00745                                 double sum_x = 0;
00746                                 double sum_y = 0;
00747                                 double sum_xx = 0;
00748                                 double sum_xy = 0;
00749 
00750                                 for (int k=max_int(i-npsm,1); k<=min_int(i+npsm,np-1); k++) {
00751                                         double y = dsnr[k];
00752                                         double x = tsnr[k];
00753 
00754                                         sum_x += x;
00755                                         sum_y += y;
00756                                         sum_xx += x * x;
00757                                         sum_xy += x * y;
00758                                         sum++;
00759                                 }
00760 
00761                                 double div = sum * sum_xx - sum_x * sum_x;
00762 //                              if (div == 0) {
00763 //                                      div = 0.0000001f;
00764 //                              }
00765 
00766         //                      *intercept = (float) ((sum_xx * sum_y - sum_x * sum_xy) / div);
00767         //                      *slope = (float) ((sum * sum_xy - sum_x * sum_y) / div);
00768 
00769                                 if (div!=0.0) r[i]=(float) ((sum * sum_xy - sum_x * sum_y) / div)*tsnr[i];
00770                                 else r[i]=0.0;
00771                                 if (r[i]<0) r[i]=0;
00772                                 
00773                                 s+=ds;
00774                         }
00775                         r[0]=0;
00776                 }
00777                 break;
00778 
00779         case CTF_WIENER_FILTER:
00780 //              if (!sf) {
00781 //                      LOGERR("CTF computation error, no SF found\n");
00782 //                      return r;
00783 //              }
00784 
00785                 for (int i = 0; i < np; i++) {
00786                         float f = s/dsbg;
00787                         int j = (int)floor(f);
00788                         float bg;
00789                         f-=j;
00790                         if (j>(int)snr.size()-2) {
00791 /*                              r[i]=snr.back();
00792                                 bg=background.back();*/
00793                                 r[i]=0;
00794                         }
00795                         else {
00796                                 r[i]=snr[j]*(1.0f-f)+snr[j+1]*f;
00797                                 bg=background[j]*(1.0f-f)+background[j+1]*f;
00798                         }
00799                         if (r[i]<0) r[i]=0;
00800                         r[i]=r[i]/(r[i]+1.0f);          // Full Wiener filter assuming perfect image with noise
00801 //                      r[i]=sqrt(r[i]/bg)/(r[i]+1.0);  // Wiener filter with 1/CTF term (sort of) to restore image then filter
00802                         s+=ds;
00803                 }
00804                 r[0]=0;
00805                 break;
00806 
00807         case CTF_TOTAL:
00808 
00809                 for (int i = 0; i < np; i++) {
00810                         float gamma = calc_gamma(g1, g2, s);
00811                         if (sf) {
00812                                 r[i] = calc_ctf1(amp1, gamma, s);
00813                                 r[i] = r[i] * r[i] * sf->get_yatx(s) + calc_noise(s);
00814                         }
00815                         else {
00816                                 r[i] = calc_ctf1(amp1, gamma, s);
00817                                 r[i] = r[i] * r[i] + calc_noise(s);
00818                         }
00819                         s += ds;
00820                 }
00821                 break;
00822         default:
00823                 break;
00824         }
00825 
00826         return r;
00827 }
00828 
00829 
00830 void EMAN2Ctf::compute_2d_real(EMData *, CtfType, XYData *)
00831 {
00832 
00833 
00834 }
00835 
00836 
00837 
00838 void EMAN2Ctf::compute_2d_complex(EMData * image, CtfType type, XYData * sf)
00839 {
00840         if (!image) {
00841                 LOGERR("image is null. cannot computer 2D complex CTF");
00842                 return;
00843         }
00844 
00845         if (image->is_complex() == false) {
00846                 LOGERR("compute_2d_complex can only work on complex images");
00847                 return;
00848         }
00849 
00850         int nx = image->get_xsize();
00851         int ny = image->get_ysize();
00852 
00853         if ((ny%2==1 && nx!=ny+1) || (ny%2==0 && nx != ny + 2)) {
00854                 printf("nx=%d  ny=%d\n",nx,ny);
00855                 LOGERR("compute_2d_complex only works on (nx, nx-2) images");
00856                 return;
00857         }
00858 
00859         float ds = 1.0f / (apix * ny);
00860         image->to_one();
00861 
00862         float *d = image->get_data();
00863         float g1 = calc_g1();
00864         float g2 = calc_g2();
00865         float amp1 = calc_amp1();
00866 
00867         if (type == CTF_BACKGROUND) {
00868                 for (int y = -ny/2; y < ny/2; y++) {
00869                         int y2=(y+ny)%ny;
00870                         int ynx = y2 * nx;
00871 
00872                         for (int x = 0; x < nx / 2; x++) {
00873                                 float s = (float) Util::hypot_fast(x, y ) * ds;
00874                                 d[x * 2 + ynx] = calc_noise(s);
00875                                 d[x * 2 + ynx + 1] = 0;                 // The phase is somewhat arbitrary
00876                         }
00877                 }
00878         }
00879         else if (type == CTF_AMP) {
00880                 for (int y = -ny/2; y < ny/2; y++) {
00881                         int y2=(y+ny)%ny;
00882                         int ynx = y2 * nx;
00883 
00884                         for (int x = 0; x < nx / 2; x++) {
00885                                 float s = (float)Util::hypot_fast(x,y ) * ds;
00886                                 float gamma = calc_gamma(g1, g2, s);
00887                                 float v = calc_ctf1(amp1, gamma, s);
00888 //                              float v = calc_amplitude(gamma);
00889                                 d[x * 2 + ynx] = v;
00890                                 d[x * 2 + ynx + 1] = 0;
00891                         }
00892                 }
00893         }
00894         else if (type == CTF_SIGN) {
00895                 for (int y = -ny/2; y < ny/2; y++) {
00896                         int y2=(y+ny)%ny;
00897                         int ynx = y2 * nx;
00898                         for (int x = 0; x < nx / 2; x++) {
00899                                 float s = (float)Util::hypot_fast(x,y ) * ds;
00900                                 float gamma = calc_gamma(g1, g2, s);
00901                                 float v = calc_amplitude(gamma);
00902                                 d[x * 2 + ynx] = v >= 0 ? 1.0f : -1.0f;
00903                                 d[x * 2 + ynx + 1] = 0;
00904                         }
00905                 }
00906         }
00907         else if (type == CTF_SNR) {
00908 
00909                 for (int y = -ny/2; y < ny/2; y++) {
00910                         int y2=(y+ny)%ny;
00911                         int ynx = y2 * nx;
00912 
00913                         for (int x = 0; x < nx / 2; x++) {
00914 
00915                                 float s = (float)Util::hypot_fast(x,y ) * ds;
00916                                 float f = s/dsbg;
00917                                 int j = (int)floor(f);
00918                                 f-=j;
00919                                 if (j>(int)snr.size()-2) d[x*2+ynx]=snr.back();
00920                                 else d[x*2+ynx]=snr[j]*(1.0f-f)+snr[j+1]*f;
00921                                 d[x * 2 + ynx + 1] = 0;
00922                         }
00923                 }
00924                 d[0]=0;
00925         }
00926         else if (type == CTF_SNR_SMOOTH) {
00927                 for (int y = -ny/2; y < ny/2; y++) {
00928                         int y2=(y+ny)%ny;
00929                         int ynx = y2 * nx;
00930 
00931                         for (int x = 0; x < nx / 2; x++) {
00932 
00933                                 float s = (float)Util::hypot_fast(x,y ) * ds;
00934                                 float f = s/dsbg;
00935                                 int j = (int)floor(f);
00936                                 f-=j;
00937                                 if (j>(int)snr.size()-2) d[x*2+ynx]=snr.back();
00938                                 else d[x*2+ynx]=snr[j]*(1.0f-f)+snr[j+1]*f;
00939                                 d[x * 2 + ynx + 1] = 0;
00940                         }
00941                 }
00942                 d[0]=0;
00943         }
00944         else if (type == CTF_WIENER_FILTER) {
00945                 if (dsbg==0) printf("Warning, DSBG set to 0\n");
00946                 for (int y = -ny/2; y < ny/2; y++) {
00947                         int y2=(y+ny)%ny;
00948                         int ynx = y2 * nx;
00949 
00950                         for (int x = 0; x < nx / 2; x++) {
00951 
00952                                 float s = (float)Util::hypot_fast(x,y ) * ds;
00953                                 float f = s/dsbg;
00954                                 int j = (int)floor(f);
00955                                 float bg,snrf;
00956                                 f-=j;
00957                                 if (j>(int)snr.size()-2) {
00958 /*                                      bg=background.back();
00959                                         d[x*2+ynx]=snr.back()/(snr.back()+1.0);*/
00960                                         d[x*2+ynx]=0;
00961                                 }
00962                                 else {
00963                                         bg=background[j]*(1.0f-f)+background[j+1]*f;
00964                                         snrf=snr[j]*(1.0f-f)+snr[j+1]*f;
00965 
00966 //                                      printf("%d\t%f\n",x,sqrt(snrf/bg)/(snrf+1.0));
00967                                         if (snrf<0) snrf=0.0;
00968 //                                      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
00969                                         d[x*2+ynx]=snrf/(snrf+1);       // This is just the simple Wiener filter
00970 
00971                                 }
00972                                 d[x * 2 + ynx + 1] = 0;
00973                         }
00974                 }
00975                 d[0]=0;
00976         }
00977         else if (type == CTF_TOTAL) {
00978                 float amp1 = calc_amp1();
00979 
00980                 for (int y = -ny/2; y < ny/2; y++) {
00981                         int y2=(y+ny)%ny;
00982                         int ynx = y2 * nx;
00983 
00984                         for (int x = 0; x < nx / 2; x++) {
00985 
00986                                 float s = (float)Util::hypot_fast(x,y ) * ds;
00987                                 float gamma = calc_gamma(g1, g2, s);
00988                                 float f = calc_ctf1(amp1, gamma, s);
00989                                 float noise = 0;
00990                                 f = f * f;
00991 
00992                                 if (sf && s) {
00993                                         f *= sf->get_yatx(s);
00994                                 }
00995                                 f+=noise;
00996 
00997                                 d[x * 2 + ynx] *= f;
00998                                 d[x * 2 + ynx + 1] = 0;
00999                         }
01000                 }
01001         }
01002         else printf("Unknown CTF image mode\n");
01003 
01004         image->update();
01005 }
01006 
01007 
01008 
01009 bool EMAN2Ctf::equal(const Ctf * ctf1) const
01010 {
01011         if (ctf1) {
01012                 const EMAN2Ctf *c = static_cast<const EMAN2Ctf *>(ctf1);
01013                 if (defocus != c->defocus ||
01014                         dfdiff != c->dfdiff ||
01015                         dfang != c->dfang ||
01016                         bfactor != c->bfactor ||
01017                         ampcont != c->ampcont ||
01018                         voltage != c->voltage ||
01019                         cs != c->cs ||
01020                         apix != c->apix ||
01021                         dsbg != c->dsbg ||
01022                         background.size() != c->background.size() ||
01023                         snr.size() != c->snr.size()
01024                         ) return false;
01025 
01026                 for (unsigned int i=0; i<background.size(); i++) {
01027                         if (background[i]!=c->background[i]) return false;
01028                 }
01029                 for (unsigned int i=0; i<snr.size(); i++) {
01030                         if (snr[i]!=c->snr[i]) return false;
01031                 }
01032                 return true;
01033         }
01034         return false;
01035 }

Generated on Thu Nov 17 12:43:44 2011 for EMAN2 by  doxygen 1.3.9.1