adjust_ctf (Adjust CTF parameters to PSD)

Collaboration diagram for adjust_ctf (Adjust CTF parameters to PSD):

Classes
class	ProgCTFEstimateFromPSD

Functions
double	evaluateIceness (const MultidimArray< double > &enhanced_ctftomodel, double Tm)
double	ROUT_Adjust_CTF (ProgCTFEstimateFromPSD &prm, CTFDescription &output_ctfmodel, bool standalone=true)

Detailed Description

Function Documentation

evaluateIceness()

double evaluateIceness	(	const MultidimArray< double > &	enhanced_ctftomodel,
		double	Tm
	)

ROUT_Adjust_CTF()

double ROUT_Adjust_CTF	(	ProgCTFEstimateFromPSD &	prm,
		CTFDescription &	output_ctfmodel,
		bool	standalone = true
	)

Core of the Adjust CTF routine. This is the routine which does everything. It returns the fitting error committed in the best fit.

Definition at line 2208 of file ctf_estimate_from_psd.cpp.

{
    DEBUG_OPEN_TEXTFILE(prm.fn_psd.removeLastExtension());
    global_prm = &prm;
    if (standalone || prm.show_optimization)
        prm.show();
    prm.produceSideInfo();
    DEBUG_TEXTFILE(formatString("After producing side info: Avg=%f",prm.ctftomodel().computeAvg()));
    DEBUG_MODEL_TEXTFILE;

    // Build initial frequency mask
    prm.value_th = -1;
    prm.min_freq_psd = prm.min_freq;
    prm.max_freq_psd = prm.max_freq;

    // Set some global variables
    prm.adjust_params = &prm.adjust;
    prm.penalize = false;
    prm.max_gauss_freq = 0;
    prm.heavy_penalization = prm.f->computeMax() * XSIZE(*prm.f) * YSIZE(*prm.f);
    prm.show_inf = 0;

    // Some variables needed by all steps
    int iter;
    double fitness;
    Matrix1D<double> steps;

    /************************************************************************/
    /* STEPs 1, 2, 3 and 4:  Find background which best fits the CTF        */
    /************************************************************************/

    prm.current_ctfmodel.enable_CTFnoise = true;
    prm.current_ctfmodel.enable_CTF = false;
    prm.evaluation_reduction = 4;

    // If initial parameters were not supplied for the gaussian curve,
    // estimate them from the CTF file
    prm.action = 0;
    if (prm.adjust(FIRST_SQRT_PARAMETER) == 0)
    {
        prm.estimate_background_sqrt_parameters();
        prm.estimate_background_gauss_parameters();
    }

    // Optimize the current background
    prm.action = 1;
    prm.penalize = true;
    prm.current_penalty = prm.penalty;
    steps.resize(BACKGROUND_CTF_PARAMETERS);
    steps.initConstant(1);
    if (!prm.modelSimplification >= 3)
        steps(7) = steps(8) = steps(10) = 0;
    powellOptimizer(*prm.adjust_params, FIRST_SQRT_PARAMETER + 1,
                    BACKGROUND_CTF_PARAMETERS, CTF_fitness, global_prm, 0.01, fitness, iter,
                    steps, prm.show_optimization);

    // Make sure that the model has physical meaning
    // (In some machines due to numerical imprecission this check is necessary
    // at the end)
    prm.current_ctfmodel.forcePhysicalMeaning();
    COPY_ctfmodel_TO_CURRENT_GUESS;

    if (prm.show_optimization)
    {
        std::cout << "Best background Fit:\n" << prm.current_ctfmodel << std::endl;
        prm.saveIntermediateResults("step01d_best_background_fit");
    }
    DEBUG_TEXTFILE(formatString("Step 4: CTF_fitness=%f",CTF_fitness));
    DEBUG_MODEL_TEXTFILE;

    /************************************************************************/
    /* STEPs 5 and 6:  Find envelope which best fits the CTF                */
    /************************************************************************/
    prm.action = 2;
    prm.current_ctfmodel.enable_CTF = true;
    prm.current_ctfmodel.phase_shift = prm.initial_ctfmodel.phase_shift;
    prm.current_ctfmodel.VPP_radius = prm.initial_ctfmodel.VPP_radius;
    if (prm.initial_ctfmodel.K == 0)
    {
        prm.current_ctfmodel.kV = prm.initial_ctfmodel.kV;
        prm.current_ctfmodel.Cs = prm.initial_ctfmodel.Cs;
        if (prm.initial_ctfmodel.Q0 != 0)
            prm.current_ctfmodel.Q0 = prm.initial_ctfmodel.Q0;
        prm.estimate_envelope_parameters();
    }
    else
    {
        prm.current_ctfmodel.K = prm.initial_ctfmodel.K;
        prm.current_ctfmodel.kV = prm.initial_ctfmodel.kV;
        prm.current_ctfmodel.DeltafU = prm.initial_ctfmodel.DeltafU;
        prm.current_ctfmodel.DeltafV = prm.initial_ctfmodel.DeltafV;
        prm.current_ctfmodel.azimuthal_angle = prm.initial_ctfmodel.azimuthal_angle;
        prm.current_ctfmodel.Cs = prm.initial_ctfmodel.Cs;
        prm.current_ctfmodel.Ca = prm.initial_ctfmodel.Ca;
        prm.current_ctfmodel.espr = prm.initial_ctfmodel.espr;
        prm.current_ctfmodel.ispr = prm.initial_ctfmodel.ispr;
        prm.current_ctfmodel.alpha = prm.initial_ctfmodel.alpha;
        prm.current_ctfmodel.DeltaF = prm.initial_ctfmodel.DeltaF;
        prm.current_ctfmodel.DeltaR = prm.initial_ctfmodel.DeltaR;
        prm.current_ctfmodel.Q0 = prm.initial_ctfmodel.Q0;
        COPY_ctfmodel_TO_CURRENT_GUESS;
    }
    DEBUG_TEXTFILE(formatString("Step 6: espr=%f",prm.current_ctfmodel.espr));
    DEBUG_MODEL_TEXTFILE;
    /************************************************************************/
    /* STEP 7:  the defocus and angular parameters                          */
    /************************************************************************/

    prm.action = 3;
    prm.evaluation_reduction = 1;
    if (prm.fastDefocusEstimate)
        prm.estimate_defoci_Zernike();
    else if (!prm.noDefocusEstimate)
        prm.estimate_defoci();
    else
    {
        prm.current_ctfmodel.DeltafU=prm.initial_ctfmodel.DeltafU;
        prm.current_ctfmodel.DeltafV=prm.initial_ctfmodel.DeltafV;
        if (prm.current_ctfmodel.DeltafV==0.0)
            prm.current_ctfmodel.DeltafV=prm.current_ctfmodel.DeltafU;
    }

    DEBUG_TEXTFILE(formatString("Step 7: DeltafU=%f",prm.current_ctfmodel.DeltafU));
    DEBUG_TEXTFILE(formatString("Step 7: DeltafV=%f",prm.current_ctfmodel.DeltafV));
    DEBUG_TEXTFILE(formatString("Step 7: azimutalAngle=%f",prm.current_ctfmodel.azimuthal_angle));
    DEBUG_MODEL_TEXTFILE;

    //This line is to test the results obtained
    //exit(1);

    /************************************************************************/
    /* STEPs 9, 10 and 11: all parameters included second Gaussian          */
    /************************************************************************/
    prm.action = 5;
    if (prm.modelSimplification < 2)
        prm.estimate_background_gauss_parameters2();

    steps.resize(ALL_CTF_PARAMETERS);
    steps.initConstant(1);
    steps(3) = 0; // kV
    steps(5) = 0; // The spherical aberration (Cs) is not optimized
    steps(37) = 0; //VPP radius not optimized
    if (prm.initial_ctfmodel.Q0 != 0)
        steps(15) = 0; // Q0
    if (prm.modelSimplification >= 3)
        steps(20) = steps(21) = steps(23) = 0;
    if (prm.modelSimplification >= 2)
        steps(24) = steps(25) = steps(26) = steps(27) = steps(28) = steps(29) = 0;
    if (prm.modelSimplification >= 1)
        steps(10) = steps(11) = 0;
    if (std::floor(prm.initial_ctfmodel.VPP_radius) == 0)
        steps(36) = 0; //VPP phase shift

    powellOptimizer(*prm.adjust_params, 0 + 1, ALL_CTF_PARAMETERS, CTF_fitness,
                    global_prm, 0.01, fitness, iter, steps, prm.show_optimization);

    prm.current_ctfmodel.forcePhysicalMeaning();
    COPY_ctfmodel_TO_CURRENT_GUESS;

    if (prm.show_optimization)
    {
        std::cout << "Best fit with Gaussian2:\n" << prm.current_ctfmodel << std::endl;
        prm.saveIntermediateResults("step04b_best_fit_with_gaussian2");
    }

    prm.evaluation_reduction = 2;
    powellOptimizer(*prm.adjust_params, 0 + 1, ALL_CTF_PARAMETERS, CTF_fitness,
                    global_prm, 0.01, fitness, iter, steps, prm.show_optimization);
    prm.current_ctfmodel.forcePhysicalMeaning();
    COPY_ctfmodel_TO_CURRENT_GUESS;

    prm.evaluation_reduction = 1;
    powellOptimizer(*prm.adjust_params, 0 + 1, ALL_CTF_PARAMETERS, CTF_fitness,
                    global_prm, 0.005, fitness, iter, steps, prm.show_optimization);
    prm.current_ctfmodel.forcePhysicalMeaning();
    COPY_ctfmodel_TO_CURRENT_GUESS;

    if (prm.refineAmplitudeContrast)
    {
        std::cout << "Refining amplitude contrast ..." << std::endl;
        steps(15) = 1; // Q0
        powellOptimizer(*prm.adjust_params, 0 + 1, ALL_CTF_PARAMETERS, CTF_fitness,
                         global_prm, 0.005, fitness, iter, steps, prm.show_optimization);
        prm.current_ctfmodel.forcePhysicalMeaning();
        COPY_ctfmodel_TO_CURRENT_GUESS;
    }

    if (prm.show_optimization)
    {
        std::cout << "Best fit:\n" << prm.current_ctfmodel << std::endl;
        prm.saveIntermediateResults("step04c_best_fit");
    }
    DEBUG_TEXTFILE(formatString("Step 11: DeltafU=%f fitness=%f",prm.current_ctfmodel.DeltafU,fitness));
    DEBUG_MODEL_TEXTFILE;

    //We adopt that always  DeltafU > DeltafV so if this is not the case we change the values and the angle
    if ( prm.current_ctfmodel.DeltafV > prm.current_ctfmodel.DeltafU)
    {
        double temp;
        temp = prm.current_ctfmodel.DeltafU;
        prm.current_ctfmodel.DeltafU = prm.current_ctfmodel.DeltafV;
        prm.current_ctfmodel.DeltafV = temp;
        prm.current_ctfmodel.azimuthal_angle -= 90;
        COPY_ctfmodel_TO_CURRENT_GUESS;
    }

    /************************************************************************/
    /* STEP 12: Produce output                                              */
    /************************************************************************/

    prm.action = 7;

    if (prm.fn_psd != "")
    {
        // Define mask between first and third zero
        prm.mask_between_zeroes.initZeros(prm.mask);
        Matrix1D<double> u(2), z1(2), z3(2);
        FOR_ALL_ELEMENTS_IN_ARRAY2D(prm.mask_between_zeroes)
        {
            VECTOR_R2(u, prm.x_digfreq(i, j), prm.y_digfreq(i, j));
            u /= u.module();
            prm.current_ctfmodel.lookFor(1, u, z1, 0);
            prm.current_ctfmodel.lookFor(3, u, z3, 0);
            if (z1.module() < prm.w_contfreq(i, j)
                && prm.w_contfreq(i, j) < z3.module())
                prm.mask_between_zeroes(i, j) = 1;
        }
        // Evaluate the correlation in this region
        CTF_fitness(prm.adjust_params->adaptForNumericalRecipes(), &prm);
        // Save results
        FileName fn_rootCTFPARAM = prm.fn_psd.withoutExtension();

        FileName fn_rootMODEL = fn_rootCTFPARAM;
        size_t atPosition=fn_rootCTFPARAM.find('@');

        if (atPosition!=std::string::npos)
        {
            fn_rootMODEL=formatString("%03d@%s",textToInteger(fn_rootCTFPARAM.substr(0, atPosition)),
                                      fn_rootCTFPARAM.substr(atPosition+1).c_str());
            fn_rootCTFPARAM=formatString("region%03d@%s",textToInteger(fn_rootCTFPARAM.substr(0, atPosition)),
                                         fn_rootCTFPARAM.substr(atPosition+1).c_str());
        }
        else
            fn_rootCTFPARAM=(String)"fullMicrograph@"+fn_rootCTFPARAM;

        prm.saveIntermediateResults(fn_rootMODEL, false);
        prm.current_ctfmodel.Tm /= prm.downsampleFactor;
        prm.current_ctfmodel.azimuthal_angle = std::fmod(prm.current_ctfmodel.azimuthal_angle,360.);
        prm.current_ctfmodel.write(fn_rootCTFPARAM + ".ctfparam_tmp");
        MetaDataVec MD;
        MD.read(fn_rootCTFPARAM + ".ctfparam_tmp");
        size_t id = MD.firstRowId();
        MD.setValue(MDL_CTF_X0, (double)output_ctfmodel.x0*prm.Tm, id);
        MD.setValue(MDL_CTF_XF, (double)output_ctfmodel.xF*prm.Tm, id);
        MD.setValue(MDL_CTF_Y0, (double)output_ctfmodel.y0*prm.Tm, id);
        MD.setValue(MDL_CTF_YF, (double)output_ctfmodel.yF*prm.Tm, id);
        MD.setValue(MDL_CTF_CRIT_FITTINGSCORE, fitness, id);
        MD.setValue(MDL_CTF_CRIT_FITTINGCORR13, prm.corr13, id);
        MD.setValue(MDL_CTF_CRIT_ICENESS, evaluateIceness(prm.psd_exp_radial, prm.Tm), id);
        MD.setValue(MDL_CTF_DOWNSAMPLE_PERFORMED, prm.downsampleFactor, id);
        MD.write(fn_rootCTFPARAM + ".ctfparam",MD_APPEND);
        fn_rootCTFPARAM = fn_rootCTFPARAM + ".ctfparam_tmp";
        fn_rootCTFPARAM.deleteFile();
    }
    output_ctfmodel = prm.current_ctfmodel;

    DEBUG_CLOSE_TEXTFILE;
    return fitness;
}