#include <common_lines.h>

Inheritance diagram for ProgCommonLine:

Collaboration diagram for ProgCommonLine:

Public Member Functions
	ProgCommonLine ()

void	readParams ()

void	defineParams ()

void	produceSideInfo ()

void	processBlock (int i, int j)

void	getAndPrepareBlock (int i, std::vector< MultidimArray< std::complex< double > > > &blockRTFs, std::vector< MultidimArray< double > > &blockRTs)

void	show ()

void	qualifyCommonLines ()

void	solveForShifts ()

void	writeResults ()

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	fn_sel
	input file More...

FileName	fn_out
	output file More...

String	outputStyle
	Output style. More...

bool	scaleDistance
	Scale output measure. More...

double	outlierFraction
	Outlier fraction. More...

double	lpf
	Low pass filter. More...

double	hpf
	High pass filter. More...

double	stepAng
	Angular sampling. More...

double	mem
	Memory limit. More...

int	Nthr
	Number of threads. More...

int	Nmpi
	Number of processors. More...

int	rank
	MPI Rank. More...

int	Nblock

MetaDataVec	SF

int	Nimg

size_t	Xdim

std::vector< CommonLine >	CLmatrix

Matrix1D< double >	shift

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

CommonLine Parameters.

Definition at line 60 of file common_lines.h.

Constructor & Destructor Documentation

◆ ProgCommonLine()

ProgCommonLine::ProgCommonLine ( )

inline

Empty constructor

Definition at line 90 of file common_lines.h.

90 : rank(0) {};

ProgCommonLine::rank

int rank

MPI Rank.

Definition: common_lines.h:87

Member Function Documentation

◆ defineParams()

void ProgCommonLine::defineParams ( )

virtual

Define parameters

Reimplemented from XmippProgram.

Definition at line 50 of file common_lines.cpp.

                                   {
     addUsageLine(
             "For every pair of images in the metadata, find the common line");
     addParamsLine("    -i <file_in>         : input selfile");
     addParamsLine("   [-o <file_out=\"commonlines.txt\">]  : Output filename");
     addParamsLine("   [--ostyle <style=matrix>]  : Output style");
     addParamsLine("          where <style>");
     addParamsLine("                matrix: Text file with the indexes of the corresponding common lines");
     addParamsLine("                line_entries: Text file with a line for each matrix entry");
     addParamsLine("   [--lpf <f=0.01>]      : low pass frequency (0<lpf<0.5)");
     addParamsLine("   [--hpf <f=0.35>]      : high pass frequency (lpf<hpf<0.5)");
     addParamsLine("                         : Set both frequencies to -1 for no filter");
     addParamsLine("   [--stepAng <s=1>]     : angular step");
     addParamsLine("   [--mem <m=1>]         : float number with the memory available in Gb");
     addParamsLine("   [--thr <thr=1>]       : number of threads for each process");
     addParamsLine("   [--scaleDistance]     : scale the output distance to 16-bit integers");
     addParamsLine("   [--outlierFraction <f=0.25>] : Fraction of expected outliers in the detection of common lines");
 }

◆ getAndPrepareBlock()

void ProgCommonLine::getAndPrepareBlock	(	int	i,
		std::vector< MultidimArray< std::complex< double > > > &	blockRTFs,
		std::vector< MultidimArray< double > > &	blockRTs
	)

Get and prepare block

Definition at line 206 of file common_lines.cpp.

                                                      {
     // Get the selfile
     MetaDataVec SFi;
     SFi.selectPart(SF, i * Nblock, Nblock);
 
     // Ask for space for all the block images
     int jmax = SFi.size();
     MultidimArray<std::complex<double> > dummyF;
     MultidimArray<double> dummy;
     for (int j = 0; j < jmax; j++)
     {
         blockRTFs.push_back(dummyF);
         blockRTs.push_back(dummy);
     }
 
     // Read and preprocess the images
     auto * th_ids = new pthread_t[Nthr];
     auto * th_args = new ThreadPrepareImages[Nthr];
     for (int nt = 0; nt < Nthr; nt++) {
         // Passing parameters to each thread
         th_args[nt].parent = this;
         th_args[nt].myThreadID = nt;
         th_args[nt].SFi = &SFi;
         th_args[nt].blockRTFs = &blockRTFs;
         th_args[nt].blockRTs = &blockRTs;
         pthread_create((th_ids + nt), nullptr, threadPrepareImages,
                 (void *) (th_args + nt));
     }
 
     // Waiting for threads to finish
     for (int nt = 0; nt < Nthr; nt++)
         pthread_join(*(th_ids + nt), nullptr);
 
     // Threads structures are not needed any more
     delete []th_ids;
     delete []th_args;
 }

