CL2DClass Class Reference

#include <mpi_classify_CL2D.h>

Collaboration diagram for CL2DClass:

Public Member Functions
	CL2DClass ()
	CL2DClass (const CL2DClass &other)
	CL2DClass (const CL2DClass &&)=delete
	~CL2DClass ()
CL2DClass &	operator= (const CL2DClass &other)
CL2DClass &	operator= (const CL2DClass &&)=delete
void	updateProjection (const MultidimArray< double > &I, const CL2DAssignment &assigned, bool force=false)
void	updateNonProjection (double corr, bool force=false)
void	transferUpdate (bool centerReference=true)
void	fitBasic (MultidimArray< double > &I, CL2DAssignment &result, bool reverse=false)
void	fit (MultidimArray< double > &I, CL2DAssignment &result)
void	lookForNeighbours (const std::vector< CL2DClass *> listP, int K)
	Look for K-nearest neighbours. More...

Public Attributes
MultidimArray< double >	P
MultidimArray< double >	Pupdate
Polar< std::complex< double > >	polarFourierP
MultidimArray< double >	rotationalCorr
Polar_fftw_plans *	plans
CorrelationAux	corrAux
RotationalCorrelationAux	rotAux
std::vector< CL2DAssignment >	currentListImg
std::vector< CL2DAssignment >	nextListImg
std::vector< double >	nextNonClassCorr
Histogram1D	histClass
Histogram1D	histNonClass
std::vector< int >	neighboursIdx

Detailed Description

CL2DClass class

Definition at line 67 of file mpi_classify_CL2D.h.

Constructor & Destructor Documentation

CL2DClass() [1/3]

CL2DClass::CL2DClass

(

)

Empty constructor

Definition at line 85 of file mpi_classify_CL2D.cpp.

{
    plans = NULL;
    P.initZeros(prm->Ydim, prm->Xdim);
    P.setXmippOrigin();
    Pupdate = P;
}

CL2DClass() [2/3]

CL2DClass::CL2DClass

(

const CL2DClass &

other

)

Copy constructor

Definition at line 93 of file mpi_classify_CL2D.cpp.

{
    *this=other;
}

CL2DClass() [3/3]

CL2DClass::CL2DClass ( const CL2DClass &&  )

delete

~CL2DClass()

CL2DClass::~CL2DClass

(

)

Destructor

Definition at line 116 of file mpi_classify_CL2D.cpp.

{
    delete plans;
}

Member Function Documentation

fit()

void CL2DClass::fit	(	MultidimArray< double > &	I,
		CL2DAssignment &	result
	)

Compute the fit of the input image with this node (check mirrors).

Definition at line 431 of file mpi_classify_CL2D.cpp.

{
    if (currentListImg.size() == 0)
        return;

    // Try this image
    MultidimArray<double> Idirect = I;
    CL2DAssignment resultDirect;
    fitBasic(Idirect, resultDirect);

    // Try its mirror
    CL2DAssignment resultMirror;
    MultidimArray<double> Imirror;
    if (prm->mirrorImages)
    {
        Imirror=I;
        fitBasic(Imirror, resultMirror, true);
    }
    else
        resultMirror.corr=-1e38;

    if (resultMirror.corr > resultDirect.corr)
    {
        I = Imirror;
        result.copyAlignment(resultMirror);
    }
    else
    {
        I = Idirect;
        result.copyAlignment(resultDirect);
    }

    result.likelihood = 0;
    if (!prm->classicalMultiref)
    {
        // Find the likelihood
        int idx;
        histClass.val2index(result.corr, idx);
        if (idx < 0)
            return;
        double likelihoodClass = 0;
        for (int i = 0; i <= idx; i++)
            likelihoodClass += DIRECT_A1D_ELEM(histClass,i);
        histNonClass.val2index(result.corr, idx);
        if (idx < 0)
            return;
        double likelihoodNonClass = 0;
        for (int i = 0; i <= idx; i++)
            likelihoodNonClass += DIRECT_A1D_ELEM(histNonClass,i);
        result.likelihood = likelihoodClass * likelihoodNonClass;
    }
}

fitBasic()

