00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
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
00067 enum CtfType
00068 {
00069 CTF_AMP,
00070 CTF_SIGN,
00071 CTF_BACKGROUND,
00072 CTF_SNR,
00073 CTF_SNR_SMOOTH,
00074 CTF_WIENER_FILTER,
00075 CTF_TOTAL
00076 };
00077 public:
00078 virtual ~ Ctf()
00079 {
00080 }
00081
00082 float defocus;
00083 float bfactor;
00084 float voltage;
00085 float cs;
00086 float apix;
00087
00088 virtual int from_string(const string & ctf) = 0;
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
00116
00117 float amplitude;
00118 float ampcont;
00119 float noise1;
00120 float noise2;
00121 float noise3;
00122 float noise4;
00123
00124
00125
00126
00127 public:
00128 EMAN1Ctf();
00129 EMAN1Ctf(const vector<float>& vf) {from_vector(vf);}
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
00226 float dfdiff;
00227 float dfang;
00228
00229 float ampcont;
00230
00231
00232
00233 float dsbg;
00234 vector<float> background;
00235 vector<float> snr;
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);}
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__