#include <ctf_estimate_from_psd_fast.h>

Inheritance diagram for ProgCTFEstimateFromPSDFast:

Collaboration diagram for ProgCTFEstimateFromPSDFast:

Public Member Functions
	ProgCTFEstimateFromPSDFast ()

void	readBasicParams (XmippProgram *program)
	Read parameters. More...

void	readParams ()
	Read parameters. More...

void	defineParams ()
	Define Parameters. More...

void	produceSideInfo ()
	Produce side information. More...

void	run ()

void	assignCTFfromParameters_fast (double *p, CTFDescription1D &ctf1Dmodel, int ia, int l)

void	assignParametersFromCTF_fast (const CTFDescription1D &ctfmodel, double *p, int ia, int l)

void	center_optimization_focus_fast (bool adjust_th, double margin)

void	generateModelSoFar_fast (MultidimArray< double > &I, bool apply_log)

void	compute_central_region_fast (double &w1, double &w2, double ang)

void	saveIntermediateResults_fast (const FileName &fn_root, bool generate_profiles)

double	CTF_fitness_object_fast (double *p)

void	estimate_background_sqrt_parameters_fast ()

void	estimate_background_gauss_parameters_fast ()

void	estimate_background_gauss_parameters2_fast ()

void	estimate_envelope_parameters_fast ()

void	showFirstDefoci_fast ()

void	estimate_defoci_fast ()

Public Member Functions inherited from ProgCTFBasicParams
	ProgCTFBasicParams ()

void	readBasicParams (XmippProgram *program)
	Read parameters. More...

void	show ()
	Show parameters. More...

void	produceSideInfo ()
	Produce side information. More...

void	generate_model_halfplane (int Ydim, int Xdim, MultidimArray< double > &model)

void	generate_model_quadrant (int Ydim, int Xdim, MultidimArray< double > &model)

Public Member Functions inherited from XmippProgram
const char *	getParam (const char *param, int arg=0)

const char *	getParam (const char param, const char subparam, int arg=0)

int	getIntParam (const char *param, int arg=0)

int	getIntParam (const char param, const char subparam, int arg=0)

double	getDoubleParam (const char *param, int arg=0)

double	getDoubleParam (const char param, const char subparam, int arg=0)

float	getFloatParam (const char *param, int arg=0)

float	getFloatParam (const char param, const char subparam, int arg=0)

void	getListParam (const char *param, StringVector &list)

int	getCountParam (const char *param)

bool	checkParam (const char *param)

bool	existsParam (const char *param)

void	addParamsLine (const String &line)

void	addParamsLine (const char *line)

ParamDef *	getParamDef (const char *param) const

virtual void	quit (int exit_code=0) const

virtual int	tryRun ()

void	initProgress (size_t total, size_t stepBin=60)

void	setProgress (size_t value=0)

void	endProgress ()

void	processDefaultComment (const char param, const char left)

void	setDefaultComment (const char param, const char comment)

virtual void	initComments ()

void	setProgramName (const char *name)

void	addUsageLine (const char *line, bool verbatim=false)

void	clearUsage ()

void	addExampleLine (const char *example, bool verbatim=true)

void	addSeeAlsoLine (const char *seeAlso)

void	addKeywords (const char *keywords)

const char *	name () const

virtual void	usage (int verb=0) const

virtual void	usage (const String &param, int verb=2)

int	version () const

virtual void	show () const

virtual void	read (int argc, const char **argv, bool reportErrors=true)

virtual void	read (int argc, char **argv, bool reportErrors=true)

void	read (const String &argumentsLine)

	XmippProgram ()

	XmippProgram (int argc, const char **argv)

virtual	~XmippProgram ()

Static Public Member Functions
static void	defineBasicParams (XmippProgram *program)
	Define basic parameters. More...

Static Public Member Functions inherited from ProgCTFBasicParams
static void	defineBasicParams (XmippProgram *program)
	Define basic parameters. More...

Public Attributes
CTFDescription1D	initial_ctfmodel

CTFDescription1D	current_ctfmodel

CTFDescription1D	ctfmodel_defoci

CTFDescription	initial_ctfmodel2D

Public Attributes inherited from ProgCTFBasicParams
FileName	fn_psd
	CTF filename. More...

double	downsampleFactor
	Downsample performed. More...

Image< double >	ctftomodel
	CTF amplitude to model. More...

Image< double >	enhanced_ctftomodel
	CTF amplitude to model. More...

Image< double >	enhanced_ctftomodel_fullsize
	CTF amplitude to model. More...

bool	show_optimization
	Show convergence values. More...

bool	selfEstimation

int	ctfmodelSize
	X dimension of particle projections (-1=the same as the psd) More...

bool	bootstrap
	Bootstrap estimation. More...

bool	refineAmplitudeContrast
	Refine amplitude contrast. More...

bool	fastDefocusEstimate
	Fast defocus estimate. More...

bool	noDefocusEstimate
	No defocus estimate. More...

double	lambdaPhase
	Regularization factor for the phase direction and unwrapping estimates (used in Zernike estimate) More...

int	sizeWindowPhase
	Size of the average window used during phase direction and unwrapping estimates (used in Zernike estimate) More...

double	min_freq
	Minimum frequency to adjust. More...

double	max_freq
	Maximum frequency to adjust. More...

double	Tm
	Sampling rate. More...

double	defocus_range
	Defocus range. More...

double	f1
	Enhancement filter low freq. More...

double	f2
	Enhancement filter high freq. More...

double	enhanced_weight
	Weight of the enhanced image. More...

Matrix1D< double >	adjust
	Set of parameters for the complete adjustment of the CTF. More...

int	modelSimplification
	Model simplification. More...

MultidimArray< double >	x_contfreq
	Frequencies in axes. More...

MultidimArray< double >	y_contfreq

MultidimArray< double >	w_contfreq

MultidimArray< double >	x_digfreq

MultidimArray< double >	y_digfreq

MultidimArray< double >	w_digfreq

MultidimArray< double >	psd_exp_radial
	PSD data. More...

MultidimArray< double >	psd_exp_enhanced_radial

MultidimArray< double >	psd_exp_enhanced_radial_2

MultidimArray< double >	psd_exp_enhanced_radial_derivative

MultidimArray< double >	psd_theo_radial_derivative

MultidimArray< double >	psd_exp_radial_derivative

MultidimArray< double >	psd_theo_radial

MultidimArray< double >	w_digfreq_r_iN

MultidimArray< double >	w_digfreq_r

MultidimArray< std::complex< double > >	psd_fft

std::vector< double >	amplitud

MultidimArray< double >	mask
	Masks. More...

MultidimArray< double >	mask_between_zeroes

MultidimArray< double >	w_count

double	value_th

double	min_freq_psd

double	max_freq_psd

double	max_gauss_freq

int	show_inf

int	action

double	corr13

bool	penalize

int	evaluation_reduction

double	heavy_penalization

double	current_penalty

MultidimArray< double > *	f

Matrix1D< double > *	adjust_params

Public Attributes inherited from XmippProgram
bool	doRun

bool	runWithoutArgs

int	verbose
	Verbosity level. More...

int	debug

Additional Inherited Members
Static Public Attributes inherited from ProgCTFBasicParams
static constexpr double	penalty = 32.0

Protected Member Functions inherited from XmippProgram
void	defineCommons ()

Protected Attributes inherited from XmippProgram
int	errorCode

ProgramDef *	progDef
	Program definition and arguments parser. More...

std::map< String, CommentList >	defaultComments

int	argc
	Original command line arguments. More...

const char **	argv

Detailed Description

Adjust CTF parameters.

Definition at line 37 of file ctf_estimate_from_psd_fast.h.

Constructor & Destructor Documentation

◆ ProgCTFEstimateFromPSDFast()

ProgCTFEstimateFromPSDFast::ProgCTFEstimateFromPSDFast ( )