void CL2DClass::fitBasic	(	MultidimArray< double > &	I,
		CL2DAssignment &	result,
		bool	reverse = false
	)

Compute the fit of the input image with this node. The input image is rotationally and traslationally aligned (2 iterations), to make it fit with the node.

Definition at line 220 of file mpi_classify_CL2D.cpp.

{
    if (reverse)
    {
        I.selfReverseX();
        I.setXmippOrigin();
    }

    Matrix2D<double> ARS, ASR, R(3, 3), INV;
    ARS.initIdentity(3);
    ASR = ARS;
    INV = ARS; // avoid creating a new transformation matrix per applyGeometry call
    MultidimArray<double> IauxSR = I, IauxRS = I;
//    Mcorr.resizeNoCopy(I);
    Polar<std::complex<double> > polarFourierI;
    corrAux.transformer1.FourierTransform(P, corrAux.FFT1, false);
#ifdef DEBUG_MORE
    Image<double> save2;
    save2()=P;
    save2.write("PPPI1.xmp");
#endif

    // Align the image with the node
    if (prm->alignImages)
    {
        double shiftXSR=INITIAL_SHIFT_THRESHOLD, shiftYSR=INITIAL_SHIFT_THRESHOLD, bestRotSR=INITIAL_ROTATE_THRESHOLD;
        double shiftXRS=INITIAL_SHIFT_THRESHOLD, shiftYRS=INITIAL_SHIFT_THRESHOLD, bestRotRS=INITIAL_ROTATE_THRESHOLD;

        for (int i = 0; i < 3; i++)
        {
            // Shift then rotate
            if (((shiftXSR > SHIFT_THRESHOLD) || (shiftXSR < (-SHIFT_THRESHOLD))) ||
                ((shiftYSR > SHIFT_THRESHOLD) || (shiftYSR < (-SHIFT_THRESHOLD))))
            {
                bestShift(P, corrAux.FFT1, IauxSR, shiftXSR, shiftYSR, corrAux);
                MAT_ELEM(ASR,0,2) += shiftXSR;
                MAT_ELEM(ASR,1,2) += shiftYSR;
                ASR.inv(INV);
                applyGeometry(xmipp_transformation::LINEAR, IauxSR, I, INV, xmipp_transformation::IS_INV, xmipp_transformation::WRAP);
            }

    #ifdef DEBUG_MORE
            save2()=IauxSR;
            save2.write("PPPIauxSR_afterShift.xmp");
            std::cout << "ASR\n" << ASR << std::endl;
    #endif

            SPEED_UP_tempsDouble;
            if (bestRotSR > ROTATE_THRESHOLD)
            {
                polarFourierTransform<true>(IauxSR, fitBasic_aux, polarFourierI, true,
                                                XSIZE(P) / 5, XSIZE(P) / 2-2, plans, 1);

                bestRotSR = best_rotation(polarFourierP, polarFourierI, rotAux);
                rotation2DMatrix(bestRotSR, R);
                M3x3_BY_M3x3(ASR,R,ASR);
                ASR.inv(INV);
                applyGeometry(xmipp_transformation::LINEAR, IauxSR, I, INV, xmipp_transformation::IS_INV, xmipp_transformation::WRAP);
            }

    #ifdef DEBUG_MORE
            save2()=IauxSR;
            save2.write("PPPIauxSR_afterShiftAndRotation.xmp");
            std::cout << "ASR\n" << ASR << std::endl;
    #endif

            // Rotate then shift
            if (bestRotRS > ROTATE_THRESHOLD)
            {
                polarFourierTransform<true>(IauxRS, fitBasic_aux, polarFourierI, true,
                                                XSIZE(P) / 5, XSIZE(P) / 2-2, plans, 1);

                bestRotRS = best_rotation(polarFourierP, polarFourierI, rotAux);
                rotation2DMatrix(bestRotRS, R);
                M3x3_BY_M3x3(ARS,R,ARS);
                ARS.inv(INV);
                applyGeometry(xmipp_transformation::LINEAR, IauxRS, I, INV, xmipp_transformation::IS_INV, xmipp_transformation::WRAP);
            }

#ifdef DEBUG_MORE
            save2()=IauxRS;
            save2.write("PPPIauxRS_afterRotation.xmp");
            std::cout << "ARS\n" << ARS << std::endl;
#endif

            if (((shiftXRS > SHIFT_THRESHOLD) || (shiftXRS < (-SHIFT_THRESHOLD))) ||
                ((shiftYRS > SHIFT_THRESHOLD) || (shiftYRS < (-SHIFT_THRESHOLD))))
            {
                bestShift(P, corrAux.FFT1, IauxRS, shiftXRS, shiftYRS, corrAux);
                MAT_ELEM(ARS,0,2) += shiftXRS;
                MAT_ELEM(ARS,1,2) += shiftYRS;
                ARS.inv(INV);
                applyGeometry(xmipp_transformation::LINEAR, IauxRS, I, INV, xmipp_transformation::IS_INV, xmipp_transformation::WRAP);
            }

    #ifdef DEBUG_MORE
            save2()=IauxRS;
            save2.write("PPPIauxRS_afterRotationAndShift.xmp");
            std::cout << "ARS\n" << ARS << std::endl;

            char c;
            std::cout << "Press any key\n";
            std::cin >> c;
    #endif
        }
    }

    // Compute the correntropy
    double corrRS=0.0, corrSR=0.0;
    MultidimArray<int> &imask = prm->mask;
    if (prm->useThresholdMask)
    {
        imask.initZeros(IauxRS);
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(IauxRS)
        if (DIRECT_MULTIDIM_ELEM(IauxRS,n)>prm->threshold)
            DIRECT_MULTIDIM_ELEM(imask,n)=1;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(IauxSR)
        if (DIRECT_MULTIDIM_ELEM(IauxSR,n)>prm->threshold)
            DIRECT_MULTIDIM_ELEM(imask,n)=1;
    }
    if (prm->useCorrelation)
    {
        corrRS = corrSR = 0;
        long N = 0;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(P)
        {
            if (DIRECT_MULTIDIM_ELEM(imask,n))
            {
                double Pval = DIRECT_MULTIDIM_ELEM(P, n);
                corrRS += Pval * DIRECT_MULTIDIM_ELEM(IauxRS, n);
                corrSR += Pval * DIRECT_MULTIDIM_ELEM(IauxSR, n);
                ++N;
            }
        }
        if (0 == N) {
            REPORT_ERROR(ERR_UNCLASSIFIED, "This should never happen. If you know how to reproduce"
                    "this issue, please contact developers.");
        }
        corrRS /= N;
        corrSR /= N;
    }
    else
    {
        corrRS = fastCorrentropy(P, IauxRS, prm->sigma,
                                 prm->gaussianInterpolator, imask);
        corrSR = fastCorrentropy(P, IauxSR, prm->sigma,
                                 prm->gaussianInterpolator, imask);
    }
#ifdef DEBUG_MORE
       std::cout << "corrRS=" << corrRS << std::endl;
       std::cout << "corrSR=" << corrSR << std::endl;
#endif


    // Prepare result
    if (reverse)
    {
        // COSS: Check that this is correct
        MAT_ELEM(ARS,0,0) *= -1;
        MAT_ELEM(ARS,1,0) *= -1;
        MAT_ELEM(ASR,0,0) *= -1;
        MAT_ELEM(ASR,1,0) *= -1;
    }

    CL2DAssignment candidateRS, candidateSR;
    candidateRS.readAlignment(ARS);
    candidateSR.readAlignment(ASR);
    if (candidateRS.shiftx * candidateRS.shiftx +
        candidateRS.shifty * candidateRS.shifty > prm->maxShift2)
        candidateRS.corr = 0;
    else
        candidateRS.corr = corrRS;
    if (candidateSR.shiftx * candidateSR.shiftx +
        candidateSR.shifty * candidateSR.shifty > prm->maxShift2)
        candidateSR.corr = 0;
    else
        candidateSR.corr = corrSR;

    if (candidateRS.corr >= candidateSR.corr)
    {
        I = IauxRS;
        result.copyAlignment(candidateRS);
    }
    else if (candidateSR.corr > candidateRS.corr)
    {
        I = IauxSR;
        result.copyAlignment(candidateSR);
    }
    else
        // This is for NaNs, Infs, ...
        result.corr = 0;

#ifdef DEBUG

    Image<double> save;
    save()=P;
    save.write("PPPI1.xmp");
    save()=I;
    save.write("PPPI2.xmp");
    save()=P-I;
    save.write("PPPdiff.xmp");
    std::cout << "sigma=" << prm->sigma << " corr=" << result.corr
    << ". Press" << std::endl;
    char c;
    std::cin >> c;
#endif
}

