without tilt

Collaboration diagram for without tilt:

Classes
class	ProgValidationNonTilt

Functions
void	ProgValidationNonTilt::readParams ()
void	ProgValidationNonTilt::defineParams ()
void	ProgValidationNonTilt::run ()
	ProgValidationNonTilt::ProgValidationNonTilt ()
void	ProgValidationNonTilt::obtainSumU (const MetaData &tempMd, std::vector< double > &sum_u, std::vector< double > &H0)
void	ProgValidationNonTilt::obtainSumU_2 (const MetaData &mdGallery, const MetaData &tempMd, std::vector< double > &sum_u, std::vector< double > &H0)
void	ProgValidationNonTilt::obtainSumW (const MetaData &tempMd, double &sum_W, std::vector< double > &sum_u, std::vector< double > &H, const double factor)
virtual void	ProgValidationNonTilt::gatherClusterability ()
	Gather alignment. More...
virtual void	ProgValidationNonTilt::synchronize ()
	Synchronize with other processors. More...

Variables
FileName	ProgValidationNonTilt::fnDir
FileName	ProgValidationNonTilt::fnSym
FileName	ProgValidationNonTilt::fnInit
FileName	ProgValidationNonTilt::fnParticles
MetaDataDb	ProgValidationNonTilt::mdPartial
size_t	ProgValidationNonTilt::rank
size_t	ProgValidationNonTilt::Nprocessors
bool	ProgValidationNonTilt::useSignificant
double	ProgValidationNonTilt::significance_noise

Detailed Description

Function Documentation

defineParams()

void ProgValidationNonTilt::defineParams ( )

virtual

Function in which the param of each Program are defined.

Reimplemented from XmippProgram.

Definition at line 49 of file validation_nontilt.cpp.

{
    addUsageLine("Validate a 3D reconstruction from its projections attending to directionality and spread of the angular assignments from a given significant value");
    addParamsLine("  [--i <md_file=\"\">]         : Metadata file with input projections");
    addParamsLine("  [--volume <md_file=\"\">]    : Volume to validate");
    addParamsLine("  [--odir <outputDir=\".\">]   : Output directory");
    addParamsLine("  [--sym <symfile=c1>]         : Enforce symmetry in projections");
    addParamsLine("  [--significance_noise<float=0.95>] : Significance of the alignment with respect the noise");
    addParamsLine("  [--useSignificant]           : Use Significant as alignment method. If not use projection matching");
}

gatherClusterability()

virtual void ProgValidationNonTilt::gatherClusterability ( )

inlinevirtual

Gather alignment.

Reimplemented in MpiProgValidationNonTilt.

Definition at line 71 of file validation_nontilt.h.

{}

obtainSumU()

void ProgValidationNonTilt::obtainSumU	(	const MetaData &	tempMd,
		std::vector< double > &	sum_u,
		std::vector< double > &	H0
	)

Definition at line 187 of file validation_nontilt.cpp.

