Xmipp  v3.23.11-Nereus
ctf_estimate_from_psd.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  *
3  * Authors: Javier Angel Velazquez Muriel (javi@cnb.csic.es)
4  * Carlos Oscar Sanchez Sorzano (coss@cnb.csic.es)
5  *
6  * Unidad de Bioinformatica of Centro Nacional de Biotecnologia , CSIC
7  *
8  * This program is free software; you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation; either version 2 of the License, or
11  * (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with this program; if not, write to the Free Software
20  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
21  * 02111-1307 USA
22  *
23  * All comments concerning this program package may be sent to the
24  * e-mail address 'xmipp@cnb.csic.es'
25  ***************************************************************************/
26 
27 //#include "ctf_estimate_psd_with_arma.h"
28 #include <fstream>
29 #include "ctf_estimate_from_psd.h"
30 #include "ctf_estimate_from_psd_fast.h"
31 #include "core/matrix2d.h"
32 #include "core/transformations.h"
33 #include "core/metadata_vec.h"
34 #include "data/numerical_tools.h"
35 #include "fringe_processing.h"
36 
37 /* prototypes */
38 double CTF_fitness(double *, void *);
39 
40 /* Number of CTF parameters */
41 constexpr int ALL_CTF_PARAMETERS = 38;
42 constexpr int CTF_PARAMETERS = 30;
43 constexpr int PARAMETRIC_CTF_PARAMETERS = 16;
44 constexpr int BACKGROUND_CTF_PARAMETERS = 14;
45 constexpr int SQRT_CTF_PARAMETERS = 8;
46 constexpr int ENVELOPE_PARAMETERS = 11;
47 constexpr int DEFOCUS_PARAMETERS = 5;
48 constexpr int FIRST_SQRT_PARAMETER = 16;
49 constexpr int FIRST_ENVELOPE_PARAMETER = 4;
50 constexpr int FIRST_DEFOCUS_PARAMETER = 0;
51 
52 //#define DEBUG_WITH_TEXTFILES
53 #ifdef DEBUG_WITH_TEXTFILES
54 std::ofstream fhDebug;
55 #define DEBUG_OPEN_TEXTFILE(fnRoot) fhDebug.open((fnRoot+"_debug.txt").c_str());
56 #define DEBUG_CLOSE_TEXTFILE fhDebug.close();
57 #define DEBUG_TEXTFILE(str) fhDebug << time (NULL) << " " << str << std::endl;
58 #define DEBUG_MODEL_TEXTFILE fhDebug << global_ctfmodel << std::endl;
59 #else
60 #define DEBUG_OPEN_TEXTFILE(fnRoot);
61 #define DEBUG_CLOSE_TEXTFILE ;
62 #define DEBUG_TEXTFILE(str);
63 #define DEBUG_MODEL_TEXTFILE
64 #endif
65 
66 /* Global variables -------------------------------------------------------- */
67 namespace AdjustCTF
68 {
69 // Some aliases
70 ProgCTFEstimateFromPSD *global_prm;
71 // CTF model and noise model
72 //CTFDescription ctfmodel_defoci;
73 
74 }
75 
76 using namespace AdjustCTF;
77 
78 #define ASSIGN_CTF_PARAM(index, paramName) if (ia <= index && l > 0) { ctfmodel.paramName = p[index]; --l; }
79 
80 /* Assign ctfmodel from a vector and viceversa ----------------------------- */
81 void ProgCTFEstimateFromPSD::assignCTFfromParameters(double *p, CTFDescription &ctfmodel, int ia,
82  int l, int modelSimplification)
83 {
84 
85  ctfmodel.Tm = Tm;
86 
87  /*
88  BE CAREFULL in add new parameters between the existing ones !!
89  below in the core of the program we set step(i) = 0 or 1
90  if you change here the value of certain param, change it everywhere!
91  */
92  ASSIGN_CTF_PARAM(0, DeltafU);
93  ASSIGN_CTF_PARAM(1, DeltafV);
94  ASSIGN_CTF_PARAM(2, azimuthal_angle);
95  ASSIGN_CTF_PARAM(3, kV);
96  ASSIGN_CTF_PARAM(4, K);
97  ASSIGN_CTF_PARAM(5, Cs);
98  ASSIGN_CTF_PARAM(6, Ca);
99  ASSIGN_CTF_PARAM(7, espr);
100  ASSIGN_CTF_PARAM(8, ispr);
101  ASSIGN_CTF_PARAM(9, alpha);
102  ASSIGN_CTF_PARAM(10, DeltaF);
103  ASSIGN_CTF_PARAM(11, DeltaR);
104  ASSIGN_CTF_PARAM(12, envR0);
105  ASSIGN_CTF_PARAM(13, envR1);
106  ASSIGN_CTF_PARAM(14, envR2);
107  ASSIGN_CTF_PARAM(15, Q0);
108  ASSIGN_CTF_PARAM(16, base_line);
109  ASSIGN_CTF_PARAM(17, sqrt_K);
110  ASSIGN_CTF_PARAM(18, sqU);
111  ASSIGN_CTF_PARAM(19, sqV);
112  ASSIGN_CTF_PARAM(20, sqrt_angle);
113  ASSIGN_CTF_PARAM(21, bgR1);
114  ASSIGN_CTF_PARAM(22, bgR2);
115  ASSIGN_CTF_PARAM(23, bgR3);
116  ASSIGN_CTF_PARAM(24, gaussian_K);
117  ASSIGN_CTF_PARAM(25, sigmaU);
118 
119  if (ia <= 26 && l > 0)
120  {
121  if (modelSimplification < 3)
122  {
123  ctfmodel.sigmaV = p[26];
124  l--;
125  } // 7 *
126  else
127  {
128  ctfmodel.sigmaV = p[25];
129  l--;
130  }
131  }
132  if (ia <= 27 && l > 0)
133  {
134  if (modelSimplification < 3)
135  {
136  ctfmodel.gaussian_angle = p[27];
137  l--;
138  } // 8 *
139  else
140  {
141  ctfmodel.gaussian_angle = 0;
142  l--;
143  }
144  }
145 
146  ASSIGN_CTF_PARAM(28, cU);
147 
148  if (ia <= 29 && l > 0)
149  {
150  if (modelSimplification < 3)
151  {
152  ctfmodel.cV = p[29];
153  l--;
154  } // 10 *
155  else
156  {
157  ctfmodel.cV = p[28];
158  l--;
159  }
160  }
161 
162  ASSIGN_CTF_PARAM(30, gaussian_K2);
163  ASSIGN_CTF_PARAM(31, sigmaU2);
164  ASSIGN_CTF_PARAM(32, sigmaV2);
165  ASSIGN_CTF_PARAM(33, gaussian_angle2);
166  ASSIGN_CTF_PARAM(34, cU2);
167  ASSIGN_CTF_PARAM(35, cV2);
168  ASSIGN_CTF_PARAM(36, phase_shift);
169  ASSIGN_CTF_PARAM(37, VPP_radius);
170 
171 }//function assignCTFfromParameters
172 
173 
174 #define ASSIGN_PARAM_CTF(index, paramName) if (ia <= index && l > 0) { p[index] = ctfmodel.paramName; --l; }
175 
176 void ProgCTFEstimateFromPSD::assignParametersFromCTF(const CTFDescription &ctfmodel, double *p, int ia,
177  int l, int modelSimplification)
178 {
179  /*
180  BE CAREFUL in add new parameters between the existing ones !!
181  below in the core of the program we set step(i) to 0 or 1
182  if you change here the value of certain param, change it everywhere!
183  */
184  ASSIGN_PARAM_CTF(0, DeltafU);
185  ASSIGN_PARAM_CTF(1, DeltafV);
186  ASSIGN_PARAM_CTF(2, azimuthal_angle);
187  ASSIGN_PARAM_CTF(3, kV);
188  ASSIGN_PARAM_CTF(4, K);
189  ASSIGN_PARAM_CTF(5, Cs);
190  ASSIGN_PARAM_CTF(6, Ca);
191  ASSIGN_PARAM_CTF(7, espr);
192  ASSIGN_PARAM_CTF(8, ispr);
193  ASSIGN_PARAM_CTF(9, alpha);
194  ASSIGN_PARAM_CTF(10, DeltaF);
195  ASSIGN_PARAM_CTF(11, DeltaR);
196  ASSIGN_PARAM_CTF(12, envR0);
197  ASSIGN_PARAM_CTF(13, envR1);
198  ASSIGN_PARAM_CTF(14, envR2);
199  ASSIGN_PARAM_CTF(15, Q0);
200  ASSIGN_PARAM_CTF(16, base_line);
201  ASSIGN_PARAM_CTF(17, sqrt_K);
202  ASSIGN_PARAM_CTF(18, sqU);
203  ASSIGN_PARAM_CTF(19, sqV);
204  ASSIGN_PARAM_CTF(20, sqrt_angle);
205  ASSIGN_PARAM_CTF(21, bgR1);
206  ASSIGN_PARAM_CTF(22, bgR2);
207  ASSIGN_PARAM_CTF(23, bgR3);
208  ASSIGN_PARAM_CTF(24, gaussian_K);
209  ASSIGN_PARAM_CTF(25, sigmaU);
210 
211  if (ia <= 26 && l > 0)
212  {
213  if (modelSimplification < 3)
214  {
215  p[26] = ctfmodel.sigmaV;
216  l--;
217  }
218  else
219  {
220  p[26] = ctfmodel.sigmaU;
221  l--;
222  }
223  }
224  if (ia <= 27 && l > 0)
225  {
226  if (modelSimplification < 3)
227  {
228  p[27] = ctfmodel.gaussian_angle;
229  l--;
230  }
231  else
232  {
233  p[27] = 0;
234  l--;
235  }
236  }
237 
238  ASSIGN_PARAM_CTF(28, cU);
239 
240  if (ia <= 29 && l > 0)
241  {
242  if (modelSimplification < 3)
243  {
244  p[29] = ctfmodel.cV;
245  l--;
246  }
247  else
248  {
249  p[29] = ctfmodel.cU;
250  l--;
251  }
252  }
253 
254  ASSIGN_PARAM_CTF(30, gaussian_K2);
255  ASSIGN_PARAM_CTF(31, sigmaU2);
256  ASSIGN_PARAM_CTF(32, sigmaV2);
257  ASSIGN_PARAM_CTF(33, gaussian_angle2);
258  ASSIGN_PARAM_CTF(34, cU2);
259  ASSIGN_PARAM_CTF(35, cV2);
260  ASSIGN_PARAM_CTF(36, phase_shift);
261  ASSIGN_PARAM_CTF(37, VPP_radius);
262 }
263 
264 #define COPY_ctfmodel_TO_CURRENT_GUESS \
265  global_prm->assignParametersFromCTF(global_prm->current_ctfmodel, \
266  MATRIX1D_ARRAY(*global_prm->adjust_params),0,ALL_CTF_PARAMETERS, \
267  global_prm->modelSimplification);
268 
269 /* Read parameters --------------------------------------------------------- */
270 void ProgCTFEstimateFromPSD::readBasicParams(XmippProgram *program)
271 {
272  ProgCTFBasicParams::readBasicParams(program);
273 
274  initial_ctfmodel.enable_CTF = initial_ctfmodel.enable_CTFnoise = true;
275  initial_ctfmodel.readParams(program);
276  if (initial_ctfmodel.DeltafU>100e3 || initial_ctfmodel.DeltafV>100e3)
277  REPORT_ERROR(ERR_ARG_INCORRECT,"Defocus cannot be larger than 10 microns (100,000 Angstroms)");
278  Tm = initial_ctfmodel.Tm;
279 
280 }
281 
282 void ProgCTFEstimateFromPSD::readParams()
283 {
284  fn_psd = getParam("--psd");
285  readBasicParams(this);
286 }
287 
288 /* Usage ------------------------------------------------------------------- */
289 void ProgCTFEstimateFromPSD::defineBasicParams(XmippProgram * program)
290 {
291 // ProgCTFBasicParams::defineBasicParams(program);
292  CTFDescription::defineParams(program);
293 
294 }
295 
296 void ProgCTFEstimateFromPSD::defineParams()
297 {
298  defineBasicParams(this);
299  ProgCTFBasicParams::defineParams();
300 }
301 
302 /* Produce side information ------------------------------------------------ */
303 void ProgCTFEstimateFromPSD::produceSideInfo()
304 {
305  adjust.resize(ALL_CTF_PARAMETERS);
306  adjust.initZeros();
307  current_ctfmodel.clear();
308  ctfmodel_defoci.clear();
309  assignParametersFromCTF(initial_ctfmodel, MATRIX1D_ARRAY(adjust), 0, ALL_CTF_PARAMETERS, true);
310 
311  // Read the CTF file, supposed to be the uncentered squared amplitudes
312  if (fn_psd != "")
313  ctftomodel.read(fn_psd);
314  f = &(ctftomodel());
315 
316  ProgCTFBasicParams::produceSideInfo();
317  current_ctfmodel.precomputeValues(x_contfreq, y_contfreq);
318 }
319 
320 /* Generate model so far ---------------------------------------------------- */
321 /* The model is taken from global_adjust and global_ctfmodel is modified */
322 void ProgCTFEstimateFromPSD::generateModelSoFar(Image<double> &I, bool apply_log = false)
323 {
324  Matrix1D<int> idx(2); // Indexes for Fourier plane
325  Matrix1D<double> freq(2); // Frequencies for Fourier plane
326  assignCTFfromParameters(MATRIX1D_ARRAY(*adjust_params), current_ctfmodel,
327  0, ALL_CTF_PARAMETERS, modelSimplification);
328  current_ctfmodel.produceSideInfo();
329  I().initZeros(*f);
330  FOR_ALL_ELEMENTS_IN_ARRAY2D(I())
331  {
332  XX(idx) = j;
333  YY(idx) = i;
334  FFT_idx2digfreq(*f, idx, freq);
335  double w=freq.module();
336  if (w>max_freq_psd)
337  continue;
338  digfreq2contfreq(freq, freq, Tm);
339 
340  // Decide what to save
341  current_ctfmodel.precomputeValues(XX(freq), YY(freq));
342  if (action <= 1)
343  I()(i, j) = current_ctfmodel.getValueNoiseAt();
344  else if (action == 2)
345  {
346  double E = current_ctfmodel.getValueDampingAt();
347  I()(i, j) = current_ctfmodel.getValueNoiseAt() + E * E;
348  }
349  else if (action >= 3 && action <= 6)
350  {
351  double ctf = current_ctfmodel.getValuePureAt();
352  I()(i, j) = current_ctfmodel.getValueNoiseAt() + ctf * ctf;
353  }
354  else
355  {
356  double ctf = current_ctfmodel.getValuePureAt();
357  I()(i, j) = ctf;
358  }
359  if (apply_log)
360  I()(i, j) = 10 * log10(I()(i, j));
361  }
362 }
363 
364 /* Save intermediate results ----------------------------------------------- */
365 /* First call to generate model so far and then save the image, and a couple
366  of cuts along X and Y.
367 
368  This function returns the fitting error.*/
369 void ProgCTFEstimateFromPSD::saveIntermediateResults(const FileName &fn_root, bool generate_profiles = true)
370 {
371  std::ofstream plotX, plotY, plot_radial;
372  Image<double> save;
373  generateModelSoFar(save, false);
374 
375  Image<double> save_ctf;
376  generate_model_halfplane(ctfmodelSize,
377  ctfmodelSize, save_ctf());
378  if (fn_root.find("@")==std::string::npos)
379  save_ctf.write(fn_root + "_ctfmodel_halfplane.xmp");
380  else
381  save_ctf.write(fn_root + "_ctfmodel_halfplane.stk");
382  generate_model_quadrant(ctfmodelSize,
383  ctfmodelSize, save_ctf());
384  if (fn_root.find("@")==std::string::npos)
385  save_ctf.write(fn_root + "_ctfmodel_quadrant.xmp");
386  else
387  save_ctf.write(fn_root + "_ctfmodel_quadrant.stk");
388 
389  if (!generate_profiles)
390  return;
391  plotX.open((fn_root + "X.txt").c_str());
392  plotY.open((fn_root + "Y.txt").c_str());
393  plot_radial.open((fn_root + "_radial.txt").c_str());
394  if (!plotX)
395  REPORT_ERROR(
396  ERR_IO_NOWRITE,
397  "save_intermediate_results::Cannot open plot file for writing\n");
398  if (!plotY)
399  REPORT_ERROR(
400  ERR_IO_NOWRITE,
401  "save_intermediate_results::Cannot open plot file for writing\n");
402  if (!plot_radial)
403  REPORT_ERROR(
404  ERR_IO_NOWRITE,
405  "save_intermediate_results::Cannot open plot file for writing\n");
406  plotX << "# freq_dig freq_angstrom model psd enhanced logModel logPsd\n";
407  plotY << "# freq_dig freq_angstrom model psd enhanced logModel logPsd\n";
408  //plot_radial << "# freq_dig freq_angstrom model psd enhanced logModel logPsd\n";
409 
410  // Generate cut along X
411  for (int i = STARTINGY(save()); i <= FINISHINGY(save()) / 2; i++)
412  {
413  if (mask(i, 0) <= 0)
414  continue;
415  plotY << w_digfreq(i, 0) << " " << w_contfreq(i, 0) << " "
416  << save()(i, 0) << " " << (*f)(i, 0) << " "
417  << enhanced_ctftomodel()(i, 0) << " "
418  << log10(save()(i, 0)) << " " << log10((*f)(i, 0))
419  << std::endl;
420  }
421 
422  // Generate cut along Y
423  for (int j = STARTINGX(save()); j <= FINISHINGX(save()) / 2; j++)
424  {
425  if (mask(0, j) <= 0)
426  continue;
427  plotX << w_digfreq(0, j) << " " << w_contfreq(0, j) << " "
428  << save()(0, j) << " " << (*f)(0, j) << " "
429  << enhanced_ctftomodel()(0, j)
430  << log10(save()(0, j)) << " " << log10((*f)(0, j))
431  << std::endl;
432  }
433 
434  // Generate radial average
435  MultidimArray<double> radial_CTFmodel_avg, radial_CTFampl_avg,
436  radial_enhanced_avg;
437  MultidimArray<int> radial_N;
438  radial_CTFmodel_avg.initZeros(YSIZE(save()) / 2);
439  radial_CTFampl_avg.initZeros(YSIZE(save()) / 2);
440  radial_enhanced_avg.initZeros(YSIZE(save()) / 2);
441  radial_N.initZeros(YSIZE(save()) / 2);
442  FOR_ALL_ELEMENTS_IN_ARRAY2D(w_digfreq)
443  {
444  if (mask(i, j) <= 0)
445  continue;
446  double model2 = save()(i, j);
447 
448  int r = w_digfreq_r(i, j);
449  radial_CTFmodel_avg(r) += model2;
450  radial_CTFampl_avg(r) += (*f)(i, j);
451  radial_enhanced_avg(r) += enhanced_ctftomodel()(i, j);
452  radial_N(r)++;
453  }
454 
455  FOR_ALL_ELEMENTS_IN_ARRAY1D(radial_CTFmodel_avg)
456  {
457  if (radial_N(i) == 0)
458  continue;
459  plot_radial << w_digfreq(i, 0) << " " << w_contfreq(i, 0)
460  << " " << radial_CTFmodel_avg(i) / radial_N(i) << " " << radial_CTFampl_avg(i) / radial_N(i) << " "
461  << radial_enhanced_avg(i) / radial_N(i)
462  << " " << log10(radial_CTFmodel_avg(i) / radial_N(i)) << " " << log10(radial_CTFampl_avg(i) / radial_N(i)) << " "
463  << std::endl;
464  }
465 
466  plotX.close();
467  plotY.close();
468  plot_radial.close();
469 }
470 
471 /* Generate model at a given size ------------------------------------------ */
472 void ProgCTFEstimateFromPSD::generate_model_quadrant(int Ydim, int Xdim,
473  MultidimArray<double> &model)
474 {
475  Matrix1D<int> idx(2); // Indexes for Fourier plane
476  Matrix1D<double> freq(2); // Frequencies for Fourier plane
477 
478  // Compute the scaled PSD
479  MultidimArray<double> enhancedPSD;
480  enhancedPSD = enhanced_ctftomodel_fullsize();
481  CenterFFT(enhancedPSD, false);
482  selfScaleToSize(xmipp_transformation::BSPLINE3, enhancedPSD, Ydim, Xdim);
483  CenterFFT(enhancedPSD, true);
484 
485  // Generate the CTF model
486  assignCTFfromParameters(MATRIX1D_ARRAY(*adjust_params), current_ctfmodel,
487  0, ALL_CTF_PARAMETERS, modelSimplification);
488  current_ctfmodel.produceSideInfo();
489 
490  // Write the two model quadrants
491  MultidimArray<int> mask_norm;
492  STARTINGX(enhancedPSD)=STARTINGY(enhancedPSD)=0;
493  mask_norm.initZeros(enhancedPSD);
494  model.resizeNoCopy(enhancedPSD);
495  FOR_ALL_ELEMENTS_IN_ARRAY2D(model)
496  {
497  if ((j >= Xdim / 2 && i >= Ydim / 2)
498  || (j < Xdim / 2 && i < Ydim / 2))
499  {
500  XX(idx) = j;
501  YY(idx) = i;
502  FFT_idx2digfreq(model, idx, freq);
503  if (fabs(XX(freq))>0.03 && fabs(YY(freq))>0.03)
504  mask_norm(i,j)=(int)mask(i,j);
505  digfreq2contfreq(freq, freq, Tm);
506 
507  current_ctfmodel.precomputeValues(XX(freq), YY(freq));
508  model(i, j) = current_ctfmodel.getValuePureAt();
509  model(i, j) *= model(i, j);
510  }
511  }
512  // Normalize the left part so that it has similar values to
513  // the enhanced PSD
514  model.rangeAdjust(enhancedPSD, &mask_norm);
515  // Copy the part of the enhancedPSD
516  FOR_ALL_ELEMENTS_IN_ARRAY2D(model)
517  {
518  if (!((j >= Xdim / 2 && i >= Ydim / 2) || (j < Xdim / 2 && i < Ydim / 2)))
519  model(i, j) = enhancedPSD(i, j);
520  else
521  {
522  XX(idx) = j;
523  YY(idx) = i;
524  FFT_idx2digfreq(model, idx, freq);
525  double w=freq.module();
526  if (w>max_freq_psd)
527  model(i,j)=0;
528 
529  }
530  }
531 
532  // Produce a centered image
533  CenterFFT(model, true);
534 }
535 
536 void ProgCTFEstimateFromPSD::generate_model_halfplane(int Ydim, int Xdim,
537  MultidimArray<double> &model)
538 {
539  Matrix1D<int> idx(2); // Indexes for Fourier plane
540  Matrix1D<double> freq(2); // Frequencies for Fourier plane
541 
542  // Compute the scaled PSD
543  MultidimArray<double> enhancedPSD;
544  enhancedPSD = enhanced_ctftomodel_fullsize();
545  CenterFFT(enhancedPSD, false);
546  selfScaleToSize(xmipp_transformation::BSPLINE3, enhancedPSD, Ydim, Xdim);
547  CenterFFT(enhancedPSD, true);
548 
549  // The left part is the CTF model
550  assignCTFfromParameters(MATRIX1D_ARRAY(*adjust_params), current_ctfmodel,
551  0, CTF_PARAMETERS, modelSimplification);
552  current_ctfmodel.produceSideInfo();
553 
554  MultidimArray<int> mask_norm;
555  STARTINGX(enhancedPSD)=STARTINGY(enhancedPSD)=0;
556  mask_norm.initZeros(enhancedPSD);
557  model.resizeNoCopy(enhancedPSD);
558  FOR_ALL_ELEMENTS_IN_ARRAY2D(model)
559  {
560  if (j <= Xdim / 2)
561  {
562  XX(idx) = j;
563  YY(idx) = i;
564  FFT_idx2digfreq(model, idx, freq);
565  double w=freq.module();
566  if (w>max_freq_psd)
567  continue;
568  if (fabs(XX(freq))>0.03 && fabs(YY(freq))>0.03)
569  mask_norm(i,j)=(int)mask(i,j);
570  digfreq2contfreq(freq, freq, Tm);
571 
572  current_ctfmodel.precomputeValues(XX(freq), YY(freq));
573  model(i, j) = current_ctfmodel.getValuePureAt();
574  model(i, j) *= model(i, j);
575  }
576  }
577  // Normalize the left part so that it has similar values to
578  // the enhanced PSD
579  model.rangeAdjust(enhancedPSD, &mask_norm);
580  // Copy the part of the enhancedPSD
581  FOR_ALL_ELEMENTS_IN_ARRAY2D(model)
582  {
583  if (!(j <= Xdim / 2))
584  model(i, j) = enhancedPSD(i, j);
585  else
586  {
587  XX(idx) = j;
588  YY(idx) = i;
589  FFT_idx2digfreq(model, idx, freq);
590  double w=freq.module();
591  if (w>max_freq_psd)
592  model(i,j)=0;
593  }
594  }
595  CenterFFT(model, true);
596 }
597 
598 /* CTF fitness ------------------------------------------------------------- */
599 /* This function measures the distance between the estimated CTF and the
600  measured CTF */
601 double ProgCTFEstimateFromPSD::CTF_fitness_object(double *p)
602 {
603  double retval;
604  // Generate CTF model
605  switch (action)
606  {
607  // Remind that p is a vector whose first element is at index 1
608  case 0:
609  assignCTFfromParameters(p - FIRST_SQRT_PARAMETER + 1,
610  current_ctfmodel, FIRST_SQRT_PARAMETER, SQRT_CTF_PARAMETERS,
611  modelSimplification);
612  if (show_inf >= 2)
613  {
614  std::cout << "Input vector:";
615  for (int i = 1; i <= SQRT_CTF_PARAMETERS; i++)
616  std::cout << p[i] << " ";
617  std::cout << std::endl;
618  }
619  break;
620  case 1:
621  assignCTFfromParameters(p - FIRST_SQRT_PARAMETER + 1,
622  current_ctfmodel, FIRST_SQRT_PARAMETER,
623  BACKGROUND_CTF_PARAMETERS, modelSimplification);
624  if (show_inf >= 2)
625  {
626  std::cout << "Input vector:";
627  for (int i = 1; i <= BACKGROUND_CTF_PARAMETERS; i++)
628  std::cout << p[i] << " ";
629  std::cout << std::endl;
630  }
631  break;
632  case 2:
633  assignCTFfromParameters(p - FIRST_ENVELOPE_PARAMETER + 1,
634  current_ctfmodel, FIRST_ENVELOPE_PARAMETER, ENVELOPE_PARAMETERS,
635  modelSimplification);
636  if (show_inf >= 2)
637  {
638  std::cout << "Input vector:";
639  for (int i = 1; i <= ENVELOPE_PARAMETERS; i++)
640  std::cout << p[i] << " ";
641  std::cout << std::endl;
642  }
643  break;
644  case 3:
645  assignCTFfromParameters(p - FIRST_DEFOCUS_PARAMETER + 1,
646  current_ctfmodel, FIRST_DEFOCUS_PARAMETER, DEFOCUS_PARAMETERS,
647  modelSimplification);
648  psd_theo_radial_derivative.initZeros();
649  if (show_inf >= 2)
650  {
651  std::cout << "Input vector:";
652  for (int i = 1; i <= DEFOCUS_PARAMETERS; i++)
653  std::cout << p[i] << " ";
654  std::cout << std::endl;
655  }
656  break;
657  case 4:
658  assignCTFfromParameters(p - 0 + 1, current_ctfmodel, 0,
659  CTF_PARAMETERS, modelSimplification);
660  psd_theo_radial.initZeros();
661  if (show_inf >= 2)
662  {
663  std::cout << "Input vector:";
664  for (int i = 1; i <= CTF_PARAMETERS; i++)
665  std::cout << p[i] << " ";
666  std::cout << std::endl;
667  }
668  break;
669  case 5:
670  case 6:
671  case 7:
672  assignCTFfromParameters(p - 0 + 1, current_ctfmodel, 0,
673  ALL_CTF_PARAMETERS, modelSimplification);
674  psd_theo_radial.initZeros();
675  if (show_inf >= 2)
676  {
677  std::cout << "Input vector:";
678  for (int i = 1; i <= ALL_CTF_PARAMETERS; i++)
679  std::cout << p[i] << " ";
680  std::cout << std::endl;
681  }
682 
683  break;
684  }
685 
686  current_ctfmodel.produceSideInfo();
687 
688  if (show_inf >= 2)
689  std::cout << "Model:\n" << current_ctfmodel << std::endl;
690  if (!current_ctfmodel.hasPhysicalMeaning())
691  {
692  if (show_inf >= 2)
693  std::cout << "Does not have physical meaning\n";
694  return heavy_penalization;
695  }
696 
697  if (action > 3 && !noDefocusEstimate
698  && (fabs(
699  (current_ctfmodel.DeltafU - ctfmodel_defoci.DeltafU)
700  / ctfmodel_defoci.DeltafU) > 0.2
701  || fabs(
702  (current_ctfmodel.DeltafV
703  - ctfmodel_defoci.DeltafV)
704  / ctfmodel_defoci.DeltafU) > 0.2))
705  {
706  if (show_inf >= 2)
707  std::cout << "Too large defocus\n";
708  return heavy_penalization;
709  }
710  /*if (initial_ctfmodel.DeltafU != 0 && action >= 3)
711  {
712  // If there is an initial model, the true solution
713  // cannot be too far
714 
715  if (fabs(initial_ctfmodel.DeltafU - current_ctfmodel.DeltafU) > defocus_range ||
716  fabs(initial_ctfmodel.DeltafV - current_ctfmodel.DeltafV) > defocus_range)
717  {std::cout << "Entra2" << std::endl;
718  if (show_inf >= 2)
719  {
720  std::cout << "Too far from hint: Initial (" << initial_ctfmodel.DeltafU << "," << initial_ctfmodel.DeltafV << ")"
721  << " current guess (" << current_ctfmodel.DeltafU << "," << current_ctfmodel.DeltafV << ") max allowed difference: "
722  << defocus_range << std::endl;
723  }
724  return heavy_penalization;
725  }
726  }*/
727 
728  /*if ((initial_ctfmodel.phase_shift != 0.0 || initial_ctfmodel.phase_shift != 1.57079) && action >= 5)
729  {
730  if (fabs(initial_ctfmodel.phase_shift - current_ctfmodel.phase_shift) > 0.05)
731  {
732  return heavy_penalization;
733  }
734  }*/
735 
736  // Now the 2D error
737  double distsum = 0;
738  int N = 0, Ncorr = 0;
739  double enhanced_avg = 0;
740  double model_avg = 0;
741  double enhanced_model = 0;
742  double enhanced2 = 0;
743  double model2 = 0;
744  double lowerLimit = 1.1 * min_freq_psd;
745  double upperLimit = 0.9 * max_freq_psd;
746  const MultidimArray<double>& local_enhanced_ctf = enhanced_ctftomodel();
747  int XdimW=XSIZE(w_digfreq);
748  int YdimW=YSIZE(w_digfreq);
749  corr13=0;
750 
751  for (int i = 0; i < YdimW; i += evaluation_reduction)
752  for (int j = 0; j < XdimW; j += evaluation_reduction)
753  {
754  if (DIRECT_A2D_ELEM(mask, i, j) <= 0)
755  continue;
756 
757  // Compute each component
758  current_ctfmodel.precomputeValues(i, j);
759  double bg = current_ctfmodel.getValueNoiseAt();
760  double envelope=0, ctf_without_damping, ctf_with_damping=0;
761  double ctf2_th=0;
762  double ctf2 = DIRECT_A2D_ELEM(*f, i, j);
763  double dist = 0;
764  double ctf_with_damping2;
765  switch (action)
766  {
767  case 0:
768  case 1:
769  ctf2_th = bg;
770  dist = fabs(ctf2 - bg);
771  if (penalize && bg > ctf2
772  && DIRECT_A2D_ELEM(w_digfreq, i, j)
773  > max_gauss_freq)
774  dist *= current_penalty;
775  break;
776  case 2:
777  envelope = current_ctfmodel.getValueDampingAt();
778  ctf2_th = bg + envelope * envelope;
779  dist = fabs(ctf2 - ctf2_th);
780  if (penalize && ctf2_th < ctf2
781  && DIRECT_A2D_ELEM(w_digfreq, i, j)
782  > max_gauss_freq)
783  dist *= current_penalty;
784  break;
785  case 3:
786  case 4:
787  case 5:
788  case 6:
789  case 7:
790  envelope = current_ctfmodel.getValueDampingAt();
791  ctf_without_damping =
792  current_ctfmodel.getValuePureWithoutDampingAt();
793  ctf_with_damping = envelope * ctf_without_damping;
794  ctf2_th = bg + ctf_with_damping * ctf_with_damping;
795 
796  if (DIRECT_A2D_ELEM(w_digfreq,i, j) < upperLimit
797  && DIRECT_A2D_ELEM(w_digfreq,i, j) > lowerLimit)
798  {
799  if (action == 3 ||
800  (action == 4 && DIRECT_A2D_ELEM(mask_between_zeroes,i,j) == 1) ||
801  (action == 7 && DIRECT_A2D_ELEM(mask_between_zeroes,i,j) == 1))
802  {
803  double enhanced_ctf =
804  DIRECT_A2D_ELEM(local_enhanced_ctf, i, j);
805  ctf_with_damping2 = ctf_with_damping * ctf_with_damping;
806  enhanced_model += enhanced_ctf * ctf_with_damping2;
807  enhanced2 += enhanced_ctf * enhanced_ctf;
808  model2 += ctf_with_damping2 * ctf_with_damping2;
809  enhanced_avg += enhanced_ctf;
810  model_avg += ctf_with_damping2;
811  Ncorr++;
812 
813  if (action==3)
814  {
815  int r = A2D_ELEM(w_digfreq_r,i, j);
816  A1D_ELEM(psd_theo_radial,r) += ctf2_th;
817  }
818  }
819  }
820  if (envelope > 1e-2)
821  dist = fabs(ctf2 - ctf2_th) / (envelope * envelope);
822  else
823  dist = fabs(ctf2 - ctf2_th);
824  // This expression comes from mapping any value so that
825  // bg becomes 0, and bg+envelope^2 becomes 1
826  // This is the transformation
827  // (x-bg) x-bg
828  // -------------=-------
829  // (bg+env^2-bg) env^2
830  // If we subtract two of this scaled values
831  // x-bg y-bg x-y
832  // ------- - ------- = -------
833  // env^2 env^2 env^2
834  break;
835  }
836 
837  distsum += dist * DIRECT_A2D_ELEM(mask,i,j);
838  N++;
839  }
840 
841  if (N > 0)
842  retval = distsum / N;
843  else
844  retval = heavy_penalization;
845 
846  if (show_inf >=2)
847  std::cout << "Fitness1=" << retval << std::endl;
848 
849  if ( (((action >= 3) && (action <= 4)) || (action == 7))
850  && (Ncorr > 0) && (enhanced_weight != 0) )
851  {
852  model_avg /= Ncorr;
853  enhanced_avg /= Ncorr;
854  double correlation_coeff = enhanced_model / Ncorr
855  - model_avg * enhanced_avg;
856  double sigma1 = sqrt(
857  fabs(enhanced2 / Ncorr - enhanced_avg * enhanced_avg));
858  double sigma2 = sqrt(fabs(model2 / Ncorr - model_avg * model_avg));
859  double maxSigma = std::max(sigma1, sigma2);
860 
861  if (sigma1 < XMIPP_EQUAL_ACCURACY || sigma2 < XMIPP_EQUAL_ACCURACY
862  || (fabs(sigma1 - sigma2) / maxSigma > 0.9 && action>=5))
863  {
864  retval = heavy_penalization;
865  if (show_inf>=2)
866  std::cout << "Fitness2=" << heavy_penalization << " sigma1=" << sigma1 << " sigma2=" << sigma2 << std::endl;
867  }
868  else
869  {
870  correlation_coeff /= sigma1 * sigma2;
871  if (action == 7)
872  corr13 = correlation_coeff;
873  else
874  retval -= enhanced_weight * correlation_coeff;
875  if (show_inf >= 2)
876  {
877  std::cout << "model_avg=" << model_avg << std::endl;
878  std::cout << "enhanced_avg=" << enhanced_avg << std::endl;
879  std::cout << "enhanced_model=" << enhanced_model / Ncorr
880  << std::endl;
881  std::cout << "sigma1=" << sigma1 << std::endl;
882  std::cout << "sigma2=" << sigma2 << std::endl;
883  std::cout << "Fitness2="
884  << -(enhanced_weight * correlation_coeff)
885  << " (" << correlation_coeff << ")" << std::endl;
886  }
887  }
888  // Correlation of the derivative of the radial profile
889  if (action==3 || evaluation_reduction==1)
890  {
891  int state=0;
892  double maxDiff=0;
893  psd_theo_radial_derivative.initZeros();
894  double lowerlimt=1.1*min_freq;
895  double upperlimit=0.9*max_freq;
896  FOR_ALL_ELEMENTS_IN_ARRAY1D(psd_theo_radial)
897  if (A1D_ELEM(w_digfreq_r_iN,i)>0)
898  {
899  A1D_ELEM(psd_theo_radial,i)*=A1D_ELEM(w_digfreq_r_iN,i);
900  double freq=A2D_ELEM(w_digfreq,i,0);
901  switch (state)
902  {
903  case 0:
904  if (freq>lowerlimt)
905  state=1;
906  break;
907  case 1:
908  if (freq>upperlimit)
909  state=2;
910  else
911  {
912  double diff=A1D_ELEM(psd_theo_radial,i)-A1D_ELEM(psd_theo_radial,i-1);
913  A1D_ELEM(psd_theo_radial_derivative,i)=diff;
914  maxDiff=std::max(maxDiff,fabs(diff));
915  }
916  break;
917  }
918  }
919 
920  double corrRadialDerivative=0,mux=0, muy=0, Ncorr=0, sigmax=0, sigmay=0;
921  double iMaxDiff=1.0/maxDiff;
922  FOR_ALL_ELEMENTS_IN_ARRAY1D(psd_theo_radial)
923  {
924  A1D_ELEM(psd_theo_radial_derivative,i)*=iMaxDiff;
925  double x=A1D_ELEM(psd_exp_enhanced_radial_derivative,i);
926  double y=A1D_ELEM(psd_theo_radial_derivative,i);
927  corrRadialDerivative+=x*y;
928  mux+=x;
929  muy+=y;
930  sigmax+=x*x;
931  sigmay+=y*y;
932  Ncorr++;
933  }
934  if (Ncorr>0)
935  {
936  double iNcorr=1.0/Ncorr;
937  corrRadialDerivative*=iNcorr;
938  mux*=iNcorr;
939  muy*=iNcorr;
940  sigmax=sqrt(fabs(sigmax*iNcorr-mux*mux));
941  sigmay=sqrt(fabs(sigmay*iNcorr-muy*muy));
942  corrRadialDerivative=(corrRadialDerivative-mux*muy)/(sigmax*sigmay);
943  }
944 
945  retval-=corrRadialDerivative;
946 
947  if (show_inf>=2)
948  {
949  std::cout << "Fitness3=" << -corrRadialDerivative << std::endl;
950  if (show_inf==3)
951  {
952  psd_exp_radial.write("PPPexpRadial.txt");
953  psd_theo_radial.write("PPPtheoRadial.txt");
954  psd_exp_radial_derivative.write("PPPexpRadialDerivative.txt");
955  psd_theo_radial_derivative.write("PPPtheoRadialDerivative.txt");
956  }
957  }
958 
959  }
960 
961  }
962 
963  // Show some debugging information
964  if (show_inf >= 2)
965  {
966  std::cout << "Fitness=" << retval << std::endl;
967  if (show_inf == 3)
968  {
969  saveIntermediateResults("PPP");
970  std::cout << "Press any key\n";
971  char c;
972  std::cin >> c;
973  }
974  }
975 
976  return retval;
977 }
978 
979 double CTF_fitness(double *p, void *vprm)
980 {
981  auto *prm=(ProgCTFEstimateFromPSD *) vprm;
982  return prm->CTF_fitness_object(p);
983 }
984 
985 /* Compute central region -------------------------------------------------- */
986 void ProgCTFEstimateFromPSD::compute_central_region(double &w1, double &w2, double ang)
987 {
988  w1 = max_freq_psd;
989  w2 = min_freq_psd;
990  Matrix1D<double> freq(2), dir(2);
991 
992  // Compute first and third zero in the given direction
993  VECTOR_R2(dir, COSD(ang), SIND(ang));
994 
995  // Detect first zero
996  current_ctfmodel.lookFor(1, dir, freq, 0);
997  if (XX(freq) == -1 && YY(freq) == -1)
998  w1 = min_freq_psd;
999  else
1000  {
1001  contfreq2digfreq(freq, freq, Tm);
1002  double w;
1003  if (XX(dir) > 0.1)
1004  w = XX(freq) / XX(dir);
1005  else
1006  w = YY(freq) / YY(dir);
1007  w1 = XMIPP_MAX(min_freq_psd, XMIPP_MIN(w1, w));
1008  }
1009 
1010  // Detect fifth zero
1011  current_ctfmodel.lookFor(5, dir, freq, 0);
1012  if (XX(freq) == -1 && YY(freq) == -1)
1013  w2 = max_freq_psd;
1014  else
1015  {
1016  double w;
1017  contfreq2digfreq(freq, freq, Tm);
1018  if (XX(dir) > 0.1)
1019  w = XX(freq) / XX(dir);
1020  else
1021  w = YY(freq) / YY(dir);
1022  w2 = XMIPP_MIN(max_freq_psd, XMIPP_MAX(w2, w));
1023  }
1024 }
1025 
1026 /* Center focus ----------------------------------------------------------- */
1027 void ProgCTFEstimateFromPSD::center_optimization_focus(bool adjust_freq, bool adjust_th, double margin = 1)
1028 {
1029  if (show_optimization)
1030  std::cout << "Freq frame before focusing=" << min_freq_psd << ","
1031  << max_freq_psd << std::endl << "Value_th before focusing="
1032  << value_th << std::endl;
1033 
1034  double w1 = min_freq_psd, w2 = max_freq_psd;
1035  if (adjust_freq)
1036  {
1037  double w1U, w2U, w1V, w2V;
1038  compute_central_region(w1U, w2U, current_ctfmodel.azimuthal_angle);
1039  compute_central_region(w1V, w2V, current_ctfmodel.azimuthal_angle + 90);
1040  w1 = XMIPP_MIN(w1U, w1V);
1041  w2 = XMIPP_MAX(w2U, w2V);
1042  min_freq_psd = XMIPP_MAX(min_freq_psd, w1 - 0.05);
1043  max_freq_psd = XMIPP_MIN(max_freq_psd, w2 + 0.01);
1044  }
1045 
1046  // Compute maximum value within central region
1047  if (adjust_th)
1048  {
1049  Image<double> save;
1050  generateModelSoFar(save);
1051  double max_val = 0;
1052  FOR_ALL_ELEMENTS_IN_ARRAY2D(w_digfreq)
1053  {
1054  double w = w_digfreq(i, j);
1055  if (w >= w1 && w <= w2)
1056  max_val = XMIPP_MAX(max_val, save()(i, j));
1057  }
1058  if (value_th != -1)
1059  value_th = XMIPP_MIN(value_th, max_val * margin);
1060  else
1061  value_th = max_val * margin;
1062  }
1063 
1064  if (show_optimization)
1065  std::cout << "Freq frame after focusing=" << min_freq_psd << ","
1066  << max_freq_psd << std::endl << "Value_th after focusing="
1067  << value_th << std::endl;
1068 }
1069 
1070 // Estimate sqrt parameters ------------------------------------------------
1071 // Results are written in global_ctfmodel
1072 void ProgCTFEstimateFromPSD::estimate_background_sqrt_parameters()
1073 {
1074  if (show_optimization)
1075  std::cout << "Computing first sqrt background ...\n";
1076 
1077  // Estimate the base line taking the value of the CTF
1078  // for the maximum X and Y frequencies
1079  double base_line = 0;
1080  int N = 0;
1081 
1082  FOR_ALL_ELEMENTS_IN_ARRAY2D(w_digfreq)
1083  if (w_digfreq(i, j) > 0.4)
1084  {
1085  N++;
1086  base_line += (*f)(i, j);
1087 
1088  }
1089  if (0 == N) {
1090  REPORT_ERROR(ERR_NUMERICAL, "N is zero (0), which would lead to division by zero");
1091  }
1092  current_ctfmodel.base_line = base_line / N;
1093  // Find the linear least squares solution for the sqrt part
1094  Matrix2D<double> A(2, 2);
1095  A.initZeros();
1096  Matrix1D<double> b(2);
1097  b.initZeros();
1098  FOR_ALL_ELEMENTS_IN_ARRAY2D(w_digfreq)
1099  {
1100  if (mask(i, j) <= 0)
1101  continue;
1102 
1103  // Compute weight for this point
1104  double weight = 1 + max_freq_psd - w_digfreq(i, j);
1105 
1106  // Compute error
1107  current_ctfmodel.precomputeValues(x_contfreq(i, j),y_contfreq(i, j));
1108  double explained = current_ctfmodel.getValueNoiseAt();
1109  double unexplained = (*f)(i, j) - explained;
1110  if (unexplained <= 0)
1111  continue;
1112  unexplained = log(unexplained);
1113 
1114  double X = -sqrt(w_contfreq(i, j));
1115  A(0, 0) += weight * X * X;
1116  A(0, 1) += weight * X;
1117  A(1, 1) += weight * 1;
1118  b(0) += X * weight * unexplained;
1119  b(1) += weight * unexplained;
1120  }
1121  A(1, 0) = A(0, 1);
1122 
1123  b = A.inv() * b;
1124 
1125  current_ctfmodel.sqU = current_ctfmodel.sqV = b(0);
1126  current_ctfmodel.sqrt_K = exp(b(1));
1127  current_ctfmodel.sqrt_angle = 0;
1128 
1129  COPY_ctfmodel_TO_CURRENT_GUESS;
1130 
1131  if (show_optimization)
1132  {
1133  std::cout << "First SQRT Fit:\n" << current_ctfmodel << std::endl;
1134  saveIntermediateResults("step01a_first_sqrt_fit");
1135  }
1136 
1137  // Now optimize .........................................................
1138  double fitness;
1139  Matrix1D<double> steps;
1140  steps.resize(SQRT_CTF_PARAMETERS);
1141  steps.initConstant(1);
1142  // Optimize without penalization
1143  if (show_optimization)
1144  std::cout << "Looking for best fitting sqrt ...\n";
1145  penalize = false;
1146  int iter;
1147  powellOptimizer(*adjust_params, FIRST_SQRT_PARAMETER + 1,
1148  SQRT_CTF_PARAMETERS, CTF_fitness, global_prm, 0.05, fitness, iter, steps,
1149  show_optimization);
1150  // Optimize with penalization
1151  if (show_optimization)
1152  std::cout << "Penalizing best fitting sqrt ...\n";
1153  penalize = true;
1154  current_penalty = 2;
1155  int imax = CEIL(log(penalty) / log(2.0));
1156  for (int i = 1; i <= imax; i++)
1157  {
1158  if (show_optimization)
1159  std::cout << " Iteration " << i << " penalty="
1160  << current_penalty << std::endl;
1161  powellOptimizer(*adjust_params, FIRST_SQRT_PARAMETER + 1,
1162  SQRT_CTF_PARAMETERS, CTF_fitness, global_prm, 0.05, fitness, iter,
1163  steps, show_optimization);
1164  current_penalty *= 2;
1165  current_penalty =
1166  XMIPP_MIN(current_penalty, penalty);
1167  }
1168  // Keep the result in global_prm->adjust
1169  current_ctfmodel.forcePhysicalMeaning();
1170  COPY_ctfmodel_TO_CURRENT_GUESS;
1171 
1172  if (show_optimization)
1173  {
1174  std::cout << "Best penalized SQRT Fit:\n" << current_ctfmodel
1175  << std::endl;
1176  saveIntermediateResults("step01b_best_penalized_sqrt_fit");
1177  }
1178 
1179  center_optimization_focus(false, true, 1.5);
1180 }
1181 
1182 // Estimate gaussian parameters --------------------------------------------
1183 //#define DEBUG
1184 void ProgCTFEstimateFromPSD::estimate_background_gauss_parameters()
1185 {
1186 
1187  if (show_optimization)
1188  std::cout << "Computing first background Gaussian parameters ...\n";
1189 
1190  // Compute radial averages
1191  MultidimArray<double> radial_CTFmodel_avg(YSIZE(*f) / 2);
1192  MultidimArray<double> radial_CTFampl_avg(YSIZE(*f) / 2);
1193  MultidimArray<int> radial_N(YSIZE(*f) / 2);
1194  double w_max_gauss = 0.25;
1195  FOR_ALL_ELEMENTS_IN_ARRAY2D(w_digfreq)
1196  {
1197  if (mask(i, j) <= 0)
1198  continue;
1199  double w = w_digfreq(i, j);
1200  if (w > w_max_gauss)
1201  continue;
1202 
1203  int r = FLOOR(w * (double)YSIZE(*f));
1204  current_ctfmodel.precomputeValues(x_contfreq(i, j),y_contfreq(i, j));
1205  radial_CTFmodel_avg(r) += current_ctfmodel.getValueNoiseAt();
1206  radial_CTFampl_avg(r) += (*f)(i, j);
1207  radial_N(r)++;
1208  }
1209 
1210  // Compute the average radial error
1211  double error2_avg = 0;
1212  int N_avg = 0;
1213  MultidimArray<double> error;
1214  error.initZeros(radial_CTFmodel_avg);
1215  FOR_ALL_ELEMENTS_IN_ARRAY1D(radial_CTFmodel_avg)
1216  {
1217  if (radial_N(i) == 0)
1218  continue;
1219  error(i) = (radial_CTFampl_avg(i) - radial_CTFmodel_avg(i))
1220  / radial_N(i);
1221  error2_avg += error(i) * error(i);
1222  N_avg++;
1223  }
1224  if (N_avg != 0)
1225  error2_avg /= N_avg;
1226 
1227 #ifdef DEBUG
1228 
1229  std::cout << "Error2 avg=" << error2_avg << std::endl;
1230 #endif
1231 
1232  // Compute the minimum radial error
1233  bool first = true, OK_to_proceed = false;
1234  double error2_min = 0, wmin=0;
1235  FOR_ALL_ELEMENTS_IN_ARRAY1D(radial_CTFmodel_avg)
1236  {
1237  if (radial_N(i) == 0)
1238  continue;
1239  double w = w_digfreq(i, 0);
1240 
1241  if (error(i) < 0 && first)
1242  continue;
1243  else if (error(i) < 0)
1244  break;
1245  double error2 = error(i) * error(i);
1246  // If the two lines cross, do not consider any error until
1247  // the cross is "old" enough
1248  if (first && error2 > 0.15 * error2_avg)
1249  OK_to_proceed = true;
1250  if (first && i > 0)
1251  OK_to_proceed &= (error(i) < error(i - 1));
1252 
1253  // If the error now is bigger than a 30% (1.69=1.3*1.3) of the error min
1254  // this must be a rebound. Stop here
1255  if (!first && error2 > 1.69 * error2_min)
1256  break;
1257  if (first && OK_to_proceed)
1258  {
1259  wmin = w;
1260  error2_min = error2;
1261  first = false;
1262  }
1263  if (!first && error2 < error2_min)
1264  {
1265  wmin = w;
1266  error2_min = error2;
1267  }
1268 #ifdef DEBUG
1269  std::cout << w << " " << error2 << " " << wmin << " " << std::endl;
1270 #endif
1271 
1272  }
1273 
1274  // Compute the frequency of the minimum error
1275  max_gauss_freq = wmin;
1276 #ifdef DEBUG
1277 
1278  std::cout << "Freq of the minimum error: " << wmin << " " << fmin << std::endl;
1279 #endif
1280 
1281  // Compute the maximum radial error
1282  first = true;
1283  double error2_max = 0, wmax=0, fmax;
1284  FOR_ALL_ELEMENTS_IN_ARRAY1D(radial_CTFmodel_avg)
1285  {
1286  if (radial_N(i) == 0)
1287  continue;
1288  double w = w_digfreq(i, 0);
1289  if (w > wmin)
1290  continue;
1291 
1292  if (error(i) < 0 && first)
1293  continue;
1294  else if (error(i) < 0)
1295  break;
1296  double error2 = error(i) * error(i);
1297  if (first)
1298  {
1299  wmax = w;
1300  error2_max = error2;
1301  first = false;
1302  }
1303  if (error2 > error2_max)
1304  {
1305  wmax = w;
1306  error2_max = error2;
1307  }
1308 #ifdef DEBUG
1309  std::cout << w << " " << error2 << " " << wmax << std::endl;
1310 #endif
1311 
1312  }
1313  fmax = current_ctfmodel.cV = current_ctfmodel.cU = wmax / Tm;
1314 #ifdef DEBUG
1315 
1316  std::cout << "Freq of the maximum error: " << wmax << " " << fmax << std::endl;
1317 #endif
1318 
1319  // Find the linear least squares solution for the gauss part
1320  Matrix2D<double> A(2, 2);
1321  A.initZeros();
1322  Matrix1D<double> b(2);
1323  b.initZeros();
1324 
1325  FOR_ALL_ELEMENTS_IN_ARRAY2D(w_digfreq)
1326  {
1327  if (mask(i, j) <= 0)
1328  continue;
1329  if (w_digfreq(i, j) > wmin)
1330  continue;
1331  double fmod = w_contfreq(i, j);
1332 
1333  // Compute weight for this point
1334  double weight = 1 + max_freq_psd - w_digfreq(i, j);
1335  // Compute error
1336  current_ctfmodel.precomputeValues(x_contfreq(i, j),y_contfreq(i, j));
1337  double explained = current_ctfmodel.getValueNoiseAt();
1338 
1339  double unexplained = (*f)(i, j) - explained;
1340  if (unexplained <= 0)
1341  continue;
1342  unexplained = log(unexplained);
1343  double F = -(fmod - fmax) * (fmod - fmax);
1344 
1345  A(0, 0) += weight * 1;
1346  A(0, 1) += weight * F;
1347  A(1, 1) += weight * F * F;
1348  b(0) += weight * unexplained;
1349  b(1) += F * weight * unexplained;
1350  }
1351  A(1, 0) = A(0, 1);
1352  if ( (A(0, 0)== 0) && (A(1, 0)== 0) && (A(1, 1)== 0))
1353  {
1354  std::cout << "A matrix es zero" << std::endl;
1355  }
1356  else
1357  {
1358  b = A.inv() * b;
1359  current_ctfmodel.sigmaU = XMIPP_MIN(fabs(b(1)), 95e3); // This value should be
1360  current_ctfmodel.sigmaV = XMIPP_MIN(fabs(b(1)), 95e3); // conformant with the physical
1361  // meaning routine in CTF.cc
1362  current_ctfmodel.gaussian_K = exp(b(0));
1363  // Store the CTF values in global_prm->adjust
1364  current_ctfmodel.forcePhysicalMeaning();
1365  COPY_ctfmodel_TO_CURRENT_GUESS;
1366 
1367  if (show_optimization)
1368  {
1369  std::cout << "First Background Fit:\n" << current_ctfmodel << std::endl;
1370  saveIntermediateResults("step01c_first_background_fit");
1371  }
1372  center_optimization_focus(false, true, 1.5);
1373  }
1374 }
1375 #undef DEBUG
1376 
1377 // Estimate second gaussian parameters -------------------------------------
1378 //#define DEBUG
1379 void ProgCTFEstimateFromPSD::estimate_background_gauss_parameters2()
1380 {
1381  if (show_optimization)
1382  std::cout << "Computing first background Gaussian2 parameters ...\n";
1383 
1384  // Compute radial averages
1385  MultidimArray<double> radial_CTFmodel_avg(YSIZE(*f) / 2);
1386  MultidimArray<double> radial_CTFampl_avg(YSIZE(*f) / 2);
1387  MultidimArray<int> radial_N(YSIZE(*f) / 2);
1388  double w_max_gauss = 0.25;
1389  FOR_ALL_ELEMENTS_IN_ARRAY2D(w_digfreq)
1390  {
1391  if (mask(i, j) <= 0)
1392  continue;
1393  double w = w_digfreq(i, j);
1394  if (w > w_max_gauss)
1395  continue;
1396 
1397  int r = FLOOR(w * (double)YSIZE(*f));
1398  double f_x = DIRECT_A2D_ELEM(x_contfreq, i, j);
1399  double f_y = DIRECT_A2D_ELEM(y_contfreq, i, j);
1400  current_ctfmodel.precomputeValues(f_x, f_y);
1401  double bg = current_ctfmodel.getValueNoiseAt();
1402  double envelope = current_ctfmodel.getValueDampingAt();
1403  double ctf_without_damping =
1404  current_ctfmodel.getValuePureWithoutDampingAt();
1405  double ctf_with_damping = envelope * ctf_without_damping;
1406  double ctf2_th = bg + ctf_with_damping * ctf_with_damping;
1407  radial_CTFmodel_avg(r) += ctf2_th;
1408  radial_CTFampl_avg(r) += (*f)(i, j);
1409  radial_N(r)++;
1410  }
1411 
1412  // Compute the average radial error
1413  MultidimArray <double> error;
1414  error.initZeros(radial_CTFmodel_avg);
1415  FOR_ALL_ELEMENTS_IN_ARRAY1D(radial_CTFmodel_avg)
1416  {
1417  if (radial_N(i) == 0)
1418  continue;
1419  error(i) = (radial_CTFampl_avg(i) - radial_CTFmodel_avg(i))
1420  / radial_N(i);
1421  }
1422 #ifdef DEBUG
1423  std::cout << "Error:\n" << error << std::endl;
1424 #endif
1425 
1426  // Compute the frequency of the minimum error
1427  double wmin = 0.15;
1428 
1429  // Compute the maximum (negative) radial error
1430  double error_max = 0, wmax=0, fmax;
1431  FOR_ALL_ELEMENTS_IN_ARRAY1D(radial_CTFmodel_avg)
1432  {
1433  if (radial_N(i) == 0)
1434  continue;
1435  double w = w_digfreq(i, 0);
1436  if (w > wmin)
1437  break;
1438  if (error(i) < error_max)
1439  {
1440  wmax = w;
1441  error_max = error(i);
1442  }
1443  }
1444  fmax = current_ctfmodel.cV2 = current_ctfmodel.cU2 = wmax / Tm;
1445 #ifdef DEBUG
1446 
1447  std::cout << "Freq of the maximum error: " << wmax << " " << fmax << std::endl;
1448 #endif
1449 
1450  // Find the linear least squares solution for the gauss part
1451  Matrix2D<double> A(2, 2);
1452  A.initZeros();
1453  Matrix1D<double> b(2);
1454  b.initZeros();
1455  int N = 0;
1456  FOR_ALL_ELEMENTS_IN_ARRAY2D(w_digfreq)
1457  {
1458  if (mask(i, j) <= 0)
1459  continue;
1460  if (w_digfreq(i, j) > wmin)
1461  continue;
1462  double fmod = w_contfreq(i, j);
1463 
1464  // Compute the zero on the direction of this point
1465  Matrix1D<double> u(2), fzero(2);
1466  XX(u) = x_contfreq(i, j) / fmod;
1467  YY(u) = y_contfreq(i, j) / fmod;
1468  current_ctfmodel.lookFor(1, u, fzero, 0);
1469  if (fmod > fzero.module())
1470  continue;
1471 
1472  // Compute weight for this point
1473  double weight = 1 + max_freq_psd - w_digfreq(i, j);
1474 
1475  // Compute error
1476  double f_x = DIRECT_A2D_ELEM(x_contfreq, i, j);
1477  double f_y = DIRECT_A2D_ELEM(y_contfreq, i, j);
1478  current_ctfmodel.precomputeValues(f_x, f_y);
1479  double bg = current_ctfmodel.getValueNoiseAt();
1480  double envelope = current_ctfmodel.getValueDampingAt();
1481  double ctf_without_damping =
1482  current_ctfmodel.getValuePureWithoutDampingAt();
1483  double ctf_with_damping = envelope * ctf_without_damping;
1484  double ctf2_th = bg + ctf_with_damping * ctf_with_damping;
1485  double explained = ctf2_th;
1486  double unexplained = explained - (*f)(i, j);
1487 
1488  if (unexplained <= 0)
1489  continue;
1490  unexplained = log(unexplained);
1491  double F = -(fmod - fmax) * (fmod - fmax);
1492  A(0, 0) += weight * 1;
1493  A(0, 1) += weight * F;
1494  A(1, 1) += weight * F * F;
1495  b(0) += weight * unexplained;
1496  b(1) += F * weight * unexplained;
1497  N++;
1498  }
1499  if (N != 0)
1500  {
1501  A(1, 0) = A(0, 1);
1502 
1503  double det=A.det();
1504  if (fabs(det)>1e-9)
1505  {
1506  b = A.inv() * b;
1507  current_ctfmodel.sigmaU2 = XMIPP_MIN(fabs(b(1)), 95e3); // This value should be
1508  current_ctfmodel.sigmaV2 = XMIPP_MIN(fabs(b(1)), 95e3); // conformant with the physical
1509  // meaning routine in CTF.cc
1510  current_ctfmodel.gaussian_K2 = exp(b(0));
1511  }
1512  else
1513  {
1514  current_ctfmodel.sigmaU2 = current_ctfmodel.sigmaV2 = 0;
1515  current_ctfmodel.gaussian_K2 = 0;
1516  }
1517  }
1518  else
1519  {
1520  current_ctfmodel.sigmaU2 = current_ctfmodel.sigmaV2 = 0;
1521  current_ctfmodel.gaussian_K2 = 0;
1522  }
1523 
1524  // Store the CTF values in global_prm->adjust
1525  current_ctfmodel.forcePhysicalMeaning();
1526  COPY_ctfmodel_TO_CURRENT_GUESS;
1527 
1528 #ifdef DEBUG
1529  // Check
1530  FOR_ALL_ELEMENTS_IN_ARRAY2D(w_digfreq)
1531  {
1532  if (mask(i, j)<=0)
1533  continue;
1534  if (w_digfreq(i, j) > wmin)
1535  continue;
1536  double fmod = w_contfreq(i, j);
1537 
1538  // Compute the zero on the direction of this point
1539  Matrix1D<double> u(2), fzero(2);
1540  XX(u) = x_contfreq(i, j) / fmod;
1541  YY(u) = y_contfreq(i, j) / fmod;
1542  global_ctfmodel.zero(1, u, fzero);
1543  if (fmod > fzero.module())
1544  continue;
1545 
1546  // Compute error
1547  double f_x = DIRECT_A2D_ELEM(x_contfreq, i, j);
1548  double f_y = DIRECT_A2D_ELEM(y_contfreq, i, j);
1549  double bg = current_ctfmodel.getValueNoiseAt(f_x, f_y);
1550  double envelope = current_ctfmodel.getValueDampingAt(f_x, f_y);
1551  double ctf_without_damping = current_ctfmodel.getValuePureWithoutDampingAt(f_x, f_y);
1552  double ctf_with_damping = envelope * ctf_without_damping;
1553  double ctf2_th = bg + ctf_with_damping * ctf_with_damping;
1554  double explained = ctf2_th;
1555  double unexplained = explained - (*f)(i, j);
1556 
1557  if (unexplained <= 0)
1558  continue;
1559  std::cout << fmod << " " << unexplained << " "
1560  << current_ctfmodel.gaussian_K2*exp(-current_ctfmodel.sigmaU2*
1561  (fmod - fmax)*(fmod - fmax)) << std::endl;
1562  }
1563 #endif
1564 
1565  if (show_optimization)
1566  {
1567  std::cout << "First Background Gaussian 2 Fit:\n" << current_ctfmodel
1568  << std::endl;
1569  saveIntermediateResults("step04a_first_background2_fit");
1570  }
1571 }
1572 #undef DEBUG
1573 
1574 // Estimate envelope parameters --------------------------------------------
1575 //#define DEBUG
1576 void ProgCTFEstimateFromPSD::estimate_envelope_parameters()
1577 {
1578  if (show_optimization)
1579  std::cout << "Looking for best fitting envelope ...\n";
1580 
1581  // Set the envelope
1582  current_ctfmodel.Ca = initial_ctfmodel.Ca;
1583  current_ctfmodel.K = 1.0;
1584  current_ctfmodel.espr = 0.0;
1585  current_ctfmodel.ispr = 0.0;
1586  current_ctfmodel.alpha = 0.0;
1587  current_ctfmodel.DeltaF = 0.0;
1588  current_ctfmodel.DeltaR = 0.0;
1589  current_ctfmodel.Q0 = initial_ctfmodel.Q0;
1590  current_ctfmodel.envR0 = 0.0;
1591  current_ctfmodel.envR1 = 0.0;
1592  COPY_ctfmodel_TO_CURRENT_GUESS;
1593 
1594  // Now optimize the envelope
1595  penalize = false;
1596  int iter;
1597  double fitness;
1598  Matrix1D<double> steps;
1599  steps.resize(ENVELOPE_PARAMETERS);
1600  steps.initConstant(1);
1601  steps(1) = 0; // Do not optimize Cs
1602  steps(5) = 0; // Do not optimize for alpha, since Ealpha depends on the
1603  // defocus
1604  if (modelSimplification >= 1)
1605  steps(6) = steps(7) = 0; // Do not optimize DeltaF and DeltaR
1606  powellOptimizer(*adjust_params, FIRST_ENVELOPE_PARAMETER + 1,
1607  ENVELOPE_PARAMETERS, CTF_fitness, global_prm, 0.05, fitness, iter, steps,
1608  show_optimization);
1609 
1610  // Keep the result in global_prm->adjust
1611  current_ctfmodel.forcePhysicalMeaning();
1612  COPY_ctfmodel_TO_CURRENT_GUESS;
1613 
1614  if (show_optimization)
1615  {
1616  std::cout << "Best envelope Fit:\n" << current_ctfmodel << std::endl;
1617  saveIntermediateResults("step02a_best_envelope_fit");
1618  }
1619 
1620  // Optimize with penalization
1621  if (show_optimization)
1622  std::cout << "Penalizing best fitting envelope ...\n";
1623  penalize = true;
1624  current_penalty = 2;
1625  int imax = CEIL(log(penalty) / log(2.0));
1626  for (int i = 1; i <= imax; i++)
1627  {
1628  if (show_optimization)
1629  std::cout << " Iteration " << i << " penalty="
1630  << current_penalty << std::endl;
1631  powellOptimizer(*adjust_params, FIRST_ENVELOPE_PARAMETER + 1,
1632  ENVELOPE_PARAMETERS, CTF_fitness, global_prm, 0.05, fitness, iter,
1633  steps, show_optimization);
1634  current_penalty *= 2;
1635  current_penalty =
1636  XMIPP_MIN(current_penalty,penalty);
1637  }
1638  // Keep the result in global_prm->adjust
1639  current_ctfmodel.forcePhysicalMeaning();
1640  COPY_ctfmodel_TO_CURRENT_GUESS;
1641 
1642  if (show_optimization)
1643  {
1644  std::cout << "Best envelope Fit:\n" << current_ctfmodel << std::endl;
1645  saveIntermediateResults("step02b_best_penalized_envelope_fit");
1646  }
1647 }
1648 #undef DEBUG
1649 
1650 // Estimate defoci ---------------------------------------------------------
1651 void ProgCTFEstimateFromPSD::showFirstDefoci()
1652 {
1653  if (show_optimization)
1654  {
1655  std::cout << "First defocus Fit:\n" << current_ctfmodel << std::endl;
1656  saveIntermediateResults("step03a_first_defocus_fit");
1657  enhanced_ctftomodel.write("step03a_enhanced_PSD.xmp");
1658  Image<double> save, save2, save3;
1659  save().resize(YSIZE(w_digfreq), XSIZE(w_digfreq));
1660  save2().resize(save());
1661  save3().resize(save());
1662  FOR_ALL_ELEMENTS_IN_ARRAY2D(save())
1663  {
1664  save()(i, j) = enhanced_ctftomodel()(i, j);
1665  double f_x = DIRECT_A2D_ELEM(x_contfreq, i, j);
1666  double f_y = DIRECT_A2D_ELEM(y_contfreq, i, j);
1667  current_ctfmodel.precomputeValues(f_x, f_y);
1668  double ctf_without_damping =
1669  current_ctfmodel.getValuePureWithoutDampingAt();
1670  save2()(i, j) = ctf_without_damping * ctf_without_damping;
1671  save3()(i, j) = -enhanced_ctftomodel()(i, j)
1672  * ctf_without_damping * ctf_without_damping;
1673  }
1674  save.write("step03a_enhanced_PSD.xmp");
1675  save2.write("step03a_fitted_CTF.xmp");
1676  save3.write("step03a_superposition.xmp");
1677  }
1678 }
1679 
1680 //#define DEBUG
1681 void ProgCTFEstimateFromPSD::estimate_defoci()
1682 {
1683  if (show_optimization)
1684  std::cout << "Looking for first defoci ...\n";
1685  double best_defocusU=0, best_defocusV=0, best_angle=0, best_K=1;
1686  double best_error = heavy_penalization * 1.1;
1687  bool first = true;
1688  int i, j;
1689  double defocusV, defocusU;
1690 
1691  double defocusV0 = 1e3, defocusU0 = 1e3;
1692  double defocusVF = 100e3, defocusUF = 100e3;
1693  double initial_defocusStep = 8e3;
1694  MultidimArray<double> error;
1695 
1696  // Check if there is no initial guess
1697  double min_allowed_defocusU = 1e3, max_allowed_defocusU = 100e3;
1698  double min_allowed_defocusV = 1e3, max_allowed_defocusV = 100e3;
1699  if (initial_ctfmodel.DeltafU != 0)
1700  {
1701  initial_defocusStep = std::min(defocus_range,20000.0);
1702  defocusU0 = std::max(
1703  1e3,
1704  initial_ctfmodel.DeltafU
1705  - defocus_range);
1706  double maxDeviation = std::max(defocus_range,
1707  0.25 * initial_ctfmodel.DeltafU);
1708  max_allowed_defocusU = std::min(100e3,
1709  initial_ctfmodel.DeltafU + maxDeviation);
1710  defocusUF = std::min(
1711  150e3,
1712  initial_ctfmodel.DeltafU
1713  + defocus_range);
1714  min_allowed_defocusU = std::max(1e3,
1715  initial_ctfmodel.DeltafU - maxDeviation);
1716  if (initial_ctfmodel.DeltafV == 0)
1717  {
1718  defocusV0 = defocusU0;
1719  min_allowed_defocusV = min_allowed_defocusU;
1720  defocusVF = defocusUF;
1721  max_allowed_defocusV = max_allowed_defocusU;
1722  }
1723  else
1724  {
1725  defocusV0 = std::max(
1726  1e3,
1727  initial_ctfmodel.DeltafV
1728  - defocus_range);
1729  max_allowed_defocusV = std::max(100e3,
1730  initial_ctfmodel.DeltafV + maxDeviation);
1731  defocusVF = std::min(
1732  150e3,
1733  initial_ctfmodel.DeltafV
1734  + defocus_range);
1735  min_allowed_defocusV = std::max(1e3,
1736  initial_ctfmodel.DeltafV - maxDeviation);
1737  }
1738  }
1739 
1740  double K_so_far = current_ctfmodel.K;
1741  Matrix1D<double> steps(DEFOCUS_PARAMETERS);
1742  steps.initConstant(1);
1743  steps(3) = 0; // Do not optimize kV
1744  steps(4) = 0; // Do not optimize K
1745  for (double defocusStep = initial_defocusStep;
1746  defocusStep >= std::min(5000., defocus_range / 2);
1747  defocusStep /= 2)
1748  {
1749  error.resize(CEIL((defocusVF - defocusV0) / defocusStep + 1),
1750  CEIL((defocusUF - defocusU0) / defocusStep + 1));
1751  error.initConstant(heavy_penalization);
1752  if (show_optimization)
1753  std::cout << "V=[" << defocusV0 << "," << defocusVF << "]\n"
1754  << "U=[" << defocusU0 << "," << defocusUF << "]\n"
1755  << "Defocus step=" << defocusStep << std::endl;
1756  for (defocusV = defocusV0, i = 0; defocusV <= defocusVF; defocusV +=
1757  defocusStep, i++)
1758  {
1759  for (defocusU = defocusU0, j = 0; defocusU <= defocusUF; defocusU +=
1760  defocusStep, j++)
1761  {
1762  bool first_angle = true;
1763  if (fabs(defocusU - defocusV) > 30e3)
1764  {
1765  error(i, j) = heavy_penalization;
1766  continue;
1767  }
1768  for (double angle = 0; angle < 180; angle += 45)
1769  {
1770  int iter;
1771  double fitness;
1772 
1773  (*adjust_params)(0) = defocusU;
1774  (*adjust_params)(1) = defocusV;
1775  (*adjust_params)(2) = angle;
1776  (*adjust_params)(4) = K_so_far;
1777 
1778  powellOptimizer(*adjust_params, FIRST_DEFOCUS_PARAMETER + 1,
1779  DEFOCUS_PARAMETERS, CTF_fitness, global_prm, 0.05,
1780  fitness, iter, steps, false);
1781 
1782  if ((first_angle || fitness < error(i, j))
1783  && (current_ctfmodel.DeltafU >= min_allowed_defocusU
1784  && current_ctfmodel.DeltafU
1785  <= max_allowed_defocusU
1786  && current_ctfmodel.DeltafV
1787  >= min_allowed_defocusV
1788  && current_ctfmodel.DeltafV
1789  <= max_allowed_defocusV))
1790  {
1791  error(i, j) = fitness;
1792  first_angle = false;
1793  if (error(i, j) < best_error || first)
1794  {
1795  best_error = error(i, j);
1796  best_defocusU = current_ctfmodel.DeltafU;
1797  best_defocusV = current_ctfmodel.DeltafV;
1798  best_angle = current_ctfmodel.azimuthal_angle;
1799  best_K = current_ctfmodel.K;
1800  first = false;
1801  if (show_optimization)
1802  {
1803  std::cout << " (DefocusU,DefocusV)=("
1804  << defocusU << "," << defocusV
1805  << "), ang=" << angle << " --> ("
1806  << current_ctfmodel.DeltafU << ","
1807  << current_ctfmodel.DeltafV << "),"
1808  << current_ctfmodel.azimuthal_angle
1809  << " K=" << current_ctfmodel.K
1810  << " error=" << error(i, j)
1811  << std::endl;
1812 #ifdef DEBUG
1813 
1814  show_inf=3;
1815  CTF_fitness(adjust_params->vdata-1,NULL);
1816  show_inf=0;
1817 
1818  Image<double> save;
1819  save() = enhanced_ctftomodel();
1820  save.write("PPPenhanced.xmp");
1821  for (int i = 0; i < YSIZE(w_digfreq); i += 1)
1822  for (int j = 0; j < XSIZE(w_digfreq); j += 1)
1823  {
1824  if (DIRECT_A2D_ELEM(mask, i, j)<=0)
1825  continue;
1826  double f_x = DIRECT_A2D_ELEM(x_contfreq, i, j);
1827  double f_y = DIRECT_A2D_ELEM(y_contfreq, i, j);
1828  double envelope = current_ctfmodel.getValueDampingAt(f_x, f_y);
1829  double ctf_without_damping = current_ctfmodel.getValuePureWithoutDampingAt(f_x, f_y);
1830  double ctf_with_damping = envelope * ctf_without_damping;
1831  double ctf_with_damping2 = ctf_with_damping * ctf_with_damping;
1832  save(i, j) = ctf_with_damping2;
1833  }
1834  save.write("PPPctf2_with_damping2.xmp");
1835  std::cout << "Press any key\n";
1836  char c;
1837  std::cin >> c;
1838 #endif
1839 
1840  }
1841  }
1842  }
1843  }
1844  }
1845  }
1846 
1847  // Compute the range of the errors
1848  double errmin = error(0, 0), errmax = error(0, 0);
1849  bool aValidErrorHasBeenFound=false;
1850  for (int ii = STARTINGY(error); ii <= FINISHINGY(error); ii++)
1851  for (int jj = STARTINGX(error); jj <= FINISHINGX(error); jj++)
1852  {
1853  if (error(ii, jj) != heavy_penalization)
1854  {
1855  aValidErrorHasBeenFound=true;
1856  if (error(ii, jj) < errmin)
1857  errmin = error(ii, jj);
1858  else if (errmax == heavy_penalization)
1859  errmax = error(ii, jj);
1860  else if (error(ii, jj) > errmax)
1861  errmax = error(ii, jj);
1862  }
1863  }
1864  if (show_optimization)
1865  std::cout << "Error matrix\n" << error << std::endl;
1866  if (!aValidErrorHasBeenFound)
1867  REPORT_ERROR(ERR_NUMERICAL,"Cannot find any good defocus within the given range");
1868 
1869  // Find those defoci which are within a 10% of the maximum
1870  if (show_inf >= 2)
1871  std::cout << "Range=" << errmax - errmin << std::endl;
1872  double best_defocusVmin = best_defocusV, best_defocusVmax =
1873  best_defocusV;
1874  double best_defocusUmin = best_defocusU, best_defocusUmax =
1875  best_defocusU;
1876  for (defocusV = defocusV0, i = 0; defocusV <= defocusVF; defocusV +=
1877  defocusStep, i++)
1878  {
1879  for (defocusU = defocusU0, j = 0; defocusU <= defocusUF; defocusU +=
1880  defocusStep, j++)
1881  {
1882  if (show_inf >= 2)
1883  std::cout << i << "," << j << " " << error(i, j) << " "
1884  << defocusU << " " << defocusV << std::endl
1885  << best_defocusUmin << " " << best_defocusUmax
1886  << std::endl << best_defocusVmin << " "
1887  << best_defocusVmax << std::endl;
1888  if (fabs((error(i, j) - errmin) / (errmax - errmin)) <= 0.1)
1889  {
1890  if (defocusV < best_defocusVmin)
1891  best_defocusVmin = defocusV;
1892  if (defocusU < best_defocusUmin)
1893  best_defocusUmin = defocusU;
1894  if (defocusV > best_defocusVmax)
1895  best_defocusVmax = defocusV;
1896  if (defocusU > best_defocusUmax)
1897  best_defocusUmax = defocusU;
1898  }
1899  }
1900  }
1901 
1902  defocusVF = std::min(max_allowed_defocusV,
1903  best_defocusVmax + defocusStep);
1904  defocusV0 = std::max(min_allowed_defocusV,
1905  best_defocusVmin - defocusStep);
1906  defocusUF = std::min(max_allowed_defocusU,
1907  best_defocusUmax + defocusStep);
1908  defocusU0 = std::max(min_allowed_defocusU,
1909  best_defocusUmin - defocusStep);
1910  i = j = 0;
1911  if (show_inf >= 2)
1912  {
1913  Image<double> save;
1914  save() = error;
1915  save.write("error.xmp");
1916  std::cout << "Press any key: Error saved\n";
1917  char c;
1918  std::cin >> c;
1919  }
1920  }
1921 
1922  current_ctfmodel.DeltafU = best_defocusU;
1923  current_ctfmodel.DeltafV = best_defocusV;
1924  current_ctfmodel.azimuthal_angle = best_angle;
1925  current_ctfmodel.K = best_K;
1926 
1927  // Keep the result in global_prm->adjust
1928  current_ctfmodel.forcePhysicalMeaning();
1929  COPY_ctfmodel_TO_CURRENT_GUESS;
1930  ctfmodel_defoci = current_ctfmodel;
1931 
1932  showFirstDefoci();
1933 }
1934 #undef DEBUG
1935 
1936 // Estimate defoci with Zernike and SPTH transform---------------------------------------------
1937 void ProgCTFEstimateFromPSD::estimate_defoci_Zernike(const MultidimArray<double> &psdToModelFullSize, double min_freq, double max_freq, double Tm,
1938  double kV, double lambdaPhase, int sizeWindowPhase,
1939  double &defocusU, double &defocusV, double &ellipseAngle, int verbose)
1940 {
1941  // Center enhanced PSD
1942  MultidimArray<double> centeredEnhancedPSD=psdToModelFullSize;
1943  CenterFFT(centeredEnhancedPSD,true);
1944 
1945  // Estimate phase, modulation and Zernikes
1946  MultidimArray<double> mod, phase;
1947  Matrix1D<double> coefs(13);
1948  coefs.initConstant(0);
1949  VEC_ELEM(coefs,0) = 1;
1950  VEC_ELEM(coefs,3) = 1;
1951  VEC_ELEM(coefs,4) = 1;
1952  VEC_ELEM(coefs,5) = 1;
1953  VEC_ELEM(coefs,12) =1;
1954 
1955  auto x=(int)((0.3*max_freq+0.7*min_freq)*std::cos(PI/4)*XSIZE(centeredEnhancedPSD)+XSIZE(centeredEnhancedPSD)/2);
1956  DEBUG_TEXTFILE(formatString("Zernike1 %d",x));
1957  DEBUG_TEXTFILE(formatString("centeredEnhancedPSD80x80 %f",centeredEnhancedPSD(80,80)));
1958  DEBUG_TEXTFILE(formatString("centeredEnhancedPSD120x120 %f",centeredEnhancedPSD(120,120)));
1959  DEBUG_TEXTFILE(formatString("centeredEnhancedPSD160x160 %f",centeredEnhancedPSD(160,160)));
1960 
1961  int numElem = 10;
1962  kV = kV*1000;
1963  double K_so_far = current_ctfmodel.K;
1964  double lambda=12.2643247/std::sqrt(kV*(1.+0.978466e-6*kV));
1965  double Z8;
1966  double Z3;
1967  double Z4;
1968  double Z5;
1969  double eAngle=0.;
1970  double deFocusAvg;
1971  double deFocusDiff;
1972 
1973  max_freq = max_freq*0.8;
1974  double fmax = max_freq;
1975 
1976  Matrix1D<double> arrayDefocusDiff(numElem);
1977  arrayDefocusDiff.initZeros();
1978 
1979  Matrix1D<double> arrayDefocusAvg(numElem);
1980  arrayDefocusAvg.initZeros();
1981 
1982  Matrix1D<double> arrayError(numElem);
1983  arrayError.initConstant(-1);
1984 
1985  int iter;
1986  double fitness;
1987  Matrix1D<double> steps(DEFOCUS_PARAMETERS);
1988  steps.initConstant(1);
1989  steps(3) = 0; // Do not optimize kV
1990  steps(4) = 0; // Do not optimize K
1991  double fmaxStep = (max_freq-min_freq)/numElem;
1992 
1993  lambdaPhase = 0.8;
1994  sizeWindowPhase = 10;
1995 
1996  Matrix1D<double> initialGlobalAdjust = *adjust_params;
1997 
1998  for (int i = 1; i < numElem; i++)
1999  {
2000  if (((fmax - min_freq)/min_freq) > 0.5)
2001  {
2002  demodulate(centeredEnhancedPSD,lambdaPhase,sizeWindowPhase,
2003  x,x,
2004  (int)(min_freq*XSIZE(centeredEnhancedPSD)),
2005  (int)(fmax*XSIZE(centeredEnhancedPSD)),
2006  phase, mod, coefs, 0);
2007 
2008  Z8=VEC_ELEM(coefs,4);
2009  Z3=VEC_ELEM(coefs,12);
2010  Z4=VEC_ELEM(coefs,3);
2011  Z5=VEC_ELEM(coefs,5);
2012 
2013  eAngle = 0.5*RAD2DEG(std::atan2(Z5,Z4))+90.0;
2014  deFocusAvg = fabs(2*Tm*Tm*(2*Z3-6*Z8)/(PI*lambda));
2015  deFocusDiff = fabs(2*Tm*Tm*(std::sqrt(Z4*Z4+Z5*Z5))/(PI*lambda));
2016 
2017  coefs.initConstant(0);
2018  VEC_ELEM(coefs,0) = 1;
2019  VEC_ELEM(coefs,3) = 1;
2020  VEC_ELEM(coefs,4) = 1;
2021  VEC_ELEM(coefs,5) = 1;
2022  VEC_ELEM(coefs,12) =1;
2023 
2024  fmax -= fmaxStep;
2025 
2026  (*adjust_params)(0) = deFocusAvg+deFocusDiff;
2027  (*adjust_params)(1) = deFocusAvg-deFocusDiff;
2028  (*adjust_params)(2) = eAngle;
2029  (*adjust_params)(4) = K_so_far;
2030  (*adjust_params)(6) = 2;
2031 
2032 
2033  fitness =0;
2034  powellOptimizer(*adjust_params, FIRST_DEFOCUS_PARAMETER + 1,
2035  DEFOCUS_PARAMETERS, CTF_fitness, global_prm, 0.05,
2036  fitness, iter, steps, false);
2037 
2038  VEC_ELEM(arrayDefocusAvg,i) = ((*adjust_params)(0) +(*adjust_params)(1))/2;
2039  VEC_ELEM(arrayDefocusDiff,i) = ((*adjust_params)(0) -(*adjust_params)(1))/2;
2040  VEC_ELEM(arrayError,i) = (-1)*fitness;
2041 
2042  }
2043  }
2044 
2045  int maxInd=arrayError.maxIndex();
2046 
2047  while ( (VEC_ELEM(arrayDefocusAvg,maxInd) < 3000) || ((VEC_ELEM(arrayDefocusAvg,maxInd) > 50000) && VEC_ELEM(arrayError,maxInd)>-1e3 ))
2048  {
2049  VEC_ELEM(arrayError,maxInd) = -1e3;
2050  VEC_ELEM(arrayDefocusAvg,maxInd) = initial_ctfmodel.DeltafU;
2051  VEC_ELEM(arrayDefocusDiff,maxInd) = initial_ctfmodel.DeltafV;
2052  maxInd=arrayError.maxIndex();
2053  }
2054  if (VEC_ELEM(arrayError,maxInd)<=-1e3)
2055  {
2056  estimate_defoci();
2057  return;
2058  }
2059 
2060  Matrix1D<double> arrayDefocusU(3);
2061  arrayDefocusU.initZeros();
2062 
2063  Matrix1D<double> arrayDefocusV(3);
2064  arrayDefocusV.initZeros();
2065 
2066  Matrix1D<double> arrayError2(3);
2067  arrayError2.initConstant(-1);
2068 
2069  //We want to take care about more parameters
2070  // We optimize for deltaU, deltaV
2071  (*adjust_params)(0) = VEC_ELEM(arrayDefocusAvg,maxInd)+VEC_ELEM(arrayDefocusDiff,maxInd);
2072  (*adjust_params)(1) = VEC_ELEM(arrayDefocusAvg,maxInd)-VEC_ELEM(arrayDefocusDiff,maxInd);
2073  (*adjust_params)(2) = eAngle;
2074  (*adjust_params)(4) = K_so_far;
2075  (*adjust_params)(6) = 2;
2076 
2077  fitness =0;
2078  powellOptimizer(*adjust_params, FIRST_DEFOCUS_PARAMETER + 1,
2079  DEFOCUS_PARAMETERS, CTF_fitness, global_prm, 0.05,
2080  fitness, iter, steps, false);
2081 
2082  VEC_ELEM(arrayDefocusU,0) = (*adjust_params)(0);
2083  VEC_ELEM(arrayDefocusV,0) = (*adjust_params)(1);
2084  VEC_ELEM(arrayError2,0) = (-1)*fitness;
2085 
2086  // We optimize for deltaU, deltaU
2087  (*adjust_params)(0) = VEC_ELEM(arrayDefocusAvg,maxInd)+VEC_ELEM(arrayDefocusDiff,maxInd);
2088  (*adjust_params)(1) = VEC_ELEM(arrayDefocusAvg,maxInd)+VEC_ELEM(arrayDefocusDiff,maxInd);
2089  (*adjust_params)(2) = eAngle;
2090  (*adjust_params)(4) = K_so_far;
2091  (*adjust_params)(6) = 2;
2092 
2093  fitness =0;
2094  powellOptimizer(*adjust_params, FIRST_DEFOCUS_PARAMETER + 1,
2095  DEFOCUS_PARAMETERS, CTF_fitness, global_prm, 0.05,
2096  fitness, iter, steps, false);
2097 
2098  VEC_ELEM(arrayDefocusU,1) = (*adjust_params)(0);
2099  VEC_ELEM(arrayDefocusV,1) = (*adjust_params)(1);
2100  VEC_ELEM(arrayError2,1) = (-1)*fitness;
2101 
2102  // We optimize for deltaV, deltaV
2103  (*adjust_params)(0) = VEC_ELEM(arrayDefocusAvg,maxInd)-VEC_ELEM(arrayDefocusDiff,maxInd);
2104  (*adjust_params)(1) = VEC_ELEM(arrayDefocusAvg,maxInd)-VEC_ELEM(arrayDefocusDiff,maxInd);
2105  (*adjust_params)(2) = eAngle;
2106  (*adjust_params)(4) = K_so_far;
2107  (*adjust_params)(6) = 2;
2108 
2109  fitness =0;
2110  powellOptimizer(*adjust_params, FIRST_DEFOCUS_PARAMETER + 1,
2111  DEFOCUS_PARAMETERS, CTF_fitness, global_prm, 0.05,
2112  fitness, iter, steps, false);
2113 
2114  VEC_ELEM(arrayDefocusU,2) = (*adjust_params)(0);
2115  VEC_ELEM(arrayDefocusV,2) = (*adjust_params)(1);
2116  VEC_ELEM(arrayError2,2) = (-1)*fitness;
2117 
2118  //Here we select the best one
2119  maxInd=arrayError2.maxIndex();
2120  defocusU = VEC_ELEM(arrayDefocusU,maxInd);
2121  defocusV = VEC_ELEM(arrayDefocusV,maxInd);
2122 
2123  (*adjust_params)(0) = defocusU;
2124  (*adjust_params)(1) = defocusV;
2125  (*adjust_params)(2) = eAngle;
2126  (*adjust_params)(4) = K_so_far;
2127  (*adjust_params)(6) = 2;
2128 
2129  while ( (0.5*(defocusU+defocusV) < 2500) || (0.5*(defocusU+defocusV) > 60000) )
2130  {
2131  VEC_ELEM(arrayError2,maxInd) = -1e3;
2132  VEC_ELEM(arrayDefocusU,maxInd) = initial_ctfmodel.DeltafU;
2133  VEC_ELEM(arrayDefocusV,maxInd) = initial_ctfmodel.DeltafV;
2134 
2135  maxInd=arrayError2.maxIndex();
2136  defocusU = VEC_ELEM(arrayDefocusU,maxInd);
2137  defocusV = VEC_ELEM(arrayDefocusV,maxInd);
2138  (*adjust_params)(0) = defocusU;
2139  (*adjust_params)(1) = defocusV;
2140  (*adjust_params)(2) = eAngle;
2141  (*adjust_params)(4) = K_so_far;
2142  (*adjust_params)(6) = 2;
2143  }
2144 
2145  if (VEC_ELEM(arrayError2,maxInd) <= 0)
2146  {
2147  COPY_ctfmodel_TO_CURRENT_GUESS;
2148  ctfmodel_defoci = current_ctfmodel;
2149 
2150  action = 5;
2151 
2152  steps.resize(ALL_CTF_PARAMETERS);
2153  steps.initConstant(1);
2154  steps(3) = 0; // kV
2155  steps(5) = 0; // The spherical aberration (Cs) is not optimized
2156  if (current_ctfmodel.Q0!=0)
2157  steps(12)=0;
2158 
2159  COPY_ctfmodel_TO_CURRENT_GUESS;
2160 
2161  evaluation_reduction = 2;
2162  powellOptimizer(*adjust_params, 0 + 1, ALL_CTF_PARAMETERS, CTF_fitness,
2163  global_prm, 0.01, fitness, iter, steps,show_optimization);
2164  COPY_ctfmodel_TO_CURRENT_GUESS;
2165 
2166  show_inf=0;
2167  action = 3;
2168  evaluation_reduction = 1;
2169 
2170  double error = -CTF_fitness(adjust_params->vdata-1,nullptr);
2171  if ( error <= -0.1)
2172  {
2173  *adjust_params = initialGlobalAdjust;
2174  COPY_ctfmodel_TO_CURRENT_GUESS;
2175  //There is nothing to do and we have to perform an exhaustive search
2176 #ifndef RELEASE_MODE
2177  std::cout << " Entering in estimate_defoci, Performing exhaustive defocus search (SLOW)" << std::endl;
2178 #endif
2179  estimate_defoci();
2180  }
2181  }
2182 }
2183 
2184 void ProgCTFEstimateFromPSD::estimate_defoci_Zernike()
2185 {
2186  if (show_optimization)
2187  std::cout << "Looking for first defoci ...\n";
2188 
2189  DEBUG_TEXTFILE("Step 6.1");
2190  DEBUG_MODEL_TEXTFILE;
2191  estimate_defoci_Zernike(enhanced_ctftomodel_fullsize(),
2192  min_freq,max_freq,Tm,
2193  initial_ctfmodel.kV,
2194  lambdaPhase,sizeWindowPhase,
2195  current_ctfmodel.DeltafU, current_ctfmodel.DeltafV, current_ctfmodel.azimuthal_angle, 0);
2196  DEBUG_TEXTFILE("Step 6.2");
2197  DEBUG_MODEL_TEXTFILE;
2198 
2199  current_ctfmodel.forcePhysicalMeaning();
2200  COPY_ctfmodel_TO_CURRENT_GUESS;
2201  ctfmodel_defoci = current_ctfmodel;
2202 
2203  showFirstDefoci();
2204 }
2205 
2206 /* Main routine ------------------------------------------------------------ */
2207 //#define DEBUG
2208 double ROUT_Adjust_CTF(ProgCTFEstimateFromPSD &prm, CTFDescription &output_ctfmodel, bool standalone)
2209 {
2210  DEBUG_OPEN_TEXTFILE(prm.fn_psd.removeLastExtension());
2211  global_prm = &prm;
2212  if (standalone || prm.show_optimization)
2213  prm.show();
2214  prm.produceSideInfo();
2215  DEBUG_TEXTFILE(formatString("After producing side info: Avg=%f",prm.ctftomodel().computeAvg()));
2216  DEBUG_MODEL_TEXTFILE;
2217 
2218  // Build initial frequency mask
2219  prm.value_th = -1;
2220  prm.min_freq_psd = prm.min_freq;
2221  prm.max_freq_psd = prm.max_freq;
2222 
2223  // Set some global variables
2224  prm.adjust_params = &prm.adjust;
2225  prm.penalize = false;
2226  prm.max_gauss_freq = 0;
2227  prm.heavy_penalization = prm.f->computeMax() * XSIZE(*prm.f) * YSIZE(*prm.f);
2228  prm.show_inf = 0;
2229 
2230  // Some variables needed by all steps
2231  int iter;
2232  double fitness;
2233  Matrix1D<double> steps;
2234 
2235  /************************************************************************/
2236  /* STEPs 1, 2, 3 and 4: Find background which best fits the CTF */
2237  /************************************************************************/
2238 
2239  prm.current_ctfmodel.enable_CTFnoise = true;
2240  prm.current_ctfmodel.enable_CTF = false;
2241  prm.evaluation_reduction = 4;
2242 
2243  // If initial parameters were not supplied for the gaussian curve,
2244  // estimate them from the CTF file
2245  prm.action = 0;
2246  if (prm.adjust(FIRST_SQRT_PARAMETER) == 0)
2247  {
2248  prm.estimate_background_sqrt_parameters();
2249  prm.estimate_background_gauss_parameters();
2250  }
2251 
2252  // Optimize the current background
2253  prm.action = 1;
2254  prm.penalize = true;
2255  prm.current_penalty = prm.penalty;
2256  steps.resize(BACKGROUND_CTF_PARAMETERS);
2257  steps.initConstant(1);
2258  if (!prm.modelSimplification >= 3)
2259  steps(7) = steps(8) = steps(10) = 0;
2260  powellOptimizer(*prm.adjust_params, FIRST_SQRT_PARAMETER + 1,
2261  BACKGROUND_CTF_PARAMETERS, CTF_fitness, global_prm, 0.01, fitness, iter,
2262  steps, prm.show_optimization);
2263 
2264  // Make sure that the model has physical meaning
2265  // (In some machines due to numerical imprecission this check is necessary
2266  // at the end)
2267  prm.current_ctfmodel.forcePhysicalMeaning();
2268  COPY_ctfmodel_TO_CURRENT_GUESS;
2269 
2270  if (prm.show_optimization)
2271  {
2272  std::cout << "Best background Fit:\n" << prm.current_ctfmodel << std::endl;
2273  prm.saveIntermediateResults("step01d_best_background_fit");
2274  }
2275  DEBUG_TEXTFILE(formatString("Step 4: CTF_fitness=%f",CTF_fitness));
2276  DEBUG_MODEL_TEXTFILE;
2277 
2278  /************************************************************************/
2279  /* STEPs 5 and 6: Find envelope which best fits the CTF */
2280  /************************************************************************/
2281  prm.action = 2;
2282  prm.current_ctfmodel.enable_CTF = true;
2283  prm.current_ctfmodel.phase_shift = prm.initial_ctfmodel.phase_shift;
2284  prm.current_ctfmodel.VPP_radius = prm.initial_ctfmodel.VPP_radius;
2285  if (prm.initial_ctfmodel.K == 0)
2286  {
2287  prm.current_ctfmodel.kV = prm.initial_ctfmodel.kV;
2288  prm.current_ctfmodel.Cs = prm.initial_ctfmodel.Cs;
2289  if (prm.initial_ctfmodel.Q0 != 0)
2290  prm.current_ctfmodel.Q0 = prm.initial_ctfmodel.Q0;
2291  prm.estimate_envelope_parameters();
2292  }
2293  else
2294  {
2295  prm.current_ctfmodel.K = prm.initial_ctfmodel.K;
2296  prm.current_ctfmodel.kV = prm.initial_ctfmodel.kV;
2297  prm.current_ctfmodel.DeltafU = prm.initial_ctfmodel.DeltafU;
2298  prm.current_ctfmodel.DeltafV = prm.initial_ctfmodel.DeltafV;
2299  prm.current_ctfmodel.azimuthal_angle = prm.initial_ctfmodel.azimuthal_angle;
2300  prm.current_ctfmodel.Cs = prm.initial_ctfmodel.Cs;
2301  prm.current_ctfmodel.Ca = prm.initial_ctfmodel.Ca;
2302  prm.current_ctfmodel.espr = prm.initial_ctfmodel.espr;
2303  prm.current_ctfmodel.ispr = prm.initial_ctfmodel.ispr;
2304  prm.current_ctfmodel.alpha = prm.initial_ctfmodel.alpha;
2305  prm.current_ctfmodel.DeltaF = prm.initial_ctfmodel.DeltaF;
2306  prm.current_ctfmodel.DeltaR = prm.initial_ctfmodel.DeltaR;
2307  prm.current_ctfmodel.Q0 = prm.initial_ctfmodel.Q0;
2308  COPY_ctfmodel_TO_CURRENT_GUESS;
2309  }
2310  DEBUG_TEXTFILE(formatString("Step 6: espr=%f",prm.current_ctfmodel.espr));
2311  DEBUG_MODEL_TEXTFILE;
2312  /************************************************************************/
2313  /* STEP 7: the defocus and angular parameters */
2314  /************************************************************************/
2315 
2316  prm.action = 3;
2317  prm.evaluation_reduction = 1;
2318  if (prm.fastDefocusEstimate)
2319  prm.estimate_defoci_Zernike();
2320  else if (!prm.noDefocusEstimate)
2321  prm.estimate_defoci();
2322  else
2323  {
2324  prm.current_ctfmodel.DeltafU=prm.initial_ctfmodel.DeltafU;
2325  prm.current_ctfmodel.DeltafV=prm.initial_ctfmodel.DeltafV;
2326  if (prm.current_ctfmodel.DeltafV==0.0)
2327  prm.current_ctfmodel.DeltafV=prm.current_ctfmodel.DeltafU;
2328  }
2329 
2330  DEBUG_TEXTFILE(formatString("Step 7: DeltafU=%f",prm.current_ctfmodel.DeltafU));
2331  DEBUG_TEXTFILE(formatString("Step 7: DeltafV=%f",prm.current_ctfmodel.DeltafV));
2332  DEBUG_TEXTFILE(formatString("Step 7: azimutalAngle=%f",prm.current_ctfmodel.azimuthal_angle));
2333  DEBUG_MODEL_TEXTFILE;
2334 
2335  //This line is to test the results obtained
2336  //exit(1);
2337 
2338  /************************************************************************/
2339  /* STEPs 9, 10 and 11: all parameters included second Gaussian */
2340  /************************************************************************/
2341  prm.action = 5;
2342  if (prm.modelSimplification < 2)
2343  prm.estimate_background_gauss_parameters2();
2344 
2345  steps.resize(ALL_CTF_PARAMETERS);
2346  steps.initConstant(1);
2347  steps(3) = 0; // kV
2348  steps(5) = 0; // The spherical aberration (Cs) is not optimized
2349  steps(37) = 0; //VPP radius not optimized
2350  if (prm.initial_ctfmodel.Q0 != 0)
2351  steps(15) = 0; // Q0
2352  if (prm.modelSimplification >= 3)
2353  steps(20) = steps(21) = steps(23) = 0;
2354  if (prm.modelSimplification >= 2)
2355  steps(24) = steps(25) = steps(26) = steps(27) = steps(28) = steps(29) = 0;
2356  if (prm.modelSimplification >= 1)
2357  steps(10) = steps(11) = 0;
2358  if (std::floor(prm.initial_ctfmodel.VPP_radius) == 0)
2359  steps(36) = 0; //VPP phase shift
2360 
2361  powellOptimizer(*prm.adjust_params, 0 + 1, ALL_CTF_PARAMETERS, CTF_fitness,
2362  global_prm, 0.01, fitness, iter, steps, prm.show_optimization);
2363 
2364  prm.current_ctfmodel.forcePhysicalMeaning();
2365  COPY_ctfmodel_TO_CURRENT_GUESS;
2366 
2367  if (prm.show_optimization)
2368  {
2369  std::cout << "Best fit with Gaussian2:\n" << prm.current_ctfmodel << std::endl;
2370  prm.saveIntermediateResults("step04b_best_fit_with_gaussian2");
2371  }
2372 
2373  prm.evaluation_reduction = 2;
2374  powellOptimizer(*prm.adjust_params, 0 + 1, ALL_CTF_PARAMETERS, CTF_fitness,
2375  global_prm, 0.01, fitness, iter, steps, prm.show_optimization);
2376  prm.current_ctfmodel.forcePhysicalMeaning();
2377  COPY_ctfmodel_TO_CURRENT_GUESS;
2378 
2379  prm.evaluation_reduction = 1;
2380  powellOptimizer(*prm.adjust_params, 0 + 1, ALL_CTF_PARAMETERS, CTF_fitness,
2381  global_prm, 0.005, fitness, iter, steps, prm.show_optimization);
2382  prm.current_ctfmodel.forcePhysicalMeaning();
2383  COPY_ctfmodel_TO_CURRENT_GUESS;
2384 
2385  if (prm.refineAmplitudeContrast)
2386  {
2387  std::cout << "Refining amplitude contrast ..." << std::endl;
2388  steps(15) = 1; // Q0
2389  powellOptimizer(*prm.adjust_params, 0 + 1, ALL_CTF_PARAMETERS, CTF_fitness,
2390  global_prm, 0.005, fitness, iter, steps, prm.show_optimization);
2391  prm.current_ctfmodel.forcePhysicalMeaning();
2392  COPY_ctfmodel_TO_CURRENT_GUESS;
2393  }
2394 
2395  if (prm.show_optimization)
2396  {
2397  std::cout << "Best fit:\n" << prm.current_ctfmodel << std::endl;
2398  prm.saveIntermediateResults("step04c_best_fit");
2399  }
2400  DEBUG_TEXTFILE(formatString("Step 11: DeltafU=%f fitness=%f",prm.current_ctfmodel.DeltafU,fitness));
2401  DEBUG_MODEL_TEXTFILE;
2402 
2403  //We adopt that always DeltafU > DeltafV so if this is not the case we change the values and the angle
2404  if ( prm.current_ctfmodel.DeltafV > prm.current_ctfmodel.DeltafU)
2405  {
2406  double temp;
2407  temp = prm.current_ctfmodel.DeltafU;
2408  prm.current_ctfmodel.DeltafU = prm.current_ctfmodel.DeltafV;
2409  prm.current_ctfmodel.DeltafV = temp;
2410  prm.current_ctfmodel.azimuthal_angle -= 90;
2411  COPY_ctfmodel_TO_CURRENT_GUESS;
2412  }
2413 
2414  /************************************************************************/
2415  /* STEP 12: Produce output */
2416  /************************************************************************/
2417 
2418  prm.action = 7;
2419 
2420  if (prm.fn_psd != "")
2421  {
2422  // Define mask between first and third zero
2423  prm.mask_between_zeroes.initZeros(prm.mask);
2424  Matrix1D<double> u(2), z1(2), z3(2);
2425  FOR_ALL_ELEMENTS_IN_ARRAY2D(prm.mask_between_zeroes)
2426  {
2427  VECTOR_R2(u, prm.x_digfreq(i, j), prm.y_digfreq(i, j));
2428  u /= u.module();
2429  prm.current_ctfmodel.lookFor(1, u, z1, 0);
2430  prm.current_ctfmodel.lookFor(3, u, z3, 0);
2431  if (z1.module() < prm.w_contfreq(i, j)
2432  && prm.w_contfreq(i, j) < z3.module())
2433  prm.mask_between_zeroes(i, j) = 1;
2434  }
2435  // Evaluate the correlation in this region
2436  CTF_fitness(prm.adjust_params->adaptForNumericalRecipes(), &prm);
2437  // Save results
2438  FileName fn_rootCTFPARAM = prm.fn_psd.withoutExtension();
2439 
2440  FileName fn_rootMODEL = fn_rootCTFPARAM;
2441  size_t atPosition=fn_rootCTFPARAM.find('@');
2442 
2443  if (atPosition!=std::string::npos)
2444  {
2445  fn_rootMODEL=formatString("%03d@%s",textToInteger(fn_rootCTFPARAM.substr(0, atPosition)),
2446  fn_rootCTFPARAM.substr(atPosition+1).c_str());
2447  fn_rootCTFPARAM=formatString("region%03d@%s",textToInteger(fn_rootCTFPARAM.substr(0, atPosition)),
2448  fn_rootCTFPARAM.substr(atPosition+1).c_str());
2449  }
2450  else
2451  fn_rootCTFPARAM=(String)"fullMicrograph@"+fn_rootCTFPARAM;
2452 
2453  prm.saveIntermediateResults(fn_rootMODEL, false);
2454  prm.current_ctfmodel.Tm /= prm.downsampleFactor;
2455  prm.current_ctfmodel.azimuthal_angle = std::fmod(prm.current_ctfmodel.azimuthal_angle,360.);
2456  prm.current_ctfmodel.write(fn_rootCTFPARAM + ".ctfparam_tmp");
2457  MetaDataVec MD;
2458  MD.read(fn_rootCTFPARAM + ".ctfparam_tmp");
2459  size_t id = MD.firstRowId();
2460  MD.setValue(MDL_CTF_X0, (double)output_ctfmodel.x0*prm.Tm, id);
2461  MD.setValue(MDL_CTF_XF, (double)output_ctfmodel.xF*prm.Tm, id);
2462  MD.setValue(MDL_CTF_Y0, (double)output_ctfmodel.y0*prm.Tm, id);
2463  MD.setValue(MDL_CTF_YF, (double)output_ctfmodel.yF*prm.Tm, id);
2464  MD.setValue(MDL_CTF_CRIT_FITTINGSCORE, fitness, id);
2465  MD.setValue(MDL_CTF_CRIT_FITTINGCORR13, prm.corr13, id);
2466  MD.setValue(MDL_CTF_CRIT_ICENESS, evaluateIceness(prm.psd_exp_radial, prm.Tm), id);
2467  MD.setValue(MDL_CTF_DOWNSAMPLE_PERFORMED, prm.downsampleFactor, id);
2468  MD.write(fn_rootCTFPARAM + ".ctfparam",MD_APPEND);
2469  fn_rootCTFPARAM = fn_rootCTFPARAM + ".ctfparam_tmp";
2470  fn_rootCTFPARAM.deleteFile();
2471  }
2472  output_ctfmodel = prm.current_ctfmodel;
2473 
2474  DEBUG_CLOSE_TEXTFILE;
2475  return fitness;
2476 }
2477 
2478 void ProgCTFEstimateFromPSD::run()
2479 {
2480  CTFDescription ctfmodel;
2481  ROUT_Adjust_CTF(*this, ctfmodel);
2482 }
2483 
2484 ProgCTFEstimateFromPSD::ProgCTFEstimateFromPSD(const ProgCTFEstimateFromPSDFast *copy)
2485 {
2486  action = copy->action;
2487  x_contfreq = copy->x_contfreq;
2488  y_contfreq = copy->y_contfreq;
2489  w_contfreq = copy->w_contfreq;
2490  x_digfreq = copy->x_digfreq;
2491  y_digfreq = copy->y_digfreq;
2492  w_digfreq = copy->w_digfreq;
2494  psd_exp_radial = copy->psd_exp_radial;
2495  psd_exp_enhanced_radial = copy->psd_exp_enhanced_radial;
2496  psd_exp_enhanced_radial_derivative = copy->psd_exp_enhanced_radial_derivative;
2497  psd_theo_radial_derivative = copy->psd_theo_radial_derivative;
2498  psd_exp_radial_derivative = copy->psd_exp_radial_derivative;
2499  psd_theo_radial = copy->psd_theo_radial;
2500  w_digfreq_r_iN = copy->w_digfreq_r_iN;
2501  w_digfreq_r = copy->w_digfreq_r;
2503  mask = copy->mask;
2504  mask_between_zeroes = copy->mask_between_zeroes;
2505  max_freq = copy->max_freq;
2506  min_freq = copy->min_freq;
2507  min_freq_psd = copy->min_freq_psd;
2508  max_freq_psd = copy->max_freq_psd;
2509  corr13 = copy->corr13;
2510 
2511  show_inf = copy->show_inf;
2512  heavy_penalization = copy->heavy_penalization;
2513  penalize = copy->penalize;
2514  evaluation_reduction = copy->evaluation_reduction;
2515  modelSimplification = copy->modelSimplification;
2516  defocus_range = copy->defocus_range;
2517  downsampleFactor = copy->downsampleFactor;
2518 
2519  enhanced_ctftomodel() = copy->enhanced_ctftomodel();
2520  enhanced_ctftomodel_fullsize() = copy->enhanced_ctftomodel_fullsize();
2521  enhanced_weight = copy->enhanced_weight;
2522 
2523  current_ctfmodel = copy->current_ctfmodel;
2524  initial_ctfmodel = copy->initial_ctfmodel;
2525  ctfmodel_defoci = copy->ctfmodel_defoci;
2526  current_ctfmodel.precomputeValues(x_contfreq,y_contfreq);
2527 
2528  Tm = copy->Tm;
2529  f = copy->f;
2530  ctfmodelSize = copy->ctfmodelSize;
2531  show_optimization = copy->show_optimization;
2532  fn_psd = copy->fn_psd;
2533 
2534 }
CTFDescription::sigmaU
double sigmaU
Spherical aberration (in milimeters). Typical value 5.6.
Definition: ctf.h:871
VECTOR_R2
#define VECTOR_R2(v, x, y)
Definition: matrix1d.h:112
ProgCTFBasicParams::defocus_range
double defocus_range
Defocus range.
Definition: ctf_estimate_from_psd_base.h:54
ENVELOPE_PARAMETERS
constexpr int ENVELOPE_PARAMETERS
Definition: ctf_estimate_from_psd.cpp:46
ProgCTFEstimateFromPSD::generate_model_quadrant
void generate_model_quadrant(int Ydim, int Xdim, MultidimArray< double > &model)
Definition: ctf_estimate_from_psd.cpp:472
numerical_tools.h
MDL_CTF_XF
The CTF is valid within (x0,y0) to (xF,yF) in the micrograph coordinates.
Definition: metadata_label.h:143
ctf_estimate_from_psd.h
YSIZE
#define YSIZE(v)
Definition: multidim_array_base.h:95
min
void min(Image< double > &op1, const Image< double > &op2)
Definition: image_operate.cpp:118
ProgCTFEstimateFromPSD::assignParametersFromCTF
void assignParametersFromCTF(const CTFDescription &ctfmodel, double *p, int ia, int l, int modelSimplification)
Definition: ctf_estimate_from_psd.cpp:176
A2D_ELEM
#define A2D_ELEM(v, i, j)
Definition: multidim_array_base.h:354
VEC_ELEM
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
ProgCTFBasicParams::enhanced_weight
double enhanced_weight
Weight of the enhanced image.
Definition: ctf_estimate_from_psd_base.h:60
MultidimArray::resize
void resize(size_t Ndim, size_t Zdim, size_t Ydim, size_t Xdim, bool copy=true)
Definition: multidim_array.cpp:983
XMIPP_MAX
#define XMIPP_MAX(x, y)
Definition: xmipp_macros.h:193
COPY_ctfmodel_TO_CURRENT_GUESS
#define COPY_ctfmodel_TO_CURRENT_GUESS
Definition: ctf_estimate_from_psd.cpp:264
ProgCTFBasicParams::f
MultidimArray< double > * f
Definition: ctf_estimate_from_psd_base.h:125
Matrix1D::module
double module() const
Definition: matrix1d.h:983
ProgCTFBasicParams::ctftomodel
Image< double > ctftomodel
CTF amplitude to model.
Definition: ctf_estimate_from_psd_base.h:25
floor
__host__ __device__ float2 floor(const float2 v)
Definition: cuda_basic_math.h:175
ProgCTFBasicParams::fastDefocusEstimate
bool fastDefocusEstimate
Fast defocus estimate.
Definition: ctf_estimate_from_psd_base.h:40
FileName::removeLastExtension
FileName removeLastExtension() const
Definition: xmipp_filename.cpp:374
MetaDataVec::read
void read(const FileName &inFile, const std::vector< MDLabel > *desiredLabels=nullptr, bool decomposeStack=true) override
Definition: metadata_vec.cpp:120
ProgCTFEstimateFromPSD::saveIntermediateResults
void saveIntermediateResults(const FileName &fn_root, bool generate_profiles)
Definition: ctf_estimate_from_psd.cpp:369
FINISHINGX
#define FINISHINGX(v)
Definition: multidim_array_base.h:55
REPORT_ERROR
#define REPORT_ERROR(nerr, ErrormMsg)
Definition: xmipp_error.h:211
c
doublereal * c
Definition: numerical_recipes.cpp:2230
MultidimArray::resizeNoCopy
void resizeNoCopy(const MultidimArray< T1 > &v)
Definition: multidim_array.h:504
metadata_vec.h
CTFDescription1D::precomputeValues
void precomputeValues(double X)
Precompute values for a given frequency.
Definition: ctf.h:376
MD_APPEND
Definition: metadata_writemode.h:34
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
FIRST_SQRT_PARAMETER
constexpr int FIRST_SQRT_PARAMETER
Definition: ctf_estimate_from_psd.cpp:48
y
static double * y
Definition: numerical_recipes.cpp:8487
ERR_IO_NOWRITE
Couldn&#39;t write to file.
Definition: xmipp_error.h:140
ProgCTFBasicParams::w_digfreq
MultidimArray< double > w_digfreq
Definition: ctf_estimate_from_psd_base.h:71
ProgCTFEstimateFromPSD::assignCTFfromParameters
void assignCTFfromParameters(double *p, CTFDescription &ctfmodel, int ia, int l, int modelSimplification)
Definition: ctf_estimate_from_psd.cpp:81
ProgCTFBasicParams::adjust_params
Matrix1D< double > * adjust_params
Definition: ctf_estimate_from_psd_base.h:126
ProgCTFBasicParams::readBasicParams
void readBasicParams(XmippProgram *program)
Read parameters.
Definition: ctf_estimate_from_psd_base.cpp:33
CTFDescription::DeltafV
double DeltafV
Defocus in V (in Angstroms). Negative values are underfocused.
Definition: ctf.h:830
ASSIGN_PARAM_CTF
#define ASSIGN_PARAM_CTF(index, paramName)
Definition: ctf_estimate_from_psd.cpp:174
DIRECT_A2D_ELEM
#define DIRECT_A2D_ELEM(v, i, j)
Definition: multidim_array_base.h:341
ImageBase::write
void write(const FileName &name="", size_t select_img=ALL_IMAGES, bool isStack=false, int mode=WRITE_OVERWRITE, CastWriteMode castMode=CW_CAST, int _swapWrite=0)
Definition: xmipp_image_base.cpp:292
CTFDescription::write
void write(const FileName &fn)
Definition: ctf.cpp:1275
ProgCTFBasicParams::y_digfreq
MultidimArray< double > y_digfreq
Definition: ctf_estimate_from_psd_base.h:70
A1D_ELEM
#define A1D_ELEM(v, i)
Definition: multidim_array_base.h:539
ProgCTFEstimateFromPSD::readBasicParams
void readBasicParams(XmippProgram *program)
Read parameters.
Definition: ctf_estimate_from_psd.cpp:270
ProgCTFBasicParams::Tm
double Tm
Sampling rate.
Definition: ctf_estimate_from_psd_base.h:52
CTFDescription1D::espr
double espr
Definition: ctf.h:251
ProgCTFBasicParams::w_digfreq_r_iN
MultidimArray< double > w_digfreq_r_iN
Definition: ctf_estimate_from_psd_base.h:80
MultidimArray::rangeAdjust
void rangeAdjust(T minF, T maxF)
Definition: multidim_array.h:1982
ProgCTFBasicParams::downsampleFactor
double downsampleFactor
Downsample performed.
Definition: ctf_estimate_from_psd_base.h:23
ProgCTFBasicParams::show
void show()
Show parameters.
Definition: ctf_estimate_from_psd_base.cpp:70
CTFDescription::cU
double cU
Gaussian center for U.
Definition: ctf.h:875
Matrix2D::inv
void inv(Matrix2D< T > &result) const
Definition: matrix2d.cpp:663
w
doublereal * w
Definition: numerical_recipes.cpp:2477
evaluateIceness
double evaluateIceness(const MultidimArray< double > &enhanced_ctftomodel, double Tm)
MultidimArray::initConstant
void initConstant(T val)
Definition: multidim_array.h:2705
ProgCTFEstimateFromPSD::produceSideInfo
void produceSideInfo()
Produce side information.
Definition: ctf_estimate_from_psd.cpp:303
MetaDataVec::write
void write(const FileName &outFile, WriteModeMetaData mode=MD_OVERWRITE) const
Definition: metadata_vec.cpp:140
CTFDescription1D::phase_shift
double phase_shift
Definition: ctf.h:305
ProgCTFBasicParams::psd_exp_enhanced_radial_derivative
MultidimArray< double > psd_exp_enhanced_radial_derivative
Definition: ctf_estimate_from_psd_base.h:76
powellOptimizer
void powellOptimizer(Matrix1D< double > &p, int i0, int n, double(*f)(double *x, void *), void *prm, double ftol, double &fret, int &iter, const Matrix1D< double > &steps, bool show)
Definition: numerical_tools.cpp:44
demodulate
void demodulate(MultidimArray< double > &im, double lambda, int size, int x, int y, int rmin, int rmax, MultidimArray< double > &phase, MultidimArray< double > &mod, Matrix1D< double > &coeffs, int verbose)
Definition: fringe_processing.cpp:671
DEFOCUS_PARAMETERS
constexpr int DEFOCUS_PARAMETERS
Definition: ctf_estimate_from_psd.cpp:47
iter
glob_prnt iter
Definition: numerical_recipes.cpp:4698
DEBUG_CLOSE_TEXTFILE
#define DEBUG_CLOSE_TEXTFILE
Definition: ctf_estimate_from_psd.cpp:61
SIND
#define SIND(x)
Definition: xmipp_macros.h:347
MDL_CTF_YF
The CTF is valid within (x0,y0) to (xF,yF) in the micrograph coordinates.
Definition: metadata_label.h:144
FFT_idx2digfreq
void FFT_idx2digfreq(T &v, const Matrix1D< int > &idx, Matrix1D< double > &freq)
Definition: xmipp_fft.h:80
ProgCTFBasicParams::fn_psd
FileName fn_psd
CTF filename.
Definition: ctf_estimate_from_psd_base.h:21
STARTINGX
#define STARTINGX(v)
Definition: multidim_array_base.h:51
ProgCTFBasicParams::w_digfreq_r
MultidimArray< double > w_digfreq_r
Definition: ctf_estimate_from_psd_base.h:81
CTFDescription1D::DeltaF
double DeltaF
Longitudinal mechanical displacement (ansgtrom). Typical value 100.
Definition: ctf.h:257
x
doublereal * x
Definition: numerical_recipes.cpp:2230
i
#define i
Definition: numerical_recipes.cpp:2493
selfScaleToSize
void selfScaleToSize(int SplineDegree, MultidimArrayGeneric &V1, int Xdim, int Ydim, int Zdim)
Definition: transformations.cpp:736
CTFDescription::defineParams
static void defineParams(XmippProgram *program)
Define parameters in the command line.
Definition: ctf.cpp:1287
ProgCTFEstimateFromPSD::CTF_fitness_object
double CTF_fitness_object(double *p)
Definition: ctf_estimate_from_psd.cpp:601
FOR_ALL_ELEMENTS_IN_ARRAY2D
#define FOR_ALL_ELEMENTS_IN_ARRAY2D(m)
Definition: multidim_array_base.h:390
Matrix1D::maxIndex
int maxIndex() const
Definition: matrix1d.cpp:610
MDL_CTF_CRIT_FITTINGCORR13
Correlation between the 1st and 3rd ring of the CTF.
Definition: metadata_label.h:198
ProgCTFBasicParams::max_freq
double max_freq
Maximum frequency to adjust.
Definition: ctf_estimate_from_psd_base.h:50
ProgCTFBasicParams::penalty
static constexpr double penalty
Definition: ctf_estimate_from_psd_base.h:123
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
CenterFFT
void CenterFFT(MultidimArray< T > &v, bool forward)
Definition: xmipp_fft.h:291
ProgCTFBasicParams::enhanced_ctftomodel_fullsize
Image< double > enhanced_ctftomodel_fullsize
CTF amplitude to model.
Definition: ctf_estimate_from_psd_base.h:29
CTFDescription::azimuthal_angle
double azimuthal_angle
Azimuthal angle (between X and U) in degrees.
Definition: ctf.h:832
ProgCTFBasicParams::ctfmodelSize
int ctfmodelSize
X dimension of particle projections (-1=the same as the psd)
Definition: ctf_estimate_from_psd_base.h:34
ProgCTFBasicParams::w_contfreq
MultidimArray< double > w_contfreq
Definition: ctf_estimate_from_psd_base.h:68
first
glob_log first
Definition: numerical_recipes.cpp:5201
ProgCTFBasicParams::psd_theo_radial
MultidimArray< double > psd_theo_radial
Definition: ctf_estimate_from_psd_base.h:79
ALL_CTF_PARAMETERS
constexpr int ALL_CTF_PARAMETERS
Definition: ctf_estimate_from_psd.cpp:41
MDL_CTF_CRIT_FITTINGSCORE
Score of the fitting.
Definition: metadata_label.h:197
b
doublereal * b
Definition: numerical_recipes.cpp:2230
ProgCTFBasicParams::psd_exp_radial
MultidimArray< double > psd_exp_radial
PSD data.
Definition: ctf_estimate_from_psd_base.h:73
ProgCTFEstimateFromPSD::center_optimization_focus
void center_optimization_focus(bool adjust_freq, bool adjust_th, double margin)
Definition: ctf_estimate_from_psd.cpp:1027
ProgCTFEstimateFromPSD::generateModelSoFar
void generateModelSoFar(Image< double > &I, bool apply_log)
Definition: ctf_estimate_from_psd.cpp:322
FLOOR
#define FLOOR(x)
Definition: xmipp_macros.h:240
log
void log(Image< double > &op)
Definition: image_operate.cpp:228
FIRST_ENVELOPE_PARAMETER
constexpr int FIRST_ENVELOPE_PARAMETER
Definition: ctf_estimate_from_psd.cpp:49
lambda
double * lambda
Definition: numerical_recipes.cpp:4414
XX
#define XX(v)
Definition: matrix1d.h:85
ProgCTFEstimateFromPSD::compute_central_region
void compute_central_region(double &w1, double &w2, double ang)
Definition: ctf_estimate_from_psd.cpp:986
contfreq2digfreq
void contfreq2digfreq(const Matrix1D< double > &contfreq, Matrix1D< double > &digfreq, double pixel_size)
Definition: xmipp_fft.h:139
MetaDataVec::setValue
bool setValue(const MDObject &mdValueIn, size_t id)
Definition: metadata_vec.cpp:282
CEIL
#define CEIL(x)
Definition: xmipp_macros.h:225
f
double * f
Definition: numerical_recipes.cpp:4414
ERR_ARG_INCORRECT
Incorrect argument received.
Definition: xmipp_error.h:113
MetaDataVec::firstRowId
size_t firstRowId() const override
Definition: metadata_vec.cpp:508
ProgCTFBasicParams::noDefocusEstimate
bool noDefocusEstimate
No defocus estimate.
Definition: ctf_estimate_from_psd_base.h:42
XMIPP_EQUAL_ACCURACY
#define XMIPP_EQUAL_ACCURACY
Definition: xmipp_macros.h:119
ctf_estimate_from_psd_fast.h
Matrix1D::resize
void resize(size_t Xdim, bool copy=true)
Definition: matrix1d.h:410
CTF_PARAMETERS
constexpr int CTF_PARAMETERS
Definition: ctf_estimate_from_psd.cpp:42
XSIZE
#define XSIZE(v)
Definition: multidim_array_base.h:91
Matrix2D::det
T det() const
Definition: matrix2d.cpp:38
CTFDescription1D::ispr
double ispr
Objective lens stability (deltaI/I) (ppm). Typical value 1.
Definition: ctf.h:253
max
void max(Image< double > &op1, const Image< double > &op2)
Definition: image_operate.cpp:128
CTFDescription1D::Tm
double Tm
Sampling rate (A/pixel)
Definition: ctf.h:240
ProgCTFBasicParams::refineAmplitudeContrast
bool refineAmplitudeContrast
Refine amplitude contrast.
Definition: ctf_estimate_from_psd_base.h:38
ERR_NUMERICAL
Error related to numerical calculation.
Definition: xmipp_error.h:179
FIRST_DEFOCUS_PARAMETER
constexpr int FIRST_DEFOCUS_PARAMETER
Definition: ctf_estimate_from_psd.cpp:50
log10
void log10(Image< double > &op)
Definition: image_operate.cpp:237
PARAMETRIC_CTF_PARAMETERS
constexpr int PARAMETRIC_CTF_PARAMETERS
Definition: ctf_estimate_from_psd.cpp:43
CTFDescription1D::Cs
double Cs
Spherical aberration (in milimeters). Typical value 5.6.
Definition: ctf.h:246
Matrix1D::initZeros
void initZeros()
Definition: matrix1d.h:592
mod
void mod(const MultidimArray< T > &x, MultidimArray< T > &m, double y)
Definition: multidim_array.h:3902
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
ProgCTFBasicParams::show_optimization
bool show_optimization
Show convergence values.
Definition: ctf_estimate_from_psd_base.h:31
ProgCTFBasicParams::psd_theo_radial_derivative
MultidimArray< double > psd_theo_radial_derivative
Definition: ctf_estimate_from_psd_base.h:77
XMIPP_MIN
#define XMIPP_MIN(x, y)
Definition: xmipp_macros.h:181
ProgCTFEstimateFromPSD::generate_model_halfplane
void generate_model_halfplane(int Ydim, int Xdim, MultidimArray< double > &model)
Definition: ctf_estimate_from_psd.cpp:536
j
#define j
Definition: numerical_recipes.cpp:2493
ProgCTFBasicParams::defineParams
void defineParams()
Define Parameters.
Definition: ctf_estimate_from_psd_base.cpp:171
steps
double steps
Definition: numerical_recipes.cpp:5166
YY
#define YY(v)
Definition: matrix1d.h:93
MDL_CTF_DOWNSAMPLE_PERFORMED
Downsampling performed to estimate the CTF.
Definition: metadata_label.h:204
ProgCTFBasicParams::x_digfreq
MultidimArray< double > x_digfreq
Definition: ctf_estimate_from_psd_base.h:69
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
error
void error(char *s)
Definition: tools.cpp:107
FINISHINGY
#define FINISHINGY(v)
Definition: multidim_array_base.h:63
FOR_ALL_ELEMENTS_IN_ARRAY1D
#define FOR_ALL_ELEMENTS_IN_ARRAY1D(v)
Definition: multidim_array_base.h:554
ProgCTFEstimateFromPSD::estimate_defoci_Zernike
void estimate_defoci_Zernike(const MultidimArray< double > &psdToModelFullSize, double min_freq, double max_freq, double Tm, double kV, double lambdaPhase, int sizeWindowPhase, double &defocusU, double &defocusV, double &ellipseAngle, int verbose)
Definition: ctf_estimate_from_psd.cpp:1937
ProgCTFBasicParams::adjust
Matrix1D< double > adjust
Set of parameters for the complete adjustment of the CTF.
Definition: ctf_estimate_from_psd_base.h:62
ProgCTFBasicParams::produceSideInfo
void produceSideInfo()
Produce side information.
Definition: ctf_estimate_from_psd_base.cpp:186
MATRIX1D_ARRAY
#define MATRIX1D_ARRAY(v)
Definition: matrix1d.h:58
ProgCTFBasicParams::modelSimplification
int modelSimplification
Model simplification.
Definition: ctf_estimate_from_psd_base.h:64
CTFDescription1D::VPP_radius
double VPP_radius
Definition: ctf.h:306
FileName::withoutExtension
FileName withoutExtension() const
Definition: xmipp_filename.cpp:345
RAD2DEG
#define RAD2DEG(r)
Definition: xmipp_macros.h:320
BACKGROUND_CTF_PARAMETERS
constexpr int BACKGROUND_CTF_PARAMETERS
Definition: ctf_estimate_from_psd.cpp:44
ProgCTFBasicParams::mask_between_zeroes
MultidimArray< double > mask_between_zeroes
Definition: ctf_estimate_from_psd_base.h:86
CTFDescription1D::Ca
double Ca
Chromatic aberration (in milimeters). Typical value 2.
Definition: ctf.h:248
MDL_CTF_CRIT_ICENESS
Iceness of the micrograph.
Definition: metadata_label.h:203
CTFDescription::lookFor
void lookFor(int n, const Matrix1D< double > &u, Matrix1D< double > &freq, int iwhat=0)
Definition: ctf.cpp:1428
ProgCTFBasicParams::min_freq
double min_freq
Minimum frequency to adjust.
Definition: ctf_estimate_from_psd_base.h:48
String
std::string String
Definition: xmipp_strings.h:34
CTFDescription::forcePhysicalMeaning
void forcePhysicalMeaning()
Definition: ctf.cpp:1781
ASSIGN_CTF_PARAM
#define ASSIGN_CTF_PARAM(index, paramName)
Definition: ctf_estimate_from_psd.cpp:78
ProgCTFBasicParams::enhanced_ctftomodel
Image< double > enhanced_ctftomodel
CTF amplitude to model.
Definition: ctf_estimate_from_psd_base.h:27
CTFDescription1D::alpha
double alpha
Convergence cone semiangle (in mrad). Typical value 0.5.
Definition: ctf.h:255
fringe_processing.h
formatString
String formatString(const char *format,...)
Definition: xmipp_strings.cpp:602
Matrix2D::initZeros
void initZeros()
Definition: matrix2d.h:626
textToInteger
int textToInteger(const char *str)
Definition: xmipp_strings.cpp:137
SQRT_CTF_PARAMETERS
constexpr int SQRT_CTF_PARAMETERS
Definition: ctf_estimate_from_psd.cpp:45
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
u
doublereal * u
Definition: numerical_recipes.cpp:2230
ProgCTFBasicParams::psd_exp_radial_derivative
MultidimArray< double > psd_exp_radial_derivative
Definition: ctf_estimate_from_psd_base.h:78
ProgCTFBasicParams::x_contfreq
MultidimArray< double > x_contfreq
Frequencies in axes.
Definition: ctf_estimate_from_psd_base.h:66
ProgCTFBasicParams::psd_exp_enhanced_radial
MultidimArray< double > psd_exp_enhanced_radial
Definition: ctf_estimate_from_psd_base.h:74
MDL_CTF_X0
The CTF is valid within (x0,y0) to (xF,yF) in the micrograph coordinates.
Definition: metadata_label.h:141
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
K
constexpr int K
Definition: integration.cpp:123
prm
ProgClassifyCL2D * prm
Definition: mpi_classify_CL2D.cpp:35
ProgCTFEstimateFromPSD::defineParams
void defineParams()
Define Parameters.
Definition: ctf_estimate_from_psd.cpp:296
CTFDescription1D::kV
double kV
Accelerating Voltage (in KiloVolts)
Definition: ctf.h:242
MultidimArray::initZeros
void initZeros(const MultidimArray< T1 > &op)
Definition: multidim_array.h:2723
ProgCTFBasicParams::mask
MultidimArray< double > mask
Masks.
Definition: ctf_estimate_from_psd_base.h:85
COSD
#define COSD(x)
Definition: xmipp_macros.h:329
DEBUG_OPEN_TEXTFILE
#define DEBUG_OPEN_TEXTFILE(fnRoot)
Definition: ctf_estimate_from_psd.cpp:60
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
PI
#define PI
Definition: tools.h:43
Matrix1D::adaptForNumericalRecipes
T * adaptForNumericalRecipes() const
Definition: matrix1d.h:844
ProgCTFEstimateFromPSD::initial_ctfmodel
CTFDescription initial_ctfmodel
CTF model.
Definition: ctf_estimate_from_psd.h:44
ProgCTFEstimateFromPSD::defineBasicParams
static void defineBasicParams(XmippProgram *program)
Define basic parameters.
Definition: ctf_estimate_from_psd.cpp:289
AdjustCTF::global_prm
ProgCTFEstimateFromPSD * global_prm
Definition: ctf_estimate_from_psd.cpp:70
ROUT_Adjust_CTF
double ROUT_Adjust_CTF(ProgCTFEstimateFromPSD &prm, CTFDescription &output_ctfmodel, bool standalone)
Definition: ctf_estimate_from_psd.cpp:2208
Matrix1D::initConstant
void initConstant(T val)
Definition: matrix1d.cpp:83
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
fitness
double fitness(double *p)
Definition: volume_align_prog.cpp:88
DEBUG_TEXTFILE
#define DEBUG_TEXTFILE(str)
Definition: ctf_estimate_from_psd.cpp:62
DEBUG_MODEL_TEXTFILE
#define DEBUG_MODEL_TEXTFILE
Definition: ctf_estimate_from_psd.cpp:63
ProgCTFBasicParams::y_contfreq
MultidimArray< double > y_contfreq
Definition: ctf_estimate_from_psd_base.h:67
CTF_fitness
double CTF_fitness(double *, void *)
Definition: ctf_estimate_from_psd.cpp:979
MultidimArray::computeMax
T computeMax() const
Definition: multidim_array.h:1537
digfreq2contfreq
void digfreq2contfreq(const Matrix1D< double > &digfreq, Matrix1D< double > &contfreq, double pixel_size)
Definition: xmipp_fft.h:125
xmipp_transformation::BSPLINE3
Definition: transformations_defines.h:32
MDL_CTF_Y0
The CTF is valid within (x0,y0) to (xF,yF) in the micrograph coordinates.
Definition: metadata_label.h:142
CTFDescription::cV
double cV
Gaussian center for V.
Definition: ctf.h:877
Generated by   doxygen 1.8.13