#include <tomo_align_tilt_series.h>

Inheritance diagram for ProgTomographAlignment:

Collaboration diagram for ProgTomographAlignment:

Public Member Functions
	~ProgTomographAlignment ()
	Destructor. More...

void	defineParams ()
	Usage. More...

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

void	show ()
	Show parameters. More...

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

void	computeAffineTransformations (bool globalAffineToUse)
	Compute affine transformations. More...

void	identifyOutliers (bool mark)
	Identify outliers. More...

void	produceInformationFromLandmarks ()
	Produce information from landmarks. More...

bool	refineLandmark (int ii, int jj, const Matrix1D< double > &rii, Matrix1D< double > &rjj, double &maxCorr, bool tryFourier) const

bool	refineLandmark (const MultidimArray< double > &pieceii, int jj, Matrix1D< double > &rjj, double actualCorrThreshold, bool reversed, double &maxCorr) const

bool	refineChain (LandmarkChain &chain, double &corrChain)

void	generateLandmarkSet ()
	Generate landmark set using a grid. More...

void	writeLandmarkSet (const FileName &fnLandmark) const
	Write landmark set. More...

void	writeTransformations (const FileName &fnTransformations) const
	Write affine transformations. More...

void	readLandmarkSet (const FileName &fnLandmark)
	Read landmark set. More...

void	removeOutlierLandmarks (const Alignment &alignment)
	Remove Outliers. More...

void	alignImages (const Alignment &alignment)
	Align images. More...

void	run ()
	Run: do the real work. More...

Public Member Functions inherited from XmippProgram
const char *	getParam (const char *param, int arg=0)

const char *	getParam (const char param, const char subparam, int arg=0)

int	getIntParam (const char *param, int arg=0)

int	getIntParam (const char param, const char subparam, int arg=0)

double	getDoubleParam (const char *param, int arg=0)

double	getDoubleParam (const char param, const char subparam, int arg=0)

float	getFloatParam (const char *param, int arg=0)

float	getFloatParam (const char param, const char subparam, int arg=0)

void	getListParam (const char *param, StringVector &list)

int	getCountParam (const char *param)

bool	checkParam (const char *param)

bool	existsParam (const char *param)

void	addParamsLine (const String &line)

void	addParamsLine (const char *line)

ParamDef *	getParamDef (const char *param) const

virtual void	quit (int exit_code=0) const

virtual int	tryRun ()

void	initProgress (size_t total, size_t stepBin=60)

void	setProgress (size_t value=0)

void	endProgress ()

void	processDefaultComment (const char param, const char left)

void	setDefaultComment (const char param, const char comment)

virtual void	initComments ()

void	setProgramName (const char *name)

void	addUsageLine (const char *line, bool verbatim=false)

void	clearUsage ()

void	addExampleLine (const char *example, bool verbatim=true)

void	addSeeAlsoLine (const char *seeAlso)

void	addKeywords (const char *keywords)

const char *	name () const

virtual void	usage (int verb=0) const

virtual void	usage (const String &param, int verb=2)

int	version () const

virtual void	show () const

virtual void	read (int argc, const char **argv, bool reportErrors=true)

virtual void	read (int argc, char **argv, bool reportErrors=true)

void	read (const String &argumentsLine)

	XmippProgram ()

	XmippProgram (int argc, const char **argv)

virtual	~XmippProgram ()

Public Attributes
FileName	fnSel
	MetaData File with all images. More...

FileName	fnSelOrig
	MetaData File with all images at the original scale. More...

FileName	fnRoot
	Output root. More...

int	numThreads
	Number of threads to use for parallel computing. More...

bool	globalAffine
	Look for local affine transformation. More...

bool	useCriticalPoints
	Use critical points. More...

size_t	Ncritical
	Number of critical points. More...

size_t	seqLength
	Sequence Length. More...

size_t	blindSeqLength

int	gridSamples
	Grid samples. More...

double	maxShiftPercentage
	Maxshift percentage. More...

int	maxIterDE
	Maximum number of iterations in Differential Evolution. More...

double	psiMax
	Maximum rotation. More...

int	maxStep
	Maximum Step for chain refinement. More...

double	deltaRot
	Delta rot. More...

double	localSize
	Local refinement size. More...

bool	optimizeTiltAngle
	Optimize tilt angle of tilt axis. More...

bool	isCapillar
	This tilt series comes from a capillar. More...

bool	dontNormalize
	Don't normalize the corrected images. More...

bool	difficult
	Difficult. More...

double	corrThreshold
	Correlation threshold for a good landmark. More...

bool	showAffine
	Show affine transformations. More...

double	thresholdAffine
	Threshold affine. More...

double	identifyOutliersZ
	Identify outlier micrographs. More...

bool	doNotIdentifyOutliers
	Do not identify outlier micrographs. More...

int	pyramidLevel
	Pyramid. More...

int	lastStep
	Last step to run (-1=run them all) More...

int	iteration

MetaDataVec	SF

MetaDataVec	SForig

std::vector< MultidimArray< unsigned char > * >	img

std::vector< MultidimArray< unsigned char > * >	maskImg

int	iMinTilt

std::vector< double >	tiltList

std::vector< std::string >	name_list

Matrix1D< double >	avgForwardPatchCorr

Matrix1D< double >	avgBackwardPatchCorr

Matrix1D< int >	isOutlier

int	Nimg

std::vector< std::vector< Matrix2D< double > > >	affineTransformations

std::vector< double >	correlationList

Matrix2D< double >	allLandmarksX

Matrix2D< double >	allLandmarksY

std::vector< std::vector< int > >	Vseti

std::vector< std::vector< int > >	Vsetj

std::vector< Matrix1D< double > >	barpi

Matrix1D< int >	ni

Alignment *	bestPreviousAlignment

bool	showRefinement

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

This is the main class

Definition at line 76 of file tomo_align_tilt_series.h.

Constructor & Destructor Documentation

◆ ~ProgTomographAlignment()

ProgTomographAlignment::~ProgTomographAlignment ( )

Destructor.

Definition at line 1621 of file tomo_align_tilt_series.cpp.

 {
     // Clear the list of images if not empty
     if (!img.empty())
     {
         for (size_t i=0; i<img.size(); i++)
             delete img[i];
         img.clear();
     }
 
     if (!maskImg.empty())
     {
         for (size_t i=0; i<maskImg.size(); i++)
             delete maskImg[i];
         maskImg.clear();
     }
 }

Member Function Documentation

◆ alignImages()

void ProgTomographAlignment::alignImages ( const Alignment & alignment )

Align images.