inline

Definition at line 44 of file ctf_estimate_from_psd_fast.h.

44 {};

Member Function Documentation

◆ assignCTFfromParameters_fast()

void ProgCTFEstimateFromPSDFast::assignCTFfromParameters_fast	(	double *	p,
		CTFDescription1D &	ctf1Dmodel,
		int	ia,
		int	l
	)

Definition at line 76 of file ctf_estimate_from_psd_fast.cpp.

 {
     ctf1Dmodel.Tm = Tm;
 
     ASSIGN_CTF_PARAM(0, Defocus);
     ASSIGN_CTF_PARAM(1, kV);
     ASSIGN_CTF_PARAM(2, K);
     ASSIGN_CTF_PARAM(3, Cs);
     ASSIGN_CTF_PARAM(4, Ca);
     ASSIGN_CTF_PARAM(5, espr);
     ASSIGN_CTF_PARAM(6, ispr); //deltaI/I
     ASSIGN_CTF_PARAM(7, alpha);
     ASSIGN_CTF_PARAM(8, DeltaF);
     ASSIGN_CTF_PARAM(9, DeltaR);
     ASSIGN_CTF_PARAM(10, envR0);
     ASSIGN_CTF_PARAM(11, envR1);
     ASSIGN_CTF_PARAM(12, envR2);
     ASSIGN_CTF_PARAM(13, Q0);
     ASSIGN_CTF_PARAM(14, base_line);
     ASSIGN_CTF_PARAM(15, sqrt_K);
     ASSIGN_CTF_PARAM(16, sq);
     ASSIGN_CTF_PARAM(17, bgR1);
     ASSIGN_CTF_PARAM(18, bgR2);
     ASSIGN_CTF_PARAM(19, bgR3);
     ASSIGN_CTF_PARAM(20, gaussian_K);
     ASSIGN_CTF_PARAM(21, sigma1);
     ASSIGN_CTF_PARAM(22, Gc1);
     ASSIGN_CTF_PARAM(23, gaussian_K2);
     ASSIGN_CTF_PARAM(24, sigma2);
     ASSIGN_CTF_PARAM(25, Gc2);
     ASSIGN_CTF_PARAM(26, phase_shift);
     ASSIGN_CTF_PARAM(27, VPP_radius);
 
 }//function assignCTFfromParameters

◆ assignParametersFromCTF_fast()

void ProgCTFEstimateFromPSDFast::assignParametersFromCTF_fast	(	const CTFDescription1D &	ctfmodel,
		double *	p,
		int	ia,
		int	l
	)

Definition at line 114 of file ctf_estimate_from_psd_fast.cpp.

 {
 
         ASSIGN_PARAM_CTF(0, Defocus);
         ASSIGN_PARAM_CTF(1, kV);
         ASSIGN_PARAM_CTF(2, K);
         ASSIGN_PARAM_CTF(3, Cs);
         ASSIGN_PARAM_CTF(4, Ca);
         ASSIGN_PARAM_CTF(5, espr);
         ASSIGN_PARAM_CTF(6, ispr); //deltaI/I
         ASSIGN_PARAM_CTF(7, alpha);
         ASSIGN_PARAM_CTF(8, DeltaF);
         ASSIGN_PARAM_CTF(9, DeltaR);
         ASSIGN_PARAM_CTF(10, envR0);
         ASSIGN_PARAM_CTF(11, envR1);
         ASSIGN_PARAM_CTF(12, envR2);
         ASSIGN_PARAM_CTF(13, Q0);
         ASSIGN_PARAM_CTF(14, base_line);
         ASSIGN_PARAM_CTF(15, sqrt_K);
         ASSIGN_PARAM_CTF(16, sq);
         ASSIGN_PARAM_CTF(17, bgR1);
         ASSIGN_PARAM_CTF(18, bgR2);
         ASSIGN_PARAM_CTF(19, bgR3);
         ASSIGN_PARAM_CTF(20, gaussian_K);
         ASSIGN_PARAM_CTF(21, sigma1);
         ASSIGN_PARAM_CTF(22, Gc1);
         ASSIGN_PARAM_CTF(23, gaussian_K2);
         ASSIGN_PARAM_CTF(24, sigma2);
         ASSIGN_PARAM_CTF(25, Gc2);
         ASSIGN_PARAM_CTF(26, phase_shift);
         ASSIGN_PARAM_CTF(27, VPP_radius);
 
 }

◆ center_optimization_focus_fast()

void ProgCTFEstimateFromPSDFast::center_optimization_focus_fast	(	bool	adjust_th,
		double	margin = `1`
	)

Definition at line 637 of file ctf_estimate_from_psd_fast.cpp.

 {
     if (show_optimization)
         std::cout << "Freq frame before focusing=" << min_freq_psd << ","
         << max_freq_psd << std::endl << "Value_th before focusing="
         << value_th << std::endl;
 
     double w1 = min_freq_psd, w2 = max_freq_psd;
 
     // Compute maximum value within central region
     if (adjust_th)
     {
         MultidimArray<double> save;
         generateModelSoFar_fast(save, false);
         double max_val = 0;
         FOR_ALL_ELEMENTS_IN_ARRAY1D(w_digfreq)
         {
             double w = w_digfreq(i);
             if (w >= w1 && w <= w2)
                 max_val = XMIPP_MAX(max_val, save(i));
         }
         if (value_th != -1)
             value_th = XMIPP_MIN(value_th, max_val * margin);
         else
             value_th = max_val * margin;
     }
 
 
     if (show_optimization)
         std::cout << "Freq frame after focusing=" << min_freq_psd << ","
         << max_freq_psd << std::endl << "Value_th after focusing="
         << value_th << std::endl;
 }

◆ compute_central_region_fast()

void ProgCTFEstimateFromPSDFast::compute_central_region_fast	(	double &	w1,
		double &	w2,
		double	ang
	)

◆ CTF_fitness_object_fast()

double ProgCTFEstimateFromPSDFast::CTF_fitness_object_fast ( double * p )

CTF fitness