{

    const size_t tempMdSz= tempMd.size();
    double xRan,yRan,zRan;
    double x,y;
    double sumWRan;
    auto * xRanArray = new double[tempMdSz];
    auto * yRanArray = new double[tempMdSz];
    auto * zRanArray  = new double[tempMdSz];
    std::vector<double> weightV;
    double a;
    std::random_device rd;
    std::mt19937 g(rd());

    for (size_t n=0; n<sum_u.size(); n++)
    {
        for (size_t nS=0; nS<tempMd.size(); nS++)
        {
             /*
                x = sin(tilt*PI/180)*cos(rot*PI/180);
                y = sin(tilt*PI/180)*sin(rot*PI/180);
                z = std::abs(cos(tilt*PI/180));
             */

            //http://mathworld.wolfram.com/SpherePointPicking.html
            x = 2*(double(std::rand())-RAND_MAX/2)/RAND_MAX;
            y = 2*(double(std::rand())-RAND_MAX/2)/RAND_MAX;

            while (x*x+y*y >= 1 )
            {
                x = 2*(std::rand()-RAND_MAX/2)/RAND_MAX;
                y = 2*(std::rand()-RAND_MAX/2)/RAND_MAX;
            }

            xRan = 2*x*std::sqrt(1-x*x-y*y);
            yRan = 2*y*std::sqrt(1-x*x-y*y);
            zRan = std::abs(1-2*(x*x+y*y));

/*          tilt = (double(std::rand())/RAND_MAX)*(PI);
            rot  = (std::rand()-RAND_MAX/2)*(2*PI/RAND_MAX);
            xRan = sin(tilt)*cos(rot);
            yRan = sin(tilt)*sin(rot);
            zRan = std::abs(cos(tilt));
*/
            //std::cout << tilt << " " << rot << std::endl;
            //std::cout << xRan << " " << yRan << " " << zRan << " " << std::endl;
            //zRan=(std::rand()-RAND_MAX/2);
            //norm = std::sqrt(xRan*xRan+yRan*yRan+zRan*zRan);
            //xRan = (xRan/norm);
            //yRan = (yRan/norm);
            //zRan = std::abs(zRan/norm);

            xRanArray[nS] = xRan;
            yRanArray[nS] = yRan;
            zRanArray[nS] = zRan;
            //tempMd.getColumnValues(MDL_WEIGHT, weightV);
            tempMd.getColumnValues(MDL_MAXCC, weightV);

            std::shuffle(weightV.begin(), weightV.end(), g);
        }

        sumWRan = 0;
        double WRan=0;
        double tempWRan, tempW1, tempW2, temp;
        for (size_t nS1=0; nS1<tempMd.size(); nS1++)
        {
            tempWRan = 1e3;
            for (size_t nS2=0; nS2<tempMd.size(); nS2++)
            {
                temp = xRanArray[nS1]*xRanArray[nS2]+yRanArray[nS1]*yRanArray[nS2]+zRanArray[nS1]*zRanArray[nS2];
                if (temp < 1)
                    a = std::abs(std::acos(temp));
                else
                    a = 0;

                if ( (a<tempWRan) && (a > 0.00001) && (temp<1) )
                {
                    tempWRan = a;
                    tempW2 = weightV[nS2];
                    tempW1 = weightV[nS1];
                    WRan = a*std::exp(std::abs(tempW1-tempW2))*std::exp(-(tempW1+tempW2));
                    if (WRan == 0)
                        WRan = a;

                }
            }
            sumWRan += WRan;
        }

        if (sumWRan == 0)
            sumWRan = 0.075*tempMd.size();

        sum_u.at(n)=sumWRan;

    }

    size_t idx = 0;
    while (idx < sum_u.size())
    {
        std::shuffle(sum_u.begin(), sum_u.end(), g);

        if(sum_u.at(0) != sum_u.at(1))
        {
            H0[idx] = sum_u.at(0)/(sum_u.at(0)+sum_u.at(1));
            idx += 1;
        }
    }

    delete[] xRanArray;
    delete[] yRanArray;
    delete[] zRanArray;

}

obtainSumU_2()

void ProgValidationNonTilt::obtainSumU_2	(	const MetaData &	mdGallery,
		const MetaData &	tempMd,
		std::vector< double > &	sum_u,
		std::vector< double > &	H0
	)

Definition at line 376 of file validation_nontilt.cpp.

{

    const size_t tempMdSz= tempMd.size();
    const size_t numGallery= mdGallery.size();

    double xRan,yRan,zRan;
    size_t indx;
    double sumWRan;
    auto * xRanArray = new double[tempMdSz];
    auto * yRanArray = new double[tempMdSz];
    auto * zRanArray  = new double[tempMdSz];
    std::vector<double> weightV;
    double a;

    double rot,tilt;
    bool mirror;
    std::random_device rd;
    std::mt19937 g(rd());
    for (size_t n=0; n<sum_u.size(); n++)
    {
        for (size_t nS=0; nS<tempMd.size(); nS++)
        {
            indx = (size_t) (double( std::rand())*numGallery)/RAND_MAX;

            while ( (indx ==0) || (indx > numGallery) )
                indx = (size_t) (double( std::rand())*numGallery )/RAND_MAX;

            mdGallery.getValue(MDL_ANGLE_ROT,rot,indx);
            mdGallery.getValue(MDL_ANGLE_TILT,tilt,indx);
            mdGallery.getValue(MDL_FLIP,mirror,indx);

            if (mirror == 1)
                tilt = tilt + 180;

            xRan = sin(tilt*PI/180.)*cos(rot*PI/180.);
            yRan = sin(tilt*PI/180.)*sin(rot*PI/180.);
            zRan = std::abs(cos(tilt*PI/180.));

            xRanArray[nS] = xRan;
            yRanArray[nS] = yRan;
            zRanArray[nS] = zRan;
            //tempMd.getColumnValues(MDL_WEIGHT, weightV);
            tempMd.getColumnValues(MDL_MAXCC, weightV);
            std::shuffle(weightV.begin(), weightV.end(), g);
        }


        sumWRan = 0;
        double WRan=0;
        double tempWRan, tempW1, tempW2, temp;
        for (size_t nS1=0; nS1<tempMd.size(); nS1++)
        {
            tempWRan = 1e3;
            for (size_t nS2=0; nS2<tempMd.size(); nS2++)
            {
                temp = xRanArray[nS1]*xRanArray[nS2]+yRanArray[nS1]*yRanArray[nS2]+zRanArray[nS1]*zRanArray[nS2];
                if (temp < 1)
                    a = std::abs(std::acos(temp));
                else
                    a = 0;

                if ( (a<tempWRan) && (a > 0.00001) && (temp<1) )
                {
                    tempWRan = a;
                    tempW2 = weightV[nS2];
                    tempW1 = weightV[nS1];
                    WRan = a*std::exp(std::abs(tempW1-tempW2))*std::exp(-(tempW1+tempW2));
                    if (WRan == 0)
                        WRan = a;
                }
            }
            sumWRan += WRan;
        }

        if (sumWRan == 0)
            sumWRan = 0.075*tempMd.size();

        sum_u.at(n)=sumWRan;

    }

    size_t idx = 0;
    while (idx < sum_u.size())
    {
        std::shuffle(sum_u.begin(), sum_u.end(), g);

        if(sum_u.at(0) != sum_u.at(1))
        {
            H0[idx] = sum_u.at(0)/(sum_u.at(0)+sum_u.at(1));
            idx += 1;
        }
    }

    delete[] xRanArray;
    delete[] yRanArray;
    delete[] zRanArray;

}