Definition at line 2245 of file tomo_align_tilt_series.cpp.

 {
     // Correct all landmarks
     Matrix1D<double> r(2);
     for (size_t i=0; i<MAT_XSIZE(allLandmarksX); i++)
     {
         Matrix2D<double> R;
         rotation2DMatrix(90-alignment.rot+alignment.psi(i),R,false);
         for (size_t j=0; j<MAT_YSIZE(allLandmarksX); j++)
             if (allLandmarksX(j,i)!=XSIZE(*img[0]))
             {
                 VECTOR_R2(r,allLandmarksX(j,i),allLandmarksY(j,i));
                 r=R*(r-alignment.di[i]+alignment.diaxis[i]);
                 allLandmarksX(j,i)=XX(r);
                 allLandmarksY(j,i)=YY(r);
             }
     }
 
     // Compute the average height of the 3D landmarks seen at 0 degrees
     Matrix1D<double> axis;
     Euler_direction(alignment.rot, alignment.tilt, 0, axis);
     double z0=0, z0N=0;
     for (size_t j=0; j<MAT_YSIZE(allLandmarksX); j++)
         if (allLandmarksX(j,iMinTilt)!=XSIZE(*img[0]))
         {
             Matrix2D<double> Raxismin, Rmin, RtiltYmin;
             rotation3DMatrix(tiltList[iMinTilt],axis,Raxismin,false);
             rotation2DMatrix(90-alignment.rot+alignment.psi(iMinTilt),Rmin);
             rotation3DMatrix(-tiltList[iMinTilt],'Y',RtiltYmin,false);
             Matrix1D<double> rjp=RtiltYmin*Rmin*Raxismin*alignment.rj[j];
             z0+=ZZ(rjp);
             z0N++;
         }
     if (z0N==0)
         REPORT_ERROR(ERR_VALUE_INCORRECT,"There is no landmark at 0 degrees");
     z0/=z0N;
     std::cout << "Average height of the landmarks at 0 degrees=" << z0 << std::endl;
     MetaDataVec DF;
 
     std::unique_ptr<MetaDataVec::id_iterator> idIter;
 
     if (!fnSelOrig.empty())
     {
         idIter = memoryUtils::make_unique<MetaDataVec::id_iterator>(SForig.ids().begin());
     }
     DF.setComment("First shift by -(shiftX,shiftY), then rotate by psi");
 
     MultidimArray<double> mask;
     MultidimArray<int> iMask;
     Image<double> I;
     Matrix2D<double> M, M1, M2, M3;
     FileName fn_corrected;
 
     for (int n=0;n<Nimg; n++)
     {
         // Align the normal image
         I.read(name_list[n]);
         mask.initZeros(I());
         FOR_ALL_ELEMENTS_IN_ARRAY2D(I())
         if (I(i,j)!=0)
             mask(i,j)=1;
         translation2DMatrix(vectorR2(-z0*sin(DEG2RAD(tiltList[n])),0),M1);
         rotation2DMatrix(90-alignment.rot+alignment.psi(n),M2);
         translation2DMatrix(-(alignment.di[n]+alignment.diaxis[n]),M3);
         M=M1*M2*M3;
         selfApplyGeometry(xmipp_transformation::BSPLINE3,I(),M,xmipp_transformation::IS_NOT_INV,xmipp_transformation::DONT_WRAP);
         selfApplyGeometry(xmipp_transformation::LINEAR,mask,M,xmipp_transformation::IS_NOT_INV,xmipp_transformation::DONT_WRAP);
         mask.binarize(0.5);
         typeCast(mask,iMask);
         double minval, maxval, avg, stddev;
         computeStats_within_binary_mask(iMask,I(),minval, maxval, avg, stddev);
         FOR_ALL_ELEMENTS_IN_ARRAY2D(iMask)
         if (iMask(i,j)==0)
             I(i,j)=0;
         else if (!dontNormalize)
             I(i,j)-=avg;
         double rot=0;
         double tilt=tiltList[n];
         double psi=0;
         I.setEulerAngles(rot, tilt, psi);
         fn_corrected.compose(n+1, fnRoot+"_corrected_", "stk");
         I.write(fn_corrected);
 
         // Align the original image
         if (fnSelOrig!="")
         {
             FileName auxFn;
             SForig.getValue( MDL_IMAGE, auxFn, **idIter);
             Image<double> Iorig;
             Iorig.read( auxFn );
             //SForig.nextObject();
             ++(*idIter);
             mask.initZeros(Iorig());
             FOR_ALL_ELEMENTS_IN_ARRAY2D(Iorig())
             if (Iorig(i,j)!=0)
                 mask(i,j)=1;
             translation2DMatrix(vectorR2(-z0*sin(DEG2RAD(tiltList[n]))*
                                          (((double)XSIZE(Iorig()))/XSIZE(I())),0),M1);
             rotation2DMatrix(90-alignment.rot+alignment.psi(n),M2);
             translation2DMatrix(-(alignment.di[n]+alignment.diaxis[n])*
                                 (((double)XSIZE(Iorig()))/XSIZE(I())),M3);
             M=M1*M2*M3;
             selfApplyGeometry(xmipp_transformation::BSPLINE3,Iorig(),M,xmipp_transformation::IS_NOT_INV,xmipp_transformation::DONT_WRAP);
             selfApplyGeometry(xmipp_transformation::LINEAR,mask,M,xmipp_transformation::IS_NOT_INV,xmipp_transformation::DONT_WRAP);
             mask.binarize(0.5);
             typeCast(mask,iMask);
             computeStats_within_binary_mask(iMask,Iorig(),minval, maxval,
                                             avg, stddev);
             FOR_ALL_ELEMENTS_IN_ARRAY2D(iMask)
             if (iMask(i,j)==0)
                 Iorig(i,j)=0;
             else if (!dontNormalize)
                 Iorig(i,j)-=avg;
             Iorig.setEulerAngles(rot, tilt, psi);
             fn_corrected.compose(n+1, fnRoot+"_corrected_originalsize_", "stk");
             Iorig.write(fn_corrected);
         }
 
         // Prepare data for the docfile
         size_t id = DF.addObject();
         DF.setValue(MDL_IMAGE, fn_corrected, id);
         DF.setValue(MDL_ANGLE_PSI, 90.-alignment.rot+alignment.psi(n), id);
         DF.setValue(MDL_SHIFT_X, XX(alignment.di[n]+alignment.diaxis[n]), id);
         DF.setValue(MDL_SHIFT_Y, YY(alignment.di[n]+alignment.diaxis[n]), id);
     }
     DF.write(fnRoot+"_correction_parameters.txt");
 #ifdef DEBUG
 
     Image<double> save;
     save().initZeros(*img[0]);
     save().setXmippOrigin();
     for (int j=0; j<YSIZE(allLandmarksX); j++)
         if (allLandmarksX(j,iMinTilt)!=XSIZE(*img[0]))
         {
             Matrix2D<double> Raxismin=rotation3DMatrix(tiltList[iMinTilt],axis);
             Raxismin.resize(3,3);
             Matrix2D<double> Rmin=rotation2DMatrix(90-alignment.rot+alignment.psi(iMinTilt));
             Matrix2D<double> RtiltYmin=rotation3DMatrix(-tiltList[iMinTilt],'Y');
             RtiltYmin.resize(3,3);
             Matrix1D<double> rjp=RtiltYmin*Rmin*Raxismin*alignment.rj[j];
             std::cout << rjp.transpose() << std::endl;
             for (int i=0; i<XSIZE(allLandmarksX); i++)
                 if (allLandmarksX(j,i)!=XSIZE(*img[0]))
                 {
                     save(allLandmarksY(j,i),allLandmarksX(j,i))=ABS(i-iMinTilt);
 
                     /*
                     Matrix2D<double> Raxis=rotation3DMatrix(tiltList[i],axis);
                     Raxis.resize(3,3);
                     Matrix2D<double> R=rotation2DMatrix(90-alignment.rot+alignment.psi(i));
                     Matrix1D<double> p=R*Raxis*alignment.rj[j];
                     if (YY(p)>=STARTINGY(save()) && YY(p)<=FINISHINGY(save()) &&
                         XX(p)>=STARTINGX(save()) && XX(p)<=FINISHINGX(save()))
                        save(YY(p),XX(p))=-ABS(i-iMinTilt);
                     */
                     Matrix2D<double> RtiltY=rotation3DMatrix(-tiltList[i],'Y');
                     RtiltY.resize(3,3);
                     Matrix1D<double> p=RtiltY*rjp;
                     if (YY(p)>=STARTINGY(save()) && YY(p)<=FINISHINGY(save()) &&
                         XX(p)>=STARTINGX(save()) && XX(p)<=FINISHINGX(save()))
                         save(YY(p),XX(p))=-ABS(i-iMinTilt);
                 }
         }
     save.write("PPPmovementsLandmarks0.xmp");
 #endif
 }

◆ computeAffineTransformations()

void ProgTomographAlignment::computeAffineTransformations ( bool globalAffineToUse )

Compute affine transformations.

Definition at line 674 of file tomo_align_tilt_series.cpp.

 {
     bool oldglobalAffine=globalAffine;
     globalAffine=globalAffineToUse;
 
     auto * th_ids = new pthread_t[numThreads];
     auto * th_args = new ThreadComputeTransformParams[numThreads];
 
     for( int nt = 0 ; nt < numThreads ; nt ++ )
     {
         // Passing parameters to each thread
         th_args[nt].parent = this;
         th_args[nt].myThreadID = nt;
         pthread_create( (th_ids+nt) , nullptr, threadComputeTransform, (void *)(th_args+nt) );
     }
 
     // Waiting for threads to finish
     for( int nt = 0 ; nt < numThreads ; nt ++ )
         pthread_join(*(th_ids+nt), nullptr);
 
     // Threads structures are not needed any more
     delete[] th_ids;
     delete[] th_args;
     globalAffine=oldglobalAffine;
 }

◆ defineParams()

void ProgTomographAlignment::defineParams ( )

virtual

Usage.

Reimplemented from XmippProgram.

Definition at line 522 of file tomo_align_tilt_series.cpp.

 {
     addUsageLine("Align a single-axis tilt series without any marker.");
     addSeeAlsoLine("tomo_align_dual_tilt_series");
     addParamsLine(" == General Options == ");
     addParamsLine("   -i <metadatafile>              : Input images");
     addParamsLine("                                  : The selfile must contain the list of micrographs");
     addParamsLine("                                  : and its tilt angles");
     addParamsLine("  [--iorig <metadatafile=\"\">]   : Metadata with images at original scale");
     addParamsLine("  [--oroot <fn_out=\"\">]         : Output alignment");
     addParamsLine("                                  : If not given, the input selfile without extension");
     addParamsLine("  [--thr <num=1>]                 : Parallel processing using \"num\" threads");
     addParamsLine("  [--lastStep+ <step=-1>]         : Last step to perform");
     addParamsLine("                                  : Step -1 -> Perform all steps");
     addParamsLine("                                  : Step  0 -> Determination of affine transformations");
     addParamsLine("                                  : Step  1 -> Determination of landmark chains");
     addParamsLine("                                  : Step  2 -> Determination of alignment parameters");
     addParamsLine("                                  : Step  3 -> Writing aligned images");
     addParamsLine(" == Step 0 (Affine alignment) Options == ");
     addParamsLine("  [--maxShiftPercentage <p=0.2>]   : Maximum shift as percentage of image size");
     addParamsLine("  [--thresholdAffine <th=0.85>]    : Threshold affine");
     addParamsLine("  [--globalAffine]                 : Look globally for affine transformations");
     addParamsLine("  [--difficult+]                   : Apply some filters before affine alignment");
     addParamsLine("  [--maxIterDE+ <n=30>]            : Maximum number of iteration in Differential Evolution");
     addParamsLine("  [--showAffine+]                  : Show affine transformations as PPP*");
     addParamsLine("  [--identifyOutliers+ <z=5>]      : Z-score to be an outlier");
     addParamsLine("  [--noOutliers+]                  : Do not identify outliers");
     addParamsLine("  [--pyramid+ <level=1>]           : Multiresolution for affine transformations");
     addParamsLine(" == Step 1 (Landmark chain) Options == ");
     addParamsLine("  [--seqLength <n=5>]              : Sequence length");
     addParamsLine("  [--localSize <size=0.04>]        : In percentage");
     addParamsLine("  [--useCriticalPoints <n=0>]      : Use critical points instead of a grid");
     addParamsLine("                                   : n is the number of critical points to choose");
     addParamsLine("                                   : in each image");
     addParamsLine("  [--threshold+ <th=-1>]           : Correlation threshold");
     addParamsLine("  [--blindSeqLength+ <n=-1>]       : Blind sequence length, -1=No blind landmarks");
     addParamsLine("  [--maxStep+ <step=4>]            : Maximum step for chain refinement");
     addParamsLine("  [--gridSamples+ <n=40>]           : Total number of samples=n*n");
     addParamsLine("  [--isCapillar+]                  : Set this flag if the tilt series is of a capillar");
     addParamsLine(" == Step 2 (Determination of alignment parameters) Options == ");
     addParamsLine("  [--psiMax+ <psi=-1>]             : Maximum psi in absolute value (degrees)");
     addParamsLine("                                  : -1 -> do not optimize for psi");
     addParamsLine("  [--deltaRot+ <rot=5>]            : In degrees. For the first optimization stage");
     addParamsLine("  [--optimizeTiltAngle+]           : Optimize tilt angle");
     addParamsLine(" == Step 3 (Produce aligned images) Options == ");
     addParamsLine("  [--dontNormalize+]               : Don't normalize the output images");
     addExampleLine("Typical run",false);
     addExampleLine("xmipp_tomo_align_tilt_series -i tiltseries.sel --thr 8");
     addExampleLine("If there are image with large shifts",false);
     addExampleLine("xmipp_tomo_align_tilt_series -i tiltseries.sel --thr 8 --globalAffine");
     addExampleLine("If there are clear landmarks that can be tracked",false);
     addExampleLine("xmipp_tomo_align_tilt_series -i tiltseries.sel --thr 8 --criticalPoints");
 }

