#include <mpi_classify_FTTRI.h>

Inheritance diagram for ProgClassifyFTTRI:

Collaboration diagram for ProgClassifyFTTRI:

Public Member Functions
	ProgClassifyFTTRI (int argc, char **argv)
	Empty constructor. More...

void	readParams ()
	Read argument from command line. More...

void	show ()
	Show. More...

void	defineParams ()
	Usage. More...

void	produceSideInfo ()
	Produce side info. More...

void	produceFTTRI ()
	Produce invariants. More...

void	estimateEpsilonInitialRange ()
	Estimate first epsilon range. More...

double	fttri_distance (const MultidimArray< double > &fttri_i, const MultidimArray< double > &fttri_j)
	Distance between two invariants. More...

void	skipRandomNumberOfUnassignedClasses (size_t &currentPointer, size_t remaining)
	Skip random number of unassigned classes. More...

size_t	wrapperFitness (double epsilon)
	Function to optimize. More...

void	epsilonClassification (double epsilon)
	Epsilon classification. More...

void	searchOptimalEpsilon ()
	Search for optimal epsilon. More...

void	removeSmallClasses ()
	Remove small classes. More...

int	findFarthestFTTRI (const MultidimArray< double > &fttriSeed, const EpsilonClass &class_i)
	Find farthest image. More...

int	findFarthest (const MultidimArray< double > &seed, const EpsilonClass &class_i, bool FTTRI)

void	splitLargeClasses (bool FTTRI)
	Split large classes. More...

void	splitLargeFTTRIClasses ()
	Split large classes. More...

void	computeClassNeighbours (bool FTTRI)
	Compute centroid neighbours. More...

size_t	reassignImagesToClasses (bool FTTRI)
	Reassign images to image classes. More...

void	writeResults (bool FTTRI)

void	alignImagesWithinClasses ()

void	computeClassCentroids (bool FTTRI)

void	run ()

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	fnIn
	Input selfile. More...

FileName	fnRoot
	Output rootname. More...

double	pad
	Padding factor. More...

double	fmax
	Maximum frequency (normalized to 0.5) More...

size_t	Rmax
	Maximum frequency in pixels. More...

double	zoom
	Zoom factor for polar conversion. More...

double	sigma1
	First weight. More...

double	sigma2
	Second weight. More...

size_t	nref
	Desired number of classes. More...

size_t	nMinImages
	Minimum number of images in a class. More...

int	Niter
	Number of iterations. More...

bool	doPhase
	Do phase optimization. More...

FileName	fnFTTRI

int	FTTRIXdim

int	FTTRIYdim

std::shared_ptr< MpiNode >	node

MpiTaskDistributor *	taskDistributor

MetaDataVec	mdIn

std::vector< size_t >	imgsId

int	padXdim

std::vector< EpsilonClass >	epsilonClasses

double	dMin

double	dMax

double	bestObjective

double	bestEpsilon

std::vector< EpsilonClass >	bestEpsilonClasses

Matrix1D< double >	classEpsilon

double	epsilonMax

Matrix1D< unsigned char >	notAssigned

Matrix1D< unsigned char >	notAssigned0

Image< double >	fttriCentroids

Image< double >	imageCentroids

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

Core analysis parameters.

Definition at line 47 of file mpi_classify_FTTRI.h.

Constructor & Destructor Documentation

◆ ProgClassifyFTTRI()

ProgClassifyFTTRI::ProgClassifyFTTRI	(	int	argc,
		char **	argv
	)

Empty constructor.

Definition at line 37 of file mpi_classify_FTTRI.cpp.

 {
     node = std::make_unique<MpiNode>(argc, argv);
     if (!node->isMaster())
         verbose=0;
     taskDistributor=nullptr;
 }

Member Function Documentation

◆ alignImagesWithinClasses()

void ProgClassifyFTTRI::alignImagesWithinClasses ( )

Align images within classes

Definition at line 1115 of file mpi_classify_FTTRI.cpp.

 {
     FileName fnCentroids, fnCandidate;
     fnCentroids=fnRoot+"_image_centroids.mrcs";
     if (node->isMaster())
     {
         std::cerr << "Aligning images within class ..." << std::endl;
         size_t Xdim, Ydim, Zdim, Ndim;
         getImageSize(mdIn, Xdim, Ydim, Zdim, Ndim);
         createEmptyFile(fnCentroids,Xdim,Ydim,1,nref,true,WRITE_OVERWRITE);
     }
     node->barrierWait();
 
     Image<double> centroid, candidate;
     if (node->isMaster())
         init_progress_bar(nref);
     Image<int> mask;
     mask.read(fnRoot+"_mask.mrc");
     centroid().resizeNoCopy(mask());
 
     MultidimArray<double> &mCentroid=centroid();
     MultidimArray<int> &mMask=mask();
     MetaDataVec MDclass, MDaux;
     Matrix2D<double> M;
     MDRowVec row;
     for (size_t i=0; i<nref; i++)
         if (((i+1)%node->size)==node->rank)
         {
             const std::vector<size_t> &class_i=bestEpsilonClasses[i].memberIdx;
             size_t nmax=class_i.size();
             mCentroid.initZeros();
 
             if (nmax!=0)
             {
                 MDclass.clear();
                 for (size_t n=0; n<nmax; n++)
                 {
                     size_t trueIdx=class_i[n];
                     size_t trueId=imgsId[trueIdx];
                     mdIn.getRow(row,trueId);
                     MDclass.addRow(row);
                 }
 
                 // Create centroid
                 alignSetOfImages(MDclass,mCentroid,5,false);
                 FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(mMask)
                 if (!DIRECT_MULTIDIM_ELEM(mMask,n))
                     DIRECT_MULTIDIM_ELEM(mCentroid,n)=0.0;
 
                 // Align images within class to centroid
                 for (size_t objId : MDclass.ids())
                 {
                     MDclass.getValue(MDL_IMAGE,fnCandidate,objId);
                     candidate.read(fnCandidate);
                     candidate().setXmippOrigin();
                     double corr=alignImages(mCentroid,candidate(),M);
                     bool flip;
                     double scale, shiftx, shifty, psi;
                     transformationMatrix2Parameters2D(M, flip, scale, shiftx, shifty, psi);
                     MDclass.setValue(MDL_SHIFT_X, shiftx, objId);
                     MDclass.setValue(MDL_SHIFT_Y, shifty, objId);
                     MDclass.setValue(MDL_ANGLE_PSI, psi, objId);
                     MDclass.setValue(MDL_MAXCC, corr, objId);
                 }
             }
 
             // Write new centroid
             centroid.write(fnCentroids,i+1,true,WRITE_REPLACE);
             MDaux.sort(MDclass,MDL_MAXCC,false);
             MDaux.write(formatString("%s_class%06d.xmd",fnRoot.c_str(),i+1));
             if (node->isMaster())
                 progress_bar(i);
         }
     node->barrierWait();
     if (node->isMaster())
     {
         progress_bar(nref);
         FileName fnClasses=fnRoot+"_classes.xmd", fnClass;
         for (size_t i=0; i<nref; i++)
         {
             fnClass=formatString("%s_class%06d.xmd",fnRoot.c_str(),i+1);
             MDclass.read(fnClass);
             MDclass.write(formatString("class%06d_images@%s",i+1,fnClasses.c_str()),MD_APPEND);
             fnClass.deleteFile();
         }
 
         MetaDataVec MDsummary;
         FileName classesBlock=(String)"classes@"+fnClasses;
         MDsummary.read(classesBlock);
         int nref=1;
         FileName fnRef;
         for (size_t objId : MDsummary.ids())
         {
             fnRef.compose(nref++,fnCentroids);
             MDsummary.setValue(MDL_IMAGE,fnRef,objId);
         }
         MDsummary.write(classesBlock,MD_APPEND);
     }
 }

