EMAN2
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                         CTF_FITREF,                     // CTF amplitude squared without B-factor and low resolution zeroed
00077                         CTF_NOISERATIO,         // 1-Noise/Total, when a particle is filtered with this it will still have noise, but the structure factor will look as if it's noise-free
00078                         CTF_INTEN                       // ctf intensity only (no background or envelope)
00079                 };
00080           public:
00081                 virtual ~ Ctf()
00082                 {
00083                 }
00084 
00085                 float defocus;                  //      Defocus in microns, note that postitive is now underfocus, whereas values in EMAN1 are negative overfocus
00086                 float bfactor;                  //      B-factor expressed using the x-ray convention (e^-B/4 s^2 in amplitude space) EMAN1 used E^-B s^2
00087                 float voltage;                  //  in kV
00088                 float cs;                               //  in mm
00089                 float apix;                             //
00090 
00091                 virtual int from_string(const string & ctf) = 0;        // The first letter of the string indicates the subclass type
00092                 virtual string to_string() const = 0;
00093 
00094                 virtual void from_dict(const Dict & dict) = 0;
00095                 virtual Dict to_dict() const = 0;
00096 
00097                 virtual void from_vector(const vector<float>& vctf) = 0;
00098                 virtual vector<float> to_vector() const = 0;
00099 
00100                 virtual float zero(int n) const = 0;
00101 
00102                 virtual vector < float >compute_1d(int size,float ds, CtfType t, XYData * struct_factor = 0) = 0;
00103                 // This computes a 1D power spectrum from an image with astigmatism compensation. This is not entirely valid.
00104                 // It will stretch/shrink the power spectrum as a function of angle to make a single coherent 1D CTF curve.
00105                 // The cost is that the shrinking means that the structure factor is being distorted to make this work. 
00106                 // Nonetheless, this is a useful tool in fitting astigmatic data
00107                 virtual vector <float>compute_1d_fromimage(int size, float ds, EMData *image) = 0;
00108                 virtual void compute_2d_real(EMData * img, CtfType t, XYData * struct_factor = 0) = 0;
00109                 virtual void compute_2d_complex(EMData * img, CtfType t, XYData * struct_factor = 0) = 0;
00110 
00111                 virtual void copy_from(const Ctf * new_ctf) = 0;
00112                 virtual bool equal(const Ctf * ctf1) const = 0;
00113 
00114           public:
00115                 enum
00116                 { CTFOS = 5 };
00117 
00118         };
00119 
00122         class EMAN1Ctf:public Ctf
00123         {
00124           public:
00125 //              float defocus;                  // 0    Defocus in microns, note that postitive is now underfocus, whereas values in EMAN1 are negative overfocus
00126 //              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
00127                 float amplitude;                // 2
00128                 float ampcont;                  // 3
00129                 float noise1;                   // 4
00130                 float noise2;                   // 5
00131                 float noise3;                   // 6
00132                 float noise4;                   // 7
00133 //              float voltage;                  // 8
00134 //              float cs;                               // 9
00135 //              float apix;                             // 10
00136 
00137           public:
00138                 EMAN1Ctf();
00139                 EMAN1Ctf(const vector<float>& vf) {from_vector(vf);}    //for unpickling
00140                 ~EMAN1Ctf();
00141 
00142                 vector < float >compute_1d(int size,float ds, CtfType type, XYData * struct_factor = 0);
00143                 vector <float>compute_1d_fromimage(int size, float ds, EMData *image);  
00144                 void compute_2d_real(EMData * image, CtfType type, XYData * struct_factor = 0);
00145                 void compute_2d_complex(EMData * image, CtfType type, XYData * struct_factor = 0);
00146 
00147                 int from_string(const string & ctf);
00148                 string to_string() const;
00149 
00150                 void from_dict(const Dict & dict);
00151                 Dict to_dict() const;
00152 
00153                 void from_vector(const vector<float>& vctf);
00154                 vector<float> to_vector() const;
00155 
00156                 void copy_from(const Ctf * new_ctf);
00157                 bool equal(const Ctf * ctf1) const;
00158 
00159                 float zero(int n) const;
00160 
00161                 float get_defocus() const
00162                 {
00163                         return defocus;
00164                 }
00165                 float get_bfactor() const
00166                 {
00167                         return bfactor;
00168                 }
00169 
00170           private:
00171                 inline float calc_amp1()
00172                 {
00173                         return (sqrt(1 - ampcont * ampcont/10000.0f));
00174                 }
00175 
00176                 inline float calc_lambda()
00177                 {
00178                         float lambda = 12.2639f / sqrt(voltage * 1000.0f + 0.97845f * voltage * voltage);
00179                         return lambda;
00180                 }
00181 
00182                 inline float calc_g1()
00183                 {
00184                         float lambda = calc_lambda();
00185                         float g1 = 2.5e6f * cs * lambda * lambda * lambda;
00186                         return g1;
00187                 }
00188 
00189                 inline float calc_g2()
00190                 {
00191                         float lambda = calc_lambda();
00192                         float g2 = 5000.0f * -defocus * lambda;
00193                         return g2;
00194                 }
00195 
00196                 inline float calc_gamma(float g1, float g2, float s)
00197                 {
00198                         float s2 = s * s;
00199                         float gamma = (float) (-2 * M_PI * (g1 * s2 * s2 + g2 * s2));
00200                         return gamma;
00201                 }
00202 
00203                 inline float calc_ctf1(float g, float gamma, float s)
00204                 {
00205                         float r = amplitude * exp(-(bfactor/4.0f * s * s)) * (g * sin(gamma) + ampcont/100.0f * cos(gamma));
00206                         return r;
00207                 }
00208 
00209                 inline float calc_amplitude(float gamma)
00210                 {
00211                         float t1 = sqrt(1.0f - ampcont * ampcont/10000.0f) * sin(gamma);
00212                         float v = amplitude * (t1 + ampcont/100.0f * cos(gamma));
00213                         return v;
00214                 }
00215 
00216                 inline float calc_noise(float s)
00217                 {
00218                         float ns = (float) M_PI / 2 * noise4 * s;
00219                         float ns2 = ns * ns;
00220                         float n = noise3 * exp(-ns2 - s * noise2 - sqrt(fabs(s)) * noise1);
00221                         return n;
00222                 }
00223 
00224                 inline float calc_snr(float g1, float gamma, float s)
00225                 {
00226                         float ctf1 = calc_ctf1(g1, gamma, s);
00227                         float ctf2 = ctf1 * ctf1 / calc_noise(s);
00228                         return ctf2;
00229                 }
00230 
00231         };
00232 
00235         class EMAN2Ctf:public Ctf
00236         {
00237           public:
00238 //              float defocus;          // defocus in microns, positive underfocus
00239                 float dfdiff;           // defocus difference for astigmatism, defocus is the major elliptical axis
00240                 float dfang;            // angle of the major elliptical axis in degrees measured counterclockwise from x
00241 //              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
00242                 float ampcont;          // amplitude contrast as a percentage ie- this should be 10 for 10% amp contrast
00243 //              float voltage;          // microscope voltage in kV
00244 //              float cs;                       // Cs in mm
00245 //              float apix;                     // A/pix value used when generating 2D results
00246                 float dsbg;                     // ds value for background and SNR
00247                 vector<float> background;       // background intensity, 1 value per radial pixel (NX/2, corners ignored)
00248                 vector<float> snr;                      // SNR, 1 value per radial pixel (NX/2, corners assumed 0)
00249 
00250                 vector<float> get_snr(){ return snr;}
00251                 void set_snr(const vector<float>& vf) {snr = vf;}
00252                 vector<float> get_background(){ return background;}
00253                 void set_background(const vector<float>& vf) {background = vf;}
00254 
00255           public:
00256                 EMAN2Ctf();
00257                 EMAN2Ctf(const vector<float>& vf) {from_vector(vf);}    //for unpickling
00258                 ~EMAN2Ctf();
00259 
00260                 vector <float> compute_1d(int size,float ds, CtfType type, XYData * struct_factor = 0);
00261                 vector <float> compute_1d_fromimage(int size, float ds, EMData *image);
00262                 void compute_2d_real(EMData * image, CtfType type, XYData * struct_factor = 0);
00263                 void compute_2d_complex(EMData * image, CtfType type, XYData * struct_factor = 0);
00264 
00265                 int from_string(const string & ctf);
00266                 string to_string() const;
00267 
00268                 void from_dict(const Dict & dict);
00269                 Dict to_dict() const;
00270 
00271                 void from_vector(const vector<float>& vctf);
00272                 vector<float> to_vector() const;
00273 
00274                 void copy_from(const Ctf * new_ctf);
00275                 bool equal(const Ctf * ctf1) const;
00276 
00277                 float zero(int n) const;
00278 
00279                 inline float stos2(float s,float dZ) {
00280                         float lmb=lambda();
00281                         return sqrt((defocus*1.0e4f-sqrt(defocus*defocus*1.0e8f -2.0e11f*cs*s*s*(defocus-dZ)*lmb*lmb+1.0e14f*cs*cs*pow(s*lmb,4.0f)))/(1.0e7f*cs*lmb*lmb));
00282                 }
00283         
00284                 private:
00285 
00286                 // Electron wavelength in A
00287                 inline float lambda() const
00288                 {
00289                         float lambda = 12.2639f / sqrt(voltage * 1000.0f + 0.97845f * voltage * voltage);
00290                         return lambda;
00291                 }
00292 
00293                 // returns defocus as a function of angle. ang in radians !
00294                 inline float df(float ang) const {
00295                         if (dfdiff==0.0) return defocus;
00296                         return defocus+dfdiff/2.0f*cos(2.0f*ang-2.0f*M_PI/180.0f*dfang);
00297                 }
00298                 
00299                 inline float calc_noise(float s) const
00300                 {
00301                         int si=(int)(s/dsbg);
00302                         if (si>(int)background.size()||si<0) return background.back();
00303                         return background[si];
00304                 }
00305 
00306         };
00307 
00308 }
00309 
00310 
00311 
00312 #endif  //eman__ctf__h__