Xmipp  v3.23.11-Nereus
ctf.h
Go to the documentation of this file.
1 /***************************************************************************
2  *
3  * Authors: Carlos Oscar S. Sorzano (coss@cnb.csic.es)
4  *
5  * Unidad de Bioinformatica of Centro Nacional de Biotecnologia , CSIC
6  *
7  * This program is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation; either version 2 of the License, or
10  * (at your option) any later version.
11  *
12  * This program is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with this program; if not, write to the Free Software
19  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
20  * 02111-1307 USA
21  *
22  * All comments concerning this program package may be sent to the
23  * e-mail address 'xmipp@cnb.csic.es'
24  ***************************************************************************/
25 
26 #ifndef _CORE_CTF_HH
27 #define _CORE_CTF_HH
28 
29 #include <complex>
30 #include "core/metadata_db.h"
31 #include "core/numerical_recipes.h"
32 #include "core/xmipp_fft.h"
33 #include "core/xmipp_macros.h"
34 
35 template<typename T>
36 class MultidimArray;
37 template<typename T>
38 class Matrix1D;
39 class XmippProgram;
40 
41 const int CTF_BASIC_LABELS_SIZE = 5;
42 const MDLabel CTF_BASIC_LABELS[] =
43  {
44  MDL_CTF_DEFOCUSU, MDL_CTF_DEFOCUSV, MDL_CTF_DEFOCUS_ANGLE, MDL_CTF_CS, MDL_CTF_Q0
45  };
46 
47 const int CTF_ALL_LABELS_SIZE = 13;
48 const MDLabel CTF_ALL_LABELS[] =
49  {
50  MDL_CTF_VOLTAGE,
51  MDL_CTF_DEFOCUSU,
52  MDL_CTF_DEFOCUSV,
53  MDL_CTF_DEFOCUS_ANGLE,
54  MDL_CTF_CS,
55  MDL_CTF_CA,
56  MDL_CTF_ENERGY_LOSS,
57  MDL_CTF_LENS_STABILITY,
58  MDL_CTF_CONVERGENCE_CONE,
59  MDL_CTF_LONGITUDINAL_DISPLACEMENT,
60  MDL_CTF_TRANSVERSAL_DISPLACEMENT,
61  MDL_CTF_Q0,
62  MDL_CTF_K
63  };
64 
66 bool containsCTFBasicLabels(const MetaData &md);
67 
76 void groupCTFMetaData(const MetaDataDb &imgMd, MetaDataDb &ctfMd, std::vector<MDLabel> &groupbyLabels);
77 
82 class PrecomputedForCTF
83 {
84 public:
85  double u_sqrt;
86  double u;
87  double u2;
88  double u3;
89  double u4;
90  double ang;
91  double deltaf;
92 };
93 
209 class CTFDescription1D
211 {
212 public:
213  // Electron wavelength (Angstroms)
214  double lambda;
215  // Squared frequency associated to the aperture
216  // double ua2;
217  // Different constants
218  double K1;
219  double K2;
220  double K3;
221  double K4;
222  double K5;
223  double K6;
224  double K7;
225  double Ksin;
226  double Kcos;
229  double D;
230  // Precomputed values
231  PrecomputedForCTF precomputed;
232  // Image of precomputed values
233  std::vector<PrecomputedForCTF> precomputedImage;
234  // Xdim size of the image
235  int precomputedImageXdim;
236 public:
238  double K;
240  double Tm;
242  double kV;
244  double Defocus;
246  double Cs;
248  double Ca;
251  double espr;
253  double ispr;
255  double alpha;
257  double DeltaF;
259  double DeltaR;
261  double Q0;
263  double x0;
265  double xF;
267  double y0;
269  double yF;
271  bool isLocalCTF;
273  bool enable_CTFnoise;
275  bool enable_CTF;
277  double base_line;
279  double gaussian_K;
281  double sigma1;
283  double Gc1;
285  double sqrt_K;
287  double sq;
289  double gaussian_K2;
291  double sigma2;
293  double Gc2;
294  // Background polynomial
295  double bgR1;
296  double bgR2;
297  double bgR3;
298  // Envelope polynomial
299  double envR0;
300  double envR1;
301  double envR2;
302  //Maximum frequency to estimate values
303  double freq_max;
304  //Extra parameters for VPP
305  double phase_shift;
306  double VPP_radius;
307 
309  CTFDescription1D()
310  {
311  clear();
312  }
313 
321  void read(const FileName &fn, bool disable_if_not_K = true);
322 
332  void readFromMetadataRow(const MetaData &MD, size_t id, bool disable_if_not_K=true);
333 
335  void readFromMdRow(const MDRow &row, bool disable_if_not_K=true);
336 
338  void setRow(MDRow &row) const;
339 
342  void write(const FileName &fn);
343 
345  static void defineParams(XmippProgram * program);
346 
348  void readParams(XmippProgram * program);
349 
351  friend std::ostream & operator << (std::ostream &out, const CTFDescription1D &ctf);
352 
354  void clear();
355 
357  void clearNoise();
358 
360  void clearPureCtf();
361 
367  inline void changeSamplingRate(double newTm)
368  {
369  Tm=newTm;
370  }
371 
373  void produceSideInfo();
374 
376  void precomputeValues(double X)
377  {
378  precomputed.ang = 0;
379  precomputed.u2 = X * X;
380  precomputed.u = sqrt(precomputed.u2);
381  precomputed.u3 = precomputed.u2 * precomputed.u;
382  precomputed.u4 = precomputed.u2 * precomputed.u2;
383  precomputed.u_sqrt = sqrt(precomputed.u);
384 
385  if (fabs(X) < XMIPP_EQUAL_ACCURACY)
386  precomputed.deltaf=0;
387  else
388  {
389 
390  /*
391  * For a derivation of this formulae confer
392  * Principles of Electron Optics page 1380
393  * in particular term defocus and twofold axial astigmatism
394  * take into account that a1 and a2 are the coefficient
395  * of the zernike polynomials difference of defocus at 0
396  * and at 45 degrees. In this case a2=0
397  */
398  precomputed.deltaf= Defocus;//(defocus_average + defocus_deviation*cos_ellipsoid_ang_2);
399 
400  }
401  }
402 
404  void precomputeValues(const MultidimArray<double> &cont_x_freq);
405 
407  void precomputeValues(int i)
408  {
409  precomputed=precomputedImage[i+precomputedImageXdim];
410  if (precomputed.deltaf==-1)
411  {
412  precomputed.deltaf= Defocus;
413  }
414  }
415 
417  double getValueAt(bool show = false) const
418  {
419  double pure_CTF = enable_CTF ? getValuePureAt(show) : 0;
420  return enable_CTFnoise ? sqrt(pure_CTF*pure_CTF + getValueNoiseAt(show)) : pure_CTF;
421  }
422 
424  inline double getValueDampingAt(bool show = false) const
425  {
426  double Eespr = exp(-K3 * precomputed.u4); // OK
427  double EdeltaF = bessj0(K5 * precomputed.u2); // OK
428  double EdeltaR = SINC(precomputed.u * DeltaR); // OK
429  double aux=K7 * precomputed.u2 * precomputed.u + precomputed.deltaf * precomputed.u;
430  double Ealpha = exp(-K6 * aux * aux);
431  double E = Eespr * EdeltaF * EdeltaR * Ealpha+envR0+envR1*precomputed.u+envR2*precomputed.u2;
432  if (E < 0)
433  E = 0;
434  if (show)
435  {
436  std::cout << " Deltaf=" << precomputed.deltaf << std::endl;
437  std::cout << " u,u2,u4=" << precomputed.u << " " << precomputed.u2
438  << " " << precomputed.u4 << std::endl;
439  std::cout << " K3,Eespr=" << K3 << " " << Eespr << std::endl;
440  std::cout << " K4,Eispr=" << K4 << " " << /*Eispr <<*/ std::endl;
441  std::cout << " K5,EdeltaF=" << K5 << " " << EdeltaF << std::endl;
442  std::cout << " EdeltaR=" << EdeltaR << std::endl;
443  std::cout << " K6,K7,Ealpha=" << K6 << " " << K7 << " " << Ealpha
444  << std::endl;
445  std::cout << " Total atenuation(E)= " << E << std::endl;
446  std::cout << " CTFdamp=" << E << std::endl;
447  }
448  return -K*E;
449  }
450 
452  inline double getValuePureAt(bool show = false) const
453  {
454  double VPP=0.0;
455  double check_VPP = round(VPP_radius*1000);
456  if(check_VPP != 0)
457  VPP = -phase_shift*(1-exp(-precomputed.u2/(2*pow(VPP_radius,2.0))));
458  else
459  VPP = 0;
460  double argument = VPP + K1 * precomputed.deltaf * precomputed.u2 + K2 *precomputed.u4;
461  double sine_part;
462  double cosine_part;
463  //sincos(argument,&sine_part, &cosine_part);
464  sine_part = sin(argument);
465  cosine_part = cos(argument);// OK
466  double Eespr = exp(-K3 * precomputed.u4); // OK
467  //CO: double Eispr=exp(-K4*u4); // OK
468  double EdeltaF = bessj0(K5 * precomputed.u2); // OK
469  double EdeltaR = SINC(precomputed.u * DeltaR); // OK
470  double aux=(K7 * precomputed.u2 * precomputed.u + precomputed.deltaf * precomputed.u);
471  double Ealpha = exp(-K6 * aux * aux); // OK
472  // CO: double E=Eespr*Eispr*EdeltaF*EdeltaR*Ealpha;
473  double E = Eespr * EdeltaF * EdeltaR * Ealpha + envR0+envR1*precomputed.u+envR2*precomputed.u2;
474  if (E < 0)
475  E = 0;
476 
477  if (show)
478  {
479  std::cout << " Deltaf=" << precomputed.deltaf << std::endl;
480  std::cout << " u,u2,u4=" << precomputed.u << " " << precomputed.u2
481  << " " << precomputed.u4 << std::endl;
482  std::cout << " K1,K2,argument=" << K1 << " " << K2 << " " << argument << std::endl;
483  std::cout << " K3,Eespr=" << K3 << " " << Eespr << std::endl;
484  //std::cout << " K4,Eispr=" << K4 << " " << /*Eispr <<*/ std::endl;
485  std::cout << " K5,EdeltaF=" << K5 << " " << EdeltaF << std::endl;
486  std::cout << " EdeltaR=" << EdeltaR << std::endl;
487  std::cout << " K6,K7,Ealpha=" << K6 << " " << K7 << " " << Ealpha
488  << std::endl;
489  std::cout << " Total atenuation(E)= " << E << std::endl;
490  std::cout << " K,Q0,base_line=" << K << "," << Q0 << "," << base_line << std::endl;
491  std::cout << " VPP=" << VPP << " VPPRadius=" << VPP_radius << " phase shift=" << phase_shift << std::endl;
492  std::cout << " CTF="
493  << -K*(Ksin*sine_part - Kcos*cosine_part)*E << std::endl;
494  }
495  return -K*(Ksin*sine_part - Kcos*cosine_part)*E;
496  }
497 
499  inline double getValuePureNoKAt() const
500  {
501  return K*getValuePureAt();
502  }
503 
505  //#define DEBUG
506  inline double getValueNoiseAt(bool show = false) const
507  {
508  double c = Gc1;
509  double sigmaG1 = sigma1;
510 
511  double c2 = Gc2;
512  double sigmaG2 = sigma2;
513 
514  if(show)
515  {
516  std::cout
517  << " sq=" << sq << "\n"
518  << " c=" << c << "\n"
519  << " sigma=" << sigmaG1 << "\n"
520  << " c2=" << c2 << "\n"
521  << " sigma2=" << sigmaG2 << "\n"
522  << " gaussian_K=" << gaussian_K << "\n"
523  << " sqrt_K=" << sqrt_K << "\n"
524  << " base_line=" << base_line << "\n";
525  std::cout << " u=" << precomputed.u << "u_sqrt=" << precomputed.u_sqrt <<") CTFnoise="
526  << base_line +
527  gaussian_K*exp(-sigmaG1*(precomputed.u - c)*(precomputed.u - c)) +
528  sqrt_K*exp(-sq*sqrt(precomputed.u)) -
529  gaussian_K2*exp(-sigmaG2*(precomputed.u - c2)*(precomputed.u - c2)) << std::endl;
530  }
531 
532  double aux=precomputed.u - c;
533  double aux2=precomputed.u - c2;
534  return base_line +
535  gaussian_K*exp(-sigmaG1*aux*aux) +
536  sqrt_K*exp(-sq*precomputed.u_sqrt) -
537  gaussian_K2*exp(-sigmaG2*aux2*aux2)+
538  bgR1*precomputed.u+bgR2*precomputed.u2+bgR3*precomputed.u3;
539  }
540 
542  double getValuePureWithoutDampingAt(bool show = false) const
543  {
544  double VPP;
545  double check_VPP = round(VPP_radius*1000);
546  if(check_VPP != 0)
547  VPP = -phase_shift*(1-exp(-precomputed.u2/(2*pow(VPP_radius,2.0))));
548  else
549  VPP = 0;
550  double argument = VPP + K1 * precomputed.deltaf * precomputed.u2 + K2 * precomputed.u4;
551  double sine_part;
552  double cosine_part;
553  //sincos(argument,&sine_part,&cosine_part);
554  sine_part = sin(argument);
555  cosine_part = cos(argument);// OK
556 
557  if (show)
558  {
559  std::cout << " Deltaf=" << precomputed.deltaf << std::endl;
560  std::cout << " u,u2,u4=" << precomputed.u << " " << precomputed.u2
561  << " " << precomputed.u4 << std::endl;
562  std::cout << " K1,K2,sin=" << K1 << " " << K2 << " "
563  << std::endl;
564  std::cout << " Q0=" << Q0 << std::endl;
565  std::cout << " VPP=" << VPP << std::endl;
566  std::cout << " CTF without damping="
567  << -(Ksin*sine_part - Kcos*cosine_part) << std::endl;
568  }
569  return -(Ksin*sine_part - Kcos*cosine_part);
570  }
571 
572  void getSineAndCosineParts(double &sine_part, double &cosine_part, double E, double u2, double deltaf, bool show) const
573  {
574  double u = sqrt(u2);
575  double u4 = u2 * u2;
576  double VPP;
577  double check_VPP = round(VPP_radius*1000);
578  if(check_VPP != 0)
579  VPP = -phase_shift*(1-exp(-u2/(2*pow(VPP_radius,2.0))));
580  else
581  VPP = 0;
582  double argument = VPP + K1 * deltaf * u2 + K2 * u4;
583  //sincos(argument,&sine_part, &cosine_part); // OK
584  sine_part = sin(argument);
585  cosine_part = cos(argument);
586  double Eespr = exp(-K3 * u4); // OK
587  //CO: double Eispr=exp(-K4*u4); // OK
588  double EdeltaF = bessj0(K5 * u2); // OK
589  double EdeltaR = sinc(u * DeltaR); // OK
590  double Ealpha = exp(-K6 * (K7 * u2 * u + deltaf * u) * (K7 * u2 * u + deltaf * u)); // OK
591  // CO: double E=Eespr*Eispr*EdeltaF*EdeltaR*Ealpha;
592  E = Eespr * EdeltaF * EdeltaR * Ealpha+envR0+envR1*precomputed.u+envR2*precomputed.u2;
593  if (E < 0)
594  E = 0;
595  if (show)
596  {
597  std::cout << " Deltaf=" << deltaf << std::endl;
598  std::cout << " u,u2,u4=" << u << " " << u2 << " " << u4 << std::endl;
599  std::cout << " K1,K2,sin=" << K1 << " " << K2 << " "
600  << sine_part << std::endl;
601  std::cout << " K3,Eespr=" << K3 << " " << Eespr << std::endl;
602  //std::cout << " K4,Eispr=" << K4 << " " << /*Eispr*/ << std::endl;
603  std::cout << " K5,EdeltaF=" << K5 << " " << EdeltaF << std::endl;
604  std::cout << " EdeltaR=" << EdeltaR << std::endl;
605  std::cout << " K6,K7,Ealpha=" << K6 << " " << K7 << " " << Ealpha
606  << std::endl;
607  std::cout << " Total atenuation(E)= " << E << std::endl;
608  std::cout << " K,Q0,base_line=" << K << "," << Q0 << "," << base_line << std::endl;
609  }
610  }
611 
613  inline double getValuePureNoPrecomputedAt(double X, bool show = false) const
614  {
615  double u2 = X * X;
616  double deltaf = Defocus;
617  double sine_part = 0;
618  double cosine_part = 0;
619  double E = 0;
620  getSineAndCosineParts(sine_part, cosine_part, E, u2, deltaf, show);
621  if (show)
622  {
623  std::cout << " (X)=(" << X << ") CTF="
624  << -K*(Ksin*sine_part - Kcos*cosine_part)*E + base_line << std::endl;
625  }
626  return -K*(Ksin*sine_part - Kcos*cosine_part)*E;
627  }
628 
629 
633  void lookFor(int n, const Matrix1D<double> &u, Matrix1D<double> &freq, int iwhat=0);
634 
636  void applyCTF(MultidimArray < std::complex<double> > &FFTI, const MultidimArray<double> &I, double Ts, bool absPhase=false);
637 
639  void applyCTF(MultidimArray <double> &I, double Ts, bool absPhase=false);
640 
642  void correctPhase(MultidimArray < std::complex<double> > &FFTI, const MultidimArray<double> &I, double Ts);
643 
645  void correctPhase(MultidimArray<double> &I, double Ts);
646 
649  template <class T1, class T2>
650  void generateCTF(const MultidimArray<T1> &sample_image, MultidimArray <T2> &CTF, double Ts=-1)
651  {
652  if ( ZSIZE(sample_image) > 1 )
653  REPORT_ERROR(ERR_MULTIDIM_DIM,"ERROR: Generate_CTF only works with 2D sample images, not 3D.");
654  generateCTF(YSIZE(sample_image), XSIZE(sample_image), CTF, Ts);
655  STARTINGX(CTF) = STARTINGX(sample_image);
656  STARTINGY(CTF) = STARTINGY(sample_image);
657  }
658 
664  void getProfile(double fmax, int nsamples, MultidimArray<double> &profiles);
665 
670  void getAverageProfile(double fmax, int nsamples, MultidimArray<double> &profiles);
671 
672  //#define DEBUG
673 
675  template <class T>
676  double initCTF(int Ydim, int Xdim, MultidimArray< T > &CTF, double Ts=-1) const
677  {
678  CTF.resizeNoCopy(Ydim, Xdim);
679  if (Ts<0)
680  Ts=Tm;
681  #ifdef DEBUG
682  std::cout << "CTF:\n" << *this << std::endl;
683  #endif
684  double iTs=1.0/Ts;
685  return iTs;
686  }
687 
689  template <class T>
690  void generateCTF(int Ydim, int Xdim, MultidimArray < T > &CTF, double Ts=-1)
691  {
692  double iTs = initCTF(Ydim, Xdim, CTF, Ts);
693  for (int i=0; i<Ydim; ++i)
694  {
695  double wy;
696  FFT_IDX2DIGFREQ(i, YSIZE(CTF), wy);
697  double fy=wy*iTs;
698  for (int j=0; j<Xdim; ++j)
699  {
700  double wx;
701  FFT_IDX2DIGFREQ(j, XSIZE(CTF), wx);
702  double fx=wx*iTs;
703  precomputeValues(fx);
704  A2D_ELEM(CTF, i, j) = (T) getValueAt();
705  #ifdef DEBUG
706  if (i == 0)
707  std::cout << i << " " << j << " " << YY(freq) << " " << XX(freq)
708  << " " << CTF(i, j) << std::endl;
709  #endif
710  }
711  }
712  }
713  #undef DEBUG
714 
715  template <class T>
716  void generateCTFWithoutDamping(int Ydim, int Xdim, MultidimArray < T > &CTF, double Ts=-1)
717  {
718  double iTs = initCTF(Ydim, Xdim, CTF, Ts=-1);
719  for (int i=0; i<Ydim; ++i)
720  {
721  double wy;
722  FFT_IDX2DIGFREQ(i, YSIZE(CTF), wy);
723  double fy=wy*iTs;
724  for (int j=0; j<Xdim; ++j)
725  {
726  double wx;
727  FFT_IDX2DIGFREQ(j, XSIZE(CTF), wx);
728  double fx=wx*iTs;
729  precomputeValues(fx);
730  A2D_ELEM(CTF, i, j) = (T) getValuePureWithoutDampingAt();
731  #ifdef DEBUG
732  if (i == 0)
733  std::cout << i << " " << j << " " << YY(freq) << " " << XX(freq)
734  << " " << CTF(i, j) << std::endl;
735  #endif
736  }
737  }
738  }
739  #undef DEBUG
740 
744  bool hasPhysicalMeaning();
745 
747  void forcePhysicalMeaning();
748 
749 };
750 
753 void generateCTFImageWith2CTFs(const MetaData &MD1, const MetaData &MD2, int Xdim, MultidimArray<double> &imgOut);
754 
756 double errorBetween2CTFs( MetaData &MD1,
757  MetaData &MD2,
758  size_t dim,
759  double minFreq=0.05,
760  double maxFreq=0.25);
764 double errorMaxFreqCTFs( MetaData &MD1,
765  double phaseRad=HALFPI);
766 
770 double errorMaxFreqCTFs2D( MetaData &MD1,
771  MetaData &MD2,
772  size_t xDim=256,
773  double phaseRad=HALFPI);
774 
777 void generatePSDCTFImage(MultidimArray<double> &img, const MetaData &MD);
778 
779 
781 
782 class CTFDescription: public CTFDescription1D
783 {
784 public:
785  // Electron wavelength (Amstrongs)
786  //double lambda;
787  // Squared frequency associated to the aperture
788  // double ua2;
789  // Different constants
790  //double K1;
791  //double K2;
792  //double K3;
793  //double K4;
794  //double K5;
795  //double K6;
796  //double K7;
797  //double Ksin;
798  //double Kcos;
799  // Azimuthal angle in radians
800  double rad_azimuth;
801  // defocus_average = -(defocus_u + defocus_v)/2
802  double defocus_average;
803  // defocus_deviation = -(defocus_u - defocus_v)/2
804  double defocus_deviation;
805  // Gaussian angle in radians
806  double rad_gaussian;
807  // Second Gaussian angle in radians
808  double rad_gaussian2;
809  // Sqrt angle in radians
810  double rad_sqrt;
813  //double D;
814  // Precomputed values
815  //PrecomputedForCTF precomputed;
816  // Image of precomputed values
817  //std::vector<PrecomputedForCTF> precomputedImage;
818  // Xdim size of the image
819  //int precomputedImageXdim;
820 public:
822  //double K;
824  //double Tm;
826  //double kV;
828  double DeltafU;
830  double DeltafV;
832  double azimuthal_angle;
833  // Radius of the aperture (in micras)
834  // double aperture;
836  //double Cs;
838  //double Ca;
841  //double espr;
843  //double ispr;
845  //double alpha;
847  //double DeltaF;
849  //double DeltaR;
851  //double Q0;
853  //double x0;
855  //double xF;
857  //double y0;
859  //double yF;
861  //bool isLocalCTF;
863  //bool enable_CTFnoise;
865  //bool enable_CTF;
867  //double base_line;
869  //double gaussian_K;
871  double sigmaU;
873  double sigmaV;
875  double cU;
877  double cV;
879  double gaussian_angle;
881  //double sqrt_K;
883  double sqU;
885  double sqV;
887  double sqrt_angle;
889  double sigmaU2;
891  double sigmaV2;
893  double cU2;
895  double cV2;
897  double gaussian_angle2;
898 
899  CTFDescription()
900  {
901  clear();
902  }
903 
904  CTFDescription(CTFDescription1D copy)
905  {
906  DeltafU = DeltafV = copy.Defocus;
907  Ca = copy.Ca;
908  Cs = copy.Cs;
909  DeltaF = copy.DeltaF;
910  DeltaR = copy.DeltaR;
911  Q0 = copy.Q0;
912  cU = cV = copy.Gc1;
913  cU2 = cV2 = copy.Gc2;
914  base_line = copy.base_line;
915  gaussian_K = copy.gaussian_K;
916  gaussian_K2 = copy.gaussian_K2;
917  sigmaU = sigmaV = copy.sigma1;
918  sigmaU2 = sigmaV2 = copy.sigma2;
919  sqU = sqV = copy.sq;
920  kV = copy.kV;
921  bgR1 = copy.bgR1;
922  bgR2 = copy.bgR2;
923  bgR3 = copy.bgR3;
924  envR0 = copy.envR0;
925  envR1 = copy.envR1;
926  envR2 = copy.envR2;
927  espr = copy.espr;
928  ispr = copy.ispr;
929  alpha = copy.alpha;
930  sqrt_K = copy.sqrt_K;
931  K = copy.K;
932  Tm = copy.Tm;
933  azimuthal_angle = 0;
934  gaussian_angle = 0;
935  gaussian_angle2 = 0;
936  sqrt_angle = 0;
937  enable_CTFnoise = copy.enable_CTFnoise;
938  phase_shift = copy.phase_shift;
939  VPP_radius = copy.VPP_radius;
940  }
941 
949  void read(const FileName &fn, bool disable_if_not_K = true);
959  void readFromMetadataRow(const MetaData &MD, size_t id, bool disable_if_not_K=true);
960 
962  void readFromMdRow(const MDRow &row, bool disable_if_not_K=true);
963 
965  void setRow(MDRow &row) const;
966 
969  void write(const FileName &fn);
970 
972  static void defineParams(XmippProgram * program);
973 
975  void readParams(XmippProgram * program);
976 
978  friend std::ostream & operator << (std::ostream &out, const CTFDescription &ctf);
979 
981  void clear();
982 
984  void clearNoise();
985 
987  void clearPureCtf();
988 
994  inline void changeSamplingRate(double newTm)
995  {
996  Tm=newTm;
997  }
999  void produceSideInfo();
1000 
1002  void precomputeValues(double X, double Y)
1003  {
1004  precomputed.ang=atan2(Y, X);
1005  precomputed.u2 = X * X + Y * Y;
1006  precomputed.u = sqrt(precomputed.u2);
1007  precomputed.u3 = precomputed.u2 * precomputed.u;
1008  precomputed.u4 = precomputed.u2 * precomputed.u2;
1009  precomputed.u_sqrt = sqrt(precomputed.u);
1010 
1011  if (fabs(X) < XMIPP_EQUAL_ACCURACY &&
1012  fabs(Y) < XMIPP_EQUAL_ACCURACY)
1013  precomputed.deltaf=0;
1014  else
1015  {
1016  double ellipsoid_ang = precomputed.ang - rad_azimuth;
1017  /*
1018  * For a derivation of this formulae confer
1019  * Principles of Electron Optics page 1380
1020  * in particular term defocus and twofold axial astigmatism
1021  * take into account that a1 and a2 are the coefficient
1022  * of the zernike polynomials difference of defocus at 0
1023  * and at 45 degrees. In this case a2=0
1024  */
1025  double cos_ellipsoid_ang_2 = cos(2*ellipsoid_ang);
1026  precomputed.deltaf= (defocus_average + defocus_deviation*cos_ellipsoid_ang_2);
1027  }
1028  }
1029 
1031  void precomputeValues(const MultidimArray<double> &cont_x_freq,
1032  const MultidimArray<double> &cont_y_freq);
1033 
1035  void precomputeValues(int i, int j)
1036  {
1037 
1038  precomputed=precomputedImage[i*precomputedImageXdim+j];
1039 
1040  if (precomputed.deltaf==-1)
1041  {
1042  double ellipsoid_ang = precomputed.ang - rad_azimuth;
1043  double cos_ellipsoid_ang_2 = cos(2*ellipsoid_ang);
1044  precomputed.deltaf= (defocus_average + defocus_deviation*cos_ellipsoid_ang_2);
1045 
1046  }
1047  }
1048 
1050  double getValueAt(bool show = false) const
1051  {
1052  double pure_CTF = enable_CTF ? getValuePureAt(show) : 0;
1053  return enable_CTFnoise ? sqrt(pure_CTF*pure_CTF + getValueNoiseAt(show)) : pure_CTF;
1054  }
1055 
1057  double getValueArgument(bool show = false) const
1058  {
1059  return K1 * precomputed.deltaf * precomputed.u2 + K2 * precomputed.u4;
1060  }
1061 
1063  inline double getPhaseAt() const
1064  {
1065  return K1 * precomputed.deltaf * precomputed.u2 + K2 *precomputed.u4;
1066  }
1067 
1069  inline double getValuePureNoPrecomputedAtxy(double X, double Y, bool show = false) const
1070  {
1071  double u2 = X * X + Y * Y;
1072  //if (u2>=ua2) return 0;
1073  double deltaf = getDeltafNoPrecomputed(X, Y);
1074  double sine_part = 0;
1075  double cosine_part = 0;
1076  double E = 0;
1077  getSineAndCosineParts(sine_part, cosine_part, E, u2, deltaf, show);
1078  if (show)
1079  {
1080  std::cout << " (X,Y)=(" << X << "," << Y << ") CTF="
1081  << -K*(Ksin*sine_part - Kcos*cosine_part)*E + base_line << std::endl;
1082  }
1083  return -K*(Ksin*sine_part - Kcos*cosine_part)*E;
1084  }
1085 
1087  inline double getValuePureNoDampingNoPrecomputedAt(double X, double Y) const
1088  {
1089  double u2 = X * X + Y * Y;
1090  //if(u2 > freq_max) return 0;
1091  double u4 = u2 * u2;
1092  double deltaf = getDeltafNoPrecomputed(X, Y);
1093  double VPP;
1094  double check_VPP = round(VPP_radius*1000);
1095  if(check_VPP != 0)
1096  VPP = -phase_shift*(1-exp(-u2/(2*pow(VPP_radius,2.0))));
1097  else
1098  VPP = 0;
1099  double argument = VPP + K1 * deltaf * u2 + K2 * u4;
1100  double sine_part;
1101  double cosine_part;
1102  //sincos(argument,&sine_part, &cosine_part); // OK
1103  sine_part = sin(argument);
1104  cosine_part = cos(argument);
1105  return -(Ksin*sine_part - Kcos*cosine_part);
1106  }
1107 
1109  double getDeltafNoPrecomputed(double X, double Y) const
1110  {
1111  if (fabs(X) < XMIPP_EQUAL_ACCURACY &&
1112  fabs(Y) < XMIPP_EQUAL_ACCURACY)
1113  return 0;
1114  double ellipsoid_ang = atan2(Y, X) - rad_azimuth;
1115  double cos_ellipsoid_ang_2 = cos(2*ellipsoid_ang);
1116  return(defocus_average + defocus_deviation*cos_ellipsoid_ang_2);
1117  }
1118 
1120  //#define DEBUG
1121  inline double getValueNoiseAt(bool show = false) const
1122  {
1123  double ellipsoid_ang = precomputed.ang - rad_gaussian;
1124  double ellipsoid_ang2 = precomputed.ang - rad_gaussian2;
1125  double ellipsoid_sqrt_ang = precomputed.ang - rad_sqrt;
1126  double cos_sqrt_ang = cos(ellipsoid_sqrt_ang);
1127  double cos_sqrt_ang_2 = cos_sqrt_ang*cos_sqrt_ang;
1128  double sin_sqrt_ang_2 = 1.0-cos_sqrt_ang_2;
1129  double sq = sqrt(sqU * sqU * cos_sqrt_ang_2 + sqV * sqV * sin_sqrt_ang_2);
1130 
1131  double cos_ang = cos(ellipsoid_ang);
1132  double cos_ang_2 = cos_ang*cos_ang;
1133  double sin_ang_2 = 1.0-cos_ang_2;
1134  double c = sqrt(cU * cU * cos_ang_2 + cV * cV * sin_ang_2);
1135  double sigma = sqrt(sigmaU * sigmaU * cos_ang_2 + sigmaV * sigmaV * sin_ang_2);
1136 
1137  double cos_ang2 = cos(ellipsoid_ang2);
1138  double cos_ang2_2 = cos_ang2*cos_ang2;
1139  double sin_ang2_2 = 1.0-cos_ang2_2;
1140  double c2 = sqrt(cU2 * cU2 * cos_ang2_2 + cV2 * cV2 * sin_ang2_2);
1141  double sigma2 = sqrt(sigmaU2 * sigmaU2 * cos_ang2_2 + sigmaV2 * sigmaV2 * sin_ang2_2);
1142 
1143  if(show)
1144  {
1145  std::cout << " ellipsoid_ang=" << RAD2DEG(ellipsoid_ang) << std::endl
1146  << " ellipsoid_sqrt_ang=" << RAD2DEG(ellipsoid_sqrt_ang) << std::endl
1147  << " sq=" << sq << "\n"
1148  << " c=" << c << "\n"
1149  << " sigma=" << sigma << "\n"
1150  << " ellipsoid_ang2=" << RAD2DEG(ellipsoid_ang2) << std::endl
1151  << " cU2=" << cU2 << " cV2=" << cV2 << "\n"
1152  << " cos_ang2=" << cos_ang2 << " sin_ang2_2=" << sin_ang2_2 << "\n"
1153  << " c2=" << c2 << "\n"
1154  << " sigmaU2=" << sigma2 << "\n"
1155  << " gaussian_K=" << gaussian_K << "\n"
1156  << " sqrt_K=" << sqrt_K << "\n"
1157  << " base_line=" << base_line << "\n";
1158  std::cout << " u=" << precomputed.u << "u_sqrt=" << precomputed.u_sqrt <<") CTFnoise="
1159  << base_line +
1160  gaussian_K*exp(-sigma*(precomputed.u - c)*(precomputed.u - c)) +
1161  sqrt_K*exp(-sq*sqrt(precomputed.u)) -
1162  gaussian_K2*exp(-sigma2*(precomputed.u - c2)*(precomputed.u - c2)) << std::endl;
1163  }
1164 
1165  double aux=precomputed.u - c;
1166  double aux2=precomputed.u - c2;
1167  return base_line +
1168  gaussian_K*exp(-sigma*aux*aux) +
1169  sqrt_K*exp(-sq*precomputed.u_sqrt) -
1170  gaussian_K2*exp(-sigma2*aux2*aux2)+
1171  bgR1*precomputed.u+bgR2*precomputed.u2+bgR3*precomputed.u3;
1172  }
1173 
1177  void lookFor(int n, const Matrix1D<double> &u, Matrix1D<double> &freq, int iwhat=0);
1178 
1180  void applyCTF(MultidimArray < std::complex<double> > &FFTI, const MultidimArray<double> &I, double Ts, bool absPhase=false);
1181 
1183  void applyCTF(MultidimArray <double> &I, double Ts, bool absPhase=false);
1184 
1186  void correctPhase(MultidimArray < std::complex<double> > &FFTI, const MultidimArray<double> &I, double Ts);
1187 
1189  void correctPhase(MultidimArray<double> &I, double Ts);
1190 
1193  template <class T1, class T2>
1194  void generateCTF(const MultidimArray<T1> &sample_image, MultidimArray <T2> &CTF, double Ts=-1)
1195  {
1196  if ( ZSIZE(sample_image) > 1 )
1197  REPORT_ERROR(ERR_MULTIDIM_DIM,"ERROR: Generate_CTF only works with 2D sample images, not 3D.");
1198  generateCTF(YSIZE(sample_image), XSIZE(sample_image), CTF, Ts);
1199  STARTINGX(CTF) = STARTINGX(sample_image);
1200  STARTINGY(CTF) = STARTINGY(sample_image);
1201  }
1202 
1208  void getProfile(double angle, double fmax, int nsamples, MultidimArray<double> &profiles);
1209 
1214  void getAverageProfile(double fmax, int nsamples, MultidimArray<double> &profiles);
1215 
1216 //#define DEBUG
1218  template <class T>
1219  void generateCTF(int Ydim, int Xdim, MultidimArray < T > &CTF, double Ts=-1)
1220  {
1221  double iTs = initCTF(Ydim, Xdim, CTF, Ts);
1222  for (int i=0; i<Ydim; ++i)
1223  {
1224  double wy;
1225  FFT_IDX2DIGFREQ(i, YSIZE(CTF), wy);
1226  double fy=wy*iTs;
1227  for (int j=0; j<Xdim; ++j)
1228  {
1229  double wx;
1230  FFT_IDX2DIGFREQ(j, XSIZE(CTF), wx);
1231  double fx=wx*iTs;
1232  precomputeValues(fx, fy);
1233  A2D_ELEM(CTF, i, j) = (T) getValueAt();
1234  #ifdef DEBUG
1235  if (i == 0)
1236  std::cout << i << " " << j << " " << fy << " " << fx
1237  << " " << CTF(i, j) << std::endl;
1238  #endif
1239  }
1240  }
1241  }
1242  #undef DEBUG
1243 
1244  template <class T>
1245  void generateCTFWithoutDamping(int Ydim, int Xdim, MultidimArray < T > &CTF, double Ts=-1)
1246  {
1247  double iTs = initCTF(Ydim, Xdim, CTF, Ts);
1248  for (int i=0; i<Ydim; ++i)
1249  {
1250  double wy;
1251  FFT_IDX2DIGFREQ(i, YSIZE(CTF), wy);
1252  double fy=wy*iTs;
1253  for (int j=0; j<Xdim; ++j)
1254  {
1255  double wx;
1256  FFT_IDX2DIGFREQ(j, XSIZE(CTF), wx);
1257  double fx=wx*iTs;
1258  precomputeValues(fx, fy);
1259  A2D_ELEM(CTF, i, j) = (T) getValuePureWithoutDampingAt();
1260  #ifdef DEBUG
1261  if (i == 0)
1262  std::cout << i << " " << j << " " << YY(freq) << " " << XX(freq)
1263  << " " << CTF(i, j) << std::endl;
1264  #endif
1265  }
1266  }
1267  }
1268  #undef DEBUG
1269 
1270  template <class T>
1271  void generateEnvelope(int Ydim, int Xdim, MultidimArray < T > &CTF, double Ts=-1)
1272  {
1273  double iTs = initCTF(Ydim, Xdim, CTF, Ts);
1274  for (int i=0; i<Ydim; ++i)
1275  {
1276  double wy;
1277  FFT_IDX2DIGFREQ(i, YSIZE(CTF), wy)
1278  double fy=wy*iTs;
1279  for (int j=0; j<Xdim; ++j)
1280  {
1281  double wx;
1282  FFT_IDX2DIGFREQ(j, XSIZE(CTF), wx)
1283  double fx=wx*iTs;
1284  precomputeValues(fx, fy);
1285  A2D_ELEM(CTF, i, j) = (T) -getValueDampingAt();
1286  #ifdef DEBUG
1287  if (i == 0)
1288  std::cout << i << " " << j << " " << YY(freq) << " " << XX(freq)
1289  << " " << CTF(i, j) << std::endl;
1290  #endif
1291  }
1292  }
1293  }
1294 
1298  bool hasPhysicalMeaning();
1299 
1301  void forcePhysicalMeaning();
1302 };
1303 
1304 #endif
CTFDescription::sigmaU
double sigmaU
Spherical aberration (in milimeters). Typical value 5.6.
Definition: ctf.h:871
CTFDescription1D::precomputedImageXdim
int precomputedImageXdim
Definition: ctf.h:235
containsCTFBasicLabels
bool containsCTFBasicLabels(const MetaData &md)
Definition: ctf.cpp:33
MDL_CTF_K
CTF gain.
Definition: metadata_label.h:163
groupCTFMetaData
void groupCTFMetaData(const MetaDataDb &imgMd, MetaDataDb &ctfMd, std::vector< MDLabel > &groupbyLabels)
Definition: ctf.cpp:41
YSIZE
#define YSIZE(v)
Definition: multidim_array_base.h:95
A2D_ELEM
#define A2D_ELEM(v, i, j)
Definition: multidim_array_base.h:354
CTFDescription1D::bgR1
double bgR1
Definition: ctf.h:295
CTF_ALL_LABELS
const MDLabel CTF_ALL_LABELS[]
Definition: ctf.h:48
MDL_CTF_DEFOCUSU
Defocus U (Angstroms)
Definition: metadata_label.h:135
CTFDescription1D::lambda
double lambda
Definition: ctf.h:214
CTFDescription1D::K1
double K1
Definition: ctf.h:218
PrecomputedForCTF::u
double u
Definition: ctf.h:86
numerical_recipes.h
CTFDescription1D::bgR2
double bgR2
Definition: ctf.h:296
CTFDescription::generateEnvelope
void generateEnvelope(int Ydim, int Xdim, MultidimArray< T > &CTF, double Ts=-1)
Definition: ctf.h:1271
CTFDescription::getDeltafNoPrecomputed
double getDeltafNoPrecomputed(double X, double Y) const
Deltaf at a given direction.
Definition: ctf.h:1109
cif::pdb::write
void write(std::ostream &os, const datablock &db)
Definition: cif2pdb.cpp:3747
CTFDescription::gaussian_angle2
double gaussian_angle2
Second Gaussian angle.
Definition: ctf.h:897
CTFDescription::rad_sqrt
double rad_sqrt
Definition: ctf.h:810
REPORT_ERROR
#define REPORT_ERROR(nerr, ErrormMsg)
Definition: xmipp_error.h:211
CTFDescription1D::envR1
double envR1
Definition: ctf.h:300
MDL_CTF_DEFOCUS_ANGLE
Defocus angle (degrees)
Definition: metadata_label.h:151
xDim
#define xDim
Definition: projection.cpp:41
CTFDescription1D::getValueAt
double getValueAt(bool show=false) const
Compute CTF at (U,V). Continuous frequencies.
Definition: ctf.h:417
c
doublereal * c
Definition: numerical_recipes.cpp:2230
MultidimArray::resizeNoCopy
void resizeNoCopy(const MultidimArray< T1 > &v)
Definition: multidim_array.h:504
CTFDescription1D::precomputeValues
void precomputeValues(double X)
Precompute values for a given frequency.
Definition: ctf.h:376
CTFDescription1D::precomputeValues
void precomputeValues(int i)
Precompute values for a given frequency.
Definition: ctf.h:407
sinc
double sinc(double x)
Definition: numerical_recipes.h:81
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
generatePSDCTFImage
void generatePSDCTFImage(MultidimArray< double > &img, const MetaData &MD)
Definition: ctf.cpp:330
CTFDescription1D::envR2
double envR2
Definition: ctf.h:301
CTFDescription::CTFDescription
CTFDescription(CTFDescription1D copy)
Definition: ctf.h:904
CTFDescription::DeltafV
double DeltafV
Defocus in V (in Angstroms). Negative values are underfocused.
Definition: ctf.h:830
metadata_db.h
CTFDescription::defocus_average
double defocus_average
Definition: ctf.h:802
CTFDescription1D::generateCTFWithoutDamping
void generateCTFWithoutDamping(int Ydim, int Xdim, MultidimArray< T > &CTF, double Ts=-1)
Definition: ctf.h:716
MDL_CTF_CONVERGENCE_CONE
Convergence cone.
Definition: metadata_label.h:159
CTFDescription1D::espr
double espr
Definition: ctf.h:251
CTFDescription::cU2
double cU2
Second Gaussian center for U.
Definition: ctf.h:893
PrecomputedForCTF::u4
double u4
Definition: ctf.h:89
CTFDescription::defocus_deviation
double defocus_deviation
Definition: ctf.h:804
CTFDescription::rad_gaussian
double rad_gaussian
Definition: ctf.h:806
CTFDescription::cU
double cU
Gaussian center for U.
Definition: ctf.h:875
CTFDescription::changeSamplingRate
void changeSamplingRate(double newTm)
Definition: ctf.h:994
CTFDescription1D::phase_shift
double phase_shift
Definition: ctf.h:305
CTFDescription::rad_azimuth
double rad_azimuth
Definition: ctf.h:800
CTFDescription1D::getValueNoiseAt
double getValueNoiseAt(bool show=false) const
Compute noise at (X,Y). Continuous frequencies, notice it is squared.
Definition: ctf.h:506
CTFDescription1D::CTFDescription1D
CTFDescription1D()
Definition: ctf.h:309
CTFDescription1D::envR0
double envR0
Definition: ctf.h:299
CTFDescription1D::gaussian_K
double gaussian_K
Gain for the gaussian term.
Definition: ctf.h:279
errorMaxFreqCTFs2D
double errorMaxFreqCTFs2D(MetaData &MD1, MetaData &MD2, size_t xDim=256, double phaseRad=HALFPI)
Definition: ctf.cpp:214
PrecomputedForCTF::u2
double u2
Definition: ctf.h:87
testing::operator<<
std::ostream & operator<<(std::ostream &os, const Message &sb)
Definition: gtest-message.h:200
CTFDescription1D::sq
double sq
Sqrt width.
Definition: ctf.h:287
MDL_CTF_CS
Spherical aberration.
Definition: metadata_label.h:152
STARTINGX
#define STARTINGX(v)
Definition: multidim_array_base.h:51
CTF_BASIC_LABELS
const MDLabel CTF_BASIC_LABELS[]
Definition: ctf.h:42
CTFDescription1D::DeltaF
double DeltaF
Longitudinal mechanical displacement (ansgtrom). Typical value 100.
Definition: ctf.h:257
CTFDescription::sqrt_angle
double sqrt_angle
Sqrt angle.
Definition: ctf.h:887
i
#define i
Definition: numerical_recipes.cpp:2493
MDL_CTF_VOLTAGE
Microscope voltage (kV)
Definition: metadata_label.h:133
CTFDescription::getValuePureNoDampingNoPrecomputedAt
double getValuePureNoDampingNoPrecomputedAt(double X, double Y) const
Compute CTF pure at (U,V). Continuous frequencies.
Definition: ctf.h:1087
CTFDescription::DeltafU
double DeltafU
Global gain. By default, 1.
Definition: ctf.h:828
CTFDescription1D::enable_CTF
bool enable_CTF
Enable CTF part.
Definition: ctf.h:275
STARTINGY
#define STARTINGY(v)
Definition: multidim_array_base.h:59
CTFDescription::azimuthal_angle
double azimuthal_angle
Azimuthal angle (between X and U) in degrees.
Definition: ctf.h:832
CTFDescription1D::Gc2
double Gc2
Second Gaussian center.
Definition: ctf.h:293
CTFDescription1D::sigma2
double sigma2
Second Gaussian width.
Definition: ctf.h:291
CTF_BASIC_LABELS_SIZE
const int CTF_BASIC_LABELS_SIZE
Definition: ctf.h:41
CTFDescription::getValueArgument
double getValueArgument(bool show=false) const
Compute pure CTF without damping at (U,V). Continuous frequencies.
Definition: ctf.h:1057
CTFDescription1D::base_line
double base_line
Global base_line.
Definition: ctf.h:277
PrecomputedForCTF::deltaf
double deltaf
Definition: ctf.h:91
CTFDescription1D::K4
double K4
Definition: ctf.h:221
CTFDescription1D::K2
double K2
Definition: ctf.h:219
CTFDescription1D::K6
double K6
Definition: ctf.h:223
CTFDescription::sqU
double sqU
Gain for the square root term.
Definition: ctf.h:883
CTFDescription1D::Gc1
double Gc1
Gaussian center.
Definition: ctf.h:283
generateCTFImageWith2CTFs
void generateCTFImageWith2CTFs(const MetaData &MD1, const MetaData &MD2, int Xdim, MultidimArray< double > &imgOut)
Definition: ctf.cpp:67
CTFDescription1D::K7
double K7
Definition: ctf.h:224
XX
#define XX(v)
Definition: matrix1d.h:85
CTFDescription::getValueAt
double getValueAt(bool show=false) const
Compute CTF at (U,V). Continuous frequencies.
Definition: ctf.h:1050
CTFDescription1D::Ksin
double Ksin
Definition: ctf.h:225
CTFDescription1D::generateCTF
void generateCTF(const MultidimArray< T1 > &sample_image, MultidimArray< T2 > &CTF, double Ts=-1)
Definition: ctf.h:650
MDL_CTF_Q0
Inelastic absorption.
Definition: metadata_label.h:162
CTFDescription1D::gaussian_K2
double gaussian_K2
Gain for the second Gaussian term.
Definition: ctf.h:289
CTF
#define CTF
Definition: fourier_filter.h:81
CTFDescription::getValueNoiseAt
double getValueNoiseAt(bool show=false) const
Compute noise at (X,Y). Continuous frequencies, notice it is squared.
Definition: ctf.h:1121
MDL_CTF_LONGITUDINAL_DISPLACEMENT
Longitudinal displacement.
Definition: metadata_label.h:160
CTFDescription1D::getValuePureNoPrecomputedAt
double getValuePureNoPrecomputedAt(double X, bool show=false) const
Compute CTF pure at (U,V). Continuous frequencies.
Definition: ctf.h:613
HALFPI
#define HALFPI
Definition: xmipp_macros.h:104
XMIPP_EQUAL_ACCURACY
#define XMIPP_EQUAL_ACCURACY
Definition: xmipp_macros.h:119
CTFDescription::CTFDescription
CTFDescription()
Definition: ctf.h:899
errorBetween2CTFs
double errorBetween2CTFs(MetaData &MD1, MetaData &MD2, size_t dim, double minFreq=0.05, double maxFreq=0.25)
Definition: ctf.cpp:107
XSIZE
#define XSIZE(v)
Definition: multidim_array_base.h:91
CTF_ALL_LABELS_SIZE
const int CTF_ALL_LABELS_SIZE
Definition: ctf.h:47
CTFDescription1D::ispr
double ispr
Objective lens stability (deltaI/I) (ppm). Typical value 1.
Definition: ctf.h:253
CTFDescription::sqV
double sqV
Sqrt width V.
Definition: ctf.h:885
CTFDescription1D::Tm
double Tm
Sampling rate (A/pixel)
Definition: ctf.h:240
CTFDescription::getValuePureNoPrecomputedAtxy
double getValuePureNoPrecomputedAtxy(double X, double Y, bool show=false) const
Compute CTF pure at (U,V). Continuous frequencies.
Definition: ctf.h:1069
ZSIZE
#define ZSIZE(v)
Definition: multidim_array_base.h:99
CTFDescription1D::bgR3
double bgR3
Definition: ctf.h:297
CTFDescription::precomputeValues
void precomputeValues(double X, double Y)
Precompute values for a given frequency.
Definition: ctf.h:1002
CTFDescription::rad_gaussian2
double rad_gaussian2
Definition: ctf.h:808
CTFDescription1D::getValueDampingAt
double getValueDampingAt(bool show=false) const
Compute CTF damping at (U,V). Continuous frequencies.
Definition: ctf.h:424
CTFDescription1D::Cs
double Cs
Spherical aberration (in milimeters). Typical value 5.6.
Definition: ctf.h:246
CTFDescription1D::DeltaR
double DeltaR
Transversal mechanical displacement (ansgtrom). Typical value 3.
Definition: ctf.h:259
CTFDescription::sigmaV
double sigmaV
Gaussian width V.
Definition: ctf.h:873
CTFDescription1D::getValuePureWithoutDampingAt
double getValuePureWithoutDampingAt(bool show=false) const
Compute pure CTF without damping at (U,V). Continuous frequencies.
Definition: ctf.h:542
FFT_IDX2DIGFREQ
__device__ float FFT_IDX2DIGFREQ(int idx, int size)
Definition: cuda_gpu_reconstruct_fourier.cpp:382
CTFDescription1D::freq_max
double freq_max
Definition: ctf.h:303
errorMaxFreqCTFs
double errorMaxFreqCTFs(MetaData &MD1, double phaseRad=HALFPI)
Definition: ctf.cpp:187
CTFDescription1D::getValuePureNoKAt
double getValuePureNoKAt() const
Compute CTF pure at (U,V). Continuous frequencies.
Definition: ctf.h:499
j
#define j
Definition: numerical_recipes.cpp:2493
CTFDescription1D::Defocus
double Defocus
Defocus (in Angstroms). Negative values are underfocused.
Definition: ctf.h:244
YY
#define YY(v)
Definition: matrix1d.h:93
CTFDescription1D::changeSamplingRate
void changeSamplingRate(double newTm)
Definition: ctf.h:367
CTFDescription::generateCTF
void generateCTF(const MultidimArray< T1 > &sample_image, MultidimArray< T2 > &CTF, double Ts=-1)
Definition: ctf.h:1194
CTFDescription1D::K
double K
Global gain. By default, 1.
Definition: ctf.h:238
CTFDescription::gaussian_angle
double gaussian_angle
Gaussian angle.
Definition: ctf.h:879
CTFDescription1D::initCTF
double initCTF(int Ydim, int Xdim, MultidimArray< T > &CTF, double Ts=-1) const
Function to initialize CTF to avoid duplicated code.
Definition: ctf.h:676
CTFDescription1D::K5
double K5
Definition: ctf.h:222
CTFDescription::generateCTF
void generateCTF(int Ydim, int Xdim, MultidimArray< T > &CTF, double Ts=-1)
Generate CTF image.
Definition: ctf.h:1219
CTFDescription1D::generateCTF
void generateCTF(int Ydim, int Xdim, MultidimArray< T > &CTF, double Ts=-1)
Generate CTF image.
Definition: ctf.h:690
CTFDescription1D::D
double D
Definition: ctf.h:229
CTFDescription1D::isLocalCTF
bool isLocalCTF
Local CTF determination.
Definition: ctf.h:271
alglib::round
int round(double x)
Definition: ap.cpp:7245
PrecomputedForCTF::u_sqrt
double u_sqrt
Definition: ctf.h:85
CTFDescription1D::VPP_radius
double VPP_radius
Definition: ctf.h:306
CTFDescription1D::precomputedImage
std::vector< PrecomputedForCTF > precomputedImage
Definition: ctf.h:233
Matrix1D
Definition: ctf.h:38
PrecomputedForCTF::ang
double ang
Definition: ctf.h:90
RAD2DEG
#define RAD2DEG(r)
Definition: xmipp_macros.h:320
CTFDescription1D::Ca
double Ca
Chromatic aberration (in milimeters). Typical value 2.
Definition: ctf.h:248
CTFDescription1D::alpha
double alpha
Convergence cone semiangle (in mrad). Typical value 0.5.
Definition: ctf.h:255
CTFDescription1D::getSineAndCosineParts
void getSineAndCosineParts(double &sine_part, double &cosine_part, double E, double u2, double deltaf, bool show) const
Definition: ctf.h:572
CTFDescription1D::y0
double y0
In the case of local CTF determination x0,xF,y0,yF determines the region where the CTF is determined...
Definition: ctf.h:267
CTFDescription1D::Q0
double Q0
Factor for the importance of the Amplitude contrast.
Definition: ctf.h:261
CTFDescription1D::x0
double x0
In the case of local CTF determination x0,xF,y0,yF determines the region where the CTF is determined...
Definition: ctf.h:263
ERR_MULTIDIM_DIM
Incorrect MultidimArray dimensions.
Definition: xmipp_error.h:173
K
constexpr int K
Definition: integration.cpp:123
MDL_CTF_ENERGY_LOSS
Energy loss.
Definition: metadata_label.h:155
CTFDescription1D::kV
double kV
Accelerating Voltage (in KiloVolts)
Definition: ctf.h:242
xmipp_fft.h
CTFDescription::cV2
double cV2
Second Gaussian center for V.
Definition: ctf.h:895
CTFDescription1D::enable_CTFnoise
bool enable_CTFnoise
Enable CTFnoise part.
Definition: ctf.h:273
CTFDescription1D::xF
double xF
In the case of local CTF determination x0,xF,y0,yF determines the region where the CTF is determined...
Definition: ctf.h:265
cif::pdb::read
file read(std::istream &is)
Definition: pdb2cif.cpp:6200
CTFDescription1D::getValuePureAt
double getValuePureAt(bool show=false) const
Compute CTF pure at (U,V). Continuous frequencies.
Definition: ctf.h:452
MDL_CTF_DEFOCUSV
Defocus V (Angstroms)
Definition: metadata_label.h:136
CTFDescription1D::sigma1
double sigma1
Gaussian width.
Definition: ctf.h:281
CTFDescription1D::sqrt_K
double sqrt_K
Gain for the square root term.
Definition: ctf.h:285
bessj0
double bessj0(double x)
Definition: numerical_recipes.cpp:250
MDL_CTF_CA
Chromatic aberration.
Definition: metadata_label.h:153
MDL_CTF_TRANSVERSAL_DISPLACEMENT
Transversal displacemente.
Definition: metadata_label.h:161
CTFDescription1D::precomputed
PrecomputedForCTF precomputed
Definition: ctf.h:231
xmipp_macros.h
fmax
double fmax
Definition: numerical_recipes.cpp:5166
CTFDescription1D::yF
double yF
In the case of local CTF determination x0,xF,y0,yF determines the region where the CTF is determined...
Definition: ctf.h:269
n
int * n
Definition: numerical_recipes.cpp:2229
CTFDescription1D::Kcos
double Kcos
Definition: ctf.h:226
MDL_CTF_LENS_STABILITY
Lens stability.
Definition: metadata_label.h:158
CTFDescription::getPhaseAt
double getPhaseAt() const
Get Phase of the CTF.
Definition: ctf.h:1063
MDLabel
MDLabel
Definition: metadata_label.h:44
CTFDescription::sigmaU2
double sigmaU2
Second Gaussian width U.
Definition: ctf.h:889
CTFDescription::precomputeValues
void precomputeValues(int i, int j)
Precompute values for a given frequency.
Definition: ctf.h:1035
SINC
#define SINC(x)
Definition: xmipp_macros.h:362
CTFDescription::sigmaV2
double sigmaV2
Second Gaussian width V.
Definition: ctf.h:891
CTFDescription::cV
double cV
Gaussian center for V.
Definition: ctf.h:877
PrecomputedForCTF::u3
double u3
Definition: ctf.h:88
CTFDescription1D::K3
double K3
Definition: ctf.h:220
CTFDescription::generateCTFWithoutDamping
void generateCTFWithoutDamping(int Ydim, int Xdim, MultidimArray< T > &CTF, double Ts=-1)
Definition: ctf.h:1245
Generated by   doxygen 1.8.13