lookForNeighbours()

void CL2DClass::lookForNeighbours	(	const std::vector< CL2DClass *>	listP,
		int	K
	)

Look for K-nearest neighbours.

Definition at line 486 of file mpi_classify_CL2D.cpp.

{
    int Q = listP.size();
    neighboursIdx.clear();
    if (K == Q)
    {
        // As many neighbours as codes
        for (int q = 0; q < Q; q++)
            neighboursIdx.push_back(q);
    }
    else
    {
        MultidimArray<double> distanceCode;
        distanceCode.initZeros(Q);
        CL2DAssignment assignment;
        MultidimArray<double> I;
        for (int q = 0; q < Q; q++)
        {
            if (listP[q] == this)
                distanceCode(q) = 1;
            else
            {
                I = listP[q]->P;
                fit(I, assignment);
                A1D_ELEM(distanceCode,q) = assignment.corr;
            }
        }

        MultidimArray<int> idx;
        distanceCode.indexSort(idx);
        for (int k = 0; k < K; k++)
            neighboursIdx.push_back(idx(Q - k - 1) - 1);
    }
#ifdef DEBUG
    Image<double> save;
    save()=P;
    save.write("PPPthis.xmp");
    for (int k=0; k<K; k++)
    {
        save()=listP[neighboursIdx[k]]->P;
        save.write((std::string)"PPPneigh"+integerToString(k,1));
    }
    std::cout << "Neighbours saved. Press any key\n";
    char c;
    std::cin >> c;
#endif
}