Definition at line 288 of file ctf_estimate_from_psd_fast.cpp.

 {
     double retval;
 
     // Generate CTF model
     switch (action)
     {
         // Remind that p is a vector whose first element is at index 1
     case 0:
         assignCTFfromParameters_fast(p - FIRST_SQRT_PARAMETER + 1,
                                 current_ctfmodel, FIRST_SQRT_PARAMETER, SQRT_CTF_PARAMETERS);
 
         if (show_inf >= 2)
         {
             std::cout << "Input vector:";
             for (int i = 1; i <= SQRT_CTF_PARAMETERS; i++)
                 std::cout << p[i] << " ";
             std::cout << std::endl;
         }
         break;
     case 1:
             assignCTFfromParameters_fast(p - FIRST_SQRT_PARAMETER + 1,
                                 current_ctfmodel, FIRST_SQRT_PARAMETER,
                                     BACKGROUND_CTF_PARAMETERS);
         if (show_inf >= 2)
         {
             std::cout << "Input vector:";
             for (int i = 1; i <= BACKGROUND_CTF_PARAMETERS; i++)
                 std::cout << p[i] << " ";
             std::cout << std::endl;
         }
         break;
     case 2:
             assignCTFfromParameters_fast(p - FIRST_ENVELOPE_PARAMETER + 1,
                                 current_ctfmodel, FIRST_ENVELOPE_PARAMETER, ENVELOPE_PARAMETERS);
         if (show_inf >= 2)
         {
             std::cout << "Input vector:";
             for (int i = 1; i <= ENVELOPE_PARAMETERS; i++)
                 std::cout << p[i] << " ";
             std::cout << std::endl;
         }
         break;
     case 3:
             assignCTFfromParameters_fast(p - FIRST_DEFOCUS_PARAMETER + 1,
                                     current_ctfmodel, FIRST_DEFOCUS_PARAMETER, DEFOCUS_PARAMETERS);
         psd_theo_radial_derivative.initZeros();
         if (show_inf >= 2)
         {
             std::cout << "Input vector:";
             for (int i = 1; i <= DEFOCUS_PARAMETERS; i++)
                 std::cout << p[i] << " ";
             std::cout << std::endl;
         }
         break;
     case 4:
             assignCTFfromParameters_fast(p - 0 + 1, current_ctfmodel, 0,
                                     CTF_PARAMETERS);
         psd_theo_radial.initZeros();
         if (show_inf >= 2)
         {
             std::cout << "Input vector:";
             for (int i = 1; i <= CTF_PARAMETERS; i++)
                 std::cout << p[i] << " ";
             std::cout << std::endl;
         }
         break;
     case 5:
     case 6:
     case 7:
         assignCTFfromParameters_fast(p - 0 + 1, current_ctfmodel, 0,
                                 ALL_CTF_PARAMETERS);
         psd_theo_radial.initZeros();
         if (show_inf >= 2)
         {
             std::cout << "Input vector:";
             for (int i = 1; i <= ALL_CTF_PARAMETERS; i++)
                 std::cout << p[i] << " ";
             std::cout << std::endl;
         }
         break;
     }
 
     current_ctfmodel.produceSideInfo();
     if (show_inf >= 2)
         std::cout << "Model:\n" << current_ctfmodel << std::endl;
     if (!current_ctfmodel.hasPhysicalMeaning())
     {
         if (show_inf >= 2)
             std::cout << "Does not have physical meaning\n";
         return heavy_penalization;
     }
     if (action > 3
         && (fabs((current_ctfmodel.Defocus - ctfmodel_defoci.Defocus)
                 / ctfmodel_defoci.Defocus) > 0.2)  && !noDefocusEstimate)
     {
         if (show_inf >= 2)
             std::cout << "Too large defocus\n";
         return heavy_penalization;
     }
     if (initial_ctfmodel.Defocus != 0 && action >= 3 && !selfEstimation && !noDefocusEstimate)
     {
         // If there is an initial model, the true solution
         // cannot be too far
         if (fabs(initial_ctfmodel.Defocus - current_ctfmodel.Defocus) > defocus_range)
         {
             if (show_inf >= 2)
             {
                 std::cout << "Too far from hint: Initial (" << initial_ctfmodel.Defocus << ")"
                 << " current guess (" << current_ctfmodel.Defocus << ") max allowed difference: "
                 << defocus_range << std::endl;
             }
             return heavy_penalization;
 
 
         }
     }
     if ((initial_ctfmodel.phase_shift != 0.0 || initial_ctfmodel.phase_shift != 1.57079) && action >= 5)
     {
         if (fabs(initial_ctfmodel.phase_shift - current_ctfmodel.phase_shift) > 0.26)
         {
             return heavy_penalization;
         }
     }
     // Now the 1D error
     double distsum = 0;
     int N = 0, Ncorr = 0;
     double enhanced_avg = 0;
     double model_avg = 0;
     double enhanced_model = 0;
     double enhanced2 = 0;
     double model2 = 0;
     double lowerLimit = 1.1 * min_freq_psd;
     double upperLimit = 0.9 * max_freq_psd;
     const MultidimArray<double>& local_enhanced_ctf = psd_exp_enhanced_radial;
     int XdimW=XSIZE(w_digfreq);
     corr13=0;
     FOR_ALL_ELEMENTS_IN_ARRAY1D(w_digfreq)
     {
         if (DIRECT_A1D_ELEM(mask, i) <= 0)
             continue;
 
         // Compute each component
         current_ctfmodel.precomputeValues(i);
         double bg = current_ctfmodel.getValueNoiseAt();
 
         double envelope=0, ctf_without_damping, ctf_with_damping=0, current_envelope = 0;
         double ctf2_th=0;
         double ctf2 = DIRECT_A1D_ELEM(psd_exp_radial, i);
         double dist = 0;
         double ctf_with_damping2;
         switch (action)
         {
         case 0:
         case 1:
             ctf2_th = bg;
             dist = fabs(ctf2 - bg);
             if (penalize && bg > ctf2 && DIRECT_A1D_ELEM(w_digfreq, i) > max_gauss_freq)
                 dist *= current_penalty;
             break;
         case 2:
             envelope = current_ctfmodel.getValueDampingAt();
             ctf2_th = bg + envelope * envelope;
             dist = fabs(ctf2 - ctf2_th);
             if (penalize && ctf2_th < ctf2 && DIRECT_A1D_ELEM(w_digfreq, i) > max_gauss_freq)
                 dist *= current_penalty;
             break;
         case 3:
         case 4:
         case 5:
         case 6:
         case 7:
             envelope = current_ctfmodel.getValueDampingAt();
             ctf_without_damping = current_ctfmodel.getValuePureWithoutDampingAt();
 
             ctf_with_damping = envelope * ctf_without_damping;
             ctf2_th = bg + ctf_with_damping * ctf_with_damping;
 
             if (DIRECT_A1D_ELEM(w_digfreq,i) < upperLimit
                             && DIRECT_A1D_ELEM(w_digfreq,i) > lowerLimit)
             {
                 if  (action == 3 ||
                      (action == 4 && DIRECT_A1D_ELEM(mask_between_zeroes,i) == 1) ||
                      (action == 7 && DIRECT_A1D_ELEM(mask_between_zeroes,i) == 1))
                 {
                     double enhanced_ctf = DIRECT_A1D_ELEM(local_enhanced_ctf, i);
                     ctf_with_damping2 = ctf_with_damping * ctf_with_damping;
                     enhanced_model += enhanced_ctf * ctf_with_damping2;
                     enhanced2 += enhanced_ctf * enhanced_ctf;
                     model2 += ctf_with_damping2 * ctf_with_damping2;
                     enhanced_avg += enhanced_ctf;
                     model_avg += ctf_with_damping2;
                     Ncorr++;
                     if (action==3)
                     {
                         int r = A1D_ELEM(w_digfreq_r,i);
                         A1D_ELEM(psd_theo_radial,r) = ctf2_th;
                     }
                 }
             }
             if (envelope > 1e-2)
                 dist = fabs(ctf2 - ctf2_th) / (envelope * envelope);
             else
                 dist = fabs(ctf2 - ctf2_th);
             // This expression comes from mapping any value so that
             // bg becomes 0, and bg+envelope^2 becomes 1
             // This is the transformation
             //        (x-bg)      x-bg
             //    -------------=-------
             //    (bg+env^2-bg)  env^2
             // If we subtract two of this scaled values
             //    x-bg      y-bg       x-y
             //   ------- - ------- = -------
             //    env^2     env^2     env^2
             break;
 
         }
 
         distsum += dist * DIRECT_A1D_ELEM(mask,i);
         N++;
     }
     if (N > 0)
     { retval = distsum/N ;
     }
     else
         retval = heavy_penalization;
     if (show_inf >=2)
         std::cout << "Fitness1=" << retval << std::endl;
     if ( (((action >= 3) && (action <= 4)) || (action == 7))
          && (Ncorr > 0) && (enhanced_weight != 0) )
     {
 
         model_avg /= Ncorr;
         enhanced_avg /= Ncorr;
         double correlation_coeff = enhanced_model/Ncorr - model_avg * enhanced_avg;
         double sigma1 = sqrt(fabs(enhanced2/Ncorr - enhanced_avg * enhanced_avg));
         double sigma2 = sqrt(fabs(model2/Ncorr - model_avg * model_avg));
         double maxSigma = std::max(sigma1, sigma2);
         if (sigma1 < XMIPP_EQUAL_ACCURACY || sigma2 < XMIPP_EQUAL_ACCURACY
             || (fabs(sigma1 - sigma2) / maxSigma > 0.9 && action>=5))
         {
             retval = heavy_penalization;
             if (show_inf>=2)
                 std::cout << "Fitness2=" << heavy_penalization << " sigma1=" << sigma1 << " sigma2=" << sigma2 << std::endl;
         }
         else
         {
             correlation_coeff /= sigma1 * sigma2;
             if (action == 7)
                 corr13 = correlation_coeff;
             else
                 retval -= enhanced_weight * correlation_coeff;
             if (show_inf >= 2)
             {
                 std::cout << "model_avg=" << model_avg << std::endl;
                 std::cout << "enhanced_avg=" << enhanced_avg << std::endl;
                 std::cout << "enhanced_model=" << enhanced_model / Ncorr
                 << std::endl;
                 std::cout << "sigma1=" << sigma1 << std::endl;
                 std::cout << "sigma2=" << sigma2 << std::endl;
                 std::cout << "Fitness2="
                 << -(enhanced_weight * correlation_coeff)
                 << " (" << correlation_coeff << ")" << std::endl;
             }
         }
 
         // Correlation of the derivative of the radial profile
         if (action==3 || evaluation_reduction==1)
         {
             int state=0;
             double maxDiff=0;
             psd_theo_radial_derivative.initZeros();
             double lowerlimt=1.1*min_freq;
             double upperlimit=0.9*max_freq;
             FOR_ALL_ELEMENTS_IN_ARRAY1D(psd_theo_radial)
             if (A1D_ELEM(w_digfreq_r_iN,i)>0)
             {
                 double freq=A1D_ELEM(w_digfreq,i);
                 switch (state)
                 {
                 case 0:
                     if (freq>lowerlimt)
                         state=1;
                     break;
                 case 1:
                     if (freq>upperlimit)
                         state=2;
                     else
                     {
                         double diff=A1D_ELEM(psd_theo_radial,i)-A1D_ELEM(psd_theo_radial,i-1);
 
                         A1D_ELEM(psd_theo_radial_derivative,i)=diff;
                         maxDiff=std::max(maxDiff,fabs(diff));
                     }
                     break;
                 }
             }
             double corrRadialDerivative=0,mux=0, muy=0, Ncorr=0, sigmax=0, sigmay=0;
             double iMaxDiff=1.0/maxDiff;
 
             FOR_ALL_ELEMENTS_IN_ARRAY1D(psd_theo_radial)
             {
                 A1D_ELEM(psd_theo_radial_derivative,i)*=iMaxDiff;
                 double x=A1D_ELEM(psd_exp_enhanced_radial_derivative,i);
                 double y=A1D_ELEM(psd_theo_radial_derivative,i);
                 corrRadialDerivative+=x*y;
                 mux+=x;
                 muy+=y;
                 sigmax+=x*x;
                 sigmay+=y*y;
                 Ncorr++;
             }
 
             if (Ncorr>0)
             {
                 double iNcorr=1.0/Ncorr;
                 corrRadialDerivative*=iNcorr;
                 mux*=iNcorr;
                 muy*=iNcorr;
                 sigmax=sqrt(fabs(sigmax*iNcorr-mux*mux));
                 sigmay=sqrt(fabs(sigmay*iNcorr-muy*muy));
                 corrRadialDerivative=(corrRadialDerivative-mux*muy)/(sigmax*sigmay);
 
             }
             retval-=corrRadialDerivative;
             if (show_inf>=2)
             {
                 std::cout << "Fitness3=" << -corrRadialDerivative << std::endl;
                 if (show_inf==3)
                 {
                     psd_exp_enhanced_radial.write("PPPexpRadial_fast.txt");
                     psd_theo_radial.write("PPPtheoRadial.txt");
                     psd_exp_enhanced_radial_derivative.write("PPPexpRadialDerivative_fast.txt");
                     psd_theo_radial_derivative.write("PPPtheoRadialDerivative.txt");
                 }
            }
         }
     }
 
     return retval;
 }