◆ computeClassCentroids()

void ProgClassifyFTTRI::computeClassCentroids ( bool FTTRI )

compute class centroids

Definition at line 715 of file mpi_classify_FTTRI.cpp.

 {
     FileName fnCentroids, fnCandidate;
     if (FTTRI)
         fnCentroids=fnRoot+"_FTTRI_centroids.mrcs";
     else
         fnCentroids=fnRoot+"_image_centroids.mrcs";
     if (node->isMaster())
     {
         std::cerr << "Computing class centroids ..." << std::endl;
         if (FTTRI)
             createEmptyFile(fnCentroids, 1+LAST_XMIPP_INDEX(FTTRIXdim), FTTRIYdim, 1, nref, true, WRITE_OVERWRITE);
         else
         {
             size_t Xdim, Ydim, Zdim, Ndim;
             getImageSize(mdIn, Xdim, Ydim, Zdim, Ndim);
             createEmptyFile(fnCentroids,Xdim,Ydim,1,nref,true,WRITE_OVERWRITE);
         }
     }
     node->barrierWait();
 
     Image<double> centroid, candidate;
     if (node->isMaster())
         init_progress_bar(nref);
     Image<int> mask;
     if (FTTRI)
     {
         centroid().resizeNoCopy(FTTRIYdim,1+LAST_XMIPP_INDEX(FTTRIXdim));
         classEpsilon.resizeNoCopy(nref);
         classEpsilon.initConstant(-1);
     }
     else
     {
         mask.read(fnRoot+"_mask.mrc");
         centroid().resizeNoCopy(mask());
     }
     MultidimArray<double> &mCentroid=centroid();
     MultidimArray<int> &mMask=mask();
     MultidimArray<double> intraclassDistance, sortedDistance;
     MetaDataVec MDclass;
     for (size_t i=0; i<nref; i++)
         if (((i+1)%node->size)==node->rank)
         {
             const std::vector<size_t> &class_i=bestEpsilonClasses[i].memberIdx;
             size_t nmax=class_i.size();
             mCentroid.initZeros();
 
             if (nmax!=0)
             {
                 if (!FTTRI)
                     MDclass.clear();
                 for (size_t n=0; n<nmax; n++)
                 {
                     size_t trueIdx=class_i[n];
                     if (FTTRI)
                     {
                         fnCandidate.compose(trueIdx+1,fnFTTRI);
                         candidate.read(fnCandidate);
                         mCentroid+=candidate();
                     }
                     else
                     {
                         mdIn.getValue(MDL_IMAGE,fnCandidate,imgsId[trueIdx]);
                         MDclass.setValue(MDL_IMAGE,fnCandidate,MDclass.addObject());
                     }
                 }
                 if (FTTRI)
                 {
                     mCentroid/=(double)nmax;
 
                     // Compute now class epsilon
                     intraclassDistance.resizeNoCopy(nmax);
                     for (size_t n=0; n<nmax; n++)
                     {
                         size_t trueIdx=class_i[n];
                         fnCandidate.compose(trueIdx+1,fnFTTRI);
                         candidate.read(fnCandidate);
                         A1D_ELEM(intraclassDistance,n)=fttri_distance(mCentroid,candidate());
                     }
                     intraclassDistance.sort(sortedDistance);
                     int idxLimit=std::min((size_t)floor(nmax*0.8),nmax-1);
                     double limit=A1D_ELEM(sortedDistance,idxLimit);
                     VEC_ELEM(classEpsilon,i)=A1D_ELEM(intraclassDistance,nmax-1);
 
                     // Robust estimate of the class centroid
                     mCentroid.initZeros();
                     double nactual=0;
                     for (size_t n=0; n<nmax; n++)
                     {
                         if (A1D_ELEM(intraclassDistance,n)>limit)
                             continue;
                         size_t trueIdx=class_i[n];
                         fnCandidate.compose(trueIdx+1,fnFTTRI);
                         candidate.read(fnCandidate);
                         mCentroid+=candidate();
                         nactual++;
                     }
                     mCentroid/=nactual;
                 }
                 else
                 {
                     alignSetOfImages(MDclass,mCentroid,5,false);
                     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(mMask)
                     if (!DIRECT_MULTIDIM_ELEM(mMask,n))
                         DIRECT_MULTIDIM_ELEM(mCentroid,n)=0.0;
                 }
             }
 
             // Write new centroid
             centroid.write(fnCentroids,i+1,true,WRITE_REPLACE);
             if (node->isMaster())
                 progress_bar(i);
         }
     if (node->isMaster())
         progress_bar(nref);
     if (FTTRI)
     {
         MPI_Allreduce(MPI_IN_PLACE, &(VEC_ELEM(classEpsilon,0)), nref, MPI_DOUBLE,
                       MPI_MAX, MPI_COMM_WORLD);
         epsilonMax=classEpsilon.computeMax();
         if (node->isMaster())
             std::cerr << "Maximum epsilon: " << epsilonMax << std::endl;
     }
 }

◆ computeClassNeighbours()

void ProgClassifyFTTRI::computeClassNeighbours ( bool FTTRI )

Compute centroid neighbours.

Definition at line 841 of file mpi_classify_FTTRI.cpp.

 {
     if (node->isMaster())
         std::cerr << "Computing class neighbours ..." << std::endl;
     MultidimArray<double> *ptrCentroids=nullptr;
     if (FTTRI)
     {
         fttriCentroids.read(fnRoot+"_FTTRI_centroids.mrcs");
         ptrCentroids=&fttriCentroids();
     }
     else
     {
         imageCentroids.read(fnRoot+"_image_centroids.mrcs");
         ptrCentroids=&imageCentroids();
     }
     const MultidimArray<double> &mCentroids=*ptrCentroids;
     for (size_t i1=0; i1<nref; ++i1)
         bestEpsilonClasses[i1].neighbours.clear();
 
     MultidimArray<double> centroid_i1, centroid_i2, centroid_i2p;
     if (FTTRI)
     {
         double limit=2*epsilonMax;
         for (size_t i1=0; i1<nref-1; ++i1)
         {
             centroid_i1.aliasImageInStack(mCentroids,i1);
             for (size_t i2=i1+1; i2<nref; ++i2)
             {
                 centroid_i2.aliasImageInStack(mCentroids,i2);
                 double d=fttri_distance(centroid_i1,centroid_i2);
                 if (d<limit)
                 {
                     bestEpsilonClasses[i1].neighbours.push_back(i2);
                     bestEpsilonClasses[i2].neighbours.push_back(i1);
                 }
             }
         }
     }
     else
     {
         MultidimArray<double> corr(nref,nref);
         Matrix2D<double> M;
         AlignmentAux aux;
         CorrelationAux aux2;
         RotationalCorrelationAux aux3;
         corr.initConstant(-1);
         int count=0;
         if (node->isMaster())
             init_progress_bar(nref);
         for (size_t i1=0; i1<nref-1; ++i1)
         {
             centroid_i1.aliasImageInStack(mCentroids,i1);
             centroid_i1.setXmippOrigin();
             for (size_t i2=i1+1; i2<nref; ++i2, ++count)
             {
                 if (((count+1)%node->size)==node->rank)
                 {
                     centroid_i2p.aliasImageInStack(mCentroids,i2);
                     centroid_i2p.setXmippOrigin();
                     centroid_i2=centroid_i2p;
                     A2D_ELEM(corr,i2,i1)=A2D_ELEM(corr,i1,i2)=alignImages(centroid_i1,centroid_i2,M,xmipp_transformation::WRAP,aux,aux2,aux3);
                 }
             }
             if (node->isMaster())
                 progress_bar(i1);
         }
         MPI_Allreduce(MPI_IN_PLACE, MULTIDIM_ARRAY(corr), nref*nref, MPI_DOUBLE,
                       MPI_MAX, MPI_COMM_WORLD);
         if (node->isMaster())
             progress_bar(nref);
 
         // For each class compute new neighbours
         MultidimArray<double> corrRow;
         MultidimArray<int> idx;
         const int Kneighbours=2;
         for (size_t i1=0; i1<nref; ++i1)
         {
             corr.getRow(i1,corrRow);
             corrRow.indexSort(idx);
             std::vector<int> &neighbours_i=bestEpsilonClasses[i1].neighbours;
             for (int k = 0; k < Kneighbours; k++)
                 neighbours_i.push_back(idx(nref - k - 1) - 1);
         }
     }
 }