obtainSumW()

void ProgValidationNonTilt::obtainSumW	(	const MetaData &	tempMd,
		double &	sum_W,
		std::vector< double > &	sum_u,
		std::vector< double > &	H,
		const double	factor
	)

Definition at line 302 of file validation_nontilt.cpp.

{
    double a;
    double rot,tilt,w;
    double x,y,z;
    double xx,yy,zz;
    bool mirror;
    double w2;
    double tempW;
    double W=0;
    double temp;
    double sumW;

    sumW = 0;
    for (size_t objId : tempMd.ids())
    {
        tempMd.getValue(MDL_ANGLE_ROT,rot,objId);
        tempMd.getValue(MDL_ANGLE_TILT,tilt,objId);
        tempMd.getValue(MDL_FLIP,mirror,objId);
        if (mirror == 1)
            tilt = tilt + 180;

        //tempMd.getValue(MDL_WEIGHT,w,objId);
        tempMd.getValue(MDL_MAXCC,w,objId);
        x = sin(tilt*PI/180.)*cos(rot*PI/180.);
        y = sin(tilt*PI/180.)*sin(rot*PI/180.);
        z = std::abs(cos(tilt*PI/180.));

        tempW = 1e3;
        for (size_t objId2 : tempMd.ids())
        {
            tempMd.getValue(MDL_ANGLE_ROT,rot,objId2);
            tempMd.getValue(MDL_ANGLE_TILT,tilt,objId2);
            //tempMd.getValue(MDL_WEIGHT,w2,objId2);
            tempMd.getValue(MDL_MAXCC,w2,objId2);
            tempMd.getValue(MDL_FLIP,mirror,objId2);
            if (mirror == 1)
                tilt = tilt + 180;

            xx = sin(tilt*PI/180.)*cos(rot*PI/180.);
            yy = sin(tilt*PI/180.)*sin(rot*PI/180.);
            zz = std::abs(cos(tilt*PI/180.));
            temp = x*xx+y*yy+z*zz;

            if (temp < 1)
                a = std::abs(std::acos(temp));
            else
                a = 0;

            if ( (a<tempW) && (a > 0.00001) && (temp<1 ))
            {
                W = a*std::exp(std::abs(w-w2))*std::exp(-(w+w2));
                tempW = a;

                if (W == 0)
                    W = a;
            }
        }
        sumW +=  W;
    }

    //Here the problem is that maybe the changes are only in psi angle. We give the best solution assuming an angular sampling of 0.5º
    if (sumW == 0)
        sumW = 0.075*tempMd.size();

    sum_W = sumW;
    for (size_t n=0; n<sum_u.size(); n++)
    {
        H[n] = sumW/(sumW+factor*sum_u.at(n));
    }
}

ProgValidationNonTilt()

ProgValidationNonTilt::ProgValidationNonTilt

(

)

Definition at line 33 of file validation_nontilt.cpp.

{
    rank=0;
    Nprocessors=1;
}

readParams()

void ProgValidationNonTilt::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 39 of file validation_nontilt.cpp.

{
    fnParticles = getParam("--i");
    fnDir = getParam("--odir");
    fnSym = getParam("--sym");
    fnInit = getParam("--volume");
    useSignificant = checkParam("--useSignificant");
    significance_noise = getDoubleParam("--significance_noise");
}

run()

void ProgValidationNonTilt::run ( )

virtual

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

Reimplemented from XmippProgram.

Definition at line 60 of file validation_nontilt.cpp.

