ProgAngularAccuracyPCA Class Reference

#include <angular_accuracy_pca.h>

Inheritance diagram for ProgAngularAccuracyPCA:

Collaboration diagram for ProgAngularAccuracyPCA:

Public Member Functions
void	readParams ()
void	defineParams ()
void	run ()
void	obtainPCAs (MetaData &SF, size_t numPCAs)
	ProgAngularAccuracyPCA ()
virtual void	gatherResults ()
	Gather alignment. More...
virtual void	synchronize ()
	Synchronize with other processors. More...
Public Member Functions inherited from XmippProgram
const char *	getParam (const char *param, int arg=0)
const char *	getParam (const char *param, const char *subparam, int arg=0)
int	getIntParam (const char *param, int arg=0)
int	getIntParam (const char *param, const char *subparam, int arg=0)
double	getDoubleParam (const char *param, int arg=0)
double	getDoubleParam (const char *param, const char *subparam, int arg=0)
float	getFloatParam (const char *param, int arg=0)
float	getFloatParam (const char *param, const char *subparam, int arg=0)
void	getListParam (const char *param, StringVector &list)
int	getCountParam (const char *param)
bool	checkParam (const char *param)
bool	existsParam (const char *param)
void	addParamsLine (const String &line)
void	addParamsLine (const char *line)
ParamDef *	getParamDef (const char *param) const
virtual void	quit (int exit_code=0) const
virtual int	tryRun ()
void	initProgress (size_t total, size_t stepBin=60)
void	setProgress (size_t value=0)
void	endProgress ()
void	processDefaultComment (const char *param, const char *left)
void	setDefaultComment (const char *param, const char *comment)
virtual void	initComments ()
void	setProgramName (const char *name)
void	addUsageLine (const char *line, bool verbatim=false)
void	clearUsage ()
void	addExampleLine (const char *example, bool verbatim=true)
void	addSeeAlsoLine (const char *seeAlso)
void	addKeywords (const char *keywords)
const char *	name () const
virtual void	usage (int verb=0) const
virtual void	usage (const String &param, int verb=2)
int	version () const
virtual void	show () const
virtual void	read (int argc, const char **argv, bool reportErrors=true)
virtual void	read (int argc, char **argv, bool reportErrors=true)
void	read (const String &argumentsLine)
	XmippProgram ()
	XmippProgram (int argc, const char **argv)
virtual	~XmippProgram ()

Public Attributes
FileName	fnPhantom
FileName	fnNeighbours
FileName	fnOut
FileName	fnOutQ
Image< double >	phantomVol
MetaDataDb	mdPartial
size_t	rank
size_t	Nprocessors
PCAMahalanobisAnalyzer	pca
int	newXdim
int	newYdim
Public Attributes inherited from XmippProgram
bool	doRun
bool	runWithoutArgs
int	verbose
	Verbosity level. More...
int	debug

Additional Inherited Members
Protected Member Functions inherited from XmippProgram
void	defineCommons ()
Protected Attributes inherited from XmippProgram
int	errorCode
ProgramDef *	progDef
	Program definition and arguments parser. More...
std::map< String, CommentList >	defaultComments
int	argc
	Original command line arguments. More...
const char **	argv

Detailed Description

Definition at line 39 of file angular_accuracy_pca.h.

Constructor & Destructor Documentation

ProgAngularAccuracyPCA()

ProgAngularAccuracyPCA::ProgAngularAccuracyPCA

(

)

Definition at line 32 of file angular_accuracy_pca.cpp.

{
    rank=0;
    Nprocessors=1;
}

Member Function Documentation

defineParams()

void ProgAngularAccuracyPCA::defineParams ( )

virtual

Function in which the param of each Program are defined.

Reimplemented from XmippProgram.

Definition at line 51 of file angular_accuracy_pca.cpp.

{
    addUsageLine("Determine the angular determination accuracy of a set of particles and a 3D reconstruction ");
    addParamsLine("  [ -i <volume_file> ]       : Voxel volume");
    addParamsLine("  [--i2 <md_file=\"\">]      : Metadata file with neighbour projections");
    addParamsLine("  [ -o <md_file=\"\">]       : Metadata file with obtained weights");
    addParamsLine("  [--dim <d=-1>]             : Scale images to this size if they are larger.");
    addParamsLine("                             : Set to -1 for no rescaling");
}