◆ processBlock()

void ProgCommonLine::processBlock	(	int	i,
		int	j
	)

Process block

Definition at line 345 of file common_lines.cpp.

 {
     if (i > j)
         return;
 
     // Preprocess each one of the selfiles
     std::vector<MultidimArray<std::complex<double> > > RTFsi, RTFsj;
     std::vector<MultidimArray<double> > RTsi, RTsj;
     getAndPrepareBlock(i, RTFsi, RTsi);
     if (i != j)
         getAndPrepareBlock(j, RTFsj, RTsj);
 
     // Compare all versus all
     // Read and preprocess the images
     auto * th_ids = new pthread_t[Nthr];
     auto * th_args = new ThreadCompareImages[Nthr];
     for (int nt = 0; nt < Nthr; nt++) {
         // Passing parameters to each thread
         th_args[nt].parent = this;
         th_args[nt].myThreadID = nt;
         th_args[nt].i = i;
         th_args[nt].j = j;
         th_args[nt].RTFsi = &RTFsi;
         th_args[nt].RTsi = &RTsi;
         if (i != j)
         {
             th_args[nt].RTFsj = &RTFsj;
             th_args[nt].RTsj = &RTsj;
         }
         else
         {
             th_args[nt].RTFsj = &RTFsi;
             th_args[nt].RTsj = &RTsi;
         }
         pthread_create((th_ids + nt), nullptr, threadCompareImages,
                 (void *) (th_args + nt));
     }
 
     // Waiting for threads to finish
     for (int nt = 0; nt < Nthr; nt++)
         pthread_join(*(th_ids + nt), nullptr);
 
     // Threads structures are not needed any more
     delete[] th_ids;
     delete[] th_args;
 }

◆ produceSideInfo()

void ProgCommonLine::produceSideInfo ( )

Produce Side Information

Definition at line 70 of file common_lines.cpp.

 {
     SF.read(fn_sel);
     Nimg = SF.size();
 
     // Compute the number of images in each block
     size_t Ydim, Zdim, Ndim;
     getImageSize(SF, Xdim, Ydim, Zdim, Ndim);
     Nblock = FLOOR(sqrt(mem*pow(2.0,30.0)/(2*Ydim*(360/stepAng)*sizeof(double))));
     Nblock = XMIPP_MIN(Nblock,CEIL(((float)Nimg)/Nmpi));
 
     // Ask for memory for the common line matrix
     CommonLine dummy;
     for (int i = 0; i < Nimg * Nimg; i++)
         CLmatrix.push_back(dummy);
 }

◆ qualifyCommonLines()

void ProgCommonLine::qualifyCommonLines ( )

Qualify common lines

Definition at line 393 of file common_lines.cpp.

 {
     std::vector<double> distance;
     size_t imax=CLmatrix.size();
     distance.reserve(imax);
     double dmax=imax;
     for (size_t i=0; i<imax; ++i)
         distance.push_back(CLmatrix[i].distanceij);
 
     std::sort(distance.begin(),distance.end());
 
     std::vector<double>::iterator dbegin=distance.begin();
     std::vector<double>::iterator low;
     for (size_t i=0; i<imax; ++i)
     {
         low=std::lower_bound(dbegin, distance.end(), CLmatrix[i].distanceij);
         CLmatrix[i].percentile=1.0-(double)(low-dbegin)/dmax;
     }
 }

◆ readParams()

void ProgCommonLine::readParams ( )

virtual

Read parameters from command line.

Reimplemented from XmippProgram.

Definition at line 35 of file common_lines.cpp.

                                 {
     fn_sel = getParam("-i");
     fn_out = getParam("-o");
     outputStyle = getParam("--ostyle");
     lpf = getDoubleParam("--lpf");
     hpf = getDoubleParam("--hpf");
     stepAng = getDoubleParam("--stepAng");
     mem = getDoubleParam("--mem");
     Nthr = getIntParam("--thr");
     scaleDistance = checkParam("--scaleDistance");
     outlierFraction = getDoubleParam("--outlierFraction");
     Nmpi = 1;
 }

