recons_misc.cpp File Reference

#include <fstream>
#include "recons_misc.h"
#include "basic_art.h"
#include "core/metadata_vec.h"
#include "core/symmetries.h"
#include "data/mask.h"
#include "data/projection.h"
#include "data/wavelet.h"
#include "symmetrize.h"

Include dependency graph for recons_misc.cpp:

Go to the source code of this file.

Macros
#define	MODE8
#define	NFEATURES   8
#define	GET_RESULTS(fh, fn, avg, cov, N, idx)
#define	POCS_N_measure   8
#define	POCS_N_use   2

Functions
void	buildReconsInfo (MetaDataVec &selfile, const FileName &fn_ctf, const SymList &SL, ReconsInfo *&IMG_Inf, bool do_not_use_symproj)
void	sortPerpendicular (int numIMG, ReconsInfo *IMG_Inf, MultidimArray< int > &ordered_list, int N)
void	noSort (int numIMG, MultidimArray< int > &ordered_list)
void	sortRandomly (int numIMG, MultidimArray< int > &ordered_list)
void	updateResidualVector (BasicARTParameters &prm, GridVolume &vol_basis, double &kappa, double &pow_residual_vol, double &pow_residual_imgs)

Macro Definition Documentation

GET_RESULTS

#define GET_RESULTS	(		fh,
			fn,
			avg,
			cov,
			N,
			idx
	)

Value:

fh.open(fn);\
    if (!fh) \
        REPORT_ERROR(ERR_IO_NOTOPEN,(std::string)"VariabilityClass::finishAnalysis: " \
                     "Cannot open "+fn+" for input"); \
    n_previous=-1; \
    while (!fh.eof()) { \
        int n; fh >> n; \
        if (n!=n_previous) { \
            n_previous=n; N++; \
            idx(n)=1; \
            avg+=v[n]; \
            aux.fromVector(v[n]); \
            cov+=aux*aux.transpose(); \
        }\
    } \
    avg/=N; \
    cov/=N; \
    aux.fromVector(avg); \
    cov-=aux*aux.transpose();

MODE8

#define MODE8

Definition at line 412 of file recons_misc.cpp.

NFEATURES

#define NFEATURES   8

POCS_N_measure

#define POCS_N_measure   8

Definition at line 696 of file recons_misc.cpp.

POCS_N_use

#define POCS_N_use   2

Definition at line 697 of file recons_misc.cpp.

Function Documentation

buildReconsInfo()

void buildReconsInfo	(	MetaDataVec &	selfile,
		const FileName &	fn_ctf,
		const SymList &	SL,
		ReconsInfo *&	IMG_Inf,
		bool	do_not_use_symproj
	)

Build from a Selection File and a Symmetry List. The result is stored in the Recons_info array which should point to NULL when it is not initialized.

Definition at line 37 of file recons_misc.cpp.