◆ defineBasicParams()

void ProgCTFEstimateFromPSDFast::defineBasicParams ( XmippProgram * program )

static

Define basic parameters.

Definition at line 175 of file ctf_estimate_from_psd_fast.cpp.

 {
     CTFDescription::defineParams(program);
 }

◆ defineParams()

void ProgCTFEstimateFromPSDFast::defineParams ( )

virtual

Define Parameters.

Reimplemented from ProgCTFBasicParams.

Definition at line 180 of file ctf_estimate_from_psd_fast.cpp.

 {
     defineBasicParams(this);
     ProgCTFBasicParams::defineParams();
 }

◆ estimate_background_gauss_parameters2_fast()

void ProgCTFEstimateFromPSDFast::estimate_background_gauss_parameters2_fast ( )

Definition at line 1262 of file ctf_estimate_from_psd_fast.cpp.

 {
     if (show_optimization)
         std::cout << "Computing first background Gaussian2 parameters ...\n";
 
     // Compute radial averages
     MultidimArray<double> radial_CTFmodel_avg(YSIZE(*f) / 2);
     MultidimArray<double> radial_CTFampl_avg(YSIZE(*f) / 2);
     MultidimArray<int> radial_N(YSIZE(*f) / 2);
     double w_max_gauss = 0.25;
     FOR_ALL_ELEMENTS_IN_ARRAY1D(w_digfreq)
     {
         if (mask(i) <= 0)
             continue;
         double w = w_digfreq(i);
         if (w > w_max_gauss)
             continue;
 
         int r = FLOOR(w * (double)YSIZE(*f));
         double f_x = DIRECT_A1D_ELEM(x_contfreq, i);
         current_ctfmodel.precomputeValues(f_x);
         double bg = current_ctfmodel.getValueNoiseAt();
         double envelope = current_ctfmodel.getValueDampingAt();
         double ctf_without_damping =
             current_ctfmodel.getValuePureWithoutDampingAt();
         double ctf_with_damping = envelope * ctf_without_damping;
         double ctf2_th = bg + ctf_with_damping * ctf_with_damping;
         radial_CTFmodel_avg(r) += ctf2_th;
         radial_CTFampl_avg(r) += psd_exp_radial(i);
         radial_N(r)++;
     }
 
     // Compute the average radial error
     MultidimArray <double> error;
     error.initZeros(radial_CTFmodel_avg);
     FOR_ALL_ELEMENTS_IN_ARRAY1D(radial_CTFmodel_avg)
     {
         if (radial_N(i) == 0)
             continue;
         error(i) = (radial_CTFampl_avg(i) - radial_CTFmodel_avg(i))/ radial_N(i);
     }
 #ifdef DEBUG
     std::cout << "Error:\n" << error << std::endl;
 #endif
 
     // Compute the frequency of the minimum error
     double wmin = 0.15;
 
     // Compute the maximum (negative) radial error
     double error_max = 0, wmax=0, fmax;
     FOR_ALL_ELEMENTS_IN_ARRAY1D(radial_CTFmodel_avg)
     {
         if (radial_N(i) == 0)
             continue;
         double w = w_digfreq(i);
         if (w > wmin)
             break;
         if (error(i) < error_max)
         {
             wmax = w;
             error_max = error(i);
         }
     }
     fmax = current_ctfmodel.Gc2 = wmax / Tm;
 #ifdef DEBUG
 
     std::cout << "Freq of the maximum error: " << wmax << " " << fmax << std::endl;
 #endif
 
     // Find the linear least squares solution for the gauss part
     Matrix2D<double> A(2, 2);
     A.initZeros();
     Matrix1D<double> b(2);
     b.initZeros();
     int N = 0;
 
     FOR_ALL_ELEMENTS_IN_ARRAY1D(w_digfreq)
     {
         if (mask(i) <= 0)
             continue;
         if (w_digfreq(i) > wmin)
             continue;
         double fmod = w_contfreq(i);
 
         // Compute the zero on the direction of this point
         Matrix1D<double> u(1), fzero(1);
         u = x_contfreq(i) / fmod;
         current_ctfmodel.lookFor(1, u, fzero, 0);
         if (fmod > fzero.module())
             continue;
 
         // Compute weight for this point
         double weight = 1 + max_freq_psd - w_digfreq(i);
         // Compute error
         double f_x = DIRECT_A1D_ELEM(x_contfreq, i);
         current_ctfmodel.precomputeValues(f_x);
         double bg = current_ctfmodel.getValueNoiseAt();
         double envelope = current_ctfmodel.getValueDampingAt();
         double ctf_without_damping =
             current_ctfmodel.getValuePureWithoutDampingAt();
         double ctf_with_damping = envelope * ctf_without_damping;
         double ctf2_th = bg + ctf_with_damping * ctf_with_damping;
         double explained = ctf2_th;
         double unexplained = explained - psd_exp_radial(i);
         if (unexplained <= 0)
             continue;
         unexplained = log(unexplained);
         double F = -(fmod - fmax) * (fmod - fmax);
         A(0, 0) += weight * 1;
         A(0, 1) += weight * F;
         A(1, 1) += weight * F * F;
         b(0) += weight * unexplained;
         b(1) += F * weight * unexplained;
         N++;
     }
 
     if (N != 0)
     {
         A(1, 0) = A(0, 1);
 
         double det=A.det();
         if (fabs(det)>1e-9)
         {
             b = A.inv() * b;
             current_ctfmodel.sigma2 = XMIPP_MIN(fabs(b(1)), 95e3); // This value should be conformant with the physical
             // meaning routine in CTF.cc
             current_ctfmodel.gaussian_K2 = exp(b(0));
         }
         else
         {
             current_ctfmodel.sigma2 = 0;
             current_ctfmodel.gaussian_K2 = 0;
         }
     }
     else
     {
         current_ctfmodel.sigma2 = 0;
         current_ctfmodel.gaussian_K2 = 0;
     }
 
     // Store the CTF values in adjust
     current_ctfmodel.forcePhysicalMeaning();
     COPY_ctfmodel_TO_CURRENT_GUESS;
 
     if (show_optimization)
     {
         std::cout << "First Background Gaussian 2 Fit:\n" << current_ctfmodel
         << std::endl;
         saveIntermediateResults_fast("step04a_first_background2_fit_fast", true);
     }
 }