◆ defineParams()

void ProgClassifyFTTRI::defineParams ( )

virtual

Usage.

Reimplemented from XmippProgram.

Definition at line 82 of file mpi_classify_FTTRI.cpp.

 {
     addUsageLine("Classify in 2D using Fourier Transform based Translational and Rotational Invariants");
     addParamsLine("    -i <infile>                : Metadata or stack with input images");
     addParamsLine("    --oroot <rootname>         : Rootname for output files");
     addParamsLine("    --nref <n>                 : Desired number of classes");
     addParamsLine("   [--padding <p=4>]           : Padding factor (it can be a non-integer factor");
     addParamsLine("   [--maxfreq <f=0.25>]        : Maximum frequency for the spectrum classification");
     addParamsLine("                               : Supply -1 for automatic estimation");
     addParamsLine("   [--zoom <f=1>]              : Perform polar transformation with this zoom factor (>=1) at low frequencies");
     addParamsLine("                               : Log-polar corresponds to an approximate factor of 2.8");
     addParamsLine("   [--nmin <n=5>]              : Minimum number of images in a class to be considered as such");
     addParamsLine("   [--iter <n=10>]             : Number of iterations for FRTTI classfication.");
     addParamsLine("                               : At each iteration, those classes whose size");
     addParamsLine("                               : is smaller than nmin are removed from the classification");
     addParamsLine("   [--sigma1+ <s=0.707>]       : First weight in the FTTRI, see paper for documentation");
     addParamsLine("   [--sigma2+ <s=1.5>]         : Second weight in the FTTRI, see paper for documentation");
     addParamsLine("   [--doPhase]                 : Do also an amplitude and phase classification");
     addExampleLine("mpirun -np 4 `which xmipp_mpi_classify_FTTRI` -i images.xmd --oroot class --nref 64");
 }

◆ epsilonClassification()

void ProgClassifyFTTRI::epsilonClassification ( double epsilon )

Epsilon classification.

Definition at line 357 of file mpi_classify_FTTRI.cpp.

 {
     epsilonClasses.clear();
     notAssigned=notAssigned0;
     auto remaining=(size_t)notAssigned.sum();
     size_t currentPointer=0;
     EpsilonClass newClass;
     FileName fnSeed, fnCandidate;
     Image<double> fttriSeed, fftriCandidate;
     std::vector<size_t> inNewClass;
 
     while (remaining>0)
     {
         // Select new class at random
         skipRandomNumberOfUnassignedClasses(currentPointer,remaining);
         fnSeed.compose(currentPointer+1,fnFTTRI);
         fttriSeed.read(fnSeed);
         MultidimArray<double> &mfftriSeed=fttriSeed();
         if (node->isMaster())
             inNewClass.push_back(currentPointer);
         VEC_ELEM(notAssigned,currentPointer)=0;
 
         // Check if any of the unassigned images belongs to this class
         FOR_ALL_ELEMENTS_IN_MATRIX1D(notAssigned)
         if (VEC_ELEM(notAssigned,i)==1 && (i+1)%node->size==node->rank)
         {
             fnCandidate.compose(i+1,fnFTTRI);
             fftriCandidate.read(fnCandidate);
             double distance=fttri_distance(mfftriSeed,fftriCandidate());
             if (distance<=epsilon)
                 inNewClass.push_back(i);
         }
 
         // Synchronize lists
         if (node->isMaster())
         {
             // Receive all new elements
             for (size_t rank=1; rank<node->size; rank++)
                 receiveListFromRank(inNewClass,rank);
             // And redistribute the new list
             for (size_t rank=1; rank<node->size; rank++)
                 sendListToRank(inNewClass,rank);
         }
         else
         {
             // Send my own new elements
             sendListToRank(inNewClass,0);
             // Receive the whole list of new elements
             receiveListFromRank(inNewClass,0);
         }
 
         // Update internal structures
         size_t imax=inNewClass.size();
         for (size_t i=0; i<imax; ++i)
             VEC_ELEM(notAssigned,inNewClass[i])=0;
         remaining=(size_t)notAssigned.sum();
         if (node->isMaster())
         {
             newClass.memberIdx=inNewClass;
             epsilonClasses.push_back(newClass);
         }
         inNewClass.clear();
     }
 }

◆ estimateEpsilonInitialRange()

void ProgClassifyFTTRI::estimateEpsilonInitialRange ( )

Estimate first epsilon range.

Definition at line 243 of file mpi_classify_FTTRI.cpp.

 {
     if (node->isMaster())
     {
         dMin=1e3;
         dMax=-1;
         MultidimArray<int> perm;
         int N=mdIn.size();
         randomPermutation(N,perm);
         N=XMIPP_MIN(N,50);
         std::vector< MultidimArray<double> > fttri;
         Image<double> aux;
         FileName fn;
         for(int i=0; i<N; i++)
         {
             fn.compose(A1D_ELEM(perm,i)+1,fnFTTRI);
             aux.read(fn);
             fttri.push_back(aux());
         }
 
         int idx=0;
         for(int i=0; i<N-1; i++)
         {
             const MultidimArray<double> &fttri_i=fttri[i];
             for(int j=i+1; j<N; j++, idx++)
             {
                 const MultidimArray<double> &fttri_j=fttri[j];
                 double d=fttri_distance(fttri_i,fttri_j);
                 dMin=std::min(d,dMin);
                 dMax=std::max(d,dMax);
             }
         }
         if (verbose)
             std::cout << "Initial epsilon range: [" << dMin << "," << dMax << "]\n";
     }
 }

◆ findFarthest()

int ProgClassifyFTTRI::findFarthest	(	const MultidimArray< double > &	seed,
		const EpsilonClass &	class_i,
		bool	FTTRI
	)

Definition at line 556 of file mpi_classify_FTTRI.cpp.

 {
     Image<double> candidate;
     FileName fnCandidate;
     int nmax=class_i.memberIdx.size();
     double maxDistance;
     if (FTTRI)
         maxDistance=-1;
     else
         maxDistance=1;
     int nMaxDistance=-1;
     const std::vector<size_t> &class_i_members=class_i.memberIdx;
     double d;
     Matrix2D<double> M;
     AlignmentAux aux;
     CorrelationAux aux2;
     RotationalCorrelationAux aux3;
     for (int n=0; n<nmax; n++)
     {
         if (FTTRI)
             fnCandidate.compose(class_i_members[n]+1,fnFTTRI);
         else
             mdIn.getValue(MDL_IMAGE,fnCandidate,imgsId[class_i_members[n]]);
         candidate.read(fnCandidate);
         candidate().setXmippOrigin();
         if (FTTRI)
             d=fttri_distance(seed,candidate());
         else
             d=alignImages(seed,candidate(),M,xmipp_transformation::WRAP,aux,aux2,aux3);
         if ((d>maxDistance && FTTRI) || (d<maxDistance && !FTTRI))
         {
             maxDistance=d;
             nMaxDistance=n;
         }
     }
     return nMaxDistance;
 }