{
    Matrix2D<double>  L(4, 4), R(4, 4);  // A matrix from the list
    FileName          fn_proj;
    FileName          fn_ctf1;
    Projection        read_proj;
    bool              is_there_ctf = false;
    bool              is_ctf_unique = false;

    int trueIMG = selfile.size();
    selfile.firstRowId();
    int numIMG;
    if (!do_not_use_symproj)
        numIMG = trueIMG * (SL.symsNo() + 1);
    else
        numIMG = trueIMG;

    // The next two ifs check whether there is a CTF file, and
    // whether it is unique
    if (fn_ctf != "")
    {
        is_there_ctf = true;
        is_ctf_unique = true;
    }
    else if (selfile.containsLabel(MDL_CTF_MODEL))
    {
        is_there_ctf = true;
        is_ctf_unique = false;
    }

    if (IMG_Inf != nullptr)
        delete [] IMG_Inf;
    if ((IMG_Inf = new ReconsInfo[numIMG]) == nullptr)
        REPORT_ERROR(ERR_MEM_NOTENOUGH, "Build_Recons_Info: No memory for the sorting");

    int i = 0; // It will account for the number of valid projections processed
    std::cout << "Reading angle information ...\n";
    init_progress_bar(trueIMG);
    for (size_t objId : selfile.ids())
    {
        ReconsInfo &imgInfo = IMG_Inf[i];
        selfile.getValue(MDL_IMAGE,fn_proj,objId);
        if (is_there_ctf && !is_ctf_unique)
            selfile.getValue(MDL_CTF_MODEL,fn_ctf1,objId);
        if (fn_proj != "")
        {
            //            read_proj.read(fn_proj, false, HEADER);
            // Filling structure
            imgInfo.fn_proj = fn_proj;
            imgInfo.row = selfile.getRowVec(objId);
            imgInfo.row.detach();
            if (is_ctf_unique)
                imgInfo.fn_ctf = fn_ctf;
            else if (is_there_ctf)
                imgInfo.fn_ctf = fn_ctf1;
            imgInfo.sym     = -1;
            imgInfo.seed    = ROUND(65535 * rnd_unif());

            imgInfo.rot = imgInfo.tilt = imgInfo.psi = 0;

            selfile.getValue(MDL_ANGLE_ROT, imgInfo.rot,objId);
            selfile.getValue(MDL_ANGLE_TILT, imgInfo.tilt,objId);
            selfile.getValue(MDL_ANGLE_PSI, imgInfo.psi,objId);
            //            read_proj.getEulerAngles(imgInfo.rot, imgInfo.tilt, imgInfo.psi);
            EULER_CLIPPING(imgInfo.rot, imgInfo.tilt, imgInfo.psi);

            // Any symmetry?
            if (SL.symsNo() > 0 && !do_not_use_symproj)
            {
                for (int j = 0; j < SL.symsNo(); j++)
                {
                    int sym_index = SYMINDEX(SL, j, i, trueIMG);
                    IMG_Inf[sym_index].fn_proj = imgInfo.fn_proj;
                    IMG_Inf[sym_index].seed   = imgInfo.seed;
                    if (is_ctf_unique)
                        IMG_Inf[sym_index].fn_ctf = fn_ctf;
                    else if (is_there_ctf)
                        IMG_Inf[sym_index].fn_ctf = imgInfo.fn_ctf;
                    IMG_Inf[sym_index].sym = j;
                    SL.getMatrices(j, L, R);
                    L.resize(3, 3); // Erase last row and column
                    R.resize(3, 3); // as only the relative orientation
                    // is useful and not the translation
                    double drot, dtilt, dpsi;
                    Euler_apply_transf(L, R,
                                       imgInfo.rot, imgInfo.tilt, imgInfo.psi,
                                       drot, dtilt, dpsi);
                    IMG_Inf[sym_index].rot = (float)drot;
                    IMG_Inf[sym_index].tilt = (float)dtilt;
                    IMG_Inf[sym_index].psi = (float)dpsi;
                }
            }
        }

        ++i; // I have processed one more image
        if (i % 25 == 0)
            progress_bar(i);
    }
    progress_bar(trueIMG);
}

noSort()

void noSort	(	int	numIMG,
		MultidimArray< int > &	ordered_list
	)

No projection sorting at all. This function directly returns the same order as in the selection file

Definition at line 222 of file recons_misc.cpp.

{
    ordered_list.initLinear(0, numIMG - 1);
}

sortPerpendicular()

void sortPerpendicular	(	int	numIMG,
		ReconsInfo *	IMG_Inf,
		MultidimArray< int > &	ordered_list,
		int	N = 2
	)

Sort projections orthogonally. This function sorts a number of images given by numIMG, whose information about their Euler angles are in IMG_inf, into an ordered list which gives the indexes. First an image is chosen randomly from the whole set. Then all images are compared to the first one, and the most perpendicular one is chosen. The remaining set of images are compared to this two images and the most perpendicular one to the former two is chosen, and so on until no image is left in the set.

If the result in ordered list is 4, 70, 54, 203, 1, 0, ... it means that the first image is the number 4, then goes the 70, then the 54, ...

If N!=-1 then the product is done only with the last N images. A very useful value is N=2

Definition at line 145 of file recons_misc.cpp.