◆ generateLandmarkSet()

void ProgTomographAlignment::generateLandmarkSet ( )

Generate landmark set using a grid.

Definition at line 1639 of file tomo_align_tilt_series.cpp.

 {
     FileName fn_tmp = fnRoot+"_landmarks.txt";
     if (!fn_tmp.exists())
     {
         auto * th_ids = new pthread_t[numThreads];
         auto * th_args=
             new ThreadGenerateLandmarkSetParams[numThreads];
         for( int nt = 0 ; nt < numThreads ; nt ++ )
         {
             th_args[nt].parent = this;
             th_args[nt].myThreadID = nt;
             if (useCriticalPoints)
                 pthread_create( (th_ids+nt) , nullptr, threadgenerateLandmarkSetCriticalPoints, (void *)(th_args+nt) );
             else
                 pthread_create( (th_ids+nt) , nullptr, threadgenerateLandmarkSetGrid, (void *)(th_args+nt) );
         }
 
         std::vector<LandmarkChain> chainList;
         int includedPoints=0;
         for( int nt = 0 ; nt < numThreads ; nt ++ )
         {
             pthread_join(*(th_ids+nt), nullptr);
             int imax=th_args[nt].chainList->size();
             for (int i=0; i<imax; i++)
             {
                 chainList.push_back((*th_args[nt].chainList)[i]);
                 includedPoints+=(*th_args[nt].chainList)[i].size();
             }
             delete th_args[nt].chainList;
         }
         delete[] th_ids;
         delete[] th_args;
 
         // Add blind landmarks
         if (blindSeqLength>0)
         {
             th_ids  = new pthread_t[numThreads];
             th_args = new ThreadGenerateLandmarkSetParams[numThreads];
             for( int nt = 0 ; nt < numThreads ; nt ++ )
             {
                 th_args[nt].parent = this;
                 th_args[nt].myThreadID = nt;
                 pthread_create( (th_ids+nt) , nullptr, threadgenerateLandmarkSetBlind, (void *)(th_args+nt) );
             }
 
             for( int nt = 0 ; nt < numThreads ; nt ++ )
             {
                 pthread_join(*(th_ids+nt), nullptr);
                 int imax=th_args[nt].chainList->size();
                 for (int i=0; i<imax; i++)
                 {
                     chainList.push_back((*th_args[nt].chainList)[i]);
                     includedPoints+=(*th_args[nt].chainList)[i].size();
                 }
                 delete th_args[nt].chainList;
             }
             delete[] th_ids;
             delete[] th_args;
         }
 
         // Generate the landmark "matrix"
         allLandmarksX.resize(chainList.size(),Nimg);
         allLandmarksY.resize(chainList.size(),Nimg);
         allLandmarksX.initConstant(XSIZE(*img[0]));
         allLandmarksY.initConstant(YSIZE(*img[0]));
         for (size_t i=0; i<chainList.size(); i++)
         {
             for (size_t j=0; j<chainList[i].size(); j++)
             {
                 int idx=chainList[i][j].imgIdx;
                 allLandmarksX(i,idx)=chainList[i][j].x;
                 allLandmarksY(i,idx)=chainList[i][j].y;
             }
         }
 
         // Write landmarks
         if (includedPoints>0)
         {
             writeLandmarkSet(fnRoot+"_landmarks.txt");
             std::cout << " Number of points="
             << includedPoints
             << " Number of chains=" << chainList.size()
             << " ( " << ((double) includedPoints)/chainList.size() << " )\n";
         }
         else
             REPORT_ERROR(ERR_VALUE_INCORRECT,"There are no landmarks meeting this threshold. Try to lower it");
     }
     else
     {
         readLandmarkSet(fnRoot+"_landmarks.txt");
     }
 }

◆ identifyOutliers()

void ProgTomographAlignment::identifyOutliers ( bool mark )

Identify outliers.

Definition at line 701 of file tomo_align_tilt_series.cpp.

 {
     isOutlier.initZeros(Nimg);
     MultidimArray<double> correlationListAux(Nimg);
     for (int i=0; i<Nimg; i++)
         correlationListAux(i)=correlationList[i];
 
     MultidimArray<double> diff;
     double medianCorr=correlationListAux.computeMedian();
     diff=correlationListAux-medianCorr;
     diff.selfABS();
     double madCorr=diff.computeMedian();
 
     std::cout << "Cost distribution= " << 1-medianCorr << " +- "
     << madCorr << std::endl;
 
     double thresholdCorr=medianCorr-identifyOutliersZ*madCorr;
     for (size_t i=1; i<VEC_XSIZE(isOutlier); i++)
     {
         bool potentialOutlier=(correlationListAux(i)<thresholdCorr);
         if (potentialOutlier)
         {
             affineTransformations[i-1][i].clear();
             if (mark)
             {
                 isOutlier(i)=true;
                 std::cout << name_list[i-1] << " [" << i
                 << "] is considered as an outlier. "
                 << "Its cost is " << 1-correlationListAux(i)
                 << std::endl;
             }
             else
                 std::cout << name_list[i-1] << " [" << i
                 << "] might be an outlier. "
                 << "Its cost is " << 1-correlationListAux(i)
                 << std::endl;
         }
         iteration++;
     }
 }

◆ produceInformationFromLandmarks()

void ProgTomographAlignment::produceInformationFromLandmarks ( )

Produce information from landmarks.

Definition at line 2414 of file tomo_align_tilt_series.cpp.

 {
     // Produce V sets
     std::vector<int> emptyVector;
 
     Vseti.clear();
     for (size_t i=0; i<MAT_XSIZE(allLandmarksX); i++)
         Vseti.push_back(emptyVector);
 
     Vsetj.clear();
     for (size_t j=0; j<MAT_YSIZE(allLandmarksX); j++)
         Vsetj.push_back(emptyVector);
 
     for (size_t j=0; j<MAT_YSIZE(allLandmarksX); j++)
         for (size_t i=0; i<MAT_XSIZE(allLandmarksX); i++)
             if (allLandmarksX(j,i)!=XSIZE(*img[0]))
             {
                 Vseti[i].push_back(j);
                 Vsetj[j].push_back(i);
             }
 
     // Count the number of landmarks per image and average projection
     // of all landmarks in a given image
     ni.initZeros(Nimg);
     barpi.clear();
     for (int i=0; i<Nimg; i++)
     {
         ni(i)=static_cast<int>(Vseti[i].size());
         Matrix1D<double> pi(2);
         for (int jj=0; jj<ni(i); jj++)
         {
             int j=Vseti[i][jj];
             XX(pi)+=allLandmarksX(j,i);
             YY(pi)+=allLandmarksY(j,i);
         }
         pi/=ni(i);
         barpi.push_back(pi);
     }
 }

◆ produceSideInfo()

void ProgTomographAlignment::produceSideInfo ( )

Produce side info.

