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

ctf.h

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 #ifndef eman__ctf__h__
00037 #define eman__ctf__h__ 1
00038 
00039 #include <cmath>
00040 
00041 #include "emobject.h"
00042 
00043 #ifdef WIN32
00044         #ifndef M_PI
00045                 #define M_PI 3.14159265358979323846f
00046         #endif  //M_PI
00047 #endif  //WIN32
00048 
00049 using std::string;
00050 using std::map;
00051 
00052 namespace EMAN
00053 {
00054         class EMData;
00055         class XYData;
00056 
00063         class Ctf
00064         {
00065           public:
00066                 // NOTE: ctf is positive for the first peak, instead of negative
00067                 enum CtfType
00068                 {
00069                         CTF_AMP,                        // ctf ampltidue only
00070                         CTF_SIGN,                       // ctf sign (+-1)
00071                         CTF_BACKGROUND,         // Background, no ctf oscillation
00072                         CTF_SNR,                        // Signal to noise ratio
00073                         CTF_SNR_SMOOTH,         // Signal to noise ratio, smoothed, algorithm may vary, but this should be more suitable for weighting
00074                         CTF_WIENER_FILTER,      // Weiner Filter = 1/(1+1/snr)
00075                         CTF_TOTAL                       // AMP*AMP+NOISE
00076                 };
00077           public:
00078                 virtual ~ Ctf()
00079                 {
00080                 }
00081 
00082                 float defocus;                  //      Defocus in microns, note that postitive is now underfocus, whereas values in EMAN1 are negative overfocus
00083                 float bfactor;                  //      B-factor expressed using the x-ray convention (e^-B/4 s^2 in amplitude space) EMAN1 used E^-B s^2
00084                 float voltage;                  //  in kV
00085                 float cs;                               //  in mm
00086                 float apix;                             //
00087 
00088                 virtual int from_string(const string & ctf) = 0;        // The first letter of the string indicates the subclass type
00089                 virtual string to_string() const = 0;
00090 
00091                 virtual void from_dict(const Dict & dict) = 0;
00092                 virtual Dict to_dict() const = 0;
00093 
00094                 virtual void from_vector(const vector<float>& vctf) = 0;
00095                 virtual vector<float> to_vector() const = 0;
00096 
00097                 virtual vector < float >compute_1d(int size,float ds, CtfType t, XYData * struct_factor = 0) = 0;
00098                 virtual void compute_2d_real(EMData * img, CtfType t, XYData * struct_factor = 0) = 0;
00099                 virtual void compute_2d_complex(EMData * img, CtfType t, XYData * struct_factor = 0) = 0;
00100 
00101                 virtual void copy_from(const Ctf * new_ctf) = 0;
00102                 virtual bool equal(const Ctf * ctf1) const = 0;
00103 
00104           public:
00105                 enum
00106                 { CTFOS = 5 };
00107 
00108         };
00109 
00112         class EMAN1Ctf:public Ctf
00113         {
00114           public:
00115 //              float defocus;                  // 0    Defocus in microns, note that postitive is now underfocus, whereas values in EMAN1 are negative overfocus
00116 //              float bfactor;                  // 1    B-factor expressed using the x-ray convention (e^-B/4 s^2 in amplitude space) EMAN1 used E^-B s^2
00117                 float amplitude;                // 2
00118                 float ampcont;                  // 3
00119                 float noise1;                   // 4
00120                 float noise2;                   // 5
00121                 float noise3;                   // 6
00122                 float noise4;                   // 7
00123 //              float voltage;                  // 8
00124 //              float cs;                               // 9
00125 //              float apix;                             // 10
00126 
00127           public:
00128                 EMAN1Ctf();
00129                 EMAN1Ctf(const vector<float>& vf) {from_vector(vf);}    //for unpickling
00130                 ~EMAN1Ctf();
00131 
00132                 vector < float >compute_1d(int size,float ds, CtfType type, XYData * struct_factor = 0);
00133                 void compute_2d_real(EMData * image, CtfType type, XYData * struct_factor = 0);
00134                 void compute_2d_complex(EMData * image, CtfType type, XYData * struct_factor = 0);
00135 
00136                 int from_string(const string & ctf);
00137                 string to_string() const;
00138 
00139                 void from_dict(const Dict & dict);
00140                 Dict to_dict() const;
00141 
00142                 void from_vector(const vector<float>& vctf);
00143                 vector<float> to_vector() const;
00144 
00145                 void copy_from(const Ctf * new_ctf);
00146                 bool equal(const Ctf * ctf1) const;
00147 
00148                 float get_defocus() const
00149                 {
00150                         return defocus;
00151                 }
00152                 float get_bfactor() const
00153                 {
00154                         return bfactor;
00155                 }
00156 
00157           private:
00158                 inline float calc_amp1()
00159                 {
00160                         return (sqrt(1 - ampcont * ampcont/10000.0f));
00161                 }
00162 
00163                 inline float calc_lambda()
00164                 {
00165                         float lambda = 12.2639f / sqrt(voltage * 1000.0f + 0.97845f * voltage * voltage);
00166                         return lambda;
00167                 }
00168 
00169                 inline float calc_g1()
00170                 {
00171                         float lambda = calc_lambda();
00172                         float g1 = 2.5e6f * cs * lambda * lambda * lambda;
00173                         return g1;
00174                 }
00175 
00176                 inline float calc_g2()
00177                 {
00178                         float lambda = calc_lambda();
00179                         float g2 = 5000.0f * -defocus * lambda;
00180                         return g2;
00181                 }
00182 
00183                 inline float calc_gamma(float g1, float g2, float s)
00184                 {
00185                         float s2 = s * s;
00186                         float gamma = (float) (-2 * M_PI * (g1 * s2 * s2 + g2 * s2));
00187                         return gamma;
00188                 }
00189 
00190                 inline float calc_ctf1(float g, float gamma, float s)
00191                 {
00192                         float r = amplitude * exp(-(bfactor/4.0f * s * s)) * (g * sin(gamma) + ampcont/100.0f * cos(gamma));
00193                         return r;
00194                 }
00195 
00196                 inline float calc_amplitude(float gamma)
00197                 {
00198                         float t1 = sqrt(1.0f - ampcont * ampcont/10000.0f) * sin(gamma);
00199                         float v = amplitude * (t1 + ampcont/100.0f * cos(gamma));
00200                         return v;
00201                 }
00202 
00203                 inline float calc_noise(float s)
00204                 {
00205                         float ns = (float) M_PI / 2 * noise4 * s;
00206                         float ns2 = ns * ns;
00207                         float n = noise3 * exp(-ns2 - s * noise2 - sqrt(fabs(s)) * noise1);
00208                         return n;
00209                 }
00210 
00211                 inline float calc_snr(float g1, float gamma, float s)
00212                 {
00213                         float ctf1 = calc_ctf1(g1, gamma, s);
00214                         float ctf2 = ctf1 * ctf1 / calc_noise(s);
00215                         return ctf2;
00216                 }
00217 
00218         };
00219 
00222         class EMAN2Ctf:public Ctf
00223         {
00224           public:
00225 //              float defocus;          // defocus in microns, positive underfocus
00226                 float dfdiff;           // defocus difference for astigmatism, defocus is the major elliptical axis
00227                 float dfang;            // angle of the major elliptical axis in degrees measured counterclockwise from x
00228 //              float bfactor;          // B-factor in 1/A^2 expressed using the x-ray convention (e^-B/4 s^2 in amplitude space) EMAN1 used E^-B s^2
00229                 float ampcont;          // amplitude contrast as a percentage ie- this should be 10 for 10% amp contrast
00230 //              float voltage;          // microscope voltage in kV
00231 //              float cs;                       // Cs in mm
00232 //              float apix;                     // A/pix value used when generating 2D results
00233                 float dsbg;                     // ds value for background and SNR
00234                 vector<float> background;       // background intensity, 1 value per radial pixel (NX/2, corners ignored)
00235                 vector<float> snr;                      // SNR, 1 value per radial pixel (NX/2, corners assumed 0)
00236 
00237                 vector<float> get_snr(){ return snr;}
00238                 void set_snr(const vector<float>& vf) {snr = vf;}
00239                 vector<float> get_background(){ return background;}
00240                 void set_background(const vector<float>& vf) {background = vf;}
00241 
00242           public:
00243                 EMAN2Ctf();
00244                 EMAN2Ctf(const vector<float>& vf) {from_vector(vf);}    //for unpickling
00245                 ~EMAN2Ctf();
00246 
00247                 vector < float >compute_1d(int size,float ds, CtfType type, XYData * struct_factor = 0);
00248                 void compute_2d_real(EMData * image, CtfType type, XYData * struct_factor = 0);
00249                 void compute_2d_complex(EMData * image, CtfType type, XYData * struct_factor = 0);
00250 
00251                 int from_string(const string & ctf);
00252                 string to_string() const;
00253 
00254                 void from_dict(const Dict & dict);
00255                 Dict to_dict() const;
00256 
00257                 void from_vector(const vector<float>& vctf);
00258                 vector<float> to_vector() const;
00259 
00260                 void copy_from(const Ctf * new_ctf);
00261                 bool equal(const Ctf * ctf1) const;
00262 
00263           private:
00264                 inline float calc_amp1()
00265                 {
00266                         return (sqrt(1 - ampcont * ampcont/10000.0f));
00267                 }
00268 
00269                 inline float calc_lambda()
00270                 {
00271                         float lambda = 12.2639f / sqrt(voltage * 1000.0f + 0.97845f * voltage * voltage);
00272                         return lambda;
00273                 }
00274 
00275                 inline float calc_g1()
00276                 {
00277                         float lambda = calc_lambda();
00278                         float g1 = 2.5e6f * cs * lambda * lambda * lambda;
00279                         return g1;
00280                 }
00281 
00282                 inline float calc_g2()
00283                 {
00284                         float lambda = calc_lambda();
00285                         float g2 = 5000.0f * -defocus * lambda;
00286                         return g2;
00287                 }
00288 
00289                 inline float calc_gamma(float g1, float g2, float s)
00290                 {
00291                         float s2 = s * s;
00292                         float gamma = (float) (-2 * M_PI * (g1 * s2 * s2 + g2 * s2));
00293                         return gamma;
00294                 }
00295 
00296                 inline float calc_ctf1(float g, float gamma, float s)
00297                 {
00298                         float r = exp(-(bfactor/4.0f * s * s)) * (g * sin(gamma) + ampcont/100.0f * cos(gamma));
00299                         return r;
00300                 }
00301 
00302                 inline float calc_amplitude(float gamma)
00303                 {
00304                         float t1 = sqrt(1.0f - ampcont * ampcont/10000.0f) * sin(gamma);
00305                         float v = (t1 + ampcont/100.0f * cos(gamma));
00306                         return v;
00307                 }
00308 
00309                 inline float calc_noise(float s)
00310                 {
00311                         int si=(int)(s/dsbg);
00312                         if (si>(int)background.size()||si<0) return background.back();
00313                         return background[si];
00314                 }
00315 
00316         };
00317 
00318 }
00319 
00320 
00321 
00322 #endif  //eman__ctf__h__

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