ml_tomo.h File Reference

#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>

Include dependency graph for ml_tomo.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
struct	structThreadExpectationSingleImage
class	ProgMLTomo
struct	ProgMLTomo::MissingInfo
struct	ProgMLTomo::AnglesInfo

Macros
#define	SIGNIFICANT_WEIGHT_LOW   1e-8
#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)

Functions
void *	threadMLTomoExpectationSingleImage (void *data)

Macro Definition Documentation

FN_ITER_MD

#define FN_ITER_MD

(

iter

)

formatString("%s_it%d.sel", fn_root.c_str(), iter)

Definition at line 53 of file ml_tomo.h.

FN_ITER_VOL

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

Definition at line 52 of file ml_tomo.h.

MLTOMO_BLOCKSIZE

#define MLTOMO_BLOCKSIZE   10

Definition at line 50 of file ml_tomo.h.

MLTOMO_DATALINELENGTH

#define MLTOMO_DATALINELENGTH   10

Definition at line 49 of file ml_tomo.h.

SIGNIFICANT_WEIGHT_LOW

#define SIGNIFICANT_WEIGHT_LOW   1e-8

Definition at line 46 of file ml_tomo.h.

SMALLANGLE

#define SMALLANGLE   2.75

Definition at line 48 of file ml_tomo.h.

Function Documentation

threadMLTomoExpectationSingleImage()

void* threadMLTomoExpectationSingleImage

(

void *

data

)

Definition at line 2702 of file ml_tomo.cpp.

{
    auto * thread_data =
        (structThreadExpectationSingleImage *) data;

    // Variables from above
    int thread_id = thread_data->thread_id;
    ProgMLTomo *prm = thread_data->prm;
    MetaDataVec *MDimg = thread_data->MDimg;
    double *wsum_sigma_noise = thread_data->wsum_sigma_noise;
    double *wsum_sigma_offset = thread_data->wsum_sigma_offset;
    double *sumfracweight = thread_data->sumfracweight;
    double *LL = thread_data->LL;
    std::vector<MultidimArray<double> > *wsumimgs = thread_data->wsumimgs;
    std::vector<MultidimArray<double> > *wsumweds = thread_data->wsumweds;
    std::vector<Image<double> > *Iref = thread_data->Iref;

    MultidimArray<double> *sumw = thread_data->sumw;
    std::vector<size_t> * imgs_id = thread_data->imgs_id;
    ThreadTaskDistributor * distributor = thread_data->distributor;

    //#define DEBUG_THREAD
#ifdef DEBUG_THREAD

    std::cerr<<"start threadMLTomoExpectationSingleImage"<<std::endl;
#endif

    // Local variables
    Image<double> img;
    FileName fn_img, fn_trans;
    std::vector<MultidimArray<double> > allref_offsets;
    Matrix1D<double> opt_offsets(3);
    float old_psi = -999.;
    double fracweight, trymindiff, dLL;
    int opt_refno, opt_angno, missno;

    //try
    {
        //Approximate number of images
        size_t myNum = MDimg->size() / prm->threads;
        //First and last image to process in each block
        size_t firstIndex, lastIndex;

        if (prm->verbose && thread_id == 0)
            init_progress_bar(myNum);

        // Loop over all images
        size_t cc = 0;
        MultidimArray<double> &docfiledata = prm->docfiledata;

        //Work while there are tasks to do
        while (distributor->getTasks(firstIndex, lastIndex))
        {
            for (size_t imgno = firstIndex; imgno <= lastIndex; ++imgno)
            {
                //TODO: Check if really needed the mutexes
                //only for read from MetaData
                pthread_mutex_lock(&mltomo_selfile_access_mutex);

                MDimg->getValue(MDL_IMAGE, fn_img, (*imgs_id)[imgno + prm->myFirstImg]);

                pthread_mutex_unlock(&mltomo_selfile_access_mutex);

                img.read(fn_img);
                img().setXmippOrigin();

                if (prm->dont_align || prm->do_only_average)
                {
                    selfTranslate(xmipp_transformation::LINEAR, img(), prm->imgs_optoffsets[imgno], xmipp_transformation::DONT_WRAP);
                }
                prm->reScaleVolume(img(), true);
                missno = (prm->do_missing) ? prm->imgs_missno[imgno] : -1;
                // These three parameters speed up expectationSingleImage
                trymindiff = prm->imgs_trymindiff[imgno];
                opt_refno = prm->imgs_optrefno[imgno];
                opt_angno = prm->imgs_optangno[imgno];
                old_psi = prm->imgs_optpsi[imgno];

                if (prm->do_ml)
                {
                    // A. Use maximum likelihood approach
                    prm->expectationSingleImage(img(), imgno, missno, old_psi, *Iref,
                                                *wsumimgs, *wsumweds, *wsum_sigma_noise, *wsum_sigma_offset,
                                                *sumw, *LL, dLL, fracweight, *sumfracweight, trymindiff,
                                                opt_refno, opt_angno, opt_offsets);
                }
                else
                {
                    // B. Use constrained correlation coefficient approach
                    prm->maxConstrainedCorrSingleImage(img(), imgno, missno, old_psi,
                                                       *Iref, *wsumimgs, *wsumweds, *sumw, fracweight, *sumfracweight,
                                                       opt_refno, opt_angno, opt_offsets);
                }
                // Store for next iteration
                prm->imgs_trymindiff[imgno] = trymindiff;
                prm->imgs_optrefno[imgno] = opt_refno;
                prm->imgs_optangno[imgno] = opt_angno;
                // Output MetaData
                //FIXME: THIS ALSO LIKE JoseMiguel did ML2D
                dAij(docfiledata, imgno, 0) = (prm->all_angle_info[opt_angno]).rot; // rot
                dAij(docfiledata, imgno, 1) = (prm->all_angle_info[opt_angno]).tilt; // tilt
                dAij(docfiledata, imgno, 2) = (prm->all_angle_info[opt_angno]).psi; // psi
                if (prm->dont_align || prm->do_only_average)
                {
                    dAij(docfiledata, imgno, 3) = prm->imgs_optoffsets[imgno](0); // Xoff
                    dAij(docfiledata, imgno, 4) = prm->imgs_optoffsets[imgno](1); // Yoff
                    dAij(docfiledata, imgno, 5) = prm->imgs_optoffsets[imgno](2); // zoff
                }
                else
                {
                    dAij(docfiledata, imgno, 3) = opt_offsets(0) / prm->scale_factor; // Xoff
                    dAij(docfiledata, imgno, 4) = opt_offsets(1) / prm->scale_factor; // Yoff
                    dAij(docfiledata, imgno, 5) = opt_offsets(2) / prm->scale_factor; // Zoff
                }
                dAij(docfiledata, imgno, 6) = (double) (opt_refno + 1); // Ref
                dAij(docfiledata, imgno, 7) = (double) (missno + 1); // Wedge number

                dAij(docfiledata, imgno, 8) = fracweight; // P_max/P_tot
                dAij(docfiledata, imgno, 9) = dLL; // log-likelihood

                if (prm->verbose && thread_id == 0)
                    progress_bar(cc);
                cc++;
            }
        }
        if (prm->verbose && thread_id == 0)
            progress_bar(myNum);
    }

#ifdef DEBUG_THREAD

    std::cerr<<"finished threadMLTomoExpectationSingleImage"<<std::endl;
#endif
    return nullptr;
}