◆ findFarthestFTTRI()

int ProgClassifyFTTRI::findFarthestFTTRI	(	const MultidimArray< double > &	fttriSeed,
		const EpsilonClass &	class_i
	)

Find farthest image.

◆ fttri_distance()

double ProgClassifyFTTRI::fttri_distance	(	const MultidimArray< double > &	fttri_i,
		const MultidimArray< double > &	fttri_j
	)

Distance between two invariants.

Definition at line 280 of file mpi_classify_FTTRI.cpp.

 {
     double retval=0;
     const size_t unroll=4;
     size_t nmax=unroll*(MULTIDIM_SIZE(fttri_i)/unroll);
     double* ptrI=nullptr;
     double* ptrJ=nullptr;
     size_t n;
     for (n=0, ptrI=MULTIDIM_ARRAY(fttri_i),ptrJ=MULTIDIM_ARRAY(fttri_j);
          n<nmax; n+=unroll, ptrI+=unroll, ptrJ+=unroll)
     {
         double diff0=*ptrI-*ptrJ;
         retval+=diff0*diff0;
         double diff1=*(ptrI+1)-*(ptrJ+1);
         retval+=diff1*diff1;
         double diff2=*(ptrI+1)-*(ptrJ+1);
         retval+=diff2*diff2;
         double diff3=*(ptrI+1)-*(ptrJ+1);
         retval+=diff3*diff3;
     }
     for (n=nmax, ptrI=MULTIDIM_ARRAY(fttri_i)+nmax, ptrJ=MULTIDIM_ARRAY(fttri_j)+nmax;
          n<MULTIDIM_SIZE(fttri_i); ++n, ++ptrI, ++ptrJ)
     {
         double diff0=*ptrI-*ptrJ;
         retval+=diff0*diff0;
     }
     return retval/(MULTIDIM_SIZE(fttri_i));
 }

◆ produceFTTRI()

void ProgClassifyFTTRI::produceFTTRI ( )

Produce invariants.

Definition at line 154 of file mpi_classify_FTTRI.cpp.

 {
     if (verbose && node->rank==0)
         std::cerr << "Computing FTTRI ...\n";
 
     // Read mask
     Image<int> mask;
     mask.read(fnRoot+"_mask.mrc");
     MultidimArray<int> &mMask=mask();
     mMask.setXmippOrigin();
 
     // Compute weights
     Matrix1D<double> R, weight1, weight2;
 
     // Process all images
     size_t Nimgs=mdIn.size();
     if (verbose && node->rank==0)
         init_progress_bar(Nimgs);
     size_t first, last;
     FileName fnImg;
     MultidimArray<double> padI(padXdim,padXdim), magFTI, filteredMagFTI, polarFilteredMagFTI, magFTpolarFilteredMagFTI;
     padI.setXmippOrigin();
     Image<double> I, avgMagFTI, centralMagFTpolarFilteredMagFTI;
     avgMagFTI().initZeros(padI);
     MultidimArray< std::complex<double> > FTI, FTpolarFilteredMagFTI;
     FourierTransformer transformer1, transformer2;
     while (taskDistributor->getTasks(first, last))
         for (size_t idx=first; idx<=last; ++idx)
         {
             if (verbose && node->rank==0)
                 progress_bar(idx);
             mdIn.getValue(MDL_IMAGE,fnImg,imgsId[idx]);
             I.read(fnImg);
             I().setXmippOrigin();
             MultidimArray<double> &mI=I();
 
             // Mask input images
             FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(mMask)
             if (!DIRECT_MULTIDIM_ELEM(mMask,n))
                 DIRECT_MULTIDIM_ELEM(mI,n)=0.0;
 
             // Pad input images
             mI.window(padI,FIRST_XMIPP_INDEX(padXdim),FIRST_XMIPP_INDEX(padXdim),
                       LAST_XMIPP_INDEX(padXdim),LAST_XMIPP_INDEX(padXdim));
 
             // Compute full Fourier transform and its magnitude
             transformer1.completeFourierTransform(padI,FTI);
             FFT_magnitude(FTI,magFTI);
             CenterFFT(magFTI,true);
             magFTI.window(filteredMagFTI,
                           FIRST_XMIPP_INDEX(Rmax),
                           FIRST_XMIPP_INDEX(Rmax),
                           LAST_XMIPP_INDEX(Rmax),
                           LAST_XMIPP_INDEX(Rmax));
 
             // Now express it in polar coordinates and apply weight
             image_convertCartesianToPolar_ZoomAtCenter(filteredMagFTI, polarFilteredMagFTI, R, zoom, 3, Rmax/2, Rmax, 0, PI, Rmax);
             if (VEC_XSIZE(weight1)==0)
             {
                 weight1.resizeNoCopy(R);
                 weight2.resizeNoCopy(R);
                 FOR_ALL_ELEMENTS_IN_MATRIX1D(R)
                 {
                     VEC_ELEM(weight1,i)=pow(VEC_ELEM(R,i),sigma1);
                     VEC_ELEM(weight2,i)=pow((Rmax-VEC_ELEM(R,i)),sigma2);
                 }
             }
             FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(polarFilteredMagFTI)
             A2D_ELEM(polarFilteredMagFTI,i,j)*=VEC_ELEM(weight1,j);
 
             // And take new Fourier transform
             transformer2.completeFourierTransform(polarFilteredMagFTI,FTpolarFilteredMagFTI);
             FFT_magnitude(FTpolarFilteredMagFTI,magFTpolarFilteredMagFTI);
             CenterFFT(magFTpolarFilteredMagFTI,true);
             FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(magFTpolarFilteredMagFTI)
             A2D_ELEM(magFTpolarFilteredMagFTI,i,j)*=VEC_ELEM(weight2,j);
             magFTpolarFilteredMagFTI.setXmippOrigin();
             magFTpolarFilteredMagFTI.window(centralMagFTpolarFilteredMagFTI(),
                                             FIRST_XMIPP_INDEX(FTTRIYdim),0,
                                             LAST_XMIPP_INDEX(FTTRIYdim),LAST_XMIPP_INDEX(FTTRIXdim));
             centralMagFTpolarFilteredMagFTI().rangeAdjust(1,255);
             centralMagFTpolarFilteredMagFTI().selfLog10();
             centralMagFTpolarFilteredMagFTI.write(fnFTTRI,idx+1,true,WRITE_REPLACE);
         }
     taskDistributor->wait();
     if (verbose && node->rank==0)
         progress_bar(Nimgs);
 }

◆ produceSideInfo()

void ProgClassifyFTTRI::produceSideInfo ( )

Produce side info.