operator=() [1/2]

CL2DClass & CL2DClass::operator=

(

const CL2DClass &

other

)

Copy assignment

Definition at line 98 of file mpi_classify_CL2D.cpp.

{
    plans = NULL;

    CL2DAssignment assignment;
    assignment.corr = 1;
    Pupdate = other.P;
    nextListImg.push_back(assignment);
    transferUpdate();

    currentListImg = other.currentListImg;
    histClass = other.histClass;
    histNonClass = other.histNonClass;
    neighboursIdx = other.neighboursIdx;

    return *this;
}

operator=() [2/2]

CL2DClass& CL2DClass::operator= ( const CL2DClass &&  )

delete

transferUpdate()

void CL2DClass::transferUpdate

(

bool

centerReference = true

)

Transfer update

Definition at line 133 of file mpi_classify_CL2D.cpp.

{
    if (nextListImg.size() > 0)
    {
        // Take from Pupdate
        double iNq = 1.0 / nextListImg.size();
        Pupdate *= iNq;
        if (prm->normalizeImages)
            Pupdate.statisticsAdjust(0., 1.);
        P = Pupdate;
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(P)
        if (!DIRECT_MULTIDIM_ELEM(prm->mask,n))
            DIRECT_MULTIDIM_ELEM(P,n) = 0;
        Pupdate.initZeros(P);

        // Make sure the image is centered
        if (centerReference && prm->alignImages)
            centerImage(P,corrAux,rotAux);
        FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(P)
        if (!DIRECT_A2D_ELEM(prm->mask,i,j))
            DIRECT_A2D_ELEM(P,i,j) = 0;

        // Compute the polar Fourier transform of the full image
        polarFourierTransform<true>(P, polarFourierP, false, XSIZE(P) / 5,
                                        XSIZE(P) / 2-2, plans, 1);
        size_t finalSize = 2 * polarFourierP.getSampleNoOuterRing() - 1;
        if (XSIZE(rotationalCorr) != finalSize)
            rotationalCorr.resize(finalSize);
        rotAux.local_transformer.setReal(rotationalCorr);

        // Take the list of images
        currentListImg = nextListImg;
        nextListImg.clear();

        // Update the histogram of corrs of elements inside and outside the class
        if (!prm->classicalMultiref)
        {
            MultidimArray<double> classCorr, nonClassCorr;

            classCorr.resizeNoCopy(currentListImg.size());
            FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY1D(classCorr)
            DIRECT_A1D_ELEM(classCorr,i) = currentListImg[i].corr;

            nonClassCorr.resizeNoCopy(nextNonClassCorr.size());
            FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY1D(nonClassCorr)
            DIRECT_A1D_ELEM(nonClassCorr,i) = nextNonClassCorr[i];
            nextNonClassCorr.clear();

            double minC=0., maxC=0.;
            classCorr.computeDoubleMinMax(minC, maxC);
            double minN=0., maxN=0.;
            nonClassCorr.computeDoubleMinMax(minN, maxN);
            double c0 = XMIPP_MIN(minC,minN);
            double cF = XMIPP_MAX(maxC,maxN);
            compute_hist(classCorr, histClass, c0, cF, 200);
            compute_hist(nonClassCorr, histNonClass, c0, cF, 200);
            histClass += 1; // Laplace correction
            histNonClass += 1;
            histClass /= histClass.sum();
            histNonClass /= histNonClass.sum();
        }

#ifdef DEBUG
        histClass.write("PPPclass.txt");
        histNonClass.write("PPPnonClass.txt");
        std::cout << "Histograms have been written. Press any key\n";
        char c;
        std::cin >> c;
#endif

    }
    else
    {
        currentListImg.clear();
        P.initZeros();
    }
}

