mlf_align2d.h

Go to the documentation of this file.

/***************************************************************************
 *
 * Authors:    Sjors Scheres           scheres@cnb.csic.es (2004)
 *
 * Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
 * 02111-1307  USA
 *
 *  All comments concerning this program package may be sent to the
 *  e-mail address 'xmipp@cnb.csic.es'
 ***************************************************************************/

#ifndef MLFALIGN2D_H
#define MLFALIGN2D_H

#include "ml2d.h"
#include <numeric>


#define FOR_ALL_MODELS() for (int refno=0;refno<model.n_ref; refno++)
#define FOR_ALL_ROTATIONS() for (size_t ipsi=0; ipsi<nr_psi; ipsi++ )
#define FOR_ALL_FLIPS() for (size_t iflip=0; iflip<nr_flip; iflip++)
#define FOR_ALL_LIMITED_TRANSLATIONS() for (size_t itrans=0; itrans<nr_trans; itrans++)
#define FOR_ALL_DEFOCUS_GROUPS() for (size_t ifocus=0; ifocus<nr_focus; ifocus++)
#define FOR_ALL_DIGITAL_FREQS() for (size_t irr = 0; irr < hdim; irr++)
#define FOR_ALL_POINTS() for (size_t ipoint = 0; ipoint < nr_points_2d; ++ipoint)

//Helper macros to define the elements of vectors:
//Vsig, Vctf and Vsnr indexed by ifocus and irr
#define VSNR_ITEM dAi(Vsnr[ifocus], irr)
#define VCTF_ITEM dAi(Vctf[ifocus], irr)
#define VDEC_ITEM dAi(Vdec[ifocus], irr)
#define VSIG_ITEM dAi(Vsig[ifocus], irr)

#define SIGNIFICANT_WEIGHT_LOW 1e-8
constexpr double  SMALLVALUE = 1e-4;
constexpr float  HISTMIN = -6.;
constexpr float  HISTMAX = 6.;
constexpr int  HISTSTEPS = 120;

#define FN_EXTRA(file) formatString("%sextra/%s", fn_root.c_str(), file)
#define FN_NOISE_IMG_MD FN_EXTRA("noise_images.xmd")
#define FN_NOISE_IMG    FN_EXTRA("noise_images.stk")
#define FN_CREF_IMG     FN_EXTRA("cref_classes.stk")
#define FN_CREF_IMG_MD  FN_EXTRA("cref_classes.xmd")

#define FN_ITER_BASE(iter) getIterExtraPath(fn_root, iter)
#define FN_REF(base, refno) formatString("%06d@%sclasses.stk", (refno), (base).c_str())
#define FN_VSIG(base, ifocus, ext) ((nr_focus > 1) ? formatString("ctf%06d@%s%s", ((ifocus) + 1), (base).c_str(), (ext)) : ((base) + "_ctf" + (ext)))


class ProgMLF2D: public ML2DBaseProgram
{
public:

    double sigma_offset;
    std::vector<double> alpha_k;
    std::vector<double> mirror_fraction;
    size_t max_nr_psi;
    bool do_variable_psi, do_variable_trans;
    std::vector< Image<double> > Ictf;
    bool limit_trans;
    size_t nr_trans;
    size_t zero_trans;
    int search_shift;
    int offsets_keepdir;

    //CTFDat ctfdat;
    bool do_ctf_correction;
    double sampling;
    std::vector<int> count_defocus;
    bool phase_flipped;
    MultidimArray<size_t> Mresol_int;
    std::vector<MultidimArray<double> > Vsig, Vctf, Vdec;
    double reduce_snr;
    size_t nr_focus;
    double lowres_limit, highres_limit, ini_highres_limit;
    bool first_iter_noctf;
    bool do_divide_ctf;
    bool do_include_allfreqs;
    double fix_high;
    std::vector<int> pointer_2d, pointer_i, pointer_j;
    size_t nr_points_prob, nr_points_2d, dnr_points_2d;
    size_t current_highres_limit;



    bool do_student;
    bool do_student_sigma_trick;


    bool do_kstest;
    int iter_write_histograms;
    Histogram1D sumhist;
    std::vector<Histogram1D > resolhist;