Definition at line 104 of file mpi_classify_FTTRI.cpp.

 {
     mdIn.read(fnIn);
 
     // Get input size
     size_t zdim, ydim, xdim, ndim;
     getImageSize(mdIn, xdim, ydim, zdim, ndim);
     if (node->isMaster())
     {
         if (zdim!=1)
             REPORT_ERROR(ERR_MULTIDIM_DIM,"This program is only intended for images, not volumes");
         if (xdim!=ydim)
             REPORT_ERROR(ERR_MULTIDIM_DIM,"This program is only intended for squared images");
     }
     padXdim=(int)(pad*xdim);
     Rmax=(size_t)floor(fmax*padXdim);
 
     size_t blockSize=mdIn.size()/(5*node->size);
     if (blockSize==0)
         blockSize=1;
     mdIn.findObjects(imgsId);
     taskDistributor = new MpiTaskDistributor(mdIn.size(), blockSize, node);
     fnFTTRI=fnRoot+"_FTTRI.mrcs";
     FTTRIXdim=(int)((Rmax+1)*0.35);
     FTTRIYdim=(int)((Rmax+1)*0.55);
     if (node->isMaster())
     {
         // Create output mask
         Image<int> mask;
         mask().initZeros(xdim,xdim);
         mask().setXmippOrigin();
         double R2=0.25*xdim*xdim;
         MultidimArray<int> &mMask=mask();
         FOR_ALL_ELEMENTS_IN_ARRAY2D(mMask)
         if (i*i+j*j<R2)
             A2D_ELEM(mMask,i,j)=1;
         mask.write(fnRoot+"_mask.mrc");
 
         // Create output FTTRI
 #ifndef FTTRI_ALREADY_COMPUTED
 
         if (fileExists(fnFTTRI))
             fnFTTRI.deleteFile();
         createEmptyFile(fnFTTRI, 1+LAST_XMIPP_INDEX(FTTRIXdim), FTTRIYdim, 1, mdIn.size(), true, WRITE_OVERWRITE);
 #endif
 
     }
     node->barrierWait();
 }

◆ readParams()

void ProgClassifyFTTRI::readParams ( )

virtual

Read argument from command line.

Reimplemented from XmippProgram.

Definition at line 46 of file mpi_classify_FTTRI.cpp.

 {
     fnIn = getParam("-i");
     fnRoot = getParam("--oroot");
     pad = getDoubleParam("--padding");
     fmax = getDoubleParam("--maxfreq");
     zoom = getDoubleParam("--zoom");
     sigma1 = getDoubleParam("--sigma1");
     sigma2 = getDoubleParam("--sigma2");
     nref = getIntParam("--nref");
     nMinImages = getIntParam("--nmin");
     Niter = getIntParam("--iter");
     doPhase  = checkParam("--doPhase");
 }

◆ reassignImagesToClasses()

size_t ProgClassifyFTTRI::reassignImagesToClasses ( bool FTTRI )

Reassign images to image classes.

Definition at line 927 of file mpi_classify_FTTRI.cpp.

 {
     Matrix1D<short int> newAssignment;
     newAssignment.resizeNoCopy(mdIn.size());
     newAssignment.initConstant(-1);
 
     if (node->isMaster())
     {
         std::cerr << "Reassigning images to classes ..." << std::endl;
         init_progress_bar(nref);
     }
 
     MultidimArray<double> *ptrCentroids=nullptr;
     if (FTTRI)
         ptrCentroids=&fttriCentroids();
     else
         ptrCentroids=&imageCentroids();
     const MultidimArray<double> &mCentroids=*ptrCentroids;
 
     Image<double> candidate;
     FileName fnCandidate;
     MultidimArray<double> own_class, neighbour, candidateCopy;
     size_t changes=0;
     AlignmentAux aux;
     CorrelationAux aux2;
     RotationalCorrelationAux aux3;
     Matrix2D<double> M;
     for (size_t i=0; i<nref; i++)
         if (((i+1)%node->size)==node->rank)
         {
             std::vector<size_t> &class_i=bestEpsilonClasses[i].memberIdx;
             std::vector<int> &neighbours_i=bestEpsilonClasses[i].neighbours;
             int nmax=class_i.size();
             int neighMax=neighbours_i.size();
             own_class.aliasImageInStack(mCentroids,i);
             own_class.setXmippOrigin();
             for (int n=0; n<nmax; ++n)
             {
                 // Get image n
                 int trueIdx=class_i[n];
                 if (FTTRI)
                     fnCandidate.compose(trueIdx+1,fnFTTRI);
                 else
                     mdIn.getValue(MDL_IMAGE,fnCandidate,imgsId[trueIdx]);
                 candidate.read(fnCandidate);
                 candidate().setXmippOrigin();
 
                 // Compare to its own class
                 const MultidimArray<double> &mCandidate=candidate();
                 double bestD;
                 if (FTTRI)
                     bestD=fttri_distance(own_class,mCandidate);
                 else
                 {
                     candidateCopy=candidate();
                     bestD=alignImages(own_class,candidateCopy,M,xmipp_transformation::WRAP,aux,aux2,aux3);
                 }
                 VEC_ELEM(newAssignment,trueIdx)=i;
 
                 // Now compare to the rest of neighbours
                 bool alreadyChanged=false;
                 for (int neigh=0; neigh<neighMax; neigh++)
                 {
                     int neighbourIdx=neighbours_i[neigh];
                     neighbour.aliasImageInStack(mCentroids,neighbourIdx);
                     neighbour.setXmippOrigin();
                     double d;
                     if (FTTRI)
                         d=fttri_distance(neighbour,mCandidate);
                     else
                     {
                         candidateCopy=candidate();
                         d=alignImages(neighbour,candidateCopy,M,xmipp_transformation::WRAP,aux,aux2,aux3);
                     }
                     if ((d<bestD && FTTRI) || (d>bestD && !FTTRI))
                     {
                         bestD=d;
                         VEC_ELEM(newAssignment,trueIdx)=neighbourIdx;
                         if (!alreadyChanged)
                         {
                             ++changes;
                             alreadyChanged=true;
                         }
                     }
                 }
             }
             if (node->isMaster())
                 progress_bar(i);
         }
 
     // Share assignments
     MPI_Allreduce(MPI_IN_PLACE, &(VEC_ELEM(newAssignment,0)), VEC_XSIZE(newAssignment), MPI_SHORT,
                   MPI_MAX, MPI_COMM_WORLD);
     MPI_Allreduce(MPI_IN_PLACE, &changes, 1, MPI_LONG, MPI_SUM, MPI_COMM_WORLD);
     if (node->isMaster())
     {
         progress_bar(nref);
         std::cerr << "Number of images changed: " << changes << std::endl;
     }
 
     // Create new best classes
     std::vector<EpsilonClass> newBestEpsilonClasses;
     newBestEpsilonClasses.resize(nref);
     FOR_ALL_ELEMENTS_IN_MATRIX1D(newAssignment)
     if (VEC_ELEM(newAssignment,i)>=0)
         newBestEpsilonClasses[VEC_ELEM(newAssignment,i)].memberIdx.push_back(i);
     bestEpsilonClasses=newBestEpsilonClasses;
     return changes;
 }

◆ removeSmallClasses()

void ProgClassifyFTTRI::removeSmallClasses ( )

Remove small classes.