{
    int   i, j;
    MultidimArray<short> chosen(numIMG);     // 1 if that image has been already
    // chosen
    double min_prod;
    int   min_prod_proj=0;
    Matrix2D<double> v(numIMG, 3);
    Matrix2D<double> euler;
    MultidimArray<double> product(numIMG);

    // Initialization
    ordered_list.resize(numIMG);
    for (i = 0; i < numIMG; i++)
    {
        Matrix1D<double> z;
        // Initially no image is chosen
        A1D_ELEM(chosen, i) = 0;

        // Compute the Euler matrix for each image and keep only
        // the third row of each one
        Euler_angles2matrix(IMG_Inf[i].rot, IMG_Inf[i].tilt, 0., euler);
        euler.getRow(2, z);
        v.setRow(i, z);
    }

    // Pick first projection as the first one to be presented
    i = 0;
    A1D_ELEM(chosen, i) = 1;
    A1D_ELEM(ordered_list, 0) = i;

    // Choose the rest of projections
    std::cout << "Sorting projections ...\n";
    init_progress_bar(numIMG - 1);
    Matrix1D<double> rowj, rowi_1, rowi_N_1;
    for (i = 1; i < numIMG; i++)
    {
        // Compute the product of not already chosen vectors with the just
        // chosen one, and select that which has minimum product
        min_prod = MAXFLOAT;
        v.getRow(A1D_ELEM(ordered_list, i - 1),rowi_1);
        if (N != -1 && i > N)
            v.getRow(A1D_ELEM(ordered_list, i - N - 1),rowi_N_1);
        for (j = 0; j < numIMG; j++)
        {
            if (!A1D_ELEM(chosen, j))
            {
                v.getRow(j,rowj);
                A1D_ELEM(product, j) += ABS(dotProduct(rowi_1,rowj));
                if (N != -1 && i > N)
                    A1D_ELEM(product, j) -= ABS(dotProduct(rowi_N_1,rowj));
                if (A1D_ELEM(product, j) < min_prod)
                {
                    min_prod = A1D_ELEM(product, j);
                    min_prod_proj = j;
                }
            }
        }

        // Store the chosen vector and mark it as chosen
        A1D_ELEM(ordered_list, i) = min_prod_proj;
        A1D_ELEM(chosen, min_prod_proj) = 1;

        // The progress bar is updated only every 10 images
        if (i % 10 == 0)
            progress_bar(i);
    }

    // A final call to progress bar to finish a possible small piece
    progress_bar(numIMG - 1);
    std::cout << std::endl;
}

sortRandomly()

void sortRandomly	(	int	numIMG,
		MultidimArray< int > &	ordered_list
	)

Randomize the projections. This function sorts randomly a number of images given by numIMG.

Definition at line 230 of file recons_misc.cpp.

{
    MultidimArray<int> chosen;

    // Initialisation
    ordered_list.resize(numIMG);
    chosen.initZeros(numIMG);

    std::cout << "Randomizing projections ...\n";
    init_progress_bar(numIMG - 1);
    int ptr = 0;
    randomize_random_generator();
    for (int i = numIMG; i > 0; i--)
    {
        // Jump a random number starting at the pointed projection
        int rnd_indx = (int) rnd_unif(0, i) + 1;
        while (rnd_indx > 0)
        {
            // Jump one not chosen image
            ptr = (ptr + 1) % numIMG;
            // Check it is not chosen, if it is, go on skipping
            while (chosen(ptr))
            {
                ptr = (ptr + 1) % numIMG;
            }
            rnd_indx--;
        }

        // Annotate this image
        A1D_ELEM(ordered_list, i - 1) = ptr;
        A1D_ELEM(chosen, ptr) = 1;

        // The progress bar is updated only every 10 images
        if (i % 10 == 0)
            progress_bar(i);
    }

    // A final call to progress bar to finish a possible small piece
    progress_bar(numIMG - 1);
    std::cout << std::endl;
}

updateResidualVector()

void updateResidualVector	(	BasicARTParameters &	prm,
		GridVolume &	vol_basis,
		double &	kappa,
		double &	pow_residual_vol,
		double &	pow_residual_imgs
	)

Update residual vector for WLS ART

Definition at line 275 of file recons_misc.cpp.