    std::vector<double> refs_avgscale;

    int rank, size, nr_vols; //mpi rank, size and 3d number of vols

    double LL, sumcorr, wsum_sigma_offset, sumw_allrefs;
    std::vector<MultidimArray<double> > wsum_Mref, wsum_ctfMref;
    std::vector< std::vector<double> > Mwsum_sigma2;
    std::vector<double> sumw, sumw2, sumwsc, sumwsc2, sumw_mirror, sumw_defocus;
    std::vector<double> Fref, Fwsum_imgs, Fwsum_ctfimgs;


    ProgMLF2D(int nr_vols = 0, int rank = 0, int size = 1);

    void defineParams();

    void readParams();


public:
    void show(bool ML3D = false);

    virtual void produceSideInfo();

    virtual void produceSideInfo2();

    virtual void defineBasicParams(XmippProgram * prog);
    virtual void defineAdditionalParams(XmippProgram * prog, const char * sectionLine);

    void estimateInitialNoiseSpectra();

    void updateWienerFilters(const MultidimArray<double> &spectral_signal,
                             std::vector<double> &sumw_defocus, int iter);

    void setCurrentSamplingRates(double current_probres_limit);

    void generateInitialReferences();


    void calculatePdfInplane();

    // Append a FourierTransform (in half format!) to a vector
    void appendFTtoVector(const MultidimArray<std::complex<double> > &Fin,
                          std::vector<double> &out);

    // get a FourierTransform (in half format!) from a vector
    void getFTfromVector(const std::vector<double> &in,
                         const int start_point,
                         MultidimArray<std::complex<double> > &out,
                         bool only_real = false);

    void rotateReference(std::vector<double> &out);

    void reverseRotateReference(const std::vector<double> &in,
                                std::vector<MultidimArray<double > > &out);

    void preselectDirections(float &phi, float &theta,
                             std::vector<double> &pdf_directions);

    void preselect_significant_model_phi(MultidimArray<double> &Mimg, std::vector<double> &offsets,
                                         std::vector <std::vector< MultidimArray<double > > > &Mref,
                                         MultidimArray<int> &Msignificant,
                                         std::vector<double > &pdf_directions);

    // Calculate the FT of a translated matrix using a phase shift in
    // Fourier space
    void fourierTranslate2D(const std::vector<double> &in,
                            MultidimArray<double> &trans,
                            std::vector<double> &out,
                            int point_start = 0);

    // If not determined yet: search optimal offsets using maxCC
    // Then for all optimal translations, calculate all translated FTs
    // for each of the flipped variants
    void calculateFourierOffsets(const MultidimArray<double> &Mimg,
                                 const std::vector<double > &offsets,
                                 std::vector<double>  &out,
                                 MultidimArray<int> &Moffsets,
                                 MultidimArray<int> &Moffsets_mirror);

    void processOneImage(const MultidimArray<double> &Mimg,
                         size_t focus, bool apply_ctf,
                         double &fracweight,  double &maxweight2, double &sum_refw2,
                         double &opt_scale, size_t &opt_refno, double &opt_psi,
                         size_t &opt_ipsi, size_t &opt_iflip,
                         MultidimArray<double> &opt_offsets,
                         std::vector<double> &opt_offsets_ref,
                         std::vector<double > &pdf_directions,
                         bool do_kstest, bool write_histograms,
                         FileName fn_img, double &KSprob);

    double performKSTest(MultidimArray<double> &Mimg,  const int focus, bool apply_ctf,
                         FileName &fn_img, bool write_histogram,
                         std::vector <std::vector< MultidimArray<std::complex<double> > > > &Fref,
                         double &opt_scale, int &opt_refno, int &opt_ipsi, int &opt_iflip,
                         MultidimArray<double> &opt_offsets);

    void iteration();
    void expectation();
    void maximization();
    void endIteration();

    virtual void writeOutputFiles(const ModelML2D &model, OutputType outputType);

    void writeNoiseFile(const FileName &fn_base, int ifocus);

    virtual void addPartialDocfileData(const MultidimArray<double> &data, size_t first, size_t last);

};

#endif