ml_tomo.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 _MLTOMO_H
#define _MLTOMO_H

#include <core/xmipp_fftw.h>
#include <core/xmipp_fft.h>
#include <core/args.h>
#include <core/xmipp_funcs.h>
#include <core/metadata_extension.h>
#include <core/xmipp_image.h>
#include <core/geometry.h>
#include <data/filters.h>
#include <data/mask.h>
#include <data/ctf.h>
#include <data/sampling.h>
#include <core/symmetries.h>
#include "symmetrize.h"
#include <core/xmipp_threads.h>
#include <vector>
#include <core/xmipp_program.h>

#define SIGNIFICANT_WEIGHT_LOW 1e-8
#undef SMALLANGLE
#define SMALLANGLE 2.75
#define MLTOMO_DATALINELENGTH 10
#define MLTOMO_BLOCKSIZE 10

#define FN_ITER_VOL(iter, base, refno) formatString("%s_it%d_%s%d.vol", fn_root.c_str(), iter, base, refno+1)
#define FN_ITER_MD(iter) formatString("%s_it%d.sel", fn_root.c_str(), iter)

class ProgMLTomo;

// Thread declaration
void * threadMLTomoExpectationSingleImage( void * data );

// This structure is needed to pass parameters to threadMLTomoExpectationSingleImage
typedef struct
{
    int thread_id;
    int thread_num;
    ProgMLTomo *prm;
    MetaDataVec *MDimg;
    int *iter;
    double *wsum_sigma_noise;
    double *wsum_sigma_offset;
    double *sumfracweight;
    double *LL;
    std::vector<MultidimArray<double> > *wsumimgs;
    std::vector<MultidimArray<double> > *wsumweds;
    std::vector<Image<double> > *Iref;
    std::vector<MultidimArray<double> > *docfiledata;
    MultidimArray<double> * sumw;
    std::vector<size_t> * imgs_id;
    ThreadTaskDistributor * distributor;
}
structThreadExpectationSingleImage ;

class ProgMLTomo: public XmippProgram
{
public:
    FileName fn_sel, fn_ref, fn_root, fn_frac, fn_sym, fn_missing, fn_doc, fn_mask;
    std::string cline;
    double sigma_noise;
    double sigma_offset;
    MultidimArray<double> alpha_k;
    bool fix_fractions;
    bool fix_sigma_offset;
    bool fix_sigma_noise;
    int istart;
    int Niter, Niter2;
    int oridim, dim, dim3, hdim;
    double ddim3;
    int nr_ref;
    bool do_keep_angles;
    size_t nr_images_global;
    size_t nr_images_local;
    size_t myFirstImg, myLastImg;
    MultidimArray<double > docfiledata;
    std::vector<double> A2, corrA2;
    double eps;
    MetaDataVec MDimg, MDref, MDmissing;
    bool mdimg_contains_angles;
    std::vector < Image<double> > Iref, Iold, Iwed;
    MultidimArray<double> P_phi, Mr2;
    bool do_generate_refs;
    std::vector <size_t> convert_refno_to_stack_position;

    std::vector<int> imgs_optrefno, imgs_optangno;
    std::vector<double> imgs_trymindiff, imgs_optpsi;
    std::vector<Matrix1D<double> > imgs_optoffsets;
    std::vector<int> imgs_missno;
    double trymindiff_factor;
    double ang_search;
    bool do_limit_psirange;
    double limit_trans;
    bool do_perturb;
    bool do_filter, do_ini_filter;
    bool do_sym;


    bool do_missing;//, do_wedge, do_pyramid, do_cone;
    bool do_impute;
    bool do_ml;
    double noimp_threshold;
    double maxres, scale_factor;
    MultidimArray<double> fourier_mask, real_mask, real_omask;
    MultidimArray<unsigned char> fourier_imask;

    bool dont_align;
    bool dont_rotate;

    bool do_mask;
    Image<double> Imask;
    double nr_mask_pixels;
    bool do_only_average;

    // Missing data information
    enum MissingType { MISSING_WEDGE_Y, MISSING_WEDGE_X, MISSING_PYRAMID, MISSING_CONE };
    struct MissingInfo
    {
        MissingType type;
        double thy0, thyF, thx0, thxF;
        double tg0_y, tgF_y, tg0_x, tgF_x;
        bool do_limit_x, do_limit_y, do_cone;
        double nr_pixels;
    };
    int nr_miss;
    std::vector<MissingInfo> all_missing_info;