Definition at line 742 of file tomo_align_tilt_series.cpp.

 {
     // Difficult images?
     if (difficult)
     {
         if (pyramidLevel==0)
             pyramidLevel=2;
         if (localSize<0.05)
             localSize=0.08;
         if (seqLength==5)
             seqLength=11;
     }
 
     bestPreviousAlignment=new Alignment(this);
     // Read input data
     SF.read(fnSel,nullptr);
     if (SF.containsLabel(MDL_ENABLED))
         SF.removeObjects(MDValueEQ(MDL_ENABLED, -1));
     Nimg=SF.size();
     if (Nimg!=0)
     {
         // Clear the list of images if not empty
         if (!img.empty())
         {
             for (size_t i=0; i<img.size(); i++)
                 delete img[i];
             img.clear();
         }
 
         // Clear the list of masks if not empty
         if (!maskImg.empty())
         {
             for (size_t i=0; i<maskImg.size(); i++)
                 delete maskImg[i];
             maskImg.clear();
         }
 
         std::cerr << "Reading input data\n";
         init_progress_bar(Nimg);
         iteration=0;
         int n=0;
 
         iMinTilt=-1;
         double minTilt=1000;
         bool nonZeroTilt=false;
         Image<double> imgaux;
         FileName fn;
         for (size_t objId : SF.ids())
         {
             SF.getValue( MDL_IMAGE, fn, objId);
             imgaux.read(fn);
             if (difficult)
             {
                 //ImageXmipp save;
                 //save()=imgaux(); save.write("PPPoriginal.xmp");
                 MultidimArray<double> Ifiltered;
                 Ifiltered=imgaux();
                 Ifiltered.setXmippOrigin();
 
                 // Reject outliers
                 reject_outliers(Ifiltered,0.5);
 
                 // Subtract background plane
                 substractBackgroundPlane(Ifiltered);
 
                 // Subtract the background (rolling ball)
                 substractBackgroundRollingBall(Ifiltered,
                                                XSIZE(Ifiltered)/10);
 
                 // Bandpass the image
                 FourierFilter FilterBP;
                 FilterBP.FilterBand=BANDPASS;
                 FilterBP.w1=1.0/XSIZE(Ifiltered);
                 FilterBP.w2=100.0/XSIZE(Ifiltered);
                 FilterBP.raised_w=1.0/XSIZE(Ifiltered);
                 FilterBP.generateMask(Ifiltered);
                 FilterBP.applyMaskSpace(Ifiltered);
 
                 // Equalize histogram
                 //histogram_equalization(Ifiltered,8);
                 imgaux()=Ifiltered;
                 //save()=Ifiltered; save.write("PPPpreprocessed.xmp");
                 //std::cout << "Press\n";
                 //char c; std::cin >> c;
             }
 
             if (!useCriticalPoints)
             {
                 auto* mask_i=new MultidimArray<unsigned char>;
                 generateMask(imgaux(),*mask_i,
                              XMIPP_MAX(ROUND(localSize*XSIZE(imgaux()))/2,5));
                 maskImg.push_back(mask_i);
             }
 
             auto* img_i=new MultidimArray<unsigned char>;
             imgaux().rangeAdjust(0,255);
             typeCast(imgaux(),*img_i);
             img_i->setXmippOrigin();
             img.push_back(img_i);
 
             double tiltAngle;
             if (SF.containsLabel(MDL_ANGLE_TILT))
                 SF.getValue(MDL_ANGLE_TILT,tiltAngle, objId);
             else
                 tiltAngle=imgaux.tilt();
             tiltList.push_back(tiltAngle);
             if (tiltAngle!=0)
                 nonZeroTilt=true;
             if (ABS(tiltAngle)<minTilt)
             {
                 iMinTilt=n;
                 minTilt=ABS(tiltAngle);
             }
             name_list.push_back(imgaux.name());
 
             progress_bar(n++);
             iteration++;
         }
         progress_bar(Nimg);
         if (!nonZeroTilt)
             REPORT_ERROR(ERR_VALUE_NOTSET,"Tilt angles have not been assigned to the input selfile");
     }
 
     // Read images at original scale
     if (!fnSelOrig.empty())
     {
         SForig.read(fnSelOrig,nullptr);
         SForig.removeDisabled();
 
         if (SForig.size()!=SF.size())
             REPORT_ERROR(ERR_MD_OBJECTNUMBER,"The number of images in both selfiles (-i and -iorig) is different");
     }
 
     // Fill the affine transformations with empty matrices
     std::vector< Matrix2D<double> > emptyRow;
     for (int i=0; i<Nimg; i++)
     {
         Matrix2D<double> A;
         emptyRow.push_back(A);
     }
     for (int i=0; i<Nimg; i++)
     {
         affineTransformations.push_back(emptyRow);
         correlationList.push_back(-1);
     }
 
     // Check if there are transformations already calculated
     FileName fn_tmp = fnRoot+"_transformations.txt";
     if (fn_tmp.exists())
     {
         std::ifstream fhIn;
         fhIn.open(fn_tmp.c_str());
         if (!fhIn)
             REPORT_ERROR(ERR_IO_NOTEXIST,(String)"Cannot open "+ fn_tmp);
         size_t linesRead=0;
         while (!fhIn.eof())
         {
             std::string line;
             getline(fhIn, line);
             if (linesRead<affineTransformations.size()-1 && line!="")
             {
                 std::vector< double > data;
                 readFloatList(line.c_str(), 8, data);
                 double tilt=data[0];
                 int i=0;
                 double bestDistance=ABS(tilt-tiltList[0]);
                 for (int j=0; j<Nimg; j++)
                 {
                     double distance=ABS(tilt-tiltList[j]);
                     if (bestDistance>distance)
                     {
                         bestDistance=distance;
                         i=j;
                     }
                 }
                 affineTransformations[i][i+1].resize(3,3);
                 affineTransformations[i][i+1].initIdentity(3);
                 affineTransformations[i][i+1](0,0)=data[2];
                 affineTransformations[i][i+1](1,0)=data[3];
                 affineTransformations[i][i+1](0,1)=data[4];
                 affineTransformations[i][i+1](1,1)=data[5];
                 affineTransformations[i][i+1](0,2)=data[6];
                 affineTransformations[i][i+1](1,2)=data[7];
 
                 affineTransformations[i+1][i]=
                     affineTransformations[i][i+1].inv();
 
                 if (showAffine)
                 {
                     MultidimArray<unsigned char>& img_i=*img[i];
                     MultidimArray<unsigned char>& img_j=*img[i+1];
                     int maxShift=FLOOR(XSIZE(*img[0])*maxShiftPercentage);
                     Matrix2D<double> Aij, Aji;
                     Aij=affineTransformations[i][i+1];
                     Aji=affineTransformations[i+1][i];
                     bool isMirror=false;
                     std::cout << "Transformation between " << i << " ("
                     << tiltList[i] << ") and " << i+1 << " ("
                     << tiltList[i+1] << std::endl;
                     double cost = computeAffineTransformation(img_i, img_j,
                                   maxShift, maxIterDE, Aij, Aji,
                                   showAffine, thresholdAffine, globalAffine,
                                   isMirror,true, pyramidLevel);
                     if (cost<0)
                     {
                         affineTransformations[i][i+1].clear();
                         affineTransformations[i+1][i].clear();
                         writeTransformations(fn_tmp);
                     }
                 }
             }
             linesRead++;
         }
         fhIn.close();
     }
 
     // Compute the rest of transformations
     computeAffineTransformations(globalAffine);
 
     // Do not show refinement
     showRefinement = false;
 
     // Check which is the distribution of correlation
     if (!useCriticalPoints)
     {
         auto X0=(int)(STARTINGX(*(img[0]))+2*localSize*XSIZE(*(img[0])));
         auto XF=(int)(FINISHINGX(*(img[0]))-2*localSize*XSIZE(*(img[0])));
         auto Y0=(int)(STARTINGY(*(img[0]))+2*localSize*YSIZE(*(img[0])));
         auto YF=(int)(FINISHINGY(*(img[0]))-2*localSize*YSIZE(*(img[0])));
         avgForwardPatchCorr.initZeros(Nimg);
         avgBackwardPatchCorr.initZeros(Nimg);
         avgForwardPatchCorr.initConstant(1);
         avgBackwardPatchCorr.initConstant(1);
 
         for (int ii=0; ii<Nimg; ++ii)
         {
             // Compute average forward patch corr
             MultidimArray<double> corrList;
             if (ii<=Nimg-2)
             {
                 corrList.initZeros(200);
                 FOR_ALL_ELEMENTS_IN_ARRAY1D(corrList)
                 {
                     Matrix1D<double> rii(3), rjj(3);
                     do
                     {
                         XX(rii)=rnd_unif(X0,XF);
                         YY(rii)=rnd_unif(Y0,YF);
                         rjj=affineTransformations[ii][ii+1]*rii;
                         refineLandmark(ii,ii+1,rii,rjj,corrList(i),false);
                     }
                     while (corrList(i)<-0.99);
                 }
                 avgForwardPatchCorr(ii)=corrList.computeAvg();
             }
 
             // Compute average backward patch corr
             if (ii>0)
             {
                 corrList.initZeros(200);
                 FOR_ALL_ELEMENTS_IN_ARRAY1D(corrList)
                 {
                     Matrix1D<double> rii(3), rjj(3);
                     do
                     {
                         XX(rii)=rnd_unif(X0,XF);
                         YY(rii)=rnd_unif(Y0,YF);
                         rjj=affineTransformations[ii][ii-1]*rii;
                         refineLandmark(ii,ii-1,rii,rjj,corrList(i),false);
                     }
                     while (corrList(i)<-0.99);
                 }
                 avgBackwardPatchCorr(ii)=corrList.computeAvg();
             }
 
             std::cout << "Image " << ii << " Average correlation forward="
             << avgForwardPatchCorr(ii)
             << " backward=" << avgBackwardPatchCorr(ii) << std::endl;
             iteration++;
         }
     }
 
     // Identify outliers
     if (!doNotIdentifyOutliers)
     {
         identifyOutliers(false);
         computeAffineTransformations(false);
         identifyOutliers(true);
     }
     else
         isOutlier.initZeros(Nimg);
     if (lastStep==0)
         exit(0);
 }

◆ readLandmarkSet()

void ProgTomographAlignment::readLandmarkSet ( const FileName & fnLandmark )

Read landmark set.

Definition at line 2170 of file tomo_align_tilt_series.cpp.

 {
     std::ifstream fhIn;
     fhIn.open(fnLandmark.c_str());
     if (!fhIn)
         REPORT_ERROR(ERR_IO_NOTEXIST,(std::string)"Cannot open "+fnLandmark+" for input");
     std::string dummyStr;
     int Nlandmark;
     fhIn >> dummyStr >> dummyStr >> dummyStr >> dummyStr >> dummyStr
     >> Nlandmark >> Nimg;
     if (Nlandmark<=0)
         REPORT_ERROR(ERR_VALUE_INCORRECT,(std::string)"No landmarks are found in "+fnLandmark);
     allLandmarksX.resize(Nlandmark,Nimg);
     allLandmarksY.resize(Nlandmark,Nimg);
     allLandmarksX.initConstant(XSIZE(*img[0]));
     allLandmarksY.initConstant(XSIZE(*img[0]));
     fhIn.exceptions ( std::ifstream::eofbit |
                       std::ifstream::failbit |
                       std::ifstream::badbit );
     while (!fhIn.eof())
     {
         try
         {
             int dummyInt, x, y, i, j;
             fhIn >> dummyInt >> x >> y >> i >> j;
             i=i-1;
             allLandmarksX(j,i)=x+STARTINGX(*img[0]);
             allLandmarksY(j,i)=y+STARTINGY(*img[0]);
         }
         catch (std::ifstream::failure e)
         {
             // Do nothing with this line
         }
     }
     fhIn.close();
     std::cout << "The file " << fnLandmark << " has been read for the landmarks\n"
     << Nlandmark << " landmarks are read\n";
 }