{
    //Clustering Tendency and Cluster Validity Stephen D. Scott
    randomize_random_generator();
    MetaDataDb md,mdGallery,mdOut,mdOut2,mdSort;
    MDRowVec row;

    FileName fnOut,fnOut2, fnGallery;
    fnOut = fnDir+"/clusteringTendency.xmd";
    fnGallery = fnDir+"/gallery.doc";
    fnOut2 = fnDir+"/validation.xmd";
    size_t nSamplesRandom = 500;

    md.read(fnParticles);
    mdGallery.read(fnGallery);
    mdSort.sort(md,MDL_IMAGE_IDX,true,-1,0);

    size_t maxNImg;
    size_t sz = md.size();

    if (useSignificant)
        mdSort.getValue(MDL_IMAGE_IDX,maxNImg,sz);
    else
    {
        mdSort.getValue(MDL_ITEM_ID,maxNImg,sz);
    }

    String expression;
    MDRowVec rowP,row2;
    SymList SL;
    int symmetry, sym_order;
    SL.readSymmetryFile(fnSym.c_str());
    SL.isSymmetryGroup(fnSym.c_str(), symmetry, sym_order);

/*
    double non_reduntant_area_of_sphere = SL.nonRedundantProjectionSphere(symmetry,sym_order);
    double area_of_sphere_no_symmetry = 4.*PI;
    double correction = std::sqrt(non_reduntant_area_of_sphere/area_of_sphere_no_symmetry);
*/
    double correction = 1;
    double validation = 0;
    double num_images = 0;

    MetaDataDb tempMd;
    std::vector<double> sum_u(nSamplesRandom);
    double sum_w=0;
    std::vector<double> H0(nSamplesRandom);
    std::vector<double> H(nSamplesRandom);
    std::vector<double> p(nSamplesRandom);

    if (rank==0)
        init_progress_bar(maxNImg);

    for (size_t idx=0; idx<=maxNImg;idx++)
    {
        if (idx%Nprocessors==rank)
        {
            if (useSignificant)
                expression = formatString("imageIndex == %lu",idx);
            else
                expression = formatString("itemId == %lu",idx);

            tempMd.importObjects(md, MDExpression(expression));


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

            //compute H_0 from noise
            obtainSumU_2(mdGallery, tempMd,sum_u,H0);
            //compute H from experimental
            obtainSumW(tempMd,sum_w,sum_u,H,correction);

            std::sort(H0.begin(),H0.end());
            std::sort(H.begin(),H.end());

            double P = 0;
            for(size_t j=0; j<sum_u.size();j++)
            {
                //P += H0.at(j)/H.at(j);
                P += H0.at(size_t((1-significance_noise)*nSamplesRandom))/H.at(j);
                p.at(j) = H0.at(j)/H.at(j);
            }

            P /= nSamplesRandom;

            if (useSignificant)
                rowP.setValue(MDL_IMAGE_IDX,idx);
            else
                rowP.setValue(MDL_ITEM_ID,idx);

            rowP.setValue(MDL_WEIGHT,P);
            mdPartial.addRow(rowP);
            tempMd.clear();

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

    if (rank==0)
        progress_bar(maxNImg);

    synchronize();
    gatherClusterability();

    if (rank == 0)
    {
        mdPartial.write(fnOut);
        std::vector<double> P;
        mdPartial.getColumnValues(MDL_WEIGHT,P);

        for (size_t idx=0; idx< P.size();idx++)
        {
            if (P[idx] > 1)
                validation += 1.;
            num_images += 1.;
        }
        validation /= num_images;

        row2.setValue(MDL_IMAGE,fnInit);
        row2.setValue(MDL_WEIGHT,validation);
        mdOut2.addRow(row2);
        mdOut2.write(fnOut2);
    }
}

synchronize()

virtual void ProgValidationNonTilt::synchronize ( )

inlinevirtual

Synchronize with other processors.

Reimplemented in MpiProgValidationNonTilt.

Definition at line 74 of file validation_nontilt.h.

{}

Variable Documentation

fnDir

FileName ProgValidationNonTilt::fnDir

Filenames

Definition at line 42 of file validation_nontilt.h.

fnInit

FileName ProgValidationNonTilt::fnInit

Definition at line 42 of file validation_nontilt.h.

fnParticles

FileName ProgValidationNonTilt::fnParticles

Definition at line 42 of file validation_nontilt.h.

fnSym

FileName ProgValidationNonTilt::fnSym

Definition at line 42 of file validation_nontilt.h.

mdPartial

MetaDataDb ProgValidationNonTilt::mdPartial

Definition at line 44 of file validation_nontilt.h.

Nprocessors

size_t ProgValidationNonTilt::Nprocessors

Definition at line 46 of file validation_nontilt.h.

rank

size_t ProgValidationNonTilt::rank

Definition at line 46 of file validation_nontilt.h.

significance_noise

double ProgValidationNonTilt::significance_noise

Definition at line 50 of file validation_nontilt.h.

useSignificant

bool ProgValidationNonTilt::useSignificant

Definition at line 48 of file validation_nontilt.h.