◆ estimate_background_gauss_parameters_fast()

void ProgCTFEstimateFromPSDFast::estimate_background_gauss_parameters_fast ( )

Definition at line 786 of file ctf_estimate_from_psd_fast.cpp.

 {
 
     if (show_optimization)
         std::cout << "Computing first background Gaussian parameters ...\n";
 
     // Compute radial averages
     MultidimArray<double> radial_CTFmodel_avg(YSIZE(*f) / 2);
     MultidimArray<double> radial_CTFampl_avg(YSIZE(*f) / 2);
     MultidimArray<int> radial_N(YSIZE(*f) / 2);
     double w_max_gauss = 0.25;
 
     FOR_ALL_ELEMENTS_IN_ARRAY1D(w_digfreq)
     {
         if (mask(i) <= 0)
             continue;
         double w = w_digfreq(i);
 
         if (w > w_max_gauss)
             continue;
 
         int r = FLOOR(w * (double)YSIZE(*f));
 
         current_ctfmodel.precomputeValues(x_contfreq(i));
         radial_CTFmodel_avg(r) += current_ctfmodel.getValueNoiseAt();
         radial_CTFampl_avg(r) += psd_exp_radial(i);
         radial_N(r)++;
 
     }
     // Compute the average radial error
     double  error2_avg = 0;
     int N_avg = 0;
     MultidimArray<double> error;
     error.initZeros(radial_CTFmodel_avg);
     FOR_ALL_ELEMENTS_IN_ARRAY1D(radial_CTFmodel_avg)
     {
         if (radial_N(i) == 0)
             continue;
         error(i) = (radial_CTFampl_avg(i) - radial_CTFmodel_avg(i))
                    / radial_N(i);
         error2_avg += error(i) * error(i);
         N_avg++;
     }
     if (N_avg != 0)
         error2_avg /= N_avg;
 
 #ifdef DEBUG
 
     std::cout << "Error2 avg=" << error2_avg << std::endl;
 #endif
 
     // Compute the minimum radial error
     bool first = true, OK_to_proceed = false;
     double error2_min = 0, wmin=0;
     FOR_ALL_ELEMENTS_IN_ARRAY1D(radial_CTFmodel_avg)
     {
         if (radial_N(i) == 0)
             continue;
 
         double w = w_digfreq(i);
 
         if (error(i) < 0 && first){
             continue;
         }
 
         else if (error(i) < 0)
         {
             break;
         }
 
         double error2 = error(i) * error(i);
         // If the two lines cross, do not consider any error until
         // the cross is "old" enough
         if (first && error2 > 0.15 * error2_avg)
             OK_to_proceed = true;
         if (first && i > 0)
             OK_to_proceed &= (error(i) < error(i - 1));
 
         // If the error now is bigger than a 30% (1.69=1.3*1.3) of the error min
         // this must be a rebound. Stop here
         if (!first && error2 > 1.69 * error2_min)
             break;
         if (first && OK_to_proceed)
         {
             wmin = w;
             error2_min = error2;
             first = false;
         }
         if (!first && error2 < error2_min)
         {
             wmin = w;
             error2_min = error2;
         }
 #ifdef DEBUG
         std::cout << w << " " << error2 << " " << wmin << " " << std::endl;
 #endif
 
     }
 
     // Compute the frequency of the minimum error
     max_gauss_freq = wmin;
 #ifdef DEBUG
 
     std::cout << "Freq of the minimum error: " << wmin << " " << fmin << std::endl;
 #endif
 
     // Compute the maximum radial error
     first = true;
     double error2_max = 0, wmax=0, fmax;
     FOR_ALL_ELEMENTS_IN_ARRAY1D(radial_CTFmodel_avg)
     {
         if (radial_N(i) == 0)
             continue;
         double w = w_digfreq(i);
         if (w > wmin)
             continue;
 
         if (error(i) < 0 && first)
             continue;
         else if (error(i) < 0)
             break;
         double error2 = error(i) * error(i);
         if (first)
         {
             wmax = w;
             error2_max = error2;
             first = false;
         }
         if (error2 > error2_max)
         {
             wmax = w;
             error2_max = error2;
         }
 #ifdef DEBUG
         std::cout << w << " " << error2 << " " << wmax << std::endl;
 #endif
 
     }
     fmax = current_ctfmodel.Gc1 = wmax / Tm;
 #ifdef DEBUG
 
     std::cout << "Freq of the maximum error: " << wmax << " " << fmax << std::endl;
 #endif
 
     // Find the linear least squares solution for the gauss part
     Matrix2D<double> A(2, 2);
     A.initZeros();
     Matrix1D<double> b(2);
     b.initZeros();
 
 
     FOR_ALL_ELEMENTS_IN_ARRAY1D(w_digfreq)
     {
 
         if (mask(i) <= 0)
             continue;
 
         if (w_digfreq(i) > wmin)
             continue;
 
         double fmod = w_contfreq(i);
 
         // Compute weight for this point
         double weight = 1 + max_freq_psd - w_digfreq(i);
 
         // Compute error
         current_ctfmodel.precomputeValues(x_contfreq(i));
         double explained = current_ctfmodel.getValueNoiseAt();
         double unexplained = psd_exp_radial(i) - explained;
         if (unexplained <= 0)
             continue;
         unexplained = log(unexplained);
         double F = -(fmod - fmax) * (fmod - fmax);
 
         A(0, 0) += weight * 1;
         A(0, 1) += weight * F;
         A(1, 1) += weight * F * F;
         b(0) += weight * unexplained;
         b(1) += F * weight * unexplained;
 
     }
     A(1, 0) = A(0, 1);
 
     if ( (A(0, 0)== 0) && (A(1, 0)== 0) && (A(1, 1)== 0))
     {
         std::cout << "Matrix A is zeros" << std::endl;
     }
     else
     {
 
         b = A.inv() * b;
         current_ctfmodel.sigma1 = XMIPP_MIN(fabs(b(1)), 95e3); // This value should beconformant with the physical
         // meaning routine in CTF.cc
         current_ctfmodel.gaussian_K = exp(b(0));
         // Store the CTF values in adjust
         current_ctfmodel.forcePhysicalMeaning();
         COPY_ctfmodel_TO_CURRENT_GUESS;
 
         if (show_optimization)
         {
             std::cout << "First Background Fit:\n" << current_ctfmodel << std::endl;
             saveIntermediateResults_fast("step01c_first_background_fit_fast", true);
         }
         center_optimization_focus_fast(true, 1.5);
 
     }
 }

