adjust_ctf (Adjust CTF parameters to PSD)

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

Classes
class	ProgCTFEstimateFromPSDFast

Functions
double	ROUT_Adjust_CTFFast (ProgCTFEstimateFromPSDFast &prm, CTFDescription1D &output_ctfmodel, bool standalone=true)

Detailed Description

Function Documentation

ROUT_Adjust_CTFFast()

double ROUT_Adjust_CTFFast	(	ProgCTFEstimateFromPSDFast &	prm,
		CTFDescription1D &	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 1417 of file ctf_estimate_from_psd_fast.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() * 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_fast();
        prm.estimate_background_gauss_parameters_fast();
    }

     // 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(1) = steps(2) = steps(4) = 0;
    powellOptimizer(*prm.adjust_params, FIRST_SQRT_PARAMETER + 1,
                    BACKGROUND_CTF_PARAMETERS, CTF_fitness_fast, 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_fast("step01d_best_background_fit_fast", true);
    }

    DEBUG_TEXTFILE(formatString("Step 4: CTF_fitness=%f",CTF_fitness_fast));
    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.kV = prm.initial_ctfmodel.kV;
    prm.current_ctfmodel.Cs = prm.initial_ctfmodel.Cs;
    prm.current_ctfmodel.phase_shift = prm.initial_ctfmodel.phase_shift;
    prm.current_ctfmodel.VPP_radius = prm.initial_ctfmodel.VPP_radius;

    if (prm.initial_ctfmodel.Q0 != 0)
        prm.current_ctfmodel.Q0 = prm.initial_ctfmodel.Q0;
    prm.estimate_envelope_parameters_fast();

    if (prm.show_optimization)
    {
        std::cout << "Best envelope Fit:\n" << prm.current_ctfmodel << std::endl;
        prm.saveIntermediateResults_fast("step02b_best_penalized_envelope_fit_fast", true);
    }

    DEBUG_TEXTFILE(formatString("Step 6: espr=%f",prm.current_ctfmodel.espr));
    DEBUG_MODEL_TEXTFILE;

    /************************************************************************/
    /* STEP 7:  defocus and angular parameters                              */
    /************************************************************************/
    prm.action = 3;
    if (!prm.noDefocusEstimate)
       prm.estimate_defoci_fast();
    else
        prm.current_ctfmodel.Defocus = prm.initial_ctfmodel.Defocus;
    DEBUG_TEXTFILE(formatString("Step 7: Defocus=%f",prm.current_ctfmodel.Defocus));
    DEBUG_MODEL_TEXTFILE;
    /************************************************************************/
    /* STEPs 9, 10 and 11: all parameters included second Gaussian          */
    /************************************************************************/
    prm.action = 5;
    if (prm.modelSimplification < 2)
        prm.estimate_background_gauss_parameters2_fast();

    steps.resize(ALL_CTF_PARAMETERS);
    steps.initConstant(1);
    steps(1) = 0; // kV
    steps(3) = 0; // The spherical aberration (Cs) is not optimized
    steps(27) = 0; //VPP radius not optimized
    if (prm.noDefocusEstimate)
        steps(0)=0; // Defocus
    if (prm.initial_ctfmodel.Q0 != 0)
        steps(13) = 0; // Q0
    if (prm.modelSimplification >= 1)
        steps(8) = steps(9) = 0;
    if (prm.initial_ctfmodel.VPP_radius == 0)
        steps(26) = 0; //VPP phase shift

    powellOptimizer(*prm.adjust_params, 0 + 1, ALL_CTF_PARAMETERS, CTF_fitness_fast,
                    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_fast("step04b_best_fit_with_gaussian2_fast", true);
    }

    /************************************************************************/
    /* STEP 12: 2D estimation parameters                                    */
    /************************************************************************/
    prm.adjust_params->resize(ALL_CTF_PARAMETERS2D);
    auto *prm2D = new ProgCTFEstimateFromPSD(&prm);
    if (prm.noDefocusEstimate)
    {
        prm2D->current_ctfmodel.DeltafU=prm.initial_ctfmodel2D.DeltafU;
        prm2D->current_ctfmodel.DeltafV=prm.initial_ctfmodel2D.DeltafV;
        prm2D->current_ctfmodel.azimuthal_angle=prm.initial_ctfmodel2D.azimuthal_angle;
        prm2D->noDefocusEstimate=true;
    }

    steps.resize(ALL_CTF_PARAMETERS2D);
    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.noDefocusEstimate)
        steps(0)=steps(1)=steps(2)=0; // Defocus and azimuthal angle
    if (prm2D->initial_ctfmodel.Q0 != 0)
        steps(15) = 0; // Q0
     if (prm2D->modelSimplification >= 3)
        steps(20) = steps(21) = steps(23) = 0;
    if (prm2D->modelSimplification >= 2)
        steps(24) = steps(25) = steps(26) = steps(27) = steps(28) = steps(29) = 0;
    if (prm2D->modelSimplification >= 1)
        steps(10) = steps(11) = 0;
    if(prm2D->initial_ctfmodel.VPP_radius == 0)
        steps(36) = 0;

    prm2D->adjust_params = prm.adjust_params;
    prm2D->adjust_params->initZeros();
    prm2D->assignParametersFromCTF(prm2D->current_ctfmodel,MATRIX1D_ARRAY(*prm2D->adjust_params),0, ALL_CTF_PARAMETERS2D, true);

    if (prm.refineAmplitudeContrast)
        steps(15) = 1; // Q0
    powellOptimizer(*prm2D->adjust_params, 0 + 1, ALL_CTF_PARAMETERS2D, CTF_fitness,
                     prm2D, 0.01, fitness, iter, steps, prm2D->show_optimization);

    prm2D->current_ctfmodel.forcePhysicalMeaning();

    //We adopt that always  DeltafU > DeltafV so if this is not the case we change the values and the angle
    if (!prm.noDefocusEstimate && prm2D->current_ctfmodel.DeltafV > prm2D->current_ctfmodel.DeltafU)
    {
        double temp;
        temp = prm2D->current_ctfmodel.DeltafU;
        prm2D->current_ctfmodel.DeltafU = prm2D->current_ctfmodel.DeltafV;
        prm2D->current_ctfmodel.DeltafV = temp;
        prm2D->current_ctfmodel.azimuthal_angle -= 90;
    }
    prm2D->assignParametersFromCTF(prm2D->current_ctfmodel,MATRIX1D_ARRAY(*prm2D->adjust_params),0, ALL_CTF_PARAMETERS2D, true);
    if (prm2D->show_optimization)
    {
        std::cout << "Best fit with 2D parameters:\n" << prm2D->current_ctfmodel << std::endl;
        prm2D->saveIntermediateResults("step05b_estimate_2D_parameters", true);
    }

    /************************************************************************/
    /* STEP 13: Produce output                                              */
    /************************************************************************/
    prm2D->action = 7;
    if (prm2D->fn_psd != "")
    {
        // Define mask between first and third zero
        prm2D->mask_between_zeroes.initZeros(prm2D->mask);
        Matrix1D<double> u(2), z1(2), z3(2);
#ifdef NEVERDEFINED
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(prm2D->mask_between_zeroes)
        {
            VECTOR_R2(u, DIRECT_MULTIDIM_ELEM(prm2D->x_digfreq, n), DIRECT_MULTIDIM_ELEM(prm2D->y_digfreq, n));
            u /= u.module();
            prm2D->current_ctfmodel.lookFor(3, u, z3, 0);
            if (DIRECT_MULTIDIM_ELEM(prm2D->w_contfreq, n) >= z3.module())
                continue;
            prm2D->current_ctfmodel.lookFor(1, u, z1, 0);
            if (z1.module() < DIRECT_MULTIDIM_ELEM(prm2D->w_contfreq, n))
                DIRECT_MULTIDIM_ELEM(prm2D->mask_between_zeroes, n) = 1;

        }
#endif
        XX(u)=1; YY(u)=0;
        prm2D->current_ctfmodel.lookFor(3, u, z3, 0);
        prm2D->current_ctfmodel.lookFor(1, u, z1, 0);
        double z1m = z1.module();
        double z3m = z3.module();
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(prm2D->mask_between_zeroes)
        {
            double wn=DIRECT_MULTIDIM_ELEM(prm2D->w_contfreq, n);
            if (wn<=z3m && wn>=z1m)
                DIRECT_MULTIDIM_ELEM(prm2D->mask_between_zeroes, n) = 1;
        }
        // Evaluate the correlation in this region
        CTF_fitness(prm2D->adjust_params->adaptForNumericalRecipes(), prm2D);
        // Save results
        FileName fn_rootCTFPARAM = prm2D->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;

        prm2D->saveIntermediateResults(fn_rootMODEL, false);
        prm2D->current_ctfmodel.Tm /= prm2D->downsampleFactor;
        prm2D->current_ctfmodel.azimuthal_angle = std::fmod(prm2D->current_ctfmodel.azimuthal_angle,360.);
        prm2D->current_ctfmodel.phase_shift = (prm2D->current_ctfmodel.phase_shift*180)/3.14;
        if(prm2D->current_ctfmodel.phase_shift<0.0)
            prm2D->current_ctfmodel.phase_shift = 0.0;

        prm2D->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*prm2D->Tm, id);
        MD.setValue(MDL_CTF_XF, (double)output_ctfmodel.xF*prm2D->Tm, id);
        MD.setValue(MDL_CTF_Y0, (double)output_ctfmodel.y0*prm2D->Tm, id);
        MD.setValue(MDL_CTF_YF, (double)output_ctfmodel.yF*prm2D->Tm, id);
        MD.setValue(MDL_CTF_CRIT_FITTINGSCORE, fitness, id);
        MD.setValue(MDL_CTF_CRIT_FITTINGCORR13, prm2D->corr13, id);
        MD.setValue(MDL_CTF_CRIT_ICENESS, evaluateIceness(prm.ctftomodel(), prm.Tm), id);
        MD.setValue(MDL_CTF_DOWNSAMPLE_PERFORMED, prm2D->downsampleFactor, id);
        MD.write(fn_rootCTFPARAM + ".ctfparam",MD_APPEND);
        fn_rootCTFPARAM = fn_rootCTFPARAM + ".ctfparam_tmp";

        fn_rootCTFPARAM.deleteFile();
    }
    output_ctfmodel = prm2D->current_ctfmodel;
    return fitness;
}