gatherResults()

virtual void ProgAngularAccuracyPCA::gatherResults ( )

inlinevirtual

Gather alignment.

Reimplemented in MpiProgAngularAccuracyPCA.

Definition at line 72 of file angular_accuracy_pca.h.

{}

obtainPCAs()

void ProgAngularAccuracyPCA::obtainPCAs	(	MetaData &	SF,
		size_t	numPCAs
	)

Definition at line 182 of file angular_accuracy_pca.cpp.

{
    size_t numIter = 200;
    pca.clear();
    size_t imgno;
    Image<double> img;
    double rot, tilt, psi;
    double shiftX, shiftY;
    bool mirror;
    size_t  Xdim, Ydim, Zdim, Ndim;
    phantomVol().getDimensions(Xdim,Ydim,Zdim,Ndim);

    if (  newXdim == -1 )
    {
        newXdim = Xdim;
        newYdim = Ydim;
    }

    Matrix2D<double> proj;;
    imgno = 0;
    Projection P;
    FileName image;
    MultidimArray<float> temp;
    MultidimArray<double> avg;
    Matrix2D<double> E, Trans(3,3);
    Matrix1D<double> opt_offsets(2);

#ifdef DEBUG
    for (size_t objId : SF.ids())
    {
        int enabled;
        SF.getValue(MDL_ENABLED, enabled, objId);
        if ( (enabled==-1)  )
        {
            imgno++;
            continue;
        }

        SF.getValue(MDL_ANGLE_ROT, rot, objId);
        SF.getValue(MDL_ANGLE_TILT, tilt, objId);
        SF.getValue(MDL_ANGLE_PSI, psi, objId);
        SF.getValue(MDL_FLIP, mirror, objId);

        if (mirror)
        {
            double newrot;
            double newtilt;
            double newpsi;
            Euler_mirrorY(rot,tilt,psi,newrot,newtilt,newpsi);
            rot = newrot;
            tilt = newtilt;
            psi = newpsi;
        }

        projectVolume(phantomVol(), P, Ydim, Xdim, rot, tilt, psi);

        Euler_angles2matrix(rot, tilt, psi, E, false);
        double angle = atan2(MAT_ELEM(E,0,1),MAT_ELEM(E,0,0));
        selfRotate(LINEAR, P(),-(angle*180)/3.14159 , WRAP);
        typeCast(P(), temp);
        selfScaleToSize(LINEAR,temp,newXdim,newYdim,1);
        temp.resize(newXdim*newYdim);
        temp.statisticsAdjust(0.f,1.f);
        pca.addVector(temp);
        imgno++;

        #ifdef DEBUG
        {
            {
                size_t val;
                SF.getValue(MDL_ITEM_ID,val, objId);
                {
                    std::cout << E << std::endl;
                    std::cout << (angle*180)/3.14159 << std::endl;
                    P.write("kk_proj.tif");
                    SF.getValue(MDL_ANGLE_PSI,psi, objId);
                    std::cout << rot << " " << tilt << " " << psi << std::endl;
                    char c;
                    std::getchar();
                }
            }
        }
        #endif

    }
    pca.subtractAvg();
    avg = pca.avg;
    pca.learnPCABasis(numPCAs,numIter);
    //pca.projectOnPCABasis(projRef);
    pca.v.clear();

#endif

    imgno = 0;
    FileName f;

    for (size_t objId : SF.ids())
    {
        int enabled;
        SF.getValue(MDL_ENABLED, enabled, objId);
        SF.getValue(MDL_SHIFT_X, shiftX, objId);
        SF.getValue(MDL_SHIFT_Y, shiftY, objId);
        SF.getValue(MDL_FLIP, mirror, objId);

        if ( enabled==-1  )
        {
            imgno++;
            continue;
        }

        //geo2TransformationMatrix(input, Trans);
        //ApplyGeoParams params;
        //params.only_apply_shifts = true;
        //img.readApplyGeo(SF,objId,params);

        SF.getValue(MDL_IMAGE, f, objId);
        img.read(f);
        SF.getValue(MDL_ANGLE_ROT, rot, objId);
        SF.getValue(MDL_ANGLE_TILT, tilt, objId);
        SF.getValue(MDL_ANGLE_PSI, psi, objId);

        if (mirror)
        {
            double newrot;
            double newtilt;
            double newpsi;
            Euler_mirrorY(rot,tilt,psi,newrot,newtilt,newpsi);
            rot = newrot;
            tilt = newtilt;
            psi = newpsi;
        }

        Euler_angles2matrix(rot, tilt, psi, E, false);
        double angle = atan2(MAT_ELEM(E,0,1),MAT_ELEM(E,0,0));
        angle=-(angle*180)/3.14159;
        rotation2DMatrix(angle, Trans, true);
        dMij(Trans, 0, 2) = shiftX;
        dMij(Trans, 1, 2) = shiftY;

        selfApplyGeometry(xmipp_transformation::BSPLINE3, img(), Trans, xmipp_transformation::IS_NOT_INV, xmipp_transformation::WRAP);

#ifdef DEBUG
        {
            std::cout <<  MAT_ELEM(Trans,0,0) << " " << MAT_ELEM(Trans,0,1) << " " << MAT_ELEM(Trans,0,2) << " " << MAT_ELEM(Trans,1,0) << " " <<  MAT_ELEM(Trans,1,1) << " " << MAT_ELEM(Trans,1,2) << " " << MAT_ELEM(Trans,2,0) << " " << MAT_ELEM(Trans,2,1) << " " << MAT_ELEM(Trans,2,2) << " " << std::endl;
            std::cout <<  MAT_ELEM(E,0,0) << " " << MAT_ELEM(E,0,1) << " " << MAT_ELEM(E,0,2) << " " << MAT_ELEM(E,1,0) << " " <<  MAT_ELEM(E,1,1) << " " << MAT_ELEM(E,1,2) << " " << MAT_ELEM(E,2,0) << " " << MAT_ELEM(E,2,1) << " " << MAT_ELEM(E,2,2) << " " << std::endl;
            size_t val;
            SF.getValue(MDL_ITEM_ID, val, objId);
            if (true)
            {
                SF.getValue(MDL_IMAGE, f, objId);
                std::cout << f << std::endl;
                img.write("kk_exp.tif");

                char c;
                std::getchar();

            }
        }
#endif

        typeCast(img(), temp);
        selfScaleToSize(xmipp_transformation::LINEAR,temp,newXdim,newYdim,1);
        temp.resize(newXdim*newYdim);
        temp.statisticsAdjust(0.f,1.f);
        pca.addVector(temp);
        imgno++;
    }

    pca.subtractAvg();
    pca.learnPCABasis(numPCAs,numIter);
    pca.projectOnPCABasis(proj);

    std::vector< MultidimArray<float> > recons(SF.size());
    for(int n=0; n<SF.size(); n++)
        recons[n] = MultidimArray<float>(newXdim*newYdim);

    pca.reconsFromPCA(proj,recons);
    pca.evaluateZScore(numPCAs,numIter, false);

    imgno = 0;
    Image<float> imgRes;
    double R2_Proj,R2_Exp;
    R2_Proj=0;
    R2_Exp=0;

    MultidimArray <int> ROI;
    ROI.resizeNoCopy(newYdim,newXdim);
    ROI.setXmippOrigin();
    FOR_ALL_ELEMENTS_IN_ARRAY2D(ROI)
    {
        double temp = std::sqrt(i*i+j*j);
        if ( temp < (newXdim/2))
            A2D_ELEM(ROI,i,j)= 1;
        else
            A2D_ELEM(ROI,i,j)= 0;
    }

    ROI.resize(newYdim*newXdim);
    ROI.setXmippOrigin();

    for (size_t objId : SF.ids())
    {
        int enabled;
        SF.getValue(MDL_ENABLED, enabled, objId);
        if ( enabled==-1  )
        {
            imgno++;
            continue;
        }

        //Projected Image
        SF.getValue(MDL_ANGLE_ROT, rot, objId);
        SF.getValue(MDL_ANGLE_TILT, tilt, objId);
        SF.getValue(MDL_ANGLE_PSI, psi, objId);
        SF.getValue(MDL_FLIP, mirror, objId);

        if (mirror)
        {
            double newrot;
            double newtilt;
            double newpsi;
            Euler_mirrorY(rot,tilt,psi,newrot,newtilt,newpsi);
            rot = newrot;
            tilt = newtilt;
            psi = newpsi;
        }

        projectVolume(phantomVol(), P, Ydim, Xdim, rot, tilt, psi);
        Euler_angles2matrix(rot, tilt, psi, E, false);
        double angle = atan2(MAT_ELEM(E,0,1),MAT_ELEM(E,0,0));
        angle = -(angle*180)/3.14159;
        selfRotate(xmipp_transformation::LINEAR, P(),angle , xmipp_transformation::WRAP);
        typeCast(P(), temp);
        selfScaleToSize(xmipp_transformation::LINEAR,temp,newXdim,newYdim,1);
        temp.resize(newXdim*newYdim);
        temp.statisticsAdjust(0.f,1.f);
        temp.setXmippOrigin();

        //Reconstructed Image
        recons[imgno].statisticsAdjust(0.f,1.f);
        recons[imgno].resize(newYdim*newXdim);
        recons[imgno].setXmippOrigin();

        R2_Proj = correlationIndex(temp,recons[imgno],&ROI);

        SF.getValue(MDL_SHIFT_X, shiftX, objId);
        SF.getValue(MDL_SHIFT_Y, shiftY, objId);

        SF.getValue(MDL_IMAGE, f, objId);
        img.read(f);

        rotation2DMatrix(angle, Trans, true);
        dMij(Trans, 0, 2) = shiftX;
        dMij(Trans, 1, 2) = shiftY;
        selfApplyGeometry(xmipp_transformation::BSPLINE3, img(), Trans, xmipp_transformation::IS_NOT_INV, xmipp_transformation::WRAP);

        typeCast(img(), temp);
        selfScaleToSize(xmipp_transformation::LINEAR,temp,newXdim,newYdim,1);
        temp.resize(newXdim*newYdim);
        temp.statisticsAdjust(0.f,1.f);
        temp.setXmippOrigin();

        R2_Exp = correlationIndex(temp,recons[imgno],&ROI);

        SF.setValue(MDL_SCORE_BY_PCA_RESIDUAL_PROJ,R2_Proj,objId);
        SF.setValue(MDL_SCORE_BY_PCA_RESIDUAL_EXP,R2_Exp,objId);
        SF.setValue(MDL_SCORE_BY_PCA_RESIDUAL,R2_Proj*R2_Exp,objId);
        SF.setValue(MDL_SCORE_BY_ZSCORE, exp(-A1D_ELEM(pca.Zscore,imgno)),objId);

        #ifdef DEBUG
        {
            size_t val;
            SF.getValue(MDL_ITEM_ID,val,objId);
            if (true)
            {
            Image<float>  imgRecons;
            Image<double> imgAvg;
            SF.getValue(MDL_IMAGE,f,objId);

            img.write("kk_exp0.tif");
            apply_binary_mask(ROI,temp,imgRecons());
            imgRecons().resize(newYdim,newXdim);
            imgRecons.write("kk_exp.tif");

            recons[imgno].statisticsAdjust(0.f,1.f);
            apply_binary_mask(ROI,recons[imgno],imgRecons());
            imgRecons().resize(newYdim,newXdim);
            imgRecons.write("kk_reconstructed.tif");

            imgAvg()=pca.avg;
            imgAvg().resize(newYdim,newXdim);
            imgAvg.write("kk_average.tif");

            //Projected Image
            SF.getValue(MDL_ANGLE_ROT,rot,objId);
            SF.getValue(MDL_ANGLE_TILT,tilt,objId);
            SF.getValue(MDL_ANGLE_PSI,psi,objId);
            SF.getValue(MDL_FLIP,mirror,objId);

            if (mirror)
            {
                double newrot;
                double newtilt;
                double newpsi;
                Euler_mirrorY(rot,tilt,psi,newrot,newtilt,newpsi);
                rot = newrot;
                tilt = newtilt;
                psi = newpsi;
            }

            projectVolume(phantomVol(), P, Ydim, Xdim, rot, tilt, psi);
            Euler_angles2matrix(rot, tilt, psi, E, false);
            double angle = atan2(MAT_ELEM(E,0,1),MAT_ELEM(E,0,0));
            selfRotate(LINEAR, P(),-(angle*180)/3.14159 , WRAP);
            P.write("kk_proj.tif");

            //std::cout <<  exp(-R2_Proj) << " " << exp(-R2_Exp) << std::endl;
            std::cout <<  R2_Proj << " " << R2_Exp <<  " " << R2_Proj*R2_Exp << std::endl;

            for(int i=0; i<numPCAs;i++)
            {
                std::cout << "proj " << MAT_ELEM(proj,i,imgno) <<  "   " <<  "projRef " <<  MAT_ELEM(proj,i,imgno) << std::endl;

            }

            char c;
            std::getchar();

            }

        }

#endif

        imgno++;
    }

    recons.clear();
    img.clear();
    temp.clear();
    imgRes.clear();
    pca.clear();

}