    // Angular sampling information
    struct AnglesInfo
    {
        size_t direction;
        double rot, rot_ori;
        double tilt, tilt_ori;
        double psi, psi_ori;
        Matrix2D<double> A, A_ori;
    };
    typedef std::vector<AnglesInfo> AnglesInfoVector;


    double angular_sampling, psi_sampling;
    AnglesInfoVector all_angle_info;
    int nr_ang;
    double pixel_size;

    double reg0, regF, reg_current;
    int reg_steps;

    bool no_SMALLANGLE;

    // sampling object
    Sampling mysampling;
    // For user-provided tilt range
    double tilt_range0, tilt_rangeF;
    // Symmetry setup
    int symmetry, sym_order;

    int threads;

    FourierTransformer transformer;

    std::vector<size_t> imgs_id;

public:
    void defineParams();

    void readParams();

    void run();

    void show();

    virtual void produceSideInfo();

    virtual void generateInitialReferences();

    virtual void setNumberOfLocalImages();

    virtual void produceSideInfo2(int nr_vols = 1);

    void perturbAngularSampling();

    void calculatePdfTranslations();

    void readMissingInfo();

    void getMissingRegion(MultidimArray<unsigned char> &Mmeasured,
                          const Matrix2D<double> &A,
                          const int missno);

    void maskSphericalAverageOutside(MultidimArray<double> &Min);

    // Resize a volume, based on the max_resol
    // if down_scale=true: go from oridim to dim
    // if down_scale=false: go from dim to oridim
    void reScaleVolume(MultidimArray<double> &Min, bool down_scale=true);

    void postProcessVolume(Image<double> &Vin, double resolution = -1.);

    void precalculateA2(std::vector< Image<double> > &Iref);

    void expectationSingleImage(MultidimArray<double> &Mimg, int imgno, const int missno, double old_rot,
                                std::vector<Image<double> > &Iref,
                                std::vector<MultidimArray<double> > &wsumimgs,
                                std::vector<MultidimArray<double> > &wsumweds,
                                double &wsum_sigma_noise, double &wsum_sigma_offset,
                                MultidimArray<double> &sumw, double &LL, double &dLL,
                                double &fracweight, double &sumfracweight, double &trymindiff,
                                int &opt_refno, int &opt_angno, Matrix1D<double> &opt_offsets);

    void maxConstrainedCorrSingleImage(MultidimArray<double> &Mimg, int imgno, int missno, double old_rot,
                                       std::vector<Image<double> > &Iref,
                                       std::vector<MultidimArray<double> > &wsumimgs,
                                       std::vector<MultidimArray<double> > &wsumweds,
                                       MultidimArray<double> &sumw, double &maxCC, double &sumCC,
                                       int &opt_refno, int &opt_angno, Matrix1D<double> &opt_offsets);

    virtual void expectation(MetaDataVec &MDimg, std::vector< Image<double> > &Iref, int iter,
                     double &LL, double &sumfracweight,
                     std::vector<MultidimArray<double> > &wsumimgs,
                     std::vector<MultidimArray<double> > &wsumweds,
                     double &wsum_sigma_noise, double &wsum_sigma_offset,
                     MultidimArray<double> &sumw);

    void maximization(std::vector<MultidimArray<double> > &wsumimgs,
                      std::vector<MultidimArray<double> > &wsumweds,
                      double &wsum_sigma_noise, double &wsum_sigma_offset,
                      MultidimArray<double> &sumw, double &sumfracweight,
                      double &sumw_allrefs, std::vector<MultidimArray<double> > &fsc, int iter);

    void calculateFsc(MultidimArray<double> &M1, MultidimArray<double> &M2,
                      MultidimArray<double> &W1, MultidimArray<double> &W2,
                      MultidimArray< double >& freq, MultidimArray< double >& fsc,
                      double &resolution);

    void regularize(int iter);

    bool checkConvergence(std::vector<double> &conv);

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

    virtual void writeOutputFiles(const int iter,
                          std::vector<MultidimArray<double> > &wsumweds,
                          double &sumw_allrefs, double &LL, double &avefracweight,
                          std::vector<double> &conv, std::vector<MultidimArray<double> > &fsc);

};
#endif