◆ readParams()

void ProgTomographAlignment::readParams ( )

virtual

Read parameters from argument line.

Reimplemented from XmippProgram.

Definition at line 452 of file tomo_align_tilt_series.cpp.

 {
     fnSel=getParam("-i");
     fnSelOrig=getParam("--iorig");
     fnRoot=getParam("--oroot");
     if (fnRoot=="")
         fnRoot=fnSel.withoutExtension();
     globalAffine=checkParam("--globalAffine");
     useCriticalPoints=checkParam("--useCriticalPoints");
     if (useCriticalPoints)
         Ncritical=getIntParam("--useCriticalPoints");
     else
         Ncritical=0;
     seqLength=getIntParam("--seqLength");
     blindSeqLength=getIntParam("--blindSeqLength");
     maxStep=getIntParam("--maxStep");
     gridSamples=getIntParam("--gridSamples");
     psiMax=getDoubleParam("--psiMax");
     deltaRot=getDoubleParam("--deltaRot");
     localSize=getDoubleParam("--localSize");
     optimizeTiltAngle=checkParam("--optimizeTiltAngle");
     isCapillar=checkParam("--isCapillar");
     dontNormalize=checkParam("--dontNormalize");
     difficult=checkParam("--difficult");
     corrThreshold=getDoubleParam("--threshold");
     maxShiftPercentage=getDoubleParam("--maxShiftPercentage");
     maxIterDE=getIntParam("--maxIterDE");
     showAffine=checkParam("--showAffine");
     thresholdAffine=getDoubleParam("--thresholdAffine");
     identifyOutliersZ = getDoubleParam("--identifyOutliers");
     doNotIdentifyOutliers = checkParam("--noOutliers");
     pyramidLevel = getIntParam("--pyramid");
     numThreads = getIntParam("--thr");
     if (numThreads<1)
         numThreads = 1;
     lastStep=getIntParam("--lastStep");
 }

◆ refineChain()

bool ProgTomographAlignment::refineChain	(	LandmarkChain &	chain,
		double &	corrChain
	)

Refine chain.

Definition at line 1997 of file tomo_align_tilt_series.cpp.

 {
 #ifdef DEBUG
     std::cout << "Chain for refinement: ";
     for (int i=0; i<chain.size(); i++)
         std::cout << chain[i].imgIdx << " ";
     std::cout << std::endl;
 #endif
 
     for (int K=0; K<2 && (chain.size()>seqLength || useCriticalPoints); K++)
     {
         sort(chain.begin(), chain.end());
         Matrix1D<double> rii(2), rjj(2), newrjj(2);
 
         if (useCriticalPoints)
         {
             // compute average piece
             int halfSize=XMIPP_MAX(ROUND(localSize*XSIZE(*img[0]))/2,5);
             int chainLength=chain.size();
             MultidimArray<double> avgPiece(2*halfSize+1,2*halfSize+1), pieceAux;
             avgPiece.setXmippOrigin();
             pieceAux=avgPiece;
             for (int n=0; n<chainLength; n++)
             {
                 int ii=chain[n].imgIdx;
                 VECTOR_R2(rii,chain[n].x,chain[n].y);
                 if (XX(rii)-halfSize<STARTINGX(*img[ii])  ||
                     XX(rii)+halfSize>FINISHINGX(*img[ii]) ||
                     YY(rii)-halfSize<STARTINGY(*img[ii])  ||
                     YY(rii)+halfSize>FINISHINGY(*img[ii]))
                 {
                     corrChain=0;
                     return false;
                 }
                 FOR_ALL_ELEMENTS_IN_ARRAY2D(avgPiece)
                 pieceAux(i,j)=(*img[ii])((int)(YY(rii)+i),(int)(XX(rii)+j));
                 if (isCapillar && ABS(chain[n].imgIdx-chain[0].imgIdx)>Nimg/2)
                     pieceAux.selfReverseY();
                 avgPiece+=pieceAux;
             }
             avgPiece/=chainLength;
 #ifdef DEBUG
 
             Image<double> save;
             save()=avgPiece;
             save.write("PPPavg.xmp");
             std::cout << "Average " << K << " -------------------  \n";
             showRefinement=true;
 #endif
 
             // Align all images with respect to this average
             corrChain=2;
             for (int j=0; j<chainLength; j++)
             {
                 int jj=chain[j].imgIdx;
                 VECTOR_R2(rjj,chain[j].x,chain[j].y);
                 double corr;
                 bool accepted=refineLandmark(avgPiece,jj,rjj,0,false,corr);
                 if (accepted)
                 {
                     chain[j].x=XX(rjj);
                     chain[j].y=YY(rjj);
                     corrChain=XMIPP_MIN(corrChain,corr);
                 }
             }
 
 #ifdef DEBUG
             showRefinement=false;
 #endif
 
         }
         else
         {
             // Refine every step
             for (int step=2; step<=maxStep; step++)
             {
                 int chainLength=chain.size();
 
                 // Refine forwards every step
                 int ileft=-1;
                 for (int i=0; i<chainLength-step; i++)
                 {
                     int ii=chain[i].imgIdx;
                     VECTOR_R2(rii,chain[i].x,chain[i].y);
                     int jj=chain[i+step].imgIdx;
                     VECTOR_R2(rjj,chain[i+step].x,chain[i+step].y);
                     newrjj=rjj;
                     double corr;
                     bool accepted=refineLandmark(ii,jj,rii,newrjj,corr,false);
                     if (((newrjj-rjj).module()<4 && accepted) || useCriticalPoints)
                     {
                         chain[i+step].x=XX(newrjj);
                         chain[i+step].y=YY(newrjj);
                     }
                     else
                         ileft=i;
                 }
 
 #ifdef DEBUG
                 // COSS showRefinement=(step==maxStep);
 #endif
                 // Refine backwards all images
                 int iright=chainLength;
                 corrChain=2;
                 for (int i=chainLength-1; i>=1; i--)
                 {
                     int ii=chain[i].imgIdx;
                     VECTOR_R2(rii,chain[i].x,chain[i].y);
                     int jj=chain[i-1].imgIdx;
                     VECTOR_R2(rjj,chain[i-1].x,chain[i-1].y);
                     newrjj=rjj;
                     double corr;
                     bool accepted=refineLandmark(ii,jj,rii,newrjj,corr,false);
                     corrChain=XMIPP_MIN(corrChain,corr);
                     if (((newrjj-rjj).module()<4 && accepted) || useCriticalPoints)
                     {
                         chain[i-1].x=XX(newrjj);
                         chain[i-1].y=YY(newrjj);
                     }
                     else
                         iright=i;
                 }
 #ifdef DEBUG
                 // COSS: showRefinement=false;
 #endif
 
                 LandmarkChain refinedChain;
                 for (int ll=ileft+1; ll<=iright-1; ll++)
                     refinedChain.push_back(chain[ll]);
                 chain=refinedChain;
                 double tilt0=tiltList[chain[0].imgIdx];
                 double tiltF=tiltList[chain[chain.size()-1].imgIdx];
                 double lengthThreshold=XMIPP_MAX(3.,FLOOR(seqLength*cos(DEG2RAD(0.5*(tilt0+tiltF)))));
                 if (chain.size()<lengthThreshold && !useCriticalPoints)
                     return false;
             }
         }
     }
 
     // Check that the chain is of the desired length
     double tilt0=tiltList[chain[0].imgIdx];
     double tiltF=tiltList[chain[chain.size()-1].imgIdx];
     double lengthThreshold=XMIPP_MAX(3,FLOOR(seqLength*cos(DEG2RAD(0.5*(tilt0+tiltF)))));
     return chain.size()>lengthThreshold;
 }

◆ refineLandmark() [1/2]

bool ProgTomographAlignment::refineLandmark	(	int	ii,
		int	jj,
		const Matrix1D< double > &	rii,
		Matrix1D< double > &	rjj,
		double &	maxCorr,
		bool	tryFourier
	)		const

Refine landmark. ii is the index of the original image, jj is the index in the image at which the landmark is being refined. rii and rjj are the corresponding landmark positions in both images.

The function returns whether the landmark is accepted or not.