◆ run()

void ProgCommonLine::run ( )

virtual

Run

Reimplemented from XmippProgram.

Definition at line 514 of file common_lines.cpp.

 {
     produceSideInfo();
 
     // Process all blocks
     int maxBlock = CEIL(((float)Nimg)/Nblock);
     if (rank == 0)
         init_progress_bar(maxBlock * maxBlock);
     for (int i = 0; i < maxBlock; i++)
         for (int j = 0; j < maxBlock; j++)
         {
             int numBlock = i * maxBlock + j;
             if ((numBlock + 1) % Nmpi == rank)
                 processBlock(i, j);
             if (rank == 0)
                 progress_bar(i * maxBlock + j);
         }
     if (rank == 0)
     {
         progress_bar(maxBlock * maxBlock);
         qualifyCommonLines();
         solveForShifts();
         writeResults();
     }
 }

◆ show()

void ProgCommonLine::show ( )

Show parameters

Definition at line 88 of file common_lines.cpp.

                           {
     std::cout << "File in:       " << fn_sel << std::endl << "File out:      "
             << fn_out << std::endl << "Lowpass:       " << lpf << std::endl
             << "Highpass:      " << hpf << std::endl << "StepAng:       "
             << stepAng << std::endl << "Memory(Gb):    " << mem << std::endl
             << "Block size:    " << Nblock << std::endl << "N.Threads:     "
             << Nthr << std::endl;
 }

◆ solveForShifts()

void ProgCommonLine::solveForShifts ( )

Solve for shifts

Definition at line 478 of file common_lines.cpp.

 {
     WeightedLeastSquaresHelper solver;
     int nmax=Nimg*(Nimg-1)/2;
     solver.A.resize(nmax,Nimg*2);
     solver.w.resize(nmax);
     solver.b.resize(nmax);
     int n=0;
     for (int j = 1; j < Nimg; j++)
         for (int i = 0; i < j; i++, n++) {
             int ii = i * Nimg + j;
             if (CLmatrix[ii].distanceij > 0) {
                 const CommonLine& cl=CLmatrix[ii];
                 double sini, cosi, sinj, cosj;
                 //sincos(DEG2RAD(cl.angi),&sini,&cosi);
                 sini = sin(DEG2RAD(cl.angi));
                 cosi = cos(DEG2RAD(cl.angi));
                 //sincos(DEG2RAD(cl.angj),&sinj,&cosj);
                 sinj = sin(DEG2RAD(cl.angj));
                 cosj = cos(DEG2RAD(cl.angj));
                 int idxi=2*i;
                 int idxj=2*j;
                 MAT_ELEM(solver.A,n,idxi)=cosi;
                 MAT_ELEM(solver.A,n,idxi+1)=sini;
                 MAT_ELEM(solver.A,n,idxj)=-cosj;
                 MAT_ELEM(solver.A,n,idxj+1)=-sinj;
                 VEC_ELEM(solver.b,n)=cl.jmax;
                 VEC_ELEM(solver.w,n)=cl.percentile;
             }
         }
 
     const double maxShiftError=1;
     ransacWeightedLeastSquares(solver,shift,maxShiftError,10000,outlierFraction,Nthr);
 }

◆ writeResults()

void ProgCommonLine::writeResults ( )

Write results