◆ estimate_background_sqrt_parameters_fast()

void ProgCTFEstimateFromPSDFast::estimate_background_sqrt_parameters_fast ( )

Definition at line 672 of file ctf_estimate_from_psd_fast.cpp.

 {
     if (show_optimization)
         std::cout << "Computing first sqrt background ...\n";
 
     // Estimate the base line taking the value of the CTF
     // for the maximum X
     double base_line = 0;
     int N = 0;
     FOR_ALL_ELEMENTS_IN_ARRAY1D(w_digfreq)
     if (w_digfreq(i) > 0.4)
     {
         N++;
         base_line += psd_exp_radial(i);
     }
     if (0 == N) {
         REPORT_ERROR(ERR_NUMERICAL, "N is zero (0), which would lead to division by zero");
     }
     current_ctfmodel.base_line = base_line / N;
 
     // Find the linear least squares solution for the sqrt part
     Matrix2D<double> A(2, 2);
     A.initZeros();
     Matrix1D<double> b(2);
     b.initZeros();
 
     FOR_ALL_ELEMENTS_IN_ARRAY1D(w_digfreq)
     {
         if (mask(i) <= 0)
             continue;
 
         // Compute weight for this point
         double weight = 1 + max_freq_psd - w_digfreq(i);
 
         // Compute error
         current_ctfmodel.precomputeValues(x_contfreq(i));
         double explained = current_ctfmodel.getValueNoiseAt();
         double unexplained = psd_exp_radial(i) - explained;
         if (unexplained <= 0)
             continue;
         unexplained = log(unexplained);
 
         double X = -sqrt(w_contfreq(i));
         A(0, 0) += weight * X * X;
         A(0, 1) += weight * X;
         A(1, 1) += weight * 1;
         b(0) += X * weight * unexplained;
         b(1) += weight * unexplained;
 
     }
     A(1, 0) = A(0, 1);
     b = A.inv() * b;
 
     current_ctfmodel.sq = b(0);
     current_ctfmodel.sqrt_K = exp(b(1));
     COPY_ctfmodel_TO_CURRENT_GUESS;
 
     if (show_optimization)
     {
         std::cout << "First SQRT Fit:\n" << current_ctfmodel << std::endl;
         saveIntermediateResults_fast("step01a_first_sqrt_fit_fast", true);
     }
 
     // Now optimize .........................................................
     double fitness;
     Matrix1D<double> steps;
     steps.resize(SQRT_CTF_PARAMETERS);
     steps.initConstant(1);
 
     // Optimize without penalization
     if (show_optimization)
         std::cout << "Looking for best fitting sqrt ...\n";
     penalize = false;
     int iter;
     powellOptimizer(*adjust_params, FIRST_SQRT_PARAMETER + 1,
                     SQRT_CTF_PARAMETERS, CTF_fitness_fast, global_prm, 0.05, fitness, iter, steps,
                     show_optimization);
 
     // Optimize with penalization
     if (show_optimization)
         std::cout << "Penalizing best fitting sqrt ...\n";
     penalize = true;
     current_penalty = 2;
     int imax = CEIL(log(penalty) / log(2.0));
 
     for (int i = 1; i <= imax; i++)
     {
         if (show_optimization)
             std::cout << "     Iteration " << i << " penalty="
             << current_penalty << std::endl;
         powellOptimizer(*adjust_params, FIRST_SQRT_PARAMETER + 1,
                         SQRT_CTF_PARAMETERS, CTF_fitness_fast, global_prm, 0.05, fitness, iter,
                         steps, show_optimization);
         current_penalty *= 2;
         current_penalty =
             XMIPP_MIN(current_penalty, penalty);
     }
     // Keep the result in adjust
     current_ctfmodel.forcePhysicalMeaning();
     COPY_ctfmodel_TO_CURRENT_GUESS;
 
     if (show_optimization)
     {
         std::cout << "Best penalized SQRT Fit:\n" << current_ctfmodel
         << std::endl;
         saveIntermediateResults_fast("step01b_best_penalized_sqrt_fit_fast", true);
     }
 
    center_optimization_focus_fast(true, 1.5);
 
 }

◆ estimate_defoci_fast()

void ProgCTFEstimateFromPSDFast::estimate_defoci_fast ( )