Definition at line 1735 of file tomo_align_tilt_series.cpp.

 {
     maxCorr=-1;
     int halfSize=XMIPP_MAX(ROUND(localSize*XSIZE(*img[ii]))/2,5);
     if (XX(rii)-halfSize<STARTINGX(*img[ii])  ||
         XX(rii)+halfSize>FINISHINGX(*img[ii]) ||
         YY(rii)-halfSize<STARTINGY(*img[ii])  ||
         YY(rii)+halfSize>FINISHINGY(*img[ii]) ||
         XX(rjj)-halfSize<STARTINGX(*img[jj])  ||
         XX(rjj)+halfSize>FINISHINGX(*img[jj]) ||
         YY(rjj)-halfSize<STARTINGY(*img[jj])  ||
         YY(rjj)+halfSize>FINISHINGY(*img[jj]))
         return 0;
 
     // Check if the two pieces are reversed
     bool reversed=isCapillar && ABS(ii-jj)>Nimg/2;
 
     // Select piece in image ii, compute its statistics and normalize
     MultidimArray<double> pieceii(2*halfSize+1,2*halfSize+1);
     pieceii.setXmippOrigin();
     const MultidimArray<unsigned char> &Iii=(*img[ii]);
     FOR_ALL_ELEMENTS_IN_ARRAY2D(pieceii)
     A2D_ELEM(pieceii,i,j)=A2D_ELEM(Iii,
                                    (int)(YY(rii)+i),(int)(XX(rii)+j));
     if (showRefinement)
     {
         Image<double> save;
         save()=pieceii;
         save.write("PPPpieceii.xmp");
         std::cout << "ii=" << ii << " jj=" << jj << std::endl;
         std::cout << "rii=" << rii.transpose() << std::endl;
     }
 
     // Choose threshold
     double actualCorrThreshold=corrThreshold;
     if (actualCorrThreshold<0 && VEC_XSIZE(avgForwardPatchCorr)>0)
     {
         if (ii<jj)
             actualCorrThreshold=XMIPP_MIN(avgForwardPatchCorr(ii),
                                           avgBackwardPatchCorr(jj));
         else
             actualCorrThreshold=XMIPP_MIN(avgBackwardPatchCorr(ii),
                                           avgForwardPatchCorr(jj));
     }
     if (showRefinement)
         std::cout << "actualCorrThreshold=" << actualCorrThreshold << std::endl;
 
     // Try Fourier
     if (tryFourier)
     {
         const MultidimArray<unsigned char> &Ijj=(*img[jj]);
         if (XX(rjj)-halfSize>=STARTINGX(Ijj) &&
             XX(rjj)+halfSize<=FINISHINGX(Ijj) &&
             YY(rjj)-halfSize>=STARTINGY(Ijj) &&
             YY(rjj)+halfSize<=FINISHINGY(Ijj))
         {
             // Take the piece at jj
             MultidimArray<double> piecejj(2*halfSize+1,2*halfSize+1);
             piecejj.setXmippOrigin();
             FOR_ALL_ELEMENTS_IN_ARRAY2D(piecejj)
             A2D_ELEM(piecejj,i,j)=A2D_ELEM(Ijj,
                                            (int)(YY(rjj)+i),(int)(XX(rjj)+j));
 
             // Compute its correlation with pieceii
             double corrOriginal=correlationIndex(pieceii,piecejj);
             if (showRefinement)
             {
                 Image<double> save;
                 save()=piecejj;
                 save.write("PPPpiecejjOriginal.xmp");
                 std::cout << "Corr original=" << corrOriginal << std::endl;
             }
 
             // Now try with the best shift
             double shiftX,shiftY;
             CorrelationAux aux;
             bestNonwrappingShift(pieceii,piecejj,shiftX,shiftY,aux);
             Matrix1D<double> fftShift(2);
             VECTOR_R2(fftShift,shiftX,shiftY);
             selfTranslate(xmipp_transformation::LINEAR,piecejj,fftShift,xmipp_transformation::WRAP);
             double corrFFT=correlationIndex(pieceii,piecejj);
             if (corrFFT>corrOriginal)
             {
                 XX(rjj)-=shiftX;
                 YY(rjj)-=shiftY;
             }
 
             if (showRefinement)
             {
                 Image<double> save;
                 save()=piecejj;
                 save.write("PPPpiecejjFFT.xmp");
                 std::cout << "FFT shift=" << fftShift.transpose() << std::endl;
                 std::cout << "Corr FFT=" << corrFFT << std::endl;
                 if (corrFFT>corrOriginal)
                 {
                     if (XX(rjj)-halfSize>=STARTINGX(Ijj) &&
                         XX(rjj)+halfSize<=FINISHINGX(Ijj) &&
                         YY(rjj)-halfSize>=STARTINGY(Ijj) &&
                         YY(rjj)+halfSize<=FINISHINGY(Ijj))
                     {
                         FOR_ALL_ELEMENTS_IN_ARRAY2D(piecejj)
                         A2D_ELEM(piecejj,i,j)=A2D_ELEM(Ijj,
                                                        (int)(YY(rjj)+i),(int)(XX(rjj)+j));
                         save()=piecejj;
                         save.write("PPPpiecejjNew.xmp");
                     }
                 }
             }
         }
     }
 
     bool retval=refineLandmark(pieceii,jj,rjj,actualCorrThreshold,
                                reversed,maxCorr);
     return retval;
 }

◆ refineLandmark() [2/2]

bool ProgTomographAlignment::refineLandmark	(	const MultidimArray< double > &	pieceii,
		int	jj,
		Matrix1D< double > &	rjj,
		double	actualCorrThreshold,
		bool	reversed,
		double &	maxCorr
	)		const

Refine landmark. The same as the previous function but an image is provided as pattern (ii) instead of an index and a position.

Definition at line 1854 of file tomo_align_tilt_series.cpp.

 {
     int halfSize=XSIZE(pieceii)/2;
 
     double mean_ii=0, stddev_ii=0;
     FOR_ALL_ELEMENTS_IN_ARRAY2D(pieceii)
     {
         mean_ii+=A2D_ELEM(pieceii,i,j);
         stddev_ii+=A2D_ELEM(pieceii,i,j)*A2D_ELEM(pieceii,i,j);
     }
     mean_ii/=MULTIDIM_SIZE(pieceii);
     stddev_ii = stddev_ii / MULTIDIM_SIZE(pieceii) - mean_ii * mean_ii;
     stddev_ii *= MULTIDIM_SIZE(pieceii) / (MULTIDIM_SIZE(pieceii) - 1);
     stddev_ii = sqrt(static_cast<double>((ABS(stddev_ii))));
     if (stddev_ii>XMIPP_EQUAL_ACCURACY)
         FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(pieceii)
         DIRECT_MULTIDIM_ELEM(pieceii,n)=
             (DIRECT_MULTIDIM_ELEM(pieceii,n)-mean_ii)/stddev_ii;
 
     // Try all possible shifts
     MultidimArray<double> corr((int)(1.5*(2*halfSize+1)),(int)(1.5*(2*halfSize+1)));
     corr.setXmippOrigin();
     corr.initConstant(-1.1);
     bool accept=false;
     double maxval=-1;
     MultidimArray<double> piecejj(2*halfSize+1,2*halfSize+1);
     piecejj.setXmippOrigin();
     if (stddev_ii>XMIPP_EQUAL_ACCURACY)
     {
         int imax=0, jmax=0;
         std::queue< Matrix1D<double> > Q;
         Q.push(vectorR2(0,0));
         while (!Q.empty())
         {
             // Get the first position to evaluate
             auto shifty=(int)YY(Q.front());
             auto shiftx=(int)XX(Q.front());
             Q.pop();
             if (!corr(shifty,shiftx)<-1)
                 continue;
 
             // Select piece in image jj, compute its statistics and normalize
             double mean_jj=0, stddev_jj=0;
             const MultidimArray<unsigned char> &Ijj=(*img[jj]);
             if (XX(rjj)+shiftx+STARTINGX(piecejj)<STARTINGX(Ijj) ||
                 YY(rjj)+shifty+STARTINGY(piecejj)<STARTINGY(Ijj) ||
                 XX(rjj)+shiftx+FINISHINGX(piecejj)>FINISHINGX(Ijj) ||
                 YY(rjj)+shifty+FINISHINGY(piecejj)>FINISHINGY(Ijj))
                 continue;
             FOR_ALL_ELEMENTS_IN_ARRAY2D(piecejj)
             {
                 double pixval=A2D_ELEM(Ijj,
                                        (int)(YY(rjj)+shifty+i),
                                        (int)(XX(rjj)+shiftx+j));
                 A2D_ELEM(piecejj,i,j)=pixval;
                 mean_jj+=pixval;
                 stddev_jj+=pixval*pixval;
             }
             mean_jj/=MULTIDIM_SIZE(piecejj);
             stddev_jj = stddev_jj / MULTIDIM_SIZE(piecejj) - mean_jj * mean_jj;
             stddev_jj *= MULTIDIM_SIZE(piecejj) / (MULTIDIM_SIZE(piecejj) - 1);
             stddev_jj = sqrt(static_cast<double>((ABS(stddev_jj))));
             if (reversed)
                 piecejj.selfReverseY();
 
             // Compute the correlation
             corr(shifty,shiftx)=0;
             double &corrRef=corr(shifty,shiftx);
             if (stddev_jj>XMIPP_EQUAL_ACCURACY)
             {
                 double istddev_jj=1.0/stddev_jj;
                 FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(piecejj)
                 corrRef+=
                     DIRECT_MULTIDIM_ELEM(pieceii,n)*
                     (DIRECT_MULTIDIM_ELEM(piecejj,n)-mean_jj)*istddev_jj;
                 corrRef/=MULTIDIM_SIZE(piecejj);
             }
 
             if (corrRef>maxval)
             {
                 maxval=corrRef;
                 imax=shifty;
                 jmax=shiftx;
                 for (int step=1; step<=5; step+=2)
                     for (int stepy=-1; stepy<=1; stepy++)
                         for (int stepx=-1; stepx<=1; stepx++)
                         {
                             int newshifty=shifty+stepy*step;
                             int newshiftx=shiftx+stepx*step;
                             if (newshifty>=STARTINGY(corr) &&
                                 newshifty<=FINISHINGY(corr) &&
                                 newshiftx>=STARTINGX(corr) &&
                                 newshiftx<=FINISHINGX(corr) &&
                                 (XX(rjj)+newshiftx-halfSize)>=STARTINGX(Ijj) &&
                                 (XX(rjj)+newshiftx+halfSize)<=FINISHINGX(Ijj) &&
                                 (YY(rjj)+newshifty-halfSize)>=STARTINGY(Ijj) &&
                                 (YY(rjj)+newshifty+halfSize)<=FINISHINGY(Ijj))
                                 if (corr(newshifty,newshiftx)<-1)
                                     Q.push(vectorR2(newshiftx,newshifty));
                         }
             }
         }
 
         if (maxval>actualCorrThreshold)
         {
             XX(rjj)+=jmax;
             if (reversed)
                 YY(rjj)-=imax;
             else
                 YY(rjj)+=imax;
             accept=true;
         }
         if (showRefinement)
         {
             Image<double> save;
             FOR_ALL_ELEMENTS_IN_ARRAY2D(piecejj)
             piecejj(i,j)=(*img[jj])((int)(YY(rjj)+i),(int)(XX(rjj)+j));
             if (reversed)
                 piecejj.selfReverseY();
             save()=piecejj;
             save.write("PPPpiecejj.xmp");
             save()=corr;
             save.write("PPPcorr.xmp");
             std::cout << "jj=" << jj << " rjj=" << rjj.transpose() << std::endl;
             std::cout << "imax=" << imax << " jmax=" << jmax << std::endl;
             std::cout << "maxval=" << maxval << std::endl;
         }
     }
     if (showRefinement)
     {
         std::cout << "Accepted=" << accept << std::endl;
         std::cout << "Press any key\n";
         char c;
         std::cin >> c;
     }
     maxCorr=maxval;
     return accept;
 }