updateNonProjection()

void CL2DClass::updateNonProjection ( double  corr, bool  force = false  )

inline

Update non-projection

Definition at line 128 of file mpi_classify_CL2D.h.

    {
        if (corr>0 || force)
            nextNonClassCorr.push_back(corr);
    }

updateProjection()

void CL2DClass::updateProjection	(	const MultidimArray< double > &	I,
		const CL2DAssignment &	assigned,
		bool	force = false
	)

Update projection.

Definition at line 121 of file mpi_classify_CL2D.cpp.

{
    if ((assigned.corr > 0 && assigned.objId != BAD_OBJID) || force)
    {
        Pupdate += I;
        nextListImg.push_back(assigned);
    }
}

Member Data Documentation

corrAux

CorrelationAux CL2DClass::corrAux

Definition at line 85 of file mpi_classify_CL2D.h.

currentListImg

std::vector<CL2DAssignment> CL2DClass::currentListImg

Definition at line 91 of file mpi_classify_CL2D.h.

histClass

Histogram1D CL2DClass::histClass

Definition at line 100 of file mpi_classify_CL2D.h.

histNonClass

Histogram1D CL2DClass::histNonClass

Definition at line 103 of file mpi_classify_CL2D.h.

neighboursIdx

std::vector<int> CL2DClass::neighboursIdx

Definition at line 106 of file mpi_classify_CL2D.h.

nextListImg

std::vector<CL2DAssignment> CL2DClass::nextListImg

Definition at line 94 of file mpi_classify_CL2D.h.

nextNonClassCorr

std::vector<double> CL2DClass::nextNonClassCorr

Definition at line 97 of file mpi_classify_CL2D.h.

P

MultidimArray<double> CL2DClass::P

Definition at line 70 of file mpi_classify_CL2D.h.

plans

Polar_fftw_plans* CL2DClass::plans

Definition at line 82 of file mpi_classify_CL2D.h.

polarFourierP

Polar<std::complex <double> > CL2DClass::polarFourierP

Definition at line 76 of file mpi_classify_CL2D.h.

Pupdate

MultidimArray<double> CL2DClass::Pupdate

Definition at line 73 of file mpi_classify_CL2D.h.

rotationalCorr

MultidimArray<double> CL2DClass::rotationalCorr

Definition at line 79 of file mpi_classify_CL2D.h.

rotAux

RotationalCorrelationAux CL2DClass::rotAux

Definition at line 88 of file mpi_classify_CL2D.h.

---

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

• xmipp/libraries/parallel/mpi_classify_CL2D.h
• xmipp/libraries/parallel/mpi_classify_CL2D.cpp