Definition at line 522 of file mpi_classify_FTTRI.cpp.

 {
     if (node->isMaster())
     {
         std::sort(bestEpsilonClasses.begin(), bestEpsilonClasses.end(), SDescendingClusterSort());
 
         std::vector< EpsilonClass > goodEpsilonClasses;
         size_t imax=bestEpsilonClasses.size();
         size_t imagesRemoved=0;
         size_t classesRemoved=0;
         imax=XMIPP_MIN(imax,nref);
         for (size_t i=0; i<imax; i++)
         {
             const std::vector<size_t> class_i=bestEpsilonClasses[i].memberIdx;
             size_t nmax=class_i.size();
             if (nmax>=nMinImages)
                 goodEpsilonClasses.push_back(bestEpsilonClasses[i]);
             else
             {
                 imagesRemoved+=nmax;
                 classesRemoved++;
                 for (size_t n=0; n<nmax; ++n)
                     VEC_ELEM(notAssigned0,class_i[n])=0;
             }
         }
         if (imagesRemoved>0)
             std::cout << "Removal of " << imagesRemoved << " images from "
             << classesRemoved << " classes for being in too small classes\n";
         bestEpsilonClasses=goodEpsilonClasses;
     }
     MPI_Bcast(MATRIX1D_ARRAY(notAssigned0), VEC_XSIZE(notAssigned0), MPI_CHAR, 0, MPI_COMM_WORLD);
 }

◆ run()

void ProgClassifyFTTRI::run ( )

virtual

Run.

Reimplemented from XmippProgram.

Definition at line 1216 of file mpi_classify_FTTRI.cpp.

 {
     show();
     produceSideInfo();
 #ifndef FTTRI_ALREADY_COMPUTED
 
     produceFTTRI();
 #endif
 
     estimateEpsilonInitialRange();
     notAssigned0.resizeNoCopy(imgsId.size());
     notAssigned0.initConstant(1);
 
     searchOptimalEpsilon();
     if (node->isMaster())
         std::cout << "Final epsilon: " << bestEpsilon << "\n";
     removeSmallClasses();
 
     // Amplitude classification
     if (node->isMaster())
         std::cout << std::endl << "Amplitude classification ..." << std::endl;
     for (int n=0; n<Niter; n++)
     {
         if (node->isMaster())
             std::cout << std::endl << "Iteration " << n << std::endl;
         splitLargeClasses(true);
         computeClassCentroids(true);
         computeClassNeighbours(true);
         reassignImagesToClasses(true);
         removeSmallClasses();
     }
 
     computeClassCentroids(true);
     if (node->isMaster())
         writeResults(true);
 
     // Amplitude and phase classification
     if (doPhase)
     {
         if (node->isMaster())
             std::cout << std::endl << "Amplitude and phase classification ..." << std::endl;
         for (int n=0; n<Niter; n++)
         {
             if (node->isMaster())
                 std::cout << std::endl << "Iteration " << n << std::endl;
             splitLargeClasses(false);
             computeClassCentroids(false);
             computeClassNeighbours(false);
             reassignImagesToClasses(false);
             removeSmallClasses();
         }
         if (node->isMaster())
             writeResults(false);
     }
 
     // Align images within classes
     node->barrierWait();
     alignImagesWithinClasses();
 
 }

◆ searchOptimalEpsilon()

void ProgClassifyFTTRI::searchOptimalEpsilon ( )

Search for optimal epsilon.