◆ removeOutlierLandmarks()

void ProgTomographAlignment::removeOutlierLandmarks ( const Alignment & alignment )

Remove Outliers.

Definition at line 2455 of file tomo_align_tilt_series.cpp.

 {
     std::cout << "Removing outliers ...\n";
 
     // Compute threshold for outliers
     Histogram1D hist;
     compute_hist(alignment.errorLandmark, hist, 100);
     // double threshold0=hist.percentil(10);
     double thresholdF=hist.percentil(90);
 
     // Identify outliers
     std::vector<bool> validLandmark;
     int invalidCounter=0;
     int Nlandmark=MAT_YSIZE(allLandmarksX);
     for (int j=0; j<Nlandmark; j++)
         if (//COSS alignment.errorLandmark(j)<threshold0 ||
             alignment.errorLandmark(j)>thresholdF)
         {
             validLandmark.push_back(false);
             invalidCounter++;
         }
         else
             validLandmark.push_back(true);
 
     // Remove outliers
     Matrix2D<double> newAllLandmarksX(Nlandmark-invalidCounter,Nimg);
     Matrix2D<double> newAllLandmarksY(Nlandmark-invalidCounter,Nimg);
     int jj=0;
     for (int j=0; j<Nlandmark; j++)
     {
         if (!validLandmark[j])
             continue;
         for (int i=0; i<Nimg; i++)
         {
             newAllLandmarksX(jj,i)=allLandmarksX(j,i);
             newAllLandmarksY(jj,i)=allLandmarksY(j,i);
         }
         jj++;
     }
     allLandmarksX=newAllLandmarksX;
     allLandmarksY=newAllLandmarksY;
 
     std::cout << invalidCounter << " out of " << Nlandmark
     << " landmarks have been removed\n";
 }

◆ run()

void ProgTomographAlignment::run ( )

virtual

Run: do the real work.

Reimplemented from XmippProgram.

Definition at line 2518 of file tomo_align_tilt_series.cpp.

 {
     produceSideInfo();
     std::cerr << "generateLandmarkSet"<< std::endl;
     generateLandmarkSet();
     std::cerr << "produceInformationFromLandmarks" << std::endl;
     produceInformationFromLandmarks();
     std::cerr << "alignment" << std::endl;
     auto *alignment=new Alignment(this);
 
     // Exhaustive search for rot
     double bestError=0, bestRot=-1;
     for (double rot=0; rot<=180-deltaRot; rot+=deltaRot)
     {
         alignment->clear();
         alignment->rot=rot;
         double error=alignment->optimizeGivenAxisDirection();
 #ifdef DEBUG
 
         std::cout << "rot= " << rot
         << " error= " << error << std::endl;
 #endif
 
         if (bestRot<0 || bestError>error)
         {
             bestRot=rot;
             bestError=error;
             *bestPreviousAlignment=*alignment;
         }
     }
     delete alignment;
     std::cout << "Best rot=" << bestRot
     << " Best error=" << bestError << std::endl;
 
     // Continuous optimization for the axis direction
     Matrix1D<double> axisAngles(2), steps(2);
     axisAngles(0)=bestRot;
     axisAngles(1)=90;
     steps.initConstant(1);
     if (!optimizeTiltAngle)
         steps(1)=0;
     double fitness;
     int iter;
     TomographAlignment::global_prm=this;
     powellOptimizer(axisAngles,1,2,&wrapperError, nullptr,
                     0.01,fitness,iter,steps,true);
 
     // Outlier removal
     for (int i=0; i<3; i++)
     {
         // Compute the best alignment
         bestPreviousAlignment->rot=axisAngles(0);
         bestPreviousAlignment->tilt=axisAngles(1);
         fitness=bestPreviousAlignment->optimizeGivenAxisDirection();
         bestPreviousAlignment->computeErrorForLandmarks();
 
         // Remove landmarks that are outliers in the current model
         removeOutlierLandmarks(*bestPreviousAlignment);
         produceInformationFromLandmarks();
         delete bestPreviousAlignment;
         bestPreviousAlignment=new Alignment(this);
         bestPreviousAlignment->rot=axisAngles(0);
         bestPreviousAlignment->tilt=axisAngles(1);
         fitness=bestPreviousAlignment->optimizeGivenAxisDirection();
 
         // Optimize again
         powellOptimizer(axisAngles,1,2,&wrapperError,nullptr,
                         0.01,fitness,iter,steps,true);
     }
     bestPreviousAlignment->rot=axisAngles(0);
     bestPreviousAlignment->tilt=axisAngles(1);
     fitness=bestPreviousAlignment->optimizeGivenAxisDirection();
 
     // Save the alignment
     writeLandmarkSet(fnRoot+"_good_landmarks.txt");
     std::ofstream fh_out;
     fh_out.open((fnRoot+"_alignment.txt").c_str());
     if (!fh_out)
         REPORT_ERROR(ERR_IO_NOWRITE,
                      (std::string)"Cannot open "+fnRoot+"_alignment.txt for output");
     fh_out << *bestPreviousAlignment;
     fh_out.close();
 
     // Correct the input images
     alignImages(*bestPreviousAlignment);
     writeLandmarkSet(fnRoot+"_good_landmarks_corrected.txt");
     delete bestPreviousAlignment;
 }

◆ show()

void ProgTomographAlignment::show ( )

Show parameters.

Definition at line 490 of file tomo_align_tilt_series.cpp.

 {
     std::cout << "Input images:       " << fnSel              << std::endl
     << "Original images:    " << fnSelOrig          << std::endl
     << "Output rootname:    " << fnRoot             << std::endl
     << "Global affine:      " << globalAffine       << std::endl
     << "Use critical points:" << useCriticalPoints  << std::endl
     << "Num critical points:" << Ncritical          << std::endl
     << "SeqLength:          " << seqLength          << std::endl
     << "BlindSeqLength:     " << blindSeqLength     << std::endl
     << "MaxStep:            " << maxStep            << std::endl
     << "Grid samples:       " << gridSamples        << std::endl
     << "Maximum psi:        " << psiMax             << std::endl
     << "Delta rot:          " << deltaRot           << std::endl
     << "Local size:         " << localSize          << std::endl
     << "Optimize tilt angle:" << optimizeTiltAngle  << std::endl
     << "isCapillar:         " << isCapillar         << std::endl
     << "DontNormalize:      " << dontNormalize      << std::endl
     << "Difficult:          " << difficult          << std::endl
     << "Threshold:          " << corrThreshold      << std::endl
     << "MaxShift Percentage:" << maxShiftPercentage << std::endl
     << "MaxIterDE:          " << maxIterDE          << std::endl
     << "Show Affine:        " << showAffine         << std::endl
     << "Threshold Affine:   " << thresholdAffine    << std::endl
     << "Identify outliers Z:" << identifyOutliersZ     << std::endl
     << "No outliers:        " << doNotIdentifyOutliers << std::endl
     << "Pyramid level:      " << pyramidLevel       << std::endl
     << "Threads to use:     " << numThreads         << std::endl
     << "Last step:          " << lastStep           << std::endl
     ;
 }

◆ writeLandmarkSet()

void ProgTomographAlignment::writeLandmarkSet ( const FileName & fnLandmark ) const

Write landmark set.

Definition at line 2146 of file tomo_align_tilt_series.cpp.

 {
     std::ofstream fhOut;
     fhOut.open(fnLandmark.c_str());
     if (!fhOut)
         REPORT_ERROR(ERR_IO_NOWRITE,(std::string)"Cannot open "+fnLandmark+" for output");
     fhOut << "Point     x       y       slice   color "
     << MAT_YSIZE(allLandmarksX) << " " << MAT_XSIZE(allLandmarksX) << std::endl;
     for (size_t i=0; i<MAT_XSIZE(allLandmarksX); i++)
     {
         int counter=0;
         for (size_t j=0; j<MAT_YSIZE(allLandmarksX); j++)
             if (allLandmarksX(j,i)!=XSIZE(*img[0]))
             {
                 fhOut << counter << " \t"
                 << ROUND(allLandmarksX(j,i)-STARTINGX(*img[0])) << " \t"
                 << ROUND(allLandmarksY(j,i)-STARTINGY(*img[0])) << " \t"
                 << i+1 << " \t" << j << std::endl;
                 counter++;
             }
     }
     fhOut.close();
 }

◆ writeTransformations()

void ProgTomographAlignment::writeTransformations ( const FileName & fnTransformations ) const

Write affine transformations.

Definition at line 2210 of file tomo_align_tilt_series.cpp.

 {
     std::ofstream fhOut;
     fhOut.open(fnTransformations.c_str());
     if (!fhOut)
         REPORT_ERROR(ERR_IO_NOWRITE,(std::string)"Cannot open "+fnTransformations+" for output");
     int imax=affineTransformations.size();
     int counter=0;
     for (int i=0; i<imax; i++)
     {
         int jmax=affineTransformations[i].size();
         for (int j=i+1; j<jmax; j++)
             if (MAT_XSIZE(affineTransformations[i][j])!=0)
             {
                 fhOut << tiltList[i] << "\t0.0\t"
                 << affineTransformations[i][j](0,0) << "\t"
                 << affineTransformations[i][j](1,0) << "\t"
                 << affineTransformations[i][j](0,1) << "\t"
                 << affineTransformations[i][j](1,1) << "\t"
                 << affineTransformations[i][j](0,2) << "\t"
                 << affineTransformations[i][j](1,2) << std::endl;
                 counter++;
             }
     }
     if (counter==imax-1)
     {
         fhOut << tiltList[tiltList.size()-1] << "\t0.0\t1.0\t0.0\t0.0\t1.0\t0.0\t0.0\n";
         fhOut << "1 0 0 1 0 0\n";
     }
     fhOut.close();
 }