readParams()

void ProgAngularAccuracyPCA::readParams ( )

virtual

Function in which each program will read parameters that it need. If some error occurs the usage will be printed out.

Reimplemented from XmippProgram.

Definition at line 38 of file angular_accuracy_pca.cpp.

{
    fnPhantom = getParam("-i");
    fnNeighbours = getParam("--i2");
    fnOut = getParam("-o");
    fnOutQ = fnOut.getDir()+"/validationAlignabilityAccuracy.xmd";
    newXdim = getIntParam("--dim");
    newYdim = newXdim;

    std::cout <<  newXdim << std::endl;

}

run()

void ProgAngularAccuracyPCA::run ( )

virtual

This function will be start running the program. it also should be implemented by derived classes.

Reimplemented from XmippProgram.

Definition at line 61 of file angular_accuracy_pca.cpp.

{
    MetaDataVec md;
    StringVector blocks;
    getBlocksInMetaDataFile(fnNeighbours, blocks);

    phantomVol.read(fnPhantom);
    phantomVol().setXmippOrigin();

    size_t numPCAs;

    if (rank==0)
        init_progress_bar(blocks.size());

    for (size_t i = 0; i < blocks.size(); ++i)
    {
        if ((i+1)%Nprocessors==rank)
        {
            md.read((String) blocks[i].c_str()+'@'+fnNeighbours);

            if (md.size() <= 1)
                continue;

            else if ( md.size() > 20 )
                numPCAs = 3;
            else if ( (md.size() >= 5) & (md.size() < 20) )
                numPCAs = 2;
            else
                numPCAs = 1;

            obtainPCAs(md,numPCAs);

            for (auto& row : md)
                mdPartial.addRow(dynamic_cast<MDRowVec&>(row));

            if (rank==0)
                progress_bar(i+1);
        }
    }

    synchronize();
    gatherResults();

    if (rank == 0)
    {
        double pcaResidualProj,pcaResidualExp,pcaResidual,Zscore,temp,qResidualProj,qResidualExp,qZscore;

        qResidualProj = 0;
        qResidualExp = 0;
        qZscore = 0;

        String expression;
        size_t maxIdx;
        MDRowSql row;
        MetaDataDb MDSort, tempMd, MDOut, MDOutQ;
        MDSort.sort(mdPartial,MDL_ITEM_ID,true,-1,0);
        MDSort.getValue(MDL_ITEM_ID,maxIdx,MDSort.lastRowId());

        for (size_t i=0; i<=maxIdx;i++)
        {
            expression = formatString("itemId == %lu",i);
            tempMd.importObjects(MDSort, MDExpression(expression));

            if (tempMd.size() <= 0)
                continue;

            pcaResidualProj = -1e3;
            pcaResidualExp = -1e3;
            pcaResidual = -1e3;
            Zscore = -1e3;

            row = tempMd.getRowSql(tempMd.firstRowId());

            for (size_t objId : tempMd.ids())
            {
                tempMd.getValue(MDL_SCORE_BY_PCA_RESIDUAL_PROJ, temp, objId);
                if (temp > pcaResidualProj)
                    pcaResidualProj=temp;

                tempMd.getValue(MDL_SCORE_BY_PCA_RESIDUAL_EXP, temp, objId);
                if (temp > pcaResidualExp)
                    pcaResidualExp=temp;

                tempMd.getValue(MDL_SCORE_BY_PCA_RESIDUAL, temp, objId);
                if (temp > pcaResidual)
                    pcaResidual=temp;

                tempMd.getValue(MDL_SCORE_BY_ZSCORE, temp, objId);
                if (temp > Zscore)
                    Zscore=temp;
            }

            qResidualProj += pcaResidualProj;
            qResidualExp  += pcaResidualExp;
            qZscore       += Zscore;

            row.setValue(MDL_SCORE_BY_PCA_RESIDUAL_PROJ,pcaResidualProj);
            row.setValue(MDL_SCORE_BY_PCA_RESIDUAL_EXP,pcaResidualExp);
            row.setValue(MDL_SCORE_BY_PCA_RESIDUAL,pcaResidual);
            row.setValue(MDL_SCORE_BY_ZSCORE,Zscore);
            MDOut.addRow(row);
            row.clear();
        }

        MDOut.write(fnOut);

        qResidualProj /= MDOut.size();
        qResidualExp /= MDOut.size();
        qZscore /= MDOut.size();

        row.setValue(MDL_IMAGE,fnPhantom);
        row.setValue(MDL_SCORE_BY_PCA_RESIDUAL_PROJ,qResidualProj);
        row.setValue(MDL_SCORE_BY_PCA_RESIDUAL_EXP,qResidualExp);
        row.setValue(MDL_SCORE_BY_ZSCORE,qZscore);

        MDOutQ.addRow(row);
        MDOutQ.write(fnOutQ);
        progress_bar(blocks.size());
    }
}