Definition at line 438 of file mpi_classify_FTTRI.cpp.

 {
     if (node->isMaster())
     {
         std::cerr << "Computing epsilon classification ..." << std::endl;
         bestEpsilon=-1;
         bestObjective=1e38;
         double epsilonLeft=dMin;
         double epsilonRight=dMax;
         size_t Nleft=wrapperFitness(epsilonLeft);
         size_t Nright=wrapperFitness(epsilonRight);
         if (Nleft!=nref && Nright!=nref)
         {
             size_t Ndivisions=0;
             do
             {
                 // Adjust margins
                 if (nref>Nleft)
                 {
                     epsilonRight=epsilonLeft;
                     Nright=Nleft;
                     epsilonLeft*=0.9;
                     Nleft=wrapperFitness(epsilonLeft);
                     if (Nleft==0)
                         break;
                 }
                 else if (nref<Nright)
                 {
                     epsilonLeft=epsilonRight;
                     Nleft=Nright;
                     epsilonRight*=1.1;
                     Nright=wrapperFitness(epsilonRight);
                     if (Nright==0)
                         break;
                 }
                 else
                 {
                     // Interval is correct, look in the middle
                     double epsilonCentral=0.5*(epsilonLeft+epsilonRight);
                     size_t Ncentral=wrapperFitness(epsilonCentral);
                     Ndivisions++;
                     if (Ncentral==nref || Ndivisions==10)
                         break;
                     if (Ncentral>nref)
                     {
                         Nleft=Ncentral;
                         epsilonLeft=epsilonCentral;
                     }
                     else
                     {
                         Nright=Ncentral;
                         epsilonRight=epsilonCentral;
                     }
                 }
             }
             while (1);
         }
         double finalize=-1e6;
         MPI_Bcast(&finalize, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
     }
     else
     {
         // Wait for orders from the master
         while (1)
         {
             double epsilon;
             MPI_Bcast(&epsilon, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
             if (epsilon>-5e5)
                 epsilonClassification(epsilon);
             else
                 break;
         }
     }
 }

◆ show()

void ProgClassifyFTTRI::show ( )

Show.

Definition at line 62 of file mpi_classify_FTTRI.cpp.

 {
     if (!verbose)
         return;
     std::cout
     << "Input file:  " << fnIn       << std::endl
     << "Output root: " << fnRoot     << std::endl
     << "Padding:     " << pad        << std::endl
     << "MaxFreq:     " << fmax       << std::endl
     << "Zoom:        " << zoom       << std::endl
     << "Sigma1:      " << sigma1     << std::endl
     << "Sigma2:      " << sigma2     << std::endl
     << "N.Classes:   " << nref       << std::endl
     << "N.Minimum:   " << nMinImages << std::endl
     << "Iterations:  " << Niter      << std::endl
     << "Do phase:    " << doPhase    << std::endl
     ;
 }

◆ skipRandomNumberOfUnassignedClasses()

void ProgClassifyFTTRI::skipRandomNumberOfUnassignedClasses	(	size_t &	currentPointer,
		size_t	remaining
	)

Skip random number of unassigned classes.

Definition at line 311 of file mpi_classify_FTTRI.cpp.

 {
     // Now the master selects the first class at random
     if (node->isMaster())
     {
         // Move current pointer to next not assigned image
         while (VEC_ELEM(notAssigned,currentPointer)==0)
             currentPointer=(currentPointer+1)%VEC_XSIZE(notAssigned);
 
         // Now skip some non empty
         auto skip=(size_t) (remaining*rnd_unif());
         while (skip>0)
         {
             currentPointer=(currentPointer+1)%VEC_XSIZE(notAssigned);
             // Adjust to next not assigned image
             while (VEC_ELEM(notAssigned,currentPointer)==0)
                 currentPointer=(currentPointer+1)%VEC_XSIZE(notAssigned);
             --skip;
         }
     }
     MPI_Bcast(&currentPointer, 1, MPI_LONG, 0, MPI_COMM_WORLD);
 }

◆ splitLargeClasses()

void ProgClassifyFTTRI::splitLargeClasses ( bool FTTRI )

Split large classes.

Definition at line 595 of file mpi_classify_FTTRI.cpp.

 {
     if (node->isMaster())
     {
         std::cerr << "Splitting large classes ..." << std::endl;
         Image<double> seed1, seed2, candidate;
         MultidimArray<double> candidateCopy;
         FileName fnSeed, fnCandidate;
         EpsilonClass c1, c2;
         std::sort(bestEpsilonClasses.begin(), bestEpsilonClasses.end(), SDescendingClusterSort());
         Matrix2D<double> M;
         AlignmentAux aux;
         CorrelationAux aux2;
         RotationalCorrelationAux aux3;
         while (bestEpsilonClasses.size()!=nref)
         {
             EpsilonClass &class_0=bestEpsilonClasses[0];
             std::vector<size_t> &class_0_members=class_0.memberIdx;
             if (class_0_members.size()<nMinImages)
                 break;
 
             // Find the image that is farthest from the center
             if (FTTRI)
                 fnSeed.compose(class_0_members[0]+1,fnFTTRI);
             else
                 mdIn.getValue(MDL_IMAGE,fnSeed,imgsId[class_0_members[0]]);
             seed1.read(fnSeed);
             seed1().setXmippOrigin();
             int n1=findFarthest(seed1(),class_0,FTTRI);
             if (FTTRI)
                 fnSeed.compose(class_0_members[n1]+1,fnFTTRI);
             else
                 mdIn.getValue(MDL_IMAGE,fnSeed,imgsId[class_0_members[n1]]);
             seed1.read(fnSeed);
             seed1().setXmippOrigin();
 
             // Now find the one that is farthest from i1
             int n2=findFarthest(seed1(),class_0,FTTRI);
             if (FTTRI)
                 fnSeed.compose(class_0_members[n2]+1,fnFTTRI);
             else
                 mdIn.getValue(MDL_IMAGE,fnSeed,imgsId[class_0_members[n2]]);
             seed2.read(fnSeed);
             seed2().setXmippOrigin();
 
             // Now split
             c1.memberIdx.clear();
             c2.memberIdx.clear();
             c1.memberIdx.push_back(class_0_members[n1]);
             c2.memberIdx.push_back(class_0_members[n2]);
             int nmax=class_0_members.size();
             const MultidimArray<double> &mSeed1=seed1();
             const MultidimArray<double> &mSeed2=seed2();
             for (int n=0; n<nmax; ++n)
             {
                 if (n==n1 || n==n2)
                     continue;
                 size_t trueIdx=class_0_members[n];
                 if (FTTRI)
                     fnCandidate.compose(trueIdx+1,fnFTTRI);
                 else
                     mdIn.getValue(MDL_IMAGE,fnCandidate,imgsId[trueIdx]);
                 candidate.read(fnCandidate);
                 if (FTTRI)
                 {
                     const MultidimArray<double> &mCandidate=candidate();
                     double d1=fttri_distance(mSeed1,mCandidate);
                     double d2=fttri_distance(mSeed2,mCandidate);
                     if (d1<d2)
                         c1.memberIdx.push_back(trueIdx);
                     else
                         c2.memberIdx.push_back(trueIdx);
                 }
                 else
                 {
                     candidate().setXmippOrigin();
                     candidateCopy=candidate();
                     double d1=alignImages(mSeed1,candidate(),M,xmipp_transformation::WRAP,aux,aux2,aux3);
                     double d2=alignImages(mSeed2,candidateCopy,M,xmipp_transformation::WRAP,aux,aux2,aux3);
                     if (d1>d2)
                         c1.memberIdx.push_back(trueIdx);
                     else
                         c2.memberIdx.push_back(trueIdx);
                 }
             }
             bestEpsilonClasses.erase(bestEpsilonClasses.begin());
             if (c1.memberIdx.size()>nMinImages)
                 bestEpsilonClasses.push_back(c1);
             if (c2.memberIdx.size()>nMinImages)
                 bestEpsilonClasses.push_back(c2);
             std::sort(bestEpsilonClasses.begin(), bestEpsilonClasses.end(), SDescendingClusterSort());
         }
         // Broadcast the best classes
         for (size_t i=0; i<nref; ++i)
         {
             int classSize=bestEpsilonClasses[i].memberIdx.size();
             MPI_Bcast(&classSize, 1, MPI_INT, 0, MPI_COMM_WORLD);
             MPI_Bcast(&(bestEpsilonClasses[i].memberIdx[0]), classSize, MPI_LONG, 0, MPI_COMM_WORLD);
         }
     }
     else
     {
         // Receive the best classes
         bestEpsilonClasses.clear();
         size_t classSize;
         size_t buffer[50000];
         EpsilonClass C;
         for (size_t i=0; i<nref; ++i)
         {
             MPI_Bcast(&classSize, 1, MPI_INT, 0, MPI_COMM_WORLD);
             MPI_Bcast(buffer, classSize, MPI_LONG, 0, MPI_COMM_WORLD);
             C.memberIdx.clear();
             for (size_t n=0; n<classSize; ++n)
                 C.memberIdx.push_back(buffer[n]);
             bestEpsilonClasses.push_back(C);
         }
     }
 }

◆ splitLargeFTTRIClasses()

void ProgClassifyFTTRI::splitLargeFTTRIClasses ( )

Split large classes.

◆ wrapperFitness()

size_t ProgClassifyFTTRI::wrapperFitness ( double epsilon )

Function to optimize.

Definition at line 422 of file mpi_classify_FTTRI.cpp.

 {
     MPI_Bcast(&epsilon, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); // Workers are waiting at searchOptimalEpsilon
     epsilonClassification(epsilon);
     size_t retval=epsilonClasses.size();
     double objective=fabs(retval-nref);
     if (objective<bestObjective)
     {
         bestObjective=objective;
         bestEpsilon=epsilon;
         bestEpsilonClasses=epsilonClasses;
     }
     std::cout << "Trying " << epsilon << " -> " << retval << std::endl;
     return retval;
 }

◆ writeResults()

void ProgClassifyFTTRI::writeResults ( bool FTTRI )

Write results

Definition at line 1038 of file mpi_classify_FTTRI.cpp.

 {
     int imax=bestEpsilonClasses.size();
     MetaDataVec MDclass, MDsummary;
     FileName fnImg, fnCandidate;
     FileName fnClasses=fnRoot+"_classes.xmd";
     if (fnClasses.exists())
         fnClasses.deleteFile();
 
     // Sort each class
     Image<double> candidate;
     MultidimArray<double> centroid;
     MultidimArray<double> distance;
     MultidimArray<int> idx;
     Matrix2D<double> M;
     AlignmentAux aux;
     CorrelationAux aux2;
     RotationalCorrelationAux aux3;
     for (int i=0; i<imax; i++)
     {
         MDclass.clear();
         const std::vector<size_t> &class_i=bestEpsilonClasses[i].memberIdx;
         size_t nmax=class_i.size();
 
         // Write this class in summary
         size_t id=MDsummary.addObject();
         MDsummary.setValue(MDL_REF, i+1, id);
         MDsummary.setValue(MDL_CLASS_COUNT,nmax,id);
         MDsummary.setValue(MDL_CLASSIFICATION_INTRACLASS_DISTANCE,classEpsilon(i),id);
 
         // Measure the distance of each class and sort
         if (FTTRI)
         {
             centroid.aliasImageInStack(fttriCentroids(),i);
             distance.resizeNoCopy(nmax);
             for (size_t n=0; n<nmax; n++)
             {
                 fnCandidate.compose(class_i[n]+1,fnFTTRI);
                 candidate.read(fnCandidate);
                 A1D_ELEM(distance,n)=fttri_distance(centroid,candidate());
             }
         }
         else
         {
             centroid.aliasImageInStack(imageCentroids(),i);
             centroid.setXmippOrigin();
             distance.resizeNoCopy(nmax);
             for (size_t n=0; n<nmax; n++)
             {
                 int trueIdx=class_i[n];
                 mdIn.getValue(MDL_IMAGE,fnCandidate,imgsId[trueIdx]);
                 candidate.read(fnCandidate);
                 candidate().setXmippOrigin();
                 A1D_ELEM(distance,n)=1-alignImages(centroid,candidate(),M,xmipp_transformation::WRAP,aux,aux2,aux3);
             }
         }
         distance.indexSort(idx);
 
         // Write the class
         MDRowVec row;
         for (size_t n=0; n<nmax; n++)
         {
             int trueIdx=A1D_ELEM(idx,n)-1;
             size_t trueId=imgsId[class_i[trueIdx]];
             mdIn.getRow(row,trueId);
             size_t idClass=MDclass.addRow(row);
             MDclass.setValue(MDL_COST,A1D_ELEM(distance,trueIdx),idClass);
             MDclass.setValue(MDL_REF,i+1,idClass);
             mdIn.setValue(MDL_REF,i+1,trueId);
         }
         MDclass.write(formatString("class%06d_images@%s",i+1,fnClasses.c_str()),MD_APPEND);
     }
     MDsummary.write(formatString("classes@%s",fnClasses.c_str()),MD_APPEND);
     mdIn.write(fnRoot+"_images.xmd");
 }

Member Data Documentation

◆ bestEpsilon

double ProgClassifyFTTRI::bestEpsilon

Definition at line 108 of file mpi_classify_FTTRI.h.

◆ bestEpsilonClasses

std::vector< EpsilonClass > ProgClassifyFTTRI::bestEpsilonClasses

Definition at line 110 of file mpi_classify_FTTRI.h.

◆ bestObjective

double ProgClassifyFTTRI::bestObjective

Definition at line 106 of file mpi_classify_FTTRI.h.

◆ classEpsilon

Matrix1D<double> ProgClassifyFTTRI::classEpsilon

Definition at line 112 of file mpi_classify_FTTRI.h.

◆ dMax

double ProgClassifyFTTRI::dMax

Definition at line 104 of file mpi_classify_FTTRI.h.

◆ dMin

double ProgClassifyFTTRI::dMin

Definition at line 104 of file mpi_classify_FTTRI.h.

◆ doPhase

bool ProgClassifyFTTRI::doPhase

Do phase optimization.

Definition at line 84 of file mpi_classify_FTTRI.h.

◆ epsilonClasses

std::vector< EpsilonClass > ProgClassifyFTTRI::epsilonClasses

Definition at line 102 of file mpi_classify_FTTRI.h.

◆ epsilonMax

double ProgClassifyFTTRI::epsilonMax

Definition at line 114 of file mpi_classify_FTTRI.h.

◆ fmax

double ProgClassifyFTTRI::fmax

Maximum frequency (normalized to 0.5)

Definition at line 60 of file mpi_classify_FTTRI.h.

◆ fnFTTRI

FileName ProgClassifyFTTRI::fnFTTRI

Definition at line 87 of file mpi_classify_FTTRI.h.

◆ fnIn

FileName ProgClassifyFTTRI::fnIn

Input selfile.

Definition at line 51 of file mpi_classify_FTTRI.h.

◆ fnRoot

FileName ProgClassifyFTTRI::fnRoot

Output rootname.

Definition at line 54 of file mpi_classify_FTTRI.h.

◆ fttriCentroids

Image<double> ProgClassifyFTTRI::fttriCentroids

Definition at line 120 of file mpi_classify_FTTRI.h.

◆ FTTRIXdim

int ProgClassifyFTTRI::FTTRIXdim

Definition at line 89 of file mpi_classify_FTTRI.h.

◆ FTTRIYdim

int ProgClassifyFTTRI::FTTRIYdim

Definition at line 90 of file mpi_classify_FTTRI.h.

◆ imageCentroids

Image<double> ProgClassifyFTTRI::imageCentroids

Definition at line 122 of file mpi_classify_FTTRI.h.

◆ imgsId

std::vector<size_t> ProgClassifyFTTRI::imgsId

Definition at line 98 of file mpi_classify_FTTRI.h.

◆ mdIn

MetaDataVec ProgClassifyFTTRI::mdIn

Definition at line 96 of file mpi_classify_FTTRI.h.

◆ Niter

int ProgClassifyFTTRI::Niter

Number of iterations.

Definition at line 81 of file mpi_classify_FTTRI.h.

◆ nMinImages

size_t ProgClassifyFTTRI::nMinImages

Minimum number of images in a class.

Definition at line 78 of file mpi_classify_FTTRI.h.

◆ node

std::shared_ptr<MpiNode> ProgClassifyFTTRI::node

Definition at line 92 of file mpi_classify_FTTRI.h.

◆ notAssigned

Matrix1D<unsigned char> ProgClassifyFTTRI::notAssigned

Definition at line 116 of file mpi_classify_FTTRI.h.

◆ notAssigned0

Matrix1D<unsigned char> ProgClassifyFTTRI::notAssigned0

Definition at line 118 of file mpi_classify_FTTRI.h.

◆ nref

size_t ProgClassifyFTTRI::nref

Desired number of classes.

Definition at line 75 of file mpi_classify_FTTRI.h.

◆ pad

double ProgClassifyFTTRI::pad

Padding factor.

Definition at line 57 of file mpi_classify_FTTRI.h.

◆ padXdim

int ProgClassifyFTTRI::padXdim

Definition at line 100 of file mpi_classify_FTTRI.h.

◆ Rmax

size_t ProgClassifyFTTRI::Rmax

Maximum frequency in pixels.

Definition at line 63 of file mpi_classify_FTTRI.h.

◆ sigma1

double ProgClassifyFTTRI::sigma1

First weight.

Definition at line 69 of file mpi_classify_FTTRI.h.

◆ sigma2

double ProgClassifyFTTRI::sigma2

Second weight.

Definition at line 72 of file mpi_classify_FTTRI.h.

◆ taskDistributor

MpiTaskDistributor* ProgClassifyFTTRI::taskDistributor

Definition at line 94 of file mpi_classify_FTTRI.h.

◆ zoom

double ProgClassifyFTTRI::zoom

Zoom factor for polar conversion.

Definition at line 66 of file mpi_classify_FTTRI.h.

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

xmipp/libraries/parallel/mpi_classify_FTTRI.h
xmipp/libraries/parallel/mpi_classify_FTTRI.cpp

Public Member Functions

Public Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ ProgClassifyFTTRI()

Member Function Documentation

◆ alignImagesWithinClasses()

◆ computeClassCentroids()

◆ computeClassNeighbours()

◆ defineParams()

◆ epsilonClassification()

◆ estimateEpsilonInitialRange()

◆ findFarthest()

◆ findFarthestFTTRI()

◆ fttri_distance()

◆ produceFTTRI()

◆ produceSideInfo()

◆ readParams()

◆ reassignImagesToClasses()

◆ removeSmallClasses()

◆ run()

◆ searchOptimalEpsilon()

◆ show()

◆ skipRandomNumberOfUnassignedClasses()

◆ splitLargeClasses()

◆ splitLargeFTTRIClasses()

◆ wrapperFitness()

◆ writeResults()

Member Data Documentation

◆ bestEpsilon

◆ bestEpsilonClasses

◆ bestObjective

◆ classEpsilon

◆ dMax

◆ dMin

◆ doPhase

◆ epsilonClasses

◆ epsilonMax

◆ fmax

◆ fnFTTRI

◆ fnIn

◆ fnRoot

◆ fttriCentroids

◆ FTTRIXdim

◆ FTTRIYdim

◆ imageCentroids

◆ imgsId

◆ mdIn

◆ Niter

◆ nMinImages

◆ node

◆ notAssigned

◆ notAssigned0

◆ nref

◆ pad

◆ padXdim

◆ Rmax

◆ sigma1

◆ sigma2

◆ taskDistributor

◆ zoom