Member Data Documentation

◆ affineTransformations

std::vector< std::vector< Matrix2D<double> > > ProgTomographAlignment::affineTransformations

Definition at line 265 of file tomo_align_tilt_series.h.

◆ allLandmarksX

Matrix2D<double> ProgTomographAlignment::allLandmarksX

Definition at line 271 of file tomo_align_tilt_series.h.

◆ allLandmarksY

Matrix2D<double> ProgTomographAlignment::allLandmarksY

Definition at line 274 of file tomo_align_tilt_series.h.

◆ avgBackwardPatchCorr

Matrix1D<double> ProgTomographAlignment::avgBackwardPatchCorr

Definition at line 255 of file tomo_align_tilt_series.h.

◆ avgForwardPatchCorr

Matrix1D<double> ProgTomographAlignment::avgForwardPatchCorr

Definition at line 252 of file tomo_align_tilt_series.h.

◆ barpi

std::vector< Matrix1D<double> > ProgTomographAlignment::barpi

Definition at line 283 of file tomo_align_tilt_series.h.

◆ bestPreviousAlignment

Alignment* ProgTomographAlignment::bestPreviousAlignment

Definition at line 289 of file tomo_align_tilt_series.h.

◆ blindSeqLength

size_t ProgTomographAlignment::blindSeqLength

Add blind landmarks. Set to -1 for no blind landmarks

Definition at line 105 of file tomo_align_tilt_series.h.

◆ correlationList

std::vector< double > ProgTomographAlignment::correlationList

Definition at line 268 of file tomo_align_tilt_series.h.

◆ corrThreshold

double ProgTomographAlignment::corrThreshold

Correlation threshold for a good landmark.

Definition at line 141 of file tomo_align_tilt_series.h.

◆ deltaRot

double ProgTomographAlignment::deltaRot

Delta rot.

Definition at line 123 of file tomo_align_tilt_series.h.

◆ difficult

bool ProgTomographAlignment::difficult

Difficult.

Definition at line 138 of file tomo_align_tilt_series.h.

◆ doNotIdentifyOutliers

bool ProgTomographAlignment::doNotIdentifyOutliers

Do not identify outlier micrographs.

Definition at line 153 of file tomo_align_tilt_series.h.

◆ dontNormalize

bool ProgTomographAlignment::dontNormalize

Don't normalize the corrected images.

Definition at line 135 of file tomo_align_tilt_series.h.

◆ fnRoot

FileName ProgTomographAlignment::fnRoot

Output root.

Definition at line 86 of file tomo_align_tilt_series.h.

◆ fnSel

FileName ProgTomographAlignment::fnSel

MetaData File with all images.

Definition at line 80 of file tomo_align_tilt_series.h.

◆ fnSelOrig

FileName ProgTomographAlignment::fnSelOrig

MetaData File with all images at the original scale.

Definition at line 83 of file tomo_align_tilt_series.h.

◆ globalAffine

bool ProgTomographAlignment::globalAffine

Look for local affine transformation.

Definition at line 92 of file tomo_align_tilt_series.h.

◆ gridSamples

int ProgTomographAlignment::gridSamples

Grid samples.

Definition at line 108 of file tomo_align_tilt_series.h.

◆ identifyOutliersZ

double ProgTomographAlignment::identifyOutliersZ

Identify outlier micrographs.

Definition at line 150 of file tomo_align_tilt_series.h.

◆ img

std::vector< MultidimArray<unsigned char> *> ProgTomographAlignment::img

Definition at line 237 of file tomo_align_tilt_series.h.

◆ iMinTilt

int ProgTomographAlignment::iMinTilt

Definition at line 243 of file tomo_align_tilt_series.h.

◆ isCapillar

bool ProgTomographAlignment::isCapillar

This tilt series comes from a capillar.

Definition at line 132 of file tomo_align_tilt_series.h.

◆ isOutlier

Matrix1D<int> ProgTomographAlignment::isOutlier

Definition at line 258 of file tomo_align_tilt_series.h.

◆ iteration

int ProgTomographAlignment::iteration

Definition at line 162 of file tomo_align_tilt_series.h.

◆ lastStep

int ProgTomographAlignment::lastStep

Last step to run (-1=run them all)

Definition at line 159 of file tomo_align_tilt_series.h.

◆ localSize

double ProgTomographAlignment::localSize

Local refinement size.

Definition at line 126 of file tomo_align_tilt_series.h.

◆ maskImg

std::vector< MultidimArray<unsigned char> *> ProgTomographAlignment::maskImg

Definition at line 240 of file tomo_align_tilt_series.h.

◆ maxIterDE

int ProgTomographAlignment::maxIterDE

Maximum number of iterations in Differential Evolution.

Definition at line 114 of file tomo_align_tilt_series.h.

◆ maxShiftPercentage

double ProgTomographAlignment::maxShiftPercentage

Maxshift percentage.

Definition at line 111 of file tomo_align_tilt_series.h.

◆ maxStep

int ProgTomographAlignment::maxStep

Maximum Step for chain refinement.

Definition at line 120 of file tomo_align_tilt_series.h.

◆ name_list

std::vector< std::string > ProgTomographAlignment::name_list

Definition at line 249 of file tomo_align_tilt_series.h.

◆ Ncritical

size_t ProgTomographAlignment::Ncritical

Number of critical points.

Definition at line 98 of file tomo_align_tilt_series.h.

◆ ni

Matrix1D<int> ProgTomographAlignment::ni

Definition at line 286 of file tomo_align_tilt_series.h.

◆ Nimg

int ProgTomographAlignment::Nimg

Definition at line 261 of file tomo_align_tilt_series.h.

◆ numThreads

int ProgTomographAlignment::numThreads

Number of threads to use for parallel computing.

Definition at line 89 of file tomo_align_tilt_series.h.

◆ optimizeTiltAngle

bool ProgTomographAlignment::optimizeTiltAngle

Optimize tilt angle of tilt axis.

Definition at line 129 of file tomo_align_tilt_series.h.

◆ psiMax

double ProgTomographAlignment::psiMax

Maximum rotation.

Definition at line 117 of file tomo_align_tilt_series.h.

◆ pyramidLevel

int ProgTomographAlignment::pyramidLevel

Pyramid.

Definition at line 156 of file tomo_align_tilt_series.h.

◆ seqLength

size_t ProgTomographAlignment::seqLength

Sequence Length.

Definition at line 101 of file tomo_align_tilt_series.h.

◆ SF

MetaDataVec ProgTomographAlignment::SF

Definition at line 231 of file tomo_align_tilt_series.h.

◆ SForig

MetaDataVec ProgTomographAlignment::SForig

Definition at line 234 of file tomo_align_tilt_series.h.

◆ showAffine

bool ProgTomographAlignment::showAffine

Show affine transformations.

Definition at line 144 of file tomo_align_tilt_series.h.

◆ showRefinement

bool ProgTomographAlignment::showRefinement

Definition at line 292 of file tomo_align_tilt_series.h.

◆ thresholdAffine

double ProgTomographAlignment::thresholdAffine

Threshold affine.

Definition at line 147 of file tomo_align_tilt_series.h.

◆ tiltList

std::vector< double > ProgTomographAlignment::tiltList

Definition at line 246 of file tomo_align_tilt_series.h.

◆ useCriticalPoints

bool ProgTomographAlignment::useCriticalPoints

Use critical points.

Definition at line 95 of file tomo_align_tilt_series.h.

◆ Vseti

std::vector< std::vector<int> > ProgTomographAlignment::Vseti

Definition at line 277 of file tomo_align_tilt_series.h.

◆ Vsetj

std::vector< std::vector<int> > ProgTomographAlignment::Vsetj

Definition at line 280 of file tomo_align_tilt_series.h.

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

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

Public Member Functions

Public Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ ~ProgTomographAlignment()

Member Function Documentation

◆ alignImages()

◆ computeAffineTransformations()

◆ defineParams()

◆ generateLandmarkSet()

◆ identifyOutliers()

◆ produceInformationFromLandmarks()

◆ produceSideInfo()

◆ readLandmarkSet()

◆ readParams()

◆ refineChain()

◆ refineLandmark() [1/2]

◆ refineLandmark() [2/2]

◆ removeOutlierLandmarks()

◆ run()

◆ show()

◆ writeLandmarkSet()

◆ writeTransformations()

Member Data Documentation

◆ affineTransformations

◆ allLandmarksX

◆ allLandmarksY

◆ avgBackwardPatchCorr

◆ avgForwardPatchCorr

◆ barpi

◆ bestPreviousAlignment

◆ blindSeqLength

◆ correlationList

◆ corrThreshold

◆ deltaRot

◆ difficult

◆ doNotIdentifyOutliers

◆ dontNormalize

◆ fnRoot

◆ fnSel

◆ fnSelOrig

◆ globalAffine

◆ gridSamples

◆ identifyOutliersZ

◆ img

◆ iMinTilt

◆ isCapillar

◆ isOutlier

◆ iteration

◆ lastStep

◆ localSize

◆ maskImg

◆ maxIterDE

◆ maxShiftPercentage

◆ maxStep

◆ name_list

◆ Ncritical

◆ ni

◆ Nimg

◆ numThreads

◆ optimizeTiltAngle

◆ psiMax

◆ pyramidLevel

◆ seqLength

◆ SF

◆ SForig

◆ showAffine

◆ showRefinement

◆ thresholdAffine

◆ tiltList

◆ useCriticalPoints

◆ Vseti

◆ Vsetj