project.cpp File Reference

#include "project.h"
#include "directions.h"
#include "project_real_shears.h"
#include <data/fourier_projection.h>
#include "core/xmipp_image_generic.h"
#include "core/metadata_sql.h"
#include <core/args.h>

Include dependency graph for project.cpp:

Go to the source code of this file.

Macros
#define	Nrot   prm.rot_range.samples
#define	Ntilt   prm.tilt_range.samples
#define	Npsi   prm.psi_range.samples
#define	proj_number(base, irot, itilt, ipsi)   base+irot*Ntilt*Npsi+itilt*Npsi+ipsi

Functions
int	translate_randomness (char *str)
void	generate_angles (int ExtProjs, const Angle_range &range, MetaDataVec &DF, char ang_name, const ParametersProjection &prm)
void	generate_even_angles (int ExtProjs, int Nrottilt, MetaDataVec &DF, const ParametersProjection &prm)
int	count_even_angles (const ParametersProjection &prm)
int	Assign_angles (MetaDataVec &DF, const ParametersProjection &prm, const FileName &fn_sym)
int	PROJECT_Effectively_project (const FileName &fnOut, bool singleProjection, projectionType projType, double sampling_rate, const ParametersProjection &prm, PROJECT_Side_Info &side, const Crystal_Projection_Parameters &prm_crystal, Projection &proj, MetaData &SF)
int	ROUT_project (ProgProject &prm, Projection &proj, MetaData &SF)

Macro Definition Documentation

Npsi

#define Npsi   prm.psi_range.samples

Definition at line 558 of file project.cpp.

Nrot

#define Nrot   prm.rot_range.samples

Definition at line 556 of file project.cpp.

Ntilt

#define Ntilt   prm.tilt_range.samples

Definition at line 557 of file project.cpp.

proj_number

#define proj_number	(		base,
			irot,
			itilt,
			ipsi
	)		base+irot*Ntilt*Npsi+itilt*Npsi+ipsi

Definition at line 559 of file project.cpp.

Function Documentation

Assign_angles()

int Assign_angles	(	MetaDataVec &	DF,
		const ParametersProjection &	prm,
		const FileName &	fn_sym
	)

Definition at line 801 of file project.cpp.

{
    int ExtProjs = 0, IntProjs = 0;    // External and internal projections
    int Nrottilt=0;                    // Number of evenly distributed

    // External generation mode
    if (prm.fn_angle != "NULL")
    {
        DF.read(prm.fn_angle);
        ExtProjs = DF.size();
    }

    // Internal generation mode
    if (prm.enable_angle_range)
    {
        randomize_random_generator();
        if (prm.rot_range.randomness != ANGLE_EVENLY)
            IntProjs = Nrot * Ntilt * Npsi;
        else
        {
            if (fn_sym == "")
            {
                Nrottilt = count_even_angles(prm);
                IntProjs = Nrottilt * Npsi;
            }
            else
                IntProjs = 0;
        }
        if (prm.rot_range.randomness != ANGLE_EVENLY)
        {
            generate_angles(ExtProjs, prm.rot_range,  DF, 'r', prm);
            generate_angles(ExtProjs, prm.tilt_range, DF, 't', prm);
            generate_angles(ExtProjs, prm.psi_range,  DF, 'p', prm);
            DF.removeLabel(MDL_ORDER);
        }
        else
        {
            if (fn_sym == "")
                generate_even_angles(ExtProjs, Nrottilt, DF, prm);
            else
            {
                SymList SL;
                if (fn_sym != "")
                    SL.readSymmetryFile(fn_sym);
                std::vector<double> rotList, tiltList;
                double rot_step_at_equator = (prm.rot_range.angF - prm.rot_range.ang0) /
                                             (double)(Nrot - 1);
                make_even_distribution(rotList, tiltList, rot_step_at_equator, SL, true);
                size_t NN=rotList.size();
                for (size_t n=0; n<NN; ++n)
                {
                    if (tiltList[n]>=prm.tilt_range.ang0 && tiltList[n]<=prm.tilt_range.angF &&
                        rotList[n]>=prm.rot_range.ang0 && rotList[n]<=prm.rot_range.angF)
                    {
                        size_t DFid=DF.addObject();
                        DF.setValue(MDL_ANGLE_ROT,rotList[n],DFid);
                        DF.setValue(MDL_ANGLE_TILT,tiltList[n],DFid);
                        DF.setValue(MDL_ANGLE_PSI,0.0,DFid);
                    }
                }
            }
        }
    }

    // Exit
    return ExtProjs + IntProjs;
}