Definition at line 1075 of file ctf_estimate_from_psd_fast.cpp.

 {
     double fitness;
     int iter;
     Matrix1D<double> steps(DEFOCUS_PARAMETERS);
     steps.initConstant(1);
     steps(1) = 0; // Do not optimize kV
     steps(2) = 0; // Do not optimize K
     if(!selfEstimation)
     {
         (*adjust_params)(0) = initial_ctfmodel.Defocus;
         (*adjust_params)(2) = current_ctfmodel.K;
         powellOptimizer(*adjust_params, FIRST_DEFOCUS_PARAMETER + 1,
                                 DEFOCUS_PARAMETERS, CTF_fitness_fast, global_prm, 0.05,
                                 fitness, iter, steps, false);
     }
     else
     {
         /*
          * Defocus estimation s^2 stigmator
          *
          * CTF = sin(PI*lambda*defocus*R^2)
          * CTF = sin(PI*lambda*defocus*w^2/Ts^2)
          * w^2 = omega
          * omega = 2*K0/Xdim   K0: position of the maximum of the psd in fourier
          *
          * Defocus = 2*K0*(2*Ts)^2/lambda
          */
         MultidimArray<double> psd_exp_radial2;
         MultidimArray<double> background;
         MultidimArray<double> psd_background;
         action = 1;
         generateModelSoFar_fast(background, false);
         action = 3;
         psd_background.initZeros(background);
 
 #ifdef DEBUG
         std::cout << "background\n" << background << std::endl;
         std::cout << "--------------------------------------" << std::endl;
 
         std::cout << "psd_exp_radial\n" << psd_exp_radial << std::endl;
         std::cout << "--------------------------------------" << std::endl;
 #endif
 
         FOR_ALL_ELEMENTS_IN_ARRAY1D(background)
         {
             if(w_digfreq(i)>min_freq && w_digfreq(i)<max_freq)
                 psd_background(i)= background(i)-psd_exp_radial(i);
         }
 
         //Uncomment to apply a high-pass FILTER to psd_background
 #ifdef FILTER
         double alpha = 0.1;
         Matrix1D<double> A_coeff(2);
         VEC_ELEM(A_coeff,0) = 1;
         VEC_ELEM(A_coeff,1) = alpha-1;
         Matrix1D<double> B_coeff(2);
         VEC_ELEM(B_coeff,0) = 1-alpha;
         VEC_ELEM(B_coeff,1) = alpha-1;
         MultidimArray<double> psd_background_filter;
         MultidimArray<double> aux_filter;
         psd_background_filter.initZeros(psd_background);
         aux_filter.initZeros(psd_background);
         matlab_filter(B_coeff,A_coeff,psd_background,psd_background_filter,aux_filter);
         psd_background = psd_background_filter;
 #endif
 
 #ifdef DEBUG
         std::cout << "psdbackground\n" << psd_background << std::endl;
         std::cout << "--------------------------------------" << std::endl;
 #endif
         psd_exp_radial2.initZeros(psd_exp_radial);
 
         double deltaW=1.0/(2*XSIZE(w_digfreq)*current_ctfmodel.Tm);
         double deltaW2=1.0/(XSIZE(w_digfreq)*pow(2*current_ctfmodel.Tm,2));
         FOR_ALL_ELEMENTS_IN_ARRAY1D(psd_exp_radial2)
         {
             double w2=i*deltaW2;
             double w=sqrt(w2);
             double widx=w/deltaW;
             size_t lowerIdx=floor(widx);
             double weight=widx-floor(widx);
             if(i < XSIZE(psd_exp_radial2)-1)
             {
                 double value =(1-weight)*A1D_ELEM(psd_background,lowerIdx)
                                                      +weight*A1D_ELEM(psd_background,lowerIdx+1);
 
                 A1D_ELEM(psd_exp_radial2,i) = value*value;
             }
         }
 #ifdef DEBUG
         std::cout << "psd_exp_radial2\n" << psd_exp_radial2 << std::endl;
         std::cout << "--------------------------------------" << std::endl;
 #endif
         FourierTransformer FourierPSD;
         FourierPSD.FourierTransform(psd_exp_radial2, psd_fft, false);
 
         const double minDefocus=2000;
         int startIndex = ceil(std::max((double)7,minDefocus*current_ctfmodel.lambda/(2*pow(2*current_ctfmodel.Tm,2))));
 #ifdef DEBUG
         std::cout << "Start index=" << minDefocus*current_ctfmodel.lambda/(2*pow(2*current_ctfmodel.Tm,2)) << std::endl;
         std::cout << "--------------------------------------" << std::endl;
 #endif
         if (current_ctfmodel.VPP_radius != 0.0)
             startIndex = 10; //avoid low frequencies
         FOR_ALL_ELEMENTS_IN_ARRAY1D(psd_fft)
         {
             if(i>=startIndex) {
                 amplitud.push_back(abs(psd_fft[i]));
 #ifdef DEBUG
         std::cout << abs(psd_fft[i]) << std::endl;
 #endif
             }
         }
 #ifdef DEBUG
         std::cout << "--------------------------------------" << std::endl;
 #endif
 
         //Uncomment to remove outliers from amplitud
 #ifdef OUTLIERS
         double totalSum = std::accumulate(amplitud.begin(), amplitud.end(), 0.0);
 
         double amplitudMean = totalSum / amplitud.size();
         double sdSum = 0;
         for(double i : amplitud)
         {
             double diff=amplitud[i]-amplitudMean;
             sdSum += diff*diff;
         }
         double amplitudSD = sqrt(sdSum/amplitud.size());
         double differenceSD = 3*amplitudSD;
         for(double i : amplitud)
         {
             if (amplitud[i]>amplitudMean+differenceSD){
                 amplitud[i] = amplitudMean+differenceSD;
             }
         }
 #endif
 
         double maxValue=-1e38;
         int finalIndex=-1;
 #ifdef DEBUG
         std::cout << "Looking for maximum ...\n";
 #endif
         for (size_t i=1; i<amplitud.size()-1; i++)
         {
 #ifdef DEBUG
             std::cout << "i=" << i << " amplitude i= " << amplitud[i] << std::endl;
 #endif
             if (amplitud[i]>maxValue && amplitud[i]>=amplitud[i+1])
             {
                 maxValue=amplitud[i];
                 finalIndex=i;
 #ifdef DEBUG
                 std::cout << "New maximum " << i << " maxValue=" << maxValue << std::endl;
 #endif
             }
         }
 
         //double maxValue = *max_element(amplitud.begin(),amplitud.end());
         //int finalIndex = distance(amplitud.begin(),max_element(amplitud.begin(),amplitud.end()));
         current_ctfmodel.Defocus = (floor(finalIndex+startIndex+1)*pow(2*current_ctfmodel.Tm,2))/
                                                 current_ctfmodel.lambda;
 
 #ifdef DEBUG
         std::cout << "maxValue: " << maxValue << std::endl;
         std::cout << "finalIndex: " << finalIndex << std::endl;
         std::cout << "defocus: " << current_ctfmodel.Defocus << std::endl;
 #endif
     }
     // Keep the result in adjust
     current_ctfmodel.forcePhysicalMeaning();
     COPY_ctfmodel_TO_CURRENT_GUESS;
     ctfmodel_defoci = current_ctfmodel;
 
     if  (show_optimization)
     {
         std::cout << "First defocus Fit:\n" << ctfmodel_defoci << std::endl;
         saveIntermediateResults_fast("step03a_first_defocus_fit_fast", true);
     }
 }

◆ estimate_envelope_parameters_fast()

void ProgCTFEstimateFromPSDFast::estimate_envelope_parameters_fast ( )

Definition at line 996 of file ctf_estimate_from_psd_fast.cpp.

 {
     if (show_optimization)
         std::cout << "Looking for best fitting envelope ...\n";
 
     // Set the envelope
     current_ctfmodel.Ca = initial_ctfmodel.Ca;
     current_ctfmodel.K = 1.0;
     current_ctfmodel.espr = 0.0;
     current_ctfmodel.ispr = 0.0;
     current_ctfmodel.alpha = 0.0;
     current_ctfmodel.DeltaF = 0.0;
     current_ctfmodel.DeltaR = 0.0;
     current_ctfmodel.Q0 = initial_ctfmodel.Q0;
     current_ctfmodel.envR0 = 0.0;
     current_ctfmodel.envR1 = 0.0;
     COPY_ctfmodel_TO_CURRENT_GUESS;
 
     // Now optimize the envelope
     penalize = false;
     int iter;
     double fitness;
     Matrix1D<double> steps;
     steps.resize(ENVELOPE_PARAMETERS);
     steps.initConstant(1);
 
     steps(1) = 0; // Do not optimize Cs
     steps(5) = 0; // Do not optimize for alpha, since Ealpha depends on the defocus
 
     if (modelSimplification >= 1)
         steps(6) = steps(7) = 0; // Do not optimize DeltaF and DeltaR
     powellOptimizer(*adjust_params, FIRST_ENVELOPE_PARAMETER + 1,
                     ENVELOPE_PARAMETERS, CTF_fitness_fast, global_prm, 0.05, fitness, iter, steps,
                     show_optimization);
 
     // Keep the result in adjust
     current_ctfmodel.forcePhysicalMeaning();
     COPY_ctfmodel_TO_CURRENT_GUESS;
 
     if (show_optimization)
     {
         std::cout << "Best envelope Fit:\n" << current_ctfmodel << std::endl;
         saveIntermediateResults_fast("step02a_best_envelope_fit_fast", true);
     }
 
     // Optimize with penalization
     if (show_optimization)
         std::cout << "Penalizing best fitting envelope ...\n";
     penalize = true;
     current_penalty = 2;
     int imax = CEIL(log(penalty) / log(2.0));
     for (int i = 1; i <= imax; i++)
     {
         if (show_optimization)
             std::cout << "     Iteration " << i << " penalty="
             << current_penalty << std::endl;
         powellOptimizer(*adjust_params, FIRST_ENVELOPE_PARAMETER + 1,
                         ENVELOPE_PARAMETERS, CTF_fitness_fast, global_prm, 0.05, fitness, iter,
                         steps, show_optimization);
         current_penalty *= 2;
         current_penalty =
             XMIPP_MIN(current_penalty, penalty);
     }
     // Keep the result in adjust
     current_ctfmodel.forcePhysicalMeaning();
     COPY_ctfmodel_TO_CURRENT_GUESS;
 
     if (show_optimization)
     {
         std::cout << "Best envelope Fit:\n" << current_ctfmodel << std::endl;
         saveIntermediateResults_fast("step02b_best_penalized_envelope_fit_fast", true);
     }
 
 }