{
    GridVolume       residual_vol;
    Projection       read_proj, dummy_proj, new_proj;
    FileName         fn_resi, fn_tmp;
    double           sqrtweight, dim2;
    Matrix2D<double> *A = nullptr;
    std::vector<MultidimArray<double> > newres_imgs;
    MultidimArray<int>    mask;

    residual_vol.resize(vol_basis);
    residual_vol.initZeros();

    // Calculate volume from all backprojected residual images
    std::cout << "Backprojection of residual images " << std::endl;
    if (!(prm.tell&TELL_SHOW_ERROR))
        init_progress_bar(prm.numIMG);

    for (int iact_proj = 0; iact_proj < prm.numIMG ; iact_proj++)
    {
        // backprojection of the weighted residual image
        sqrtweight = sqrt(prm.residual_imgs[iact_proj].weight() / prm.sum_weight);
        read_proj = prm.residual_imgs[iact_proj];
        read_proj() *= sqrtweight;
        dummy_proj().resize(read_proj());

        dummy_proj.setAngles(prm.IMG_Inf[iact_proj].rot,
                              prm.IMG_Inf[iact_proj].tilt,
                              prm.IMG_Inf[iact_proj].psi);

        project_GridVolume(residual_vol, prm.basis, dummy_proj,
                           read_proj, YSIZE(read_proj()), XSIZE(read_proj()),
                           prm.IMG_Inf[iact_proj].rot,
                           prm.IMG_Inf[iact_proj].tilt,
                           prm.IMG_Inf[iact_proj].psi, BACKWARD, prm.eq_mode,
                           prm.GVNeq, nullptr, nullptr, prm.ray_length, prm.threads);

        if (!(prm.tell&TELL_SHOW_ERROR))
            if (iact_proj % XMIPP_MAX(1, prm.numIMG / 60) == 0)
                progress_bar(iact_proj);
    }
    if (!(prm.tell&TELL_SHOW_ERROR))
        progress_bar(prm.numIMG);

    // Convert to voxels: solely for output of power of residual volume
    Image<double>      residual_vox;
    int Xoutput_volume_size = (prm.Xoutput_volume_size == 0) ?
                              prm.projXdim : prm.Xoutput_volume_size;
    int Youtput_volume_size = (prm.Youtput_volume_size == 0) ?
                              prm.projYdim : prm.Youtput_volume_size;
    int Zoutput_volume_size = (prm.Zoutput_volume_size == 0) ?
                              prm.projXdim : prm.Zoutput_volume_size;
    prm.basis.changeToVoxels(residual_vol, &(residual_vox()),
                             Zoutput_volume_size, Youtput_volume_size, Xoutput_volume_size);
    pow_residual_vol = residual_vox().sum2() / MULTIDIM_SIZE(residual_vox());
    residual_vox.clear();

    std::cout << "Projection of residual volume; kappa = " << kappa << std::endl;
    if (!(prm.tell&TELL_SHOW_ERROR))
        init_progress_bar(prm.numIMG);

    // Now that we have the residual volume: project in all directions
    pow_residual_imgs = 0.;
    new_proj().resize(read_proj());
    mask.resize(read_proj());
    BinaryCircularMask(mask, YSIZE(read_proj()) / 2, INNER_MASK);

    dim2 = (double)YSIZE(read_proj()) * XSIZE(read_proj());
    for (int iact_proj = 0; iact_proj < prm.numIMG ; iact_proj++)
    {
        project_GridVolume(residual_vol, prm.basis, new_proj,
                           dummy_proj, YSIZE(read_proj()), XSIZE(read_proj()),
                           prm.IMG_Inf[iact_proj].rot,
                           prm.IMG_Inf[iact_proj].tilt,
                           prm.IMG_Inf[iact_proj].psi, FORWARD, prm.eq_mode,
                           prm.GVNeq, A, nullptr, prm.ray_length, prm.threads);

        sqrtweight = sqrt(prm.residual_imgs[iact_proj].weight() / prm.sum_weight);

        // Next lines like normalization in [EHL] (2.18)?
        FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(new_proj())
        {
            dAij(dummy_proj(), i, j) = XMIPP_MAX(1., dAij(dummy_proj(), i, j)); // to avoid division by zero
            dAij(new_proj(), i, j) /= dAij(dummy_proj(), i, j);
        }
        new_proj() *= sqrtweight * kappa;

        /*
        fn_tmp="residual_"+integerToString(iact_proj);
        dummy_proj()=1000*prm.residual_imgs[iact_proj]();
        dummy_proj.write(fn_tmp+".old");
        */

        prm.residual_imgs[iact_proj]() -= new_proj();
        pow_residual_imgs += prm.residual_imgs[iact_proj]().sum2();

        // Mask out edges of the images
        apply_binary_mask(mask, prm.residual_imgs[iact_proj](), prm.residual_imgs[iact_proj](), 0.);

        /*
        dummy_proj()=1000*new_proj();
        dummy_proj.write(fn_tmp+".change");
        dummy_proj()=1000*prm.residual_imgs[iact_proj]();
        dummy_proj.write(fn_tmp+".new");
        */

        if (!(prm.tell&TELL_SHOW_ERROR))
            if (iact_proj % XMIPP_MAX(1, prm.numIMG / 60) == 0)
                progress_bar(iact_proj);
    }

    pow_residual_imgs /= dim2;
    newres_imgs.clear();

    if (!(prm.tell&TELL_SHOW_ERROR))
        progress_bar(prm.numIMG);

}