count_even_angles()

int count_even_angles

(

const ParametersProjection &

prm

)

Definition at line 775 of file project.cpp.

{
    int N = 0;
    int limit = prm.tilt_range.samples;
    double rot_step_at_equator = (prm.rot_range.angF - prm.rot_range.ang0) /
                                 (double)(Nrot - 1);
    for (int i = 0; i < limit; i++)
    {
        double tilt = prm.tilt_range.ang0 +
                      (prm.tilt_range.angF - prm.tilt_range.ang0) /
                      (double)(Ntilt - 1) * i;
        double rot_step;
        if (tilt != 0 && tilt != 180)
            rot_step = rot_step_at_equator / sin(DEG2RAD(tilt));
        else
            rot_step = prm.rot_range.angF - prm.rot_range.ang0 + 1;
        int iterRot = static_cast<int>((prm.rot_range.angF - prm.rot_range.ang0) / rot_step);
        N+=iterRot+1;
    }
    N++; // This shouldn't be necessary but some GCC optimization
    // sometimes doesn't do well its work. For instance if we
    // add std::cout << N after N++ in the loop, then it works perfectly
    return N;
}

generate_angles()

void generate_angles	(	int	ExtProjs,
		const Angle_range &	range,
		MetaDataVec &	DF,
		char	ang_name,
		const ParametersProjection &	prm
	)

Definition at line 560 of file project.cpp.

{
    double ang;
    int   N1=0, N2=0;
    int   i, j, k;
    size_t iproj, idx=0;
    int   limit=0;
    MetaDataVec DFaux=DF;

    // Select loop limit ....................................................
    switch (range.randomness)
    {
    case ANGLE_RANGE_DETERMINISTIC:
        limit = range.samples;
        break;
    case ANGLE_RANGE_RANDOM_GROUPS:
        limit = range.samples;
        break;
    case ANGLE_RANGE_RANDOM       :
        limit = Nrot * Ntilt * Npsi;
        break;
    }

    // Which column to write in the document file ...........................
    switch (ang_name)
    {
    case 'r':
        idx = 0;
        break;
    case 't':
        idx = 1;
        break;
    case 'p':
        idx = 2;
        break;
    }

    double unif_min = cos(DEG2RAD(range.angF));
    double unif_max = cos(DEG2RAD(range.ang0));
    for (i = 0; i < limit; i++)
    {
        // Select angle .....................................................
        if (range.randomness == ANGLE_RANGE_DETERMINISTIC)
        {
            if (range.samples > 1)
                ang = range.ang0 +
                      (range.angF - range.ang0) / (double)(range.samples - 1) * i;
            else
                ang = range.ang0;
        }
        else
        {
            switch (ang_name)
            {
            case 'r':
            case 'p':
                ang = rnd_unif(range.ang0, range.angF);
                break;
            case 't':
                ang = RAD2DEG(acos(rnd_unif(unif_min, unif_max)));
                break;
            }
        }

        // Copy this angle to those projections belonging to this group .....
        // If there is any group
        if (range.randomness != ANGLE_RANGE_RANDOM)
        {
            switch (ang_name)
            {
            case 'r':
                N1 = Ntilt;
                N2 = Npsi;
                break;
            case 't':
                N1 = Nrot;
                N2 = Npsi;
                break;
            case 'p':
                N1 = Nrot;
                N2 = Ntilt;
                break;
            }
            for (j = 0; j < N1; j++)
                for (k = 0; k < N2; k++)
                {
                    switch (ang_name)
                    {
                    case 'r':
                        iproj = proj_number(ExtProjs, i, j, k);
                        break;
                    case 't':
                        iproj = proj_number(ExtProjs, j, i, k);
                        break;
                    case 'p':
                        iproj = proj_number(ExtProjs, j, k, i);
                        break;
                    }
                    size_t idx_tmp=DFaux.firstObject(MDValueEQ(MDL_ORDER,iproj));
                    if (idx_tmp==BAD_OBJID)
                    {
                        idx_tmp=DFaux.addObject();
                        DFaux.setValue(MDL_ORDER,iproj,idx_tmp);
                    }
                    switch (idx)
                    {
                    case 0:
                        DFaux.setValue(MDL_ANGLE_ROT,ang,idx_tmp);
                        break;
                    case 1:
                        DFaux.setValue(MDL_ANGLE_TILT,ang,idx_tmp);
                        break;
                    case 2:
                        DFaux.setValue(MDL_ANGLE_PSI,ang,idx_tmp);
                        break;
                    }
                }
        }
        else
        {
            size_t iproj=ExtProjs + i;
            size_t dfidx=DFaux.firstObject(MDValueEQ(MDL_ORDER,iproj));
            if (dfidx==BAD_OBJID)
            {
                dfidx=DFaux.addObject();
                DFaux.setValue(MDL_ORDER,iproj,dfidx);
            }
            switch (idx)
            {
            case 0:
                DFaux.setValue(MDL_ANGLE_ROT,ang,dfidx);
                break;
            case 1:
                DFaux.setValue(MDL_ANGLE_TILT,ang,dfidx);
                break;
            case 2:
                DFaux.setValue(MDL_ANGLE_PSI,ang,dfidx);
                break;
            }
        }
    }
    DF.sort(DFaux,MDL_ORDER);
}