◆ generateModelSoFar_fast()

void ProgCTFEstimateFromPSDFast::generateModelSoFar_fast	(	MultidimArray< double > &	I,
		bool	apply_log = `false`
	)

Definition at line 203 of file ctf_estimate_from_psd_fast.cpp.

 {
     Matrix1D<int> idx(1); // Indexes for Fourier plane
     Matrix1D<double> freq(1); // Frequencies for Fourier plane
 
     assignCTFfromParameters_fast(MATRIX1D_ARRAY(*adjust_params), current_ctfmodel,
                             0, ALL_CTF_PARAMETERS);
     current_ctfmodel.produceSideInfo();
 
     I.initZeros(psd_exp_radial);
     FOR_ALL_ELEMENTS_IN_ARRAY1D(I)
     {
         XX(idx) = i;
         FFT_idx2digfreq(*f, idx, freq);
         double w=freq.module();
         if (w>max_freq_psd)
             continue;
         digfreq2contfreq(freq, freq, Tm);
         // Decide what to save
         current_ctfmodel.precomputeValues(XX(freq));
         if (action <= 1)
             I(i) = current_ctfmodel.getValueNoiseAt();
         else if (action == 2)
         {
             double E = current_ctfmodel.getValueDampingAt();
             I(i) = current_ctfmodel.getValueNoiseAt() + E * E;
 
         }
         else if (action >= 3 && action <= 6)
         {
             double ctf = current_ctfmodel.getValuePureAt();
             I(i) = current_ctfmodel.getValueNoiseAt() + ctf * ctf;
         }
         else
         {
             double ctf = current_ctfmodel.getValuePureAt();
             I(i) = ctf;
         }
         if (apply_log)
             I(i) = 10 * log10(I(i));
 
     }
 }

◆ produceSideInfo()

void ProgCTFEstimateFromPSDFast::produceSideInfo ( )

Produce side information.

Definition at line 187 of file ctf_estimate_from_psd_fast.cpp.

 {
     adjust.resize(ALL_CTF_PARAMETERS);
     adjust.initZeros();
     current_ctfmodel.clear();
     ctfmodel_defoci.clear();
     assignParametersFromCTF_fast(initial_ctfmodel, MATRIX1D_ARRAY(adjust), 0, ALL_CTF_PARAMETERS);
     // Read the CTF file, supposed to be the uncentered squared amplitudes
     if (fn_psd != "")
         ctftomodel.read(fn_psd);
     f = &(ctftomodel());
 
     ProgCTFBasicParams::produceSideInfo();
     current_ctfmodel.precomputeValues(x_contfreq);
 }

◆ readBasicParams()

void ProgCTFEstimateFromPSDFast::readBasicParams ( XmippProgram * program )

Read parameters.

Definition at line 155 of file ctf_estimate_from_psd_fast.cpp.

 {
     ProgCTFBasicParams::readBasicParams(program);
 
     initial_ctfmodel.enable_CTF = initial_ctfmodel.enable_CTFnoise = true;
     initial_ctfmodel.readParams(program);
     initial_ctfmodel2D.readParams(program);
 
     if (initial_ctfmodel.Defocus>100e3)
         REPORT_ERROR(ERR_ARG_INCORRECT,"Defocus cannot be larger than 10 microns (100,000 Angstroms)");
     Tm = initial_ctfmodel.Tm;
 }

◆ readParams()

void ProgCTFEstimateFromPSDFast::readParams ( )

virtual

Read parameters.

Reimplemented from ProgCTFBasicParams.

Definition at line 168 of file ctf_estimate_from_psd_fast.cpp.

 {
     fn_psd = getParam("--psd");
     readBasicParams(this);
 }

◆ run()

void ProgCTFEstimateFromPSDFast::run ( )

virtual

Run

Reimplemented from XmippProgram.

Definition at line 1692 of file ctf_estimate_from_psd_fast.cpp.

 {
         CTFDescription1D ctf1Dmodel;
         ROUT_Adjust_CTFFast(*this, ctf1Dmodel);
 }

◆ saveIntermediateResults_fast()

void ProgCTFEstimateFromPSDFast::saveIntermediateResults_fast	(	const FileName &	fn_root,
		bool	generate_profiles = `true`
	)

Definition at line 247 of file ctf_estimate_from_psd_fast.cpp.

 {
     std::ofstream plot_radial;
     MultidimArray<double> save;
     generateModelSoFar_fast(save, false);
     if (!generate_profiles)
         return;
     plot_radial.open((fn_root + "_radial.txt").c_str());
     if (!plot_radial)
         REPORT_ERROR(
             ERR_IO_NOWRITE,
             "save_intermediate_results::Cannot open plot file for writing\n");
 
     //plot_radial << "# freq_dig freq_angstrom model psd enhanced logModel logPsd\n";
 
     // Generate radial average
     MultidimArray<double> radial_CTFmodel_avg;
     radial_CTFmodel_avg.initZeros(psd_exp_enhanced_radial);
 
     FOR_ALL_ELEMENTS_IN_ARRAY1D(w_digfreq)
     {
         if (mask(i) <= 0)
             continue;
         double model2 = save(i);
 
         int r = w_digfreq_r(i);
         radial_CTFmodel_avg(r) = model2;
 
         plot_radial << w_digfreq(r, 0) << " " << w_contfreq(r, 0)
         << " " << radial_CTFmodel_avg(r) << " "
         << psd_exp_enhanced_radial(r) << " " << psd_exp_radial(r)
         << " " << log10(radial_CTFmodel_avg(r)) << " "
         << std::endl;
     }
 
     plot_radial.close();
 }

◆ showFirstDefoci_fast()

void ProgCTFEstimateFromPSDFast::showFirstDefoci_fast ( )

Member Data Documentation

◆ ctfmodel_defoci

CTFDescription1D ProgCTFEstimateFromPSDFast::ctfmodel_defoci

Definition at line 41 of file ctf_estimate_from_psd_fast.h.

◆ current_ctfmodel

CTFDescription1D ProgCTFEstimateFromPSDFast::current_ctfmodel

Definition at line 41 of file ctf_estimate_from_psd_fast.h.

◆ initial_ctfmodel

CTFDescription1D ProgCTFEstimateFromPSDFast::initial_ctfmodel

Definition at line 41 of file ctf_estimate_from_psd_fast.h.

◆ initial_ctfmodel2D

CTFDescription ProgCTFEstimateFromPSDFast::initial_ctfmodel2D

Definition at line 42 of file ctf_estimate_from_psd_fast.h.

The documentation for this class was generated from the following files:

xmipp/libraries/reconstruction/ctf_estimate_from_psd_fast.h
xmipp/libraries/reconstruction/ctf_estimate_from_psd_fast.cpp

Public Member Functions

Static Public Member Functions

Public Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ ProgCTFEstimateFromPSDFast()

Member Function Documentation

◆ assignCTFfromParameters_fast()

◆ assignParametersFromCTF_fast()

◆ center_optimization_focus_fast()

◆ compute_central_region_fast()

◆ CTF_fitness_object_fast()

◆ defineBasicParams()

◆ defineParams()

◆ estimate_background_gauss_parameters2_fast()

◆ estimate_background_gauss_parameters_fast()

◆ estimate_background_sqrt_parameters_fast()

◆ estimate_defoci_fast()

◆ estimate_envelope_parameters_fast()

◆ generateModelSoFar_fast()

◆ produceSideInfo()

◆ readBasicParams()

◆ readParams()

◆ run()

◆ saveIntermediateResults_fast()

◆ showFirstDefoci_fast()

Member Data Documentation

◆ ctfmodel_defoci

◆ current_ctfmodel

◆ initial_ctfmodel

◆ initial_ctfmodel2D