Definition at line 414 of file common_lines.cpp.

 {
     // Look for the minimum and maximum of the common line matrix
     double minVal = 2, maxVal = -2;
     if (scaleDistance)
     {
         for (int i = 0; i < Nimg; i++)
             for (int j = 0; j < Nimg; j++)
             {
                 double val = CLmatrix[i * Nimg + j].distanceij;
                 if (val > 0)
                 {
                     minVal = XMIPP_MIN(minVal,val);
                     maxVal = XMIPP_MAX(maxVal,val);
                 }
             }
     }
 
     // Write the common line matrix
     if (outputStyle == "matrix") {
         Matrix2D<int> CL(Nimg,Nimg);
         CL.initConstant(-1);
         for (int j = 1; j < Nimg; j++)
             for (int i = 0; i < j; i++) {
                 int ii = i * Nimg + j;
                 MAT_ELEM(CL,i,j)= (int)round(CLmatrix[ii].angi/stepAng);
                 MAT_ELEM(CL,j,i)= (int)round(CLmatrix[ii].angj/stepAng);
             }
         CL.write(fn_out);
     } else {
         std::ofstream fh_out;
         fh_out.open(fn_out.c_str());
         if (!fh_out)
             REPORT_ERROR(ERR_IO_NOWRITE, fn_out);
         for (int j = 1; j < Nimg; j++)
             for (int i = 0; i < j; i++) {
                 int ii = i * Nimg + j;
                 if (CLmatrix[ii].distanceij > 0) {
                     fh_out << j << " " << i << " ";
                     if (scaleDistance)
                         fh_out << round(65535*(CLmatrix[ii].distanceij-minVal)/
                                 (maxVal-minVal)) << " ";
                     else
                         fh_out << CLmatrix[ii].distanceij << " ";
                     fh_out << round(CLmatrix[ii].angi/stepAng) << " "
                            << round(CLmatrix[ii].angj/stepAng) << " "
                            << CLmatrix[ii].jmax << " "
                            << CLmatrix[ii].percentile << std::endl;
                 }
             }
         fh_out.close();
     }
 
     // Write the aligned images
     int idx=0;
     for (size_t objId : SF.ids())
     {
         SF.setValue(MDL_SHIFT_X,-shift(idx++)/2, objId); // *** FIXME: COSS Why /2?
         SF.setValue(MDL_SHIFT_Y,-shift(idx++)/2, objId);
     }
     SF.write(fn_out.insertBeforeExtension("_aligned_images"));
 }

Member Data Documentation

◆ CLmatrix

std::vector<CommonLine> ProgCommonLine::CLmatrix

Definition at line 137 of file common_lines.h.

◆ fn_out

FileName ProgCommonLine::fn_out

output file

Definition at line 66 of file common_lines.h.

◆ fn_sel

FileName ProgCommonLine::fn_sel

input file

Definition at line 64 of file common_lines.h.

◆ hpf

double ProgCommonLine::hpf

High pass filter.

Definition at line 76 of file common_lines.h.

◆ lpf

double ProgCommonLine::lpf

Low pass filter.

Definition at line 74 of file common_lines.h.

◆ mem

double ProgCommonLine::mem

Memory limit.

Definition at line 81 of file common_lines.h.

◆ Nblock

int ProgCommonLine::Nblock

Definition at line 125 of file common_lines.h.

◆ Nimg

int ProgCommonLine::Nimg

Definition at line 131 of file common_lines.h.

◆ Nmpi

int ProgCommonLine::Nmpi

Number of processors.

Definition at line 85 of file common_lines.h.

◆ Nthr

int ProgCommonLine::Nthr

Number of threads.

Definition at line 83 of file common_lines.h.

◆ outlierFraction

double ProgCommonLine::outlierFraction

Outlier fraction.

Definition at line 72 of file common_lines.h.

◆ outputStyle

String ProgCommonLine::outputStyle

Output style.

Definition at line 68 of file common_lines.h.

◆ rank

int ProgCommonLine::rank

MPI Rank.

Definition at line 87 of file common_lines.h.

◆ scaleDistance

bool ProgCommonLine::scaleDistance

Scale output measure.

Definition at line 70 of file common_lines.h.

◆ SF

MetaDataVec ProgCommonLine::SF

Definition at line 128 of file common_lines.h.

◆ shift

Matrix1D<double> ProgCommonLine::shift

Definition at line 140 of file common_lines.h.

◆ stepAng

double ProgCommonLine::stepAng

Angular sampling.

Definition at line 78 of file common_lines.h.

◆ Xdim

size_t ProgCommonLine::Xdim

Definition at line 134 of file common_lines.h.

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

xmipp/legacy/libraries/reconstruction/common_lines.h
xmipp/legacy/libraries/reconstruction/common_lines.cpp

Public Member Functions

Public Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ ProgCommonLine()

Member Function Documentation

◆ defineParams()

◆ getAndPrepareBlock()

◆ processBlock()

◆ produceSideInfo()

◆ qualifyCommonLines()

◆ readParams()

◆ run()

◆ show()

◆ solveForShifts()

◆ writeResults()

Member Data Documentation

◆ CLmatrix

◆ fn_out

◆ fn_sel

◆ hpf

◆ lpf

◆ mem

◆ Nblock

◆ Nimg

◆ Nmpi

◆ Nthr

◆ outlierFraction

◆ outputStyle

◆ rank

◆ scaleDistance

◆ SF

◆ shift

◆ stepAng

◆ Xdim