generate_even_angles()

void generate_even_angles	(	int	ExtProjs,
		int	Nrottilt,
		MetaDataVec &	DF,
		const ParametersProjection &	prm
	)

Definition at line 706 of file project.cpp.

{
    // We will run over the tilt angle in a deterministic way
    // then for every tilt angle, a rot_step is computed so that
    // it keeps the same distance in the circle generated by tilt
    // as the sample distance at the equator (tilt=90).
    MetaDataVec DFaux;
    int N = 0;
    int limit = prm.tilt_range.samples;
    double rot_step_at_equator = (prm.rot_range.angF - prm.rot_range.ang0) /
                                 (double)(Nrot - 1);
    for (int i = 0; i < limit; i++)
    {
        // Compute the corresponding deterministic tilt
        double tilt = prm.tilt_range.ang0 +
                      (prm.tilt_range.angF - prm.tilt_range.ang0) /
                      (double)(Ntilt - 1) * i;
        // Now compute the corresponding rotational angles
        double rot_step;
        if (tilt != 0 && tilt != 180)
            rot_step = rot_step_at_equator / sin(DEG2RAD(tilt));
        else
            rot_step = prm.rot_range.angF - prm.rot_range.ang0 + 1;
        int iterRot = static_cast<int>((prm.rot_range.angF - prm.rot_range.ang0) / rot_step);
        for (int irot = 0; irot <= iterRot; irot ++)
        {
            double rot = prm.rot_range.ang0+irot*rot_step;
            // Copy this angle to those projections belonging to this group .....
            // If there is any group
            for (int k = 0; k < Npsi; k++)
            {
                // Select psi
                double psi;
                if (prm.psi_range.randomness == ANGLE_RANGE_DETERMINISTIC)
                {
                    if (prm.psi_range.samples > 1)
                    {
                        psi = prm.psi_range.ang0 +
                              (prm.psi_range.angF - prm.psi_range.ang0) /
                              (double)(prm.psi_range.samples - 1) * k;
                    }
                    else
                        psi = prm.psi_range.ang0;
                }
                else
                    psi = prm.psi_range.ang0 +
                          (prm.psi_range.angF - prm.psi_range.ang0) /
                          (double)(Npsi - 1) * k;

                size_t iproj = ExtProjs + N + Nrottilt * k;
                size_t idx_tmp=DFaux.firstObject(MDValueEQ(MDL_ORDER,iproj));
                if (idx_tmp==BAD_OBJID)
                {
                    idx_tmp=DFaux.addObject();
                    DFaux.setValue(MDL_ORDER,iproj,idx_tmp);
                }
                DFaux.setValue(MDL_ANGLE_ROT,rot,idx_tmp);
                DFaux.setValue(MDL_ANGLE_TILT,tilt,idx_tmp);
                DFaux.setValue(MDL_ANGLE_PSI,psi,idx_tmp);
            }
            N++;
        }
    }
    DF.sort(DFaux,MDL_ORDER);
    DF.removeLabel(MDL_ORDER);
}

PROJECT_Effectively_project()

int PROJECT_Effectively_project	(	const FileName &	fnOut,
		bool	singleProjection,
		projectionType	projType,
		double	sampling_rate,
		const ParametersProjection &	prm,
		PROJECT_Side_Info &	side,
		const Crystal_Projection_Parameters &	prm_crystal,
		Projection &	proj,
		MetaData &	SF
	)

1 .5 NEAREST

Definition at line 939 of file project.cpp.