synchronize()

virtual void ProgAngularAccuracyPCA::synchronize ( )

inlinevirtual

Synchronize with other processors.

Reimplemented in MpiProgAngularAccuracyPCA.

Definition at line 75 of file angular_accuracy_pca.h.

{}

Member Data Documentation

fnNeighbours

FileName ProgAngularAccuracyPCA::fnNeighbours

Definition at line 43 of file angular_accuracy_pca.h.

fnOut

FileName ProgAngularAccuracyPCA::fnOut

Definition at line 43 of file angular_accuracy_pca.h.

fnOutQ

FileName ProgAngularAccuracyPCA::fnOutQ

Definition at line 43 of file angular_accuracy_pca.h.

fnPhantom

FileName ProgAngularAccuracyPCA::fnPhantom

Filenames

Definition at line 43 of file angular_accuracy_pca.h.

mdPartial

MetaDataDb ProgAngularAccuracyPCA::mdPartial

Definition at line 47 of file angular_accuracy_pca.h.

newXdim

int ProgAngularAccuracyPCA::newXdim

Definition at line 53 of file angular_accuracy_pca.h.

newYdim

int ProgAngularAccuracyPCA::newYdim

Definition at line 55 of file angular_accuracy_pca.h.

Nprocessors

size_t ProgAngularAccuracyPCA::Nprocessors

Definition at line 49 of file angular_accuracy_pca.h.

pca

PCAMahalanobisAnalyzer ProgAngularAccuracyPCA::pca

Definition at line 51 of file angular_accuracy_pca.h.

phantomVol

Image<double> ProgAngularAccuracyPCA::phantomVol

Definition at line 45 of file angular_accuracy_pca.h.

rank

size_t ProgAngularAccuracyPCA::rank

Definition at line 49 of file angular_accuracy_pca.h.

---

The documentation for this class was generated from the following files:

• xmipp/libraries/reconstruction/angular_accuracy_pca.h
• xmipp/libraries/reconstruction/angular_accuracy_pca.cpp