{
    int NumProjs = 0;
    SF.clear();
    if (!fnOut.isInStack())
        fnOut.deleteFile();
    std::cerr << "Projecting ...\n";
    init_progress_bar(side.DF.size());
    SF.setComment("First set of angles=actual angles; Second set of angles=noisy angles");

    //#define DEBUG
#ifdef DEBUG

    MetaDataVec mdShifts;
    MetaDataVec mdRotations;
#endif

    bool existFlip = false;
    bool hasCTF=false;

    int projIdx=FIRST_IMAGE;
    FileName fn_proj;              // Projection name
    RealShearsInfo *Vshears=nullptr;
    FourierProjector *Vfourier=nullptr;
    if (projType == SHEARS && side.phantomMode==PROJECT_Side_Info::VOXEL)
        Vshears=new RealShearsInfo(side.phantomVol());
    if (projType == FOURIER && side.phantomMode==PROJECT_Side_Info::VOXEL)
        Vfourier=new FourierProjector(side.phantomVol(),side.paddFactor,side.maxFrequency,side.BSplineDeg);
    fn_proj=fnOut;
    if (side.doCrystal)
        createEmptyFile(fn_proj, prm_crystal.crystal_Xdim, prm_crystal.crystal_Ydim,
                        1, side.DF.size(), true,WRITE_OVERWRITE);
    else
        createEmptyFile(fn_proj, prm.proj_Xdim, prm.proj_Ydim,
                        1, side.DF.size(), true,WRITE_OVERWRITE);

    if (side.DF.containsLabel(MDL_FLIP))
        existFlip = true;


    for (size_t objId : side.DF.ids())
    {
        size_t DFmov_objId=SF.addObject();
        if (!singleProjection)
            fn_proj.compose(projIdx,fnOut);
        SF.setValue(MDL_IMAGE,fn_proj,DFmov_objId);
        SF.setValue(MDL_ENABLED,1,DFmov_objId);

        // Choose angles .....................................................
        double rot, tilt, psi, x=0, y=0;    // Actual projecting angles
        bool flip=false;
        side.DF.getValue(MDL_ANGLE_ROT,rot,objId);
        side.DF.getValue(MDL_ANGLE_TILT,tilt,objId);
        side.DF.getValue(MDL_ANGLE_PSI,psi,objId);

        if (prm.applyShift)
        {
            if (side.DF.containsLabel(MDL_SHIFT_X))
                side.DF.getValue(MDL_SHIFT_X,x,objId);
            if (side.DF.containsLabel(MDL_SHIFT_Y))
                side.DF.getValue(MDL_SHIFT_Y,y,objId);
        }
        realWRAP(rot, 0, 360);
        realWRAP(tilt, 0, 360);
        realWRAP(psi, 0, 360);

        if (existFlip)
            side.DF.getValue(MDL_FLIP,flip,objId);


        SF.setValue(MDL_ANGLE_ROT,rot,DFmov_objId);
        SF.setValue(MDL_ANGLE_TILT,tilt,DFmov_objId);
        SF.setValue(MDL_ANGLE_PSI,psi,DFmov_objId);

        if ((NumProjs % XMIPP_MAX(1, side.DF.size() / 60)) == 0)
            progress_bar(NumProjs);

        // Choose Center displacement ........................................
        double shiftX = rnd_gaus(prm.Ncenter_avg, prm.Ncenter_dev)+x;
        double shiftY = rnd_gaus(prm.Ncenter_avg, prm.Ncenter_dev)+y;
        SF.setValue(MDL_SHIFT_X,shiftX,DFmov_objId);
        SF.setValue(MDL_SHIFT_Y,shiftY,DFmov_objId);
        //ROB: totally anoying behaviour in xmipp
        shiftX=-shiftX;
        shiftY=-shiftY;

        //If there is CTF information use it
        CTFDescription ctf;
        if ( (side.DF.containsLabel(MDL_CTF_DEFOCUSU) || side.DF.containsLabel(MDL_CTF_MODEL)) && (prm.doCTFCorrection) )
        {
            //JV he tenido que tocar esta funcion para poder acceder al sampling rate
            hasCTF = true;
            MDRowVec row;
            side.DF.getRow(row,objId);
            ctf.readFromMdRow(row);
            ctf.produceSideInfo();
        }

#ifdef DEBUG

        mdShifts.addObject();
        mdShifts.setValue(MDL_IMAGE,fn_proj,DFmov_objId);
        mdShifts.setValue(MDL_SHIFT_X,-shiftX,DFmov_objId);
        mdShifts.setValue(MDL_SHIFT_Y,-shiftY,DFmov_objId);
        mdShifts.setValue(MDL_ENABLED,1,DFmov_objId);


        mdRotations.addObject();
        mdRotations.setValue(MDL_IMAGE,fn_proj,DFmov_objId);
        mdRotations.setValue(MDL_ANGLE_ROT,-rot,DFmov_objId);
        mdRotations.setValue(MDL_ANGLE_TILT,-tilt,DFmov_objId);
        mdRotations.setValue(MDL_ANGLE_PSI,-psi,DFmov_objId);


#endif
        // Really project ....................................................
        if (side.phantomMode==PROJECT_Side_Info::VOXEL)
        {
            if (projType == SHEARS)
                projectVolume(*Vshears, proj, prm.proj_Ydim, prm.proj_Xdim,
                              rot, tilt, psi);
            else if (projType == FOURIER)
                projectVolume(*Vfourier, proj, prm.proj_Ydim, prm.proj_Xdim,
                              rot, tilt, psi);
            else if (projType == REALSPACE)
                projectVolume(side.phantomVol(), proj, prm.proj_Ydim, prm.proj_Xdim,
                              rot, tilt, psi);

            if (hasCTF)
                ctf.applyCTF(proj(),sampling_rate, prm.doPhaseFlip);

            hasCTF = false;

            Matrix1D<double> shifts(2);
            XX(shifts) = shiftX;
            YY(shifts) = shiftY;
            selfTranslate(xmipp_transformation::LINEAR,IMGMATRIX(proj), shifts);
        }
        else if (side.phantomMode==PROJECT_Side_Info::PDB)
        {
            PDBPhantom aux;
            aux=side.phantomPDB;
            aux.shift(shiftX,shiftY,0);
            projectPDB(side.phantomPDB, side.interpolator,
                       proj, prm.proj_Ydim, prm.proj_Xdim,
                       rot, tilt, psi);
        }
        else if (side.phantomMode==PROJECT_Side_Info::XMIPP)
        {
            Phantom aux;
            aux = side.phantomDescr;
            aux.shift(shiftX, shiftY, 0);
            if (prm_crystal.crystal_Xdim == 0)
            {
                // Project a single mathematical volume
                aux.project_to(proj,
                               prm.proj_Ydim, prm.proj_Xdim,
                               rot, tilt, psi);
            }
            else
            { // Project mathematical volume as a crystal
                project_crystal(aux, proj, prm, side, prm_crystal, rot, tilt, psi);
            }
        }
        if (flip)
            proj().selfReverseX();

        // Add noise in angles and voxels ....................................
        SF.setValue(MDL_ANGLE_ROT2,rot,DFmov_objId);
        SF.setValue(MDL_ANGLE_TILT2,tilt,DFmov_objId);
        SF.setValue(MDL_ANGLE_PSI2,psi,DFmov_objId);
        rot  += rnd_gaus(prm.rot_range.Navg,  prm.rot_range.Ndev);
        tilt += rnd_gaus(prm.tilt_range.Navg, prm.tilt_range.Ndev);
        psi  += rnd_gaus(prm.psi_range.Navg,  prm.psi_range.Ndev);


        SF.setValue(MDL_ANGLE_ROT,realWRAP(rot, 0, 360),DFmov_objId);
        SF.setValue(MDL_ANGLE_TILT,realWRAP(tilt, 0, 360),DFmov_objId);
        SF.setValue(MDL_ANGLE_PSI,realWRAP(psi, 0, 360),DFmov_objId);

        IMGMATRIX(proj).addNoise(prm.Npixel_avg, prm.Npixel_dev, "gaussian");

        // Save ..............................................................
        if (singleProjection)
            proj.write(fn_proj);
        else
        {
            proj.write(fn_proj,projIdx,true,WRITE_OVERWRITE);
        }
        projIdx++;
        NumProjs++;
    }
    progress_bar(side.DF.size());

#ifdef DEBUG

    mdShifts.write("shifts.xmd");
    mdRotations.write("rotations.xmd");
#endif
#undef DEBUG
    // release memory
    delete Vshears;
    delete Vfourier;

    return NumProjs;
}

translate_randomness()

int translate_randomness

(

char *

str

)

Definition at line 189 of file project.cpp.

{
    if (str == nullptr)
        return ANGLE_RANGE_DETERMINISTIC;
    if (strcmp(str, "random_group") == 0)
        return ANGLE_RANGE_RANDOM_GROUPS;
    if (strcmp(str, "random") == 0)
        return ANGLE_RANGE_RANDOM;
    if (strcmp(str, "even") == 0)
        return ANGLE_EVENLY;
    REPORT_ERROR(ERR_PARAM_INCORRECT,
                 (String)"Prog_Project_Parameters::read: Not recognized randomness: "
                 + str);
}