#include <angular_projection_matching.h>

Inheritance diagram for ProgAngularProjectionMatching:

Collaboration diagram for ProgAngularProjectionMatching:

Public Member Functions
virtual void	readParams ()
	Read arguments from command line. More...

virtual void	defineParams ()
	Define arguments accepted. More...

void	show ()

void	run ()

virtual void	produceSideInfo ()

void	rotationallyAlignOneImage (Matrix2D< double > &img, int imgno, int &opt_samplenr, double &opt_psi, bool &opt_flip, double &maxcorr)

void	translationallyAlignOneImage (MultidimArray< double > &img, const int &samplenr, const double &psi, const bool &opt_flip, double &opt_xoff, double &opt_yoff, double &maxcorr)

void	scaleAlignOneImage (MultidimArray< double > &img, const int &samplenr, const double &psi, const bool &opt_flip, const double &opt_xoff, const double &opt_yoff, const double &old_scale, double &opt_scale, double &maxcorr)

void	getCurrentReference (int refno, Polar_fftw_plans &local_plans)

virtual void	processAllImages ()

void	processSomeImages (const std::vector< size_t > &imagesToProcess)

void	getCurrentImage (size_t imgid, Image< double > &img)

virtual void	writeOutputFiles ()

void	destroyAndClean ()

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_exp

FileName	fn_ref

FileName	fn_out

FileName	fn_ctf

MetaDataDb	DFexp

MetaDataDb	DFo

size_t	dim

size_t	paddim

double	pad

double	max_shift

int	Ri

int	Ro

double	avail_memory

int	max_nr_refs_in_memory

int	max_nr_imgs_in_memory

int	total_nr_refs

int	counter_refs_in_memory

std::vector< int >	pointer_allrefs2refsinmem

std::vector< int >	pointer_refsinmem2allrefs

std::vector< size_t >	convert_refno_to_stack_position

std::vector< size_t >	ids

Polar< std::complex< double > > *	fP_ref

Polar< std::complex< double > > *	fP_img

Polar< std::complex< double > > *	fPm_img

MultidimArray< double > *	proj_ref

Polar_fftw_plans	global_plans

double *	stddev_ref

double *	stddev_img

Sampling	mysampling

bool	loop_forward_refs

int	search5d_shift

int	search5d_step

std::vector< int >	search5d_xoff

std::vector< int >	search5d_yoff

MultidimArray< double >	Mctf

bool	phase_flipped

int	threads

int	numOrientations

size_t	nr_trans

barrier_t	thread_barrier

bool	do_scale

WriteModeMetaData	do_overwrite

double	scale_step

double	scale_nsteps

int	progress_bar_step

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

projection_matching parameters.

Definition at line 64 of file angular_projection_matching.h.

Member Function Documentation

◆ defineParams()

void ProgAngularProjectionMatching::defineParams ( )

virtual

Define arguments accepted.

Reimplemented from XmippProgram.

Reimplemented in MpiProgAngularProjectionMatching.

Definition at line 83 of file angular_projection_matching.cpp.

 {
     addUsageLine("Perform a discrete angular assignment using projection matching in real space.");
     addUsageLine("This program is relatively fast, using polar coordinates for the in-plane ");
     addUsageLine("angular searches and the 5-dimensional search of rotation angles and origin ");
     addUsageLine("offsets is broken in two: first the angles are search in a 3D-search; then, ");
     addUsageLine("for the optimal orientation the origin offsets are searched (2D).");
     addUsageLine(" ");
     addUsageLine("The output of the program consists of a document file with all assigned angles");
     addUsageLine("and rotations. This file also contains a column for the maximum cross-correlation ");
     addUsageLine("coefficient. Note that the program does not alter the image headers. ");
     addUsageLine("The recommended use of this program is within the python script of the ");
     addUsageLine("xmipp_protocol_projmatch.py");
     addSeeAlsoLine("angular_discrete_assign, angular_continuous_assign, angular_project_library");
     addExampleLine("Example of use: Sample at 2 pixel step size for 5D shift search",false);
     addExampleLine("xmipp_angular_projection_matching -i experimental.doc -o assigned_angles.doc --ref reference.stk --search5d_step 2");
     addParamsLine("   -i <doc_file>                : Docfile with input images");
     addParamsLine("   -o <output_filename>         : Output filename");
     addParamsLine("   -r <stackFile>               : Reference projections");
     addParamsLine("     alias --ref;");
     addParamsLine("  [--search5d_shift <s5dshift=0>]: Search range (in +/- pix) for 5D shift search");
     addParamsLine("  [--search5d_step <s5dstep=2>]  : Step size for 5D shift search (in pix)");
     addParamsLine("  [--Ri <ri=1>]               : Inner radius to limit rotational search");
     addParamsLine("  [--Ro <ro=-1>]              : Outer radius to limit rotational search");
     addParamsLine("                        : ro = -1 -> dim/2-1");
     addParamsLine("  [-s <step=1> <n_steps=3>]    : scale step factor (1 means 0.01 in/de-crements) and number of steps around 1.");
     addParamsLine("                               : with default values: 1 0.01 | 0.02 | 0.03");
     addParamsLine("    alias --scale;");
     addParamsLine("==+Extra parameters==");
     addParamsLine("  [--mem <mem=1>]             : Available memory for reference library (Gb)");
     addParamsLine("  [--max_shift <max_shift=-1>]   : Max. change in origin offset (+/- pixels; neg= no limit)");
     addParamsLine("  [--ctf <filename>]            : CTF to apply to the reference projections, either a");
     addParamsLine("                     : CTF parameter file or a 2D image with the CTF amplitudes");
     addParamsLine("  [--pad <pad=1>]             : Padding factor (for CTF correction only)");
     addParamsLine("  [--phase_flipped]            : Use this if the experimental images have been phase flipped");
     addParamsLine("  [--thr <threads=1>]           : Number of concurrent threads");
     addParamsLine("  [--number_orientations <numOrientations=1>]  : Number of possible orientations for each experimental image");
     addParamsLine("  [--append]                : Append (versus overwrite) data to the output file");
 }

◆ destroyAndClean()

void ProgAngularProjectionMatching::destroyAndClean ( )

destroy and clean

Definition at line 197 of file angular_projection_matching.cpp.

 {
     delete [] fP_ref;
     delete [] proj_ref;
     delete [] fP_img;
     delete [] fPm_img;
     delete [] stddev_ref;
     delete [] stddev_img;
 
 }

◆ getCurrentImage()

void ProgAngularProjectionMatching::getCurrentImage	(	size_t	imgid,
		Image< double > &	img
	)

Read current image into memory and translate according to previous optimal Xoff and Yoff

Definition at line 1194 of file angular_projection_matching.cpp.

 {
     FileName fn_img;
     Matrix2D<double> A;
     //init A just in case
     A.initIdentity(3);
 
     // jump to line imgno+1 in DFexp, get data and filename
     DFexp.getValue(MDL_IMAGE,fn_img, imgid);
 
     // Read actual image
     img.read(fn_img);
     img().setXmippOrigin();
 
     // Store translation in header and apply it to the actual image
     //No need to get initial angles since those came with the reference projection
     double shiftX=0., shiftY=0.;
     DFexp.getValue(MDL_SHIFT_X,shiftX, imgid);
     DFexp.getValue(MDL_SHIFT_Y,shiftY, imgid);
 
     img.setShifts(shiftX,shiftY);
 
     //FIXME ROB
     //img.setEulerAngles(0.,0.,psi);
     img.setEulerAngles(0.,0.,0);
     //
     img.setFlip(0.);
 
     double scale;
     scale = 1.;
     if(DFexp.containsLabel(MDL_SCALE))
         DFexp.getValue(MDL_SCALE, scale, imgid);
     img.setScale(scale);
 
     img.getTransformationMatrix(A,true);
     A=A.inv();
     //std::cerr << "DEBUG_ROB, A:" << A << std::endl;
     //img.write("before.spi");
     if (!A.isIdentity())
         selfApplyGeometry(xmipp_transformation::BSPLINE3, img(), A, xmipp_transformation::IS_INV, xmipp_transformation::WRAP);
     //img.write("after.spi");
     //exit(0);
 }

◆ getCurrentReference()

void ProgAngularProjectionMatching::getCurrentReference	(	int	refno,
		Polar_fftw_plans &	local_plans
	)

Get pointer to the current reference image If this image wasn't stored in memory yet, read it from disc and store FT of the polar transform as well as the original image

a delete _DATA_ALL

Definition at line 408 of file angular_projection_matching.cpp.

 {
     FileName                      fnt;
     Image<double>                 img;
     double                        mean,stddev;
     MultidimArray<double>         Maux;
     Polar<double>                 P;
     Polar<std::complex <double> > fP;
     FourierTransformer                     local_transformer;
 
     // Image was not stored yet: read it from disc and store
     //    std::vector<size_t>::const_iterator found =
     //        std::find((mysampling.no_redundant_sampling_points_index).begin(),
     //                  (mysampling.no_redundant_sampling_points_index).end(),
     //                  (size_t)refno
     //                 );
     //    //found = found - (mysampling.no_redundant_sampling_points_index).begin();
     //    _pointer = found - (mysampling.no_redundant_sampling_points_index).begin();
     //    if (found == (mysampling.no_redundant_sampling_points_index).end())
     //        REPORT_ERROR(ERR_VALUE_INCORRECT, "Wrong reference number");
     //    fnt.compose(_pointer + FIRST_IMAGE, fn_ref);
     fnt.compose(convert_refno_to_stack_position[refno] + FIRST_IMAGE, fn_ref);
     img.read(fnt, _DATA_ALL);
     img().setXmippOrigin();
 
     //#define DEBUG
 #ifdef DEBUG
 
     double rot_tmp,tilt_tmp,psi_tmp;
     img.getEulerAngles(rot_tmp,tilt_tmp,psi_tmp);
 
     {
         std::cerr << "index_found: " << convert_refno_to_stack_position[refno] << std::endl;
         std::cerr << "refno: " << refno<<std::endl;
         std::cerr << "reading image " << fnt << std::endl;
         std::cerr << "rot_tmp,tilt_tmp,psi_tmp: " << rot_tmp<< " "<< tilt_tmp<< " "<<psi_tmp<< std::endl;
         //        std::cerr << "XXXXno_redundant_sampling_points_indexXXXXXX" <<std::endl;
         //        for (std::vector<size_t>::iterator i =
         //                    mysampling.my_neighbors.begin();
         //                i != mysampling.my_neighbors.end();
         //                ++i)
         //            std::cerr << *i << " ";
         std::cerr << std::endl;
     }
 #endif
 #undef DEBUG
     // Apply CTF
     if (fn_ctf!="")
     {
         MultidimArray<std::complex<double> > Faux;
         if (paddim > dim)
         {
             // pad real-space image
             int x0 = FIRST_XMIPP_INDEX(paddim);
             int xF = LAST_XMIPP_INDEX(paddim);
             img().selfWindow(x0, x0, xF,xF);
         }
         local_transformer.FourierTransform(img(),Faux);
         FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(Faux)
         {
             dAij(Faux,i,j) *= dAij(Mctf,i,j);
         }
         local_transformer.inverseFourierTransform(Faux,img());
 
         if (paddim > dim)
         {
             // de-pad real-space image
             int x0 = FIRST_XMIPP_INDEX(dim);
             int xF = LAST_XMIPP_INDEX(dim);
             img().selfWindow(x0, x0, xF,xF);
         }
     }
 
     // Calculate FTs of polar rings and its stddev
     produceSplineCoefficients(xmipp_transformation::BSPLINE3,Maux,img());
     P.getPolarFromCartesianBSpline(Maux,Ri,Ro);
     P.computeAverageAndStddev(mean,stddev);
     P -= mean;
     fourierTransformRings(P,fP,local_plans,true);
 
     pthread_mutex_lock(  &update_refs_in_memory_mutex );
 
     int counter = counter_refs_in_memory % max_nr_refs_in_memory;
     pointer_allrefs2refsinmem[refno] = counter;
     if (pointer_refsinmem2allrefs[counter] != -1)
     {
         // This position was in use already
         // Images will be overwritten, so reset the
         // pointer_allrefs2refsinmem of the old images to -1
         pointer_allrefs2refsinmem[pointer_refsinmem2allrefs[counter]] = -1;
     }
     pointer_refsinmem2allrefs[counter] = refno;
     fP_ref[counter] = fP;
     stddev_ref[counter] = stddev;
     proj_ref[counter] = img();
     //#define DEBUG
 #ifdef DEBUG
 
     std::cerr<<"counter= "<<counter<<"refno= "<<refno<<" stddev = "<<stddev;
     std::cerr<<" refsinmem2allrefs= "<<pointer_refsinmem2allrefs[counter];
     std::cerr<<" allrefs2refsinmem= "<<pointer_allrefs2refsinmem[pointer_refsinmem2allrefs[counter]] <<std::endl;
     std::cerr << "pointer_allrefs2refsinmem" <<std::endl;
     for (std::vector<int>::iterator i = pointer_allrefs2refsinmem.begin();
          i != pointer_allrefs2refsinmem.end();
          ++i)
         std::cerr << *i << " ";
     std::cerr <<std::endl;
     std::cerr << "pointer_refsinmem2allrefs" <<std::endl;
     for (std::vector<int>::iterator i = pointer_refsinmem2allrefs.begin();
          i != pointer_refsinmem2allrefs.end();
          ++i)
         std::cerr << *i << " ";
     std::cerr <<std::endl;
 #endif
 
     counter_refs_in_memory++;
     pthread_mutex_unlock(  &update_refs_in_memory_mutex );
     //    local_transformer.cleanup();
 }

◆ processAllImages()

void ProgAngularProjectionMatching::processAllImages ( )

virtual

Get images to process. This function will return the id's of images to process. It will be specially useful for MPI case when images will be distributed by the master node. Return false if there aren't more images to process

Reimplemented in MpiProgAngularProjectionMatching.

Definition at line 977 of file angular_projection_matching.cpp.

 {
     //Init progress bar
     size_t total_number_of_images = DFexp.size();
     if (verbose)
     {
         progress_bar_step = XMIPP_MAX(1, total_number_of_images / 80);
         init_progress_bar(total_number_of_images);
     }
     processSomeImages(ids);
     if (verbose)
         progress_bar(total_number_of_images);
 }

◆ processSomeImages()

void ProgAngularProjectionMatching::processSomeImages ( const std::vector< size_t > & imagesToProcess )

Loop over all images

a

Definition at line 991 of file angular_projection_matching.cpp.

 {
     Image<double> img;
     //double
     //opt_rot=0.,
             //opt_tilt=0.;
                      //opt_psi=0.,
                              //opt_xoff=0.,
                                       //opt_yoff=0.,
                                                //opt_scale=0.;
                                                          //maxcorr=-99.e99;
                                                          
     auto opt_flip  = std::vector<bool> (numOrientations, false);
     auto opt_refno = std::vector<int> (numOrientations, 0);
     auto opt_psi   = std::vector<double> (numOrientations, 0);
     auto maxcorr   = std::vector<double> (numOrientations, 0);
     auto opt_xoff  = std::vector<double> (numOrientations, 0);
     auto opt_yoff  = std::vector<double> (numOrientations, 0);
     auto opt_scale = std::vector<double> (numOrientations, 0);
     auto opt_rot   = std::vector<double> (numOrientations, 0);
     auto opt_tilt  = std::vector<double> (numOrientations, 0);
 
     size_t itemId=0;
     size_t nr_images = imagesToProcess.size();
     size_t idNew, imgid;
     FileName fn;
 
     // Call threads to calculate the rotational alignment of each image in the selfile
     auto * th_ids = (pthread_t *)malloc( threads * sizeof( pthread_t));
 
     // Allocate threads.
     auto * threads_d = (structThreadRotationallyAlignOneImage *)
             malloc ( threads * sizeof( structThreadRotationallyAlignOneImage ) );
 
     // Allocate threads vectors.
     for( int c = 0 ; c < threads ; c++ )
     {
         threads_d[c].opt_refno = (int *)malloc (sizeof(int)*numOrientations);
         threads_d[c].opt_psi =  (double *)malloc (sizeof(double)*numOrientations);
         threads_d[c].opt_flip =  (bool *) malloc (sizeof(bool)*numOrientations);
         threads_d[c].maxcorr =  (double *)malloc (sizeof(double)*numOrientations);
         threads_d[c].numOrientations = numOrientations;
     }
 
     for (size_t imgno = 0; imgno < nr_images; imgno++)
     {
         imgid = imagesToProcess[imgno];
         
         FileName pp;
         DFexp.getValue(MDL_IMAGE,pp, imgid);
 
         getCurrentImage(imgid, img);
         for( int c = 0 ; c < threads ; c++ )
         {
             threads_d[c].thread_id = c;
             threads_d[c].prm = this;
             threads_d[c].img = &img();
             threads_d[c].this_image = imgid;
             for (size_t i=0; i<numOrientations ;i++)
             {
                 threads_d[c].maxcorr[i] = -99.e99;
                 threads_d[c].opt_refno[i] = -1;
                 opt_scale[i] = 1;
                 opt_refno[i] = -1;
             }
             pthread_create( (th_ids+c), nullptr, threadRotationallyAlignOneImage, (void *)(threads_d+c) );
         }
         // Wait for threads to finish
         for( int c = 0 ; c < threads ; c++ )
             pthread_join(*(th_ids+c),nullptr);
 
         //Get optimal refno, psi, flip and maxcorr
         auto * indexThreads = new int[threads](); // init to zero
         size_t bestThreadCorr = 0;
         double tempCorr;
         size_t counterValidCorrs = 0;
         bool validCorr = false;
         for( int n = 0 ; n < numOrientations ; n++ )
         {
             tempCorr = -99.e99;
             validCorr = false;
             for( int c = 0 ; c < threads ; c++ )
             {
                 if (threads_d[c].maxcorr[indexThreads[c]] > tempCorr)
                 {
                     if (not validCorr)
                     {
                         counterValidCorrs++;
                         validCorr = true;
                     }
 
                     bestThreadCorr = c;
                     tempCorr = threads_d[c].maxcorr[indexThreads[c]];
                 }
 
             }
 
             if (not validCorr)
                 break;
 
             opt_refno[n] = threads_d[bestThreadCorr].opt_refno[indexThreads[bestThreadCorr]];
             opt_psi[n] = threads_d[bestThreadCorr].opt_psi[indexThreads[bestThreadCorr]];
             opt_flip[n] = threads_d[bestThreadCorr].opt_flip[indexThreads[bestThreadCorr]];
             maxcorr[n] = threads_d[bestThreadCorr].maxcorr[indexThreads[bestThreadCorr]];
             indexThreads[bestThreadCorr]++;
 
 
             if (opt_refno[n]>=0)
             {
                 opt_rot[n]  = XX(mysampling.no_redundant_sampling_points_angles[convert_refno_to_stack_position[opt_refno[n]]]);
                 opt_tilt[n] = YY(mysampling.no_redundant_sampling_points_angles[convert_refno_to_stack_position[opt_refno[n]]]);
             }
             else
             {
                 opt_rot[n]=0;
                 opt_tilt[n]=0;
             }
         }
 
         delete[] indexThreads;
         // Flip order to loop through references
         loop_forward_refs = !loop_forward_refs;
 
         //#define DEBUG
 #ifdef DEBUG
         std::cerr << "DEBUG_ROB, img.name():" << img.name() << std::endl;
 #endif
 #undef DEBUG
 
         for(int n=0; n < counterValidCorrs; n++)
         {
             //std::cout << " " << std::endl;
             //std::cout << maxcorr[n] << std::endl;
             translationallyAlignOneImage(img(),
                                          opt_refno[n],
                                          opt_psi[n],
                                          opt_flip[n],
                                          opt_xoff[n],
                                          opt_yoff[n],
                                          maxcorr[n]);
 
             //std::cout << maxcorr[n] << std::endl;
             if(do_scale && scale_nsteps > 0)
             {
                 // Compute a better scale (scale_min -> scale_max)
                 scaleAlignOneImage(img(), opt_refno[n], opt_psi[n], opt_flip[n], opt_xoff[n], opt_yoff[n], img.scale(), opt_scale[n], maxcorr[n]);
             }
 
             //Add the previously applied scale to the newly found one
             opt_scale[n] *= img.scale();
             // Divide opt_shifts by old_scale
 
             // Add previously applied translation to the newly found one
             opt_xoff[n] += img.Xoff();
             opt_yoff[n] += img.Yoff();
 
             // Output
             DFexp.getValue(MDL_IMAGE, fn, imgid);
             DFexp.getValue(MDL_ITEM_ID, itemId, imgid);
 
             idNew = DFo.addObject();
             DFo.setValue(MDL_ITEM_ID,itemId,idNew);
             DFo.setValue(MDL_IMAGE, fn,idNew);
             DFo.setValue(MDL_ANGLE_ROT, opt_rot[n],idNew);
             DFo.setValue(MDL_ANGLE_TILT,opt_tilt[n],idNew);
             DFo.setValue(MDL_ANGLE_PSI, opt_psi[n],idNew);
             
             //opt_xoff=0;
             DFo.setValue(MDL_SHIFT_X,   opt_xoff[n],idNew);
             DFo.setValue(MDL_SHIFT_Y,   opt_yoff[n],idNew);
             DFo.setValue(MDL_REF,      opt_refno[n] /*+ FIRST_IMAGE*/,idNew);
             DFo.setValue(MDL_FLIP,     opt_flip[n],idNew);
             DFo.setValue(MDL_SCALE,    opt_scale[n],idNew);
             DFo.setValue(MDL_MAXCC,    maxcorr[n],idNew);
 
         }
 
         /* FIXME ROB     */
         //opt_psi += img.psi();
         /*         */
 
         //opt_scale=1.0;
 
         if (verbose && imgno % progress_bar_step == 0)
             progress_bar(imgno);
         //#endif
         //DFo.write("/dev/stderr");
     }//loop over images
     free(th_ids);
 
     for( int c = 0 ; c < threads ; c++ )
     {
         free(threads_d[c].opt_refno);
         free(threads_d[c].opt_psi);
         free(threads_d[c].opt_flip);
         free(threads_d[c].maxcorr);
     }
 
     free(threads_d);
 
 }//function processSomeImages

◆ produceSideInfo()

void ProgAngularProjectionMatching::produceSideInfo ( )

virtual

Make shiftmask and calculate nr_psi

Reimplemented in MpiProgAngularProjectionMatching.

Definition at line 209 of file angular_projection_matching.cpp.

 {
 
     Image<double>    img,empty;
     Projection       proj;
     MetaDataVec      DF;
     MetaDataVec      SFr,emptySF;
     SymList          SL;
     FileName         fn_img;
     MultidimArray<double> Maux;
     MultidimArray<double> dataline(3);
     Polar<double>    P;
     Polar<std::complex <double> > fP;
 
     // Read Selfile and get dimensions
     MetaDataDb MetaDataLeft;
     MetaDataDb MetaDataRight;
     String blockName = "all_exp_images";
 
     FileName inputFn = fn_exp.removeBlockName();
     if (existsBlockInMetaDataFile(inputFn, blockName))
     {
         FileName rightFn;
         rightFn.compose(blockName, inputFn);
         MetaDataRight.read(rightFn);
         MetaDataLeft.read(fn_exp);
         //Join with angles and shifhts
         DFexp.join1(MetaDataLeft, MetaDataRight, MDL_IMAGE,LEFT);
     }
     else
         DFexp.read(fn_exp);
     //DFexp.removeDisabled();
 
 
     // Thread barrier
     barrier_init(&thread_barrier, threads);
 
     // Read one image to get dim
     DFexp.getValue(MDL_IMAGE,fn_img,DFexp.firstRowId());
     img.read(fn_img);
     dim = XSIZE(img());
 
     // Check that the reference and the experimental images are of the same
     // size
     FileName fnt;
     fnt.compose(FIRST_IMAGE, fn_ref);
     Image<double> imgRef;
     imgRef.read(fnt,HEADER);
 
     if (!imgRef().sameShape(img()))
         REPORT_ERROR(ERR_MULTIDIM_SIZE,
                      "Check that the reference volume and the experimental images are of the same size");
 
     // Set padding dimension
     paddim=ROUND(pad*dim);
 
     // Set max_shift
     if (max_shift<0)
         max_shift = dim/2;
 
     // Set ring defaults
     if (Ri<1)
         Ri=1;
     if (Ro<0)
         Ro=(dim/2)-1;
 
     // Calculate necessary memory per image
     produceSplineCoefficients(xmipp_transformation::BSPLINE3,Maux,img());
     P.getPolarFromCartesianBSpline(Maux,Ri,Ro);
     P.calculateFftwPlans(global_plans);
     fourierTransformRings(P,fP,global_plans,false);
     double memory_per_ref = 0.;
     for (int i = 0; i < fP.getRingNo(); i++)
     {
         memory_per_ref += (double) fP.getSampleNo(i) * 2 * sizeof(double);
     }
     memory_per_ref += dim * dim * sizeof(double);
     max_nr_imgs_in_memory = ROUND( 1024 * 1024 * 1024 * avail_memory / memory_per_ref);
 
     // Set up angular sampling
     mysampling.readSamplingFile(fn_ref.removeAllExtensions(),false);
     total_nr_refs = mysampling.no_redundant_sampling_points_angles.size();
 
     //#define DEBUG
 #ifdef DEBUG
 
     std::cerr << "XXXXmysampling.no_redundant_sampling_points_indexXXXXXX" <<std::endl;
     for (std::vector<size_t>::iterator i =
              mysampling.no_redundant_sampling_points_index.begin();
          i != mysampling.no_redundant_sampling_points_index.end();
          ++i)
         std::cerr << *i << " ";
     std::cerr << std::endl;
     std::cerr << "exit" <<std::endl;
 #endif
 #undef DEBUG
 
     convert_refno_to_stack_position.resize(mysampling.numberSamplesAsymmetricUnit, -1);
     for (size_t i = 0; i < mysampling.no_redundant_sampling_points_index.size(); i++)
         convert_refno_to_stack_position[mysampling.no_redundant_sampling_points_index[i]] = i;
 
     // Don't reserve more memory than necessary
     max_nr_refs_in_memory = XMIPP_MIN(max_nr_imgs_in_memory, total_nr_refs);
 
     // Initialize pointers for reference retrieval
     pointer_allrefs2refsinmem.resize(mysampling.numberSamplesAsymmetricUnit,-1);
     pointer_refsinmem2allrefs.resize(max_nr_refs_in_memory,-1);
     counter_refs_in_memory = 0;
     loop_forward_refs=true;
 
     // Initialize 5D search vectors
     search5d_xoff.clear();
     search5d_yoff.clear();
     // Make sure origin is included
     if(search5d_step == 0)
     {
         printf("***********************************************************\n");
         printf("*   ERROR: search step should be different from 0\n");
         printf("*   search step set to 1 \n");
         printf("***********************************************************\n");
 
         search5d_step = 1;
     }
     int myfinal=search5d_shift + search5d_shift%search5d_step;
     nr_trans = 0;//number translations in 5D
     for (int xoff = -myfinal; xoff <= myfinal; xoff+= search5d_step)
     {
         for (int yoff = -myfinal; yoff <= myfinal; yoff+= search5d_step)
         {
             // Only take a circle (not a square)
             if ( xoff*xoff + yoff*yoff <= search5d_shift*search5d_shift)
             {
                 search5d_xoff.push_back(xoff);
                 search5d_yoff.push_back(yoff);
                 nr_trans++;
             }
         }
     }
 
     // Initialize all arrays
     try
     {
         fP_ref = new Polar<std::complex<double> >[max_nr_refs_in_memory];
         proj_ref = new MultidimArray<double>[max_nr_refs_in_memory];
         fP_img = new Polar<std::complex<double> >[nr_trans];
         fPm_img = new Polar<std::complex<double> >[nr_trans];
 
         stddev_ref = new double[max_nr_refs_in_memory];
         stddev_img = new double[nr_trans];
     }
     catch (std::bad_alloc&)
     {
         REPORT_ERROR(ERR_MEM_BADREQUEST,"Error allocating memory in produceSideInfo");
     }
 
     // CTF stuff
     if (fn_ctf != "")
     {
         if (!fn_ctf.isMetaData())
         {
             Image<double> img;
             img.read(fn_ctf);
             Mctf=img();
             if (XSIZE(Mctf) != paddim)
             {
                 std::cerr<<"image size= "<<dim<<" padding factor= "<<pad<<" padded image size= "<<paddim<<" Wiener filter size= "<<XSIZE(Mctf)<<std::endl;
                 REPORT_ERROR(ERR_VALUE_INCORRECT,
                              "Incompatible padding factor for this CTF filter");
             }
         }
         else
         {
             CTFDescription ctf;
             MultidimArray<std::complex<double> >  ctfmask;
             ctf.read(fn_ctf);
             if (ABS(ctf.DeltafV - ctf.DeltafU) >1.)
             {
                 REPORT_ERROR(ERR_VALUE_INCORRECT,
                              "ERROR!! Only non-astigmatic CTFs are allowed!");
             }
             ctf.enable_CTF = true;
             ctf.produceSideInfo();
             ctf.generateCTF(paddim, paddim, ctfmask);
             Mctf.resize(paddim,paddim);
             FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(Mctf)
             {
                 if (phase_flipped)
                     dAij(Mctf, i, j) = fabs(dAij(ctfmask, i, j).real());
                 else
                     dAij(Mctf, i, j) = dAij(ctfmask, i, j).real();
             }
         }
     }
 
     //Store the id's of each experimental image from metadata
     DFexp.findObjects(ids);
 }

◆ readParams()

void ProgAngularProjectionMatching::readParams ( )

virtual

Read arguments from command line.

Reimplemented from XmippProgram.

Reimplemented in MpiProgAngularProjectionMatching.

Definition at line 44 of file angular_projection_matching.cpp.

 {
     fn_exp  = getParam("-i");
     fn_out  = getParam("-o");
     fn_ref  = getParam("--ref");
 
     // Additional commands
     pad=XMIPP_MAX(1.,getDoubleParam("--pad"));
     Ri=getIntParam("--Ri");
     Ro=getIntParam("--Ro");
     search5d_shift  = getIntParam("--search5d_shift");
     search5d_step = getIntParam("--search5d_step");
     max_shift = getDoubleParam("--max_shift");
     numOrientations = getIntParam("--number_orientations");
 
     avail_memory = getDoubleParam("--mem");
     if (checkParam("--ctf"))
         fn_ctf  = getParam("--ctf");
     phase_flipped = checkParam("--phase_flipped");
     threads = getIntParam("--thr");
 
     do_scale = checkParam("--scale");
     if (checkParam("--append"))
         do_overwrite = MD_APPEND;
     else
         do_overwrite = MD_OVERWRITE;
 
     if(do_scale)
     {
         scale_step = getDoubleParam("--scale",0);
         //Scale Step can't be 0
         if (scale_step == 0)
             scale_step = 1;
         scale_nsteps = getDoubleParam("--scale",1);
     }
 
 }

◆ rotationallyAlignOneImage()

void ProgAngularProjectionMatching::rotationallyAlignOneImage	(	Matrix2D< double > &	img,
		int	imgno,
		int &	opt_samplenr,
		double &	opt_psi,
		bool &	opt_flip,
		double &	maxcorr
	)

Rotational alignment using polar coordinates The input image is assumed to be in FTs of polar rings

◆ run()

void ProgAngularProjectionMatching::run ( )

virtual

Run.

Reimplemented from XmippProgram.

Definition at line 182 of file angular_projection_matching.cpp.

 {
     produceSideInfo();
 
     show();
 
     processAllImages();
 
     writeOutputFiles();
 
     destroyAndClean();
     if (verbose)
         std::cout << "done!"<<std::endl;
 }

◆ scaleAlignOneImage()

void ProgAngularProjectionMatching::scaleAlignOneImage	(	MultidimArray< double > &	img,
		const int &	samplenr,
		const double &	psi,
		const bool &	opt_flip,
		const double &	opt_xoff,
		const double &	opt_yoff,
		const double &	old_scale,
		double &	opt_scale,
		double &	maxcorr
	)

Translational alignment using cartesian coordinates The optimal direction is re-projected from the volume

a

Definition at line 871 of file angular_projection_matching.cpp.

 {
     MultidimArray<double> Mscale,Mtrans,Mref,Mimg;
     int refno;
 
     Mscale.setXmippOrigin();
     Mtrans.setXmippOrigin();
     Mref.setXmippOrigin();
     Mimg.setXmippOrigin();
 
 #ifdef TIMING
 
     TimeStamp t0,t1,t2;
     time_config();
     annotate_time(&t0);
 #endif
 
     // ROTATE + SHIFT
     // Transformation matrix
     Matrix2D<double> A(3,3);
     A.initIdentity();
     double ang, cosine, sine;
     ang = DEG2RAD(opt_psi);
     cosine = cos(ang);
     sine = sin(ang);
 
     // Rotation
     MAT_ELEM(A,0, 0) = cosine;
     MAT_ELEM(A,0, 1) = sine;
     MAT_ELEM(A,1, 0) = -sine;
     MAT_ELEM(A,1, 1) = cosine;
 
     // Shift
     MAT_ELEM(A,0, 2) = -opt_xoff;
     MAT_ELEM(A,1, 2) = -opt_yoff;
 
     // Multiply shifts by old_scale
 
     if (opt_flip)
     {
         MAT_ELEM(A,0, 2) = opt_xoff;
         MAT_ELEM(A,0, 0) *= -1.;
         MAT_ELEM(A,0, 1) *= -1.;
     }
 
     // Get pointer to the correct reference image in memory
     refno = pointer_allrefs2refsinmem[opt_refno];
     if (refno == -1)
     {
         // Reference is not stored in memory (anymore): (re-)read from disc
         getCurrentReference(opt_refno,global_plans);
         refno = pointer_allrefs2refsinmem[opt_refno];
     }
     applyGeometry(xmipp_transformation::LINEAR, Mref, proj_ref[refno], A, xmipp_transformation::IS_NOT_INV, xmipp_transformation::DONT_WRAP);
 
 
     Mtrans = img;
 
     // Scale search
     double corr;
     opt_scale = 1;
     //maxcorr = -999;
 
     // 1 (0.01 * scale_step * scale_nsteps)
     double scale = 1 - 0.01 * scale_step * scale_nsteps;
     while(scale <= 1 + 0.01 * scale_step * scale_nsteps)
     {
         if(scale == 1.)
             continue;
 
         // apply current scale
         A.initIdentity();
         A /= scale;
         applyGeometry(xmipp_transformation::LINEAR, Mscale, Mtrans, A, xmipp_transformation::IS_INV, xmipp_transformation::DONT_WRAP);
 
         //Image spread correction (if scale != 1) for scale search
         Mscale = Mscale / (scale * scale * old_scale * old_scale);
 
         corr = correlationIndex(Mref,Mscale);
 
         // best scale update
         if(corr > maxcorr)
         {
             opt_scale = scale;
             maxcorr = corr;
         }
         scale += 0.01 * scale_step;
     }
 
 #ifdef TIMING
 
     float total_trans = elapsed_time(t0);
     std::cerr<<" trans%% "<<total_trans <<std::endl;
 #endif
 
 }

◆ show()

void ProgAngularProjectionMatching::show ( )

Show.

Definition at line 124 of file angular_projection_matching.cpp.

 {
     if (!verbose)
         return;
 
     std::cout
     << "  Input images            : "<< fn_exp << std::endl
     << "  Output file             : "<< fn_out << std::endl
     ;
     if (Ri>0)
         std::cout << "  Inner radius rot-search : " << Ri<< std::endl;
     if (Ro>0)
         std::cout << "  Outer radius rot-search : " << Ro << std::endl;
     if (max_nr_refs_in_memory<total_nr_refs)
     {
         std::cout << "  Number of references    : " << total_nr_refs << std::endl
         << "  Nr. refs in memory      : " << max_nr_refs_in_memory << " (using " << avail_memory <<" Gb)" << std::endl
         ;
     }
     else
     {
         std::cout << "  Number of references    : " << total_nr_refs << " (all stored in memory)" << std::endl;
     }
     std::cout << "  Max. allowed shift      : +/- " <<max_shift<<" pixels"<<std::endl;
     if (search5d_shift > 0)
     {
         std::cout << "  5D-search shift range   : "<<search5d_shift<<" pixels (sampled "<<nr_trans<<" times)"<<std::endl;
     }
     if (fn_ctf!="")
     {
         if (!fn_ctf.isMetaData())
         {
             std::cout << "  CTF image               :  " <<fn_ctf<<std::endl;
             if (pad > 1.)
                 std::cout << "  Padding factor          : "<< pad << std::endl;
         }
         else
         {
             std::cout << "  CTF parameter file      :  " <<fn_ctf<<std::endl;
             if (phase_flipped)
                 std::cout << "    + Assuming images have been phase flipped " << std::endl;
             else
                 std::cout << "    + Assuming images have not been phase flipped " << std::endl;
         }
     }
     if (threads>1)
     {
         std::cout << "  -> Using "<<threads<<" parallel threads"<<std::endl;
     }
 
     if (numOrientations != 1)
          std::cout << "  -> Using "<<numOrientations<<" possible orientations for each particle"<<std::endl;
 
 
     std::cout << " ================================================================="<<std::endl;
 }

◆ translationallyAlignOneImage()

void ProgAngularProjectionMatching::translationallyAlignOneImage	(	MultidimArray< double > &	img,
		const int &	samplenr,
		const double &	psi,
		const bool &	opt_flip,
		double &	opt_xoff,
		double &	opt_yoff,
		double &	maxcorr
	)

Translational alignment using cartesian coordinates The optimal direction is re-projected from the volume

Definition at line 776 of file angular_projection_matching.cpp.

 {
 
     MultidimArray<double> Mtrans,Mimg,Mref;
     int refno;
     Mtrans.setXmippOrigin();
     Mimg.setXmippOrigin();
     Mref.setXmippOrigin();
 #ifdef TIMING
 
     TimeStamp t0,t1,t2;
     time_config();
     annotate_time(&t0);
 #endif
 
 #ifdef DEBUG
 
     std::cerr<<"start trans: opt_refno= "<<opt_refno<<" pointer= "<<pointer_allrefs2refsinmem[opt_refno]<<" opt_psi= "<<opt_psi<<"opt_flip= "<<opt_flip<<std::endl;
 #endif
 
     // Get pointer to the correct reference image in memory
     refno = pointer_allrefs2refsinmem[opt_refno];
     if (refno == -1)
     {
         // Reference is not stored in memory (anymore): (re-)read from disc
         getCurrentReference(opt_refno,global_plans);
         refno = pointer_allrefs2refsinmem[opt_refno];
     }
 
     // Rotate stored reference projection by phi degrees
     rotate(xmipp_transformation::BSPLINE3,Mref,proj_ref[refno],opt_psi,xmipp_transformation::DONT_WRAP);
     //rotate(BSPLINE3,Mref,proj_ref[refno],-opt_psi,DONT_WRAP);
 
 #ifdef DEBUG
 
     std::cerr<<"rotated ref "<<std::endl;
 #endif
 
     if (opt_flip)
     {
         // Flip experimental image
         Matrix2D<double> A(3,3);
         A.initIdentity();
         MAT_ELEM(A,0, 0) = -1.;
         //MAT_ELEM(A,0, 1) *= -1.;
         applyGeometry(xmipp_transformation::LINEAR, Mimg, img, A, xmipp_transformation::IS_INV, xmipp_transformation::DONT_WRAP);
     }
     else
         Mimg = img;
 
 
     // Perform the actual search for the optimal shift
     if (max_shift>0)
     {
         CorrelationAux aux;
         bestShift(Mref,Mimg,opt_xoff,opt_yoff,aux);
     }
     else
         opt_xoff = opt_yoff = 0.;
     if (opt_xoff * opt_xoff + opt_yoff * opt_yoff > max_shift * max_shift)
         opt_xoff = opt_yoff = 0.;
     //#define DEBUG
 #ifdef DEBUG
 
     std::cerr<<"optimal shift "<<opt_xoff<<" "<<opt_yoff<<std::endl;
 #endif
 
     // Calculate standard cross-correlation coefficient
     translate(xmipp_transformation::LINEAR,Mtrans,Mimg,vectorR2(opt_xoff,opt_yoff),true);
     maxcorr = correlationIndex(Mref,Mtrans);
 
 #ifdef DEBUG
 
     std::cerr<<"optimal shift corr "<<maxcorr<<std::endl;
 #endif
 #undef DEBUG
     // Correct X-shift for mirrored images
     if (opt_flip)
         opt_xoff *= -1.;
 
 #ifdef TIMING
 
     float total_trans = elapsed_time(t0);
     std::cerr<<" trans%% "<<total_trans <<std::endl;
 #endif
 
 }

◆ writeOutputFiles()

void ProgAngularProjectionMatching::writeOutputFiles ( )

virtual

Write out results to disk This function should be override in MPI class, only master should write.

Reimplemented in MpiProgAngularProjectionMatching.

Definition at line 1238 of file angular_projection_matching.cpp.

 {
     DFo.write(fn_out,do_overwrite);
 }

Member Data Documentation

◆ avail_memory

double ProgAngularProjectionMatching::avail_memory

Available memory for storage of all references (in Gb)

Definition at line 84 of file angular_projection_matching.h.

◆ convert_refno_to_stack_position

std::vector<size_t> ProgAngularProjectionMatching::convert_refno_to_stack_position

Vector to assign reference number to stack positions

Definition at line 100 of file angular_projection_matching.h.

◆ counter_refs_in_memory

int ProgAngularProjectionMatching::counter_refs_in_memory

Counter for current filling of memory with references

Definition at line 95 of file angular_projection_matching.h.

◆ DFexp

MetaDataDb ProgAngularProjectionMatching::DFexp

Docfile with experimental images

Definition at line 72 of file angular_projection_matching.h.

◆ DFo

MetaDataDb ProgAngularProjectionMatching::DFo

Docfile with results

Definition at line 74 of file angular_projection_matching.h.

◆ dim

size_t ProgAngularProjectionMatching::dim

dimension of the images and padded images

Definition at line 76 of file angular_projection_matching.h.

◆ do_overwrite

WriteModeMetaData ProgAngularProjectionMatching::do_overwrite

Definition at line 137 of file angular_projection_matching.h.

◆ do_scale

bool ProgAngularProjectionMatching::do_scale

scale params

Definition at line 136 of file angular_projection_matching.h.

◆ fn_ctf

FileName ProgAngularProjectionMatching::fn_ctf

Definition at line 70 of file angular_projection_matching.h.

◆ fn_exp

FileName ProgAngularProjectionMatching::fn_exp

Filenames

Definition at line 70 of file angular_projection_matching.h.

◆ fn_out

FileName ProgAngularProjectionMatching::fn_out

Definition at line 70 of file angular_projection_matching.h.

◆ fn_ref

FileName ProgAngularProjectionMatching::fn_ref

Definition at line 70 of file angular_projection_matching.h.

◆ fP_img

Polar<std::complex<double> > * ProgAngularProjectionMatching::fP_img

Definition at line 104 of file angular_projection_matching.h.

◆ fP_ref

Polar<std::complex<double> >* ProgAngularProjectionMatching::fP_ref

Array with Polars of references and of translated images and their mirrors

Definition at line 104 of file angular_projection_matching.h.

◆ fPm_img

Polar<std::complex<double> > * ProgAngularProjectionMatching::fPm_img

Definition at line 104 of file angular_projection_matching.h.

◆ global_plans

Polar_fftw_plans ProgAngularProjectionMatching::global_plans

Global plans for fftw transformers of all polar rings

Definition at line 108 of file angular_projection_matching.h.

◆ ids

std::vector<size_t> ProgAngularProjectionMatching::ids

Array containing the images ids in metadata

Definition at line 102 of file angular_projection_matching.h.

◆ loop_forward_refs

bool ProgAngularProjectionMatching::loop_forward_refs

Flag whether to loop from low to high or from high to low through the references

Definition at line 115 of file angular_projection_matching.h.

◆ max_nr_imgs_in_memory

int ProgAngularProjectionMatching::max_nr_imgs_in_memory

Maximum number of references that can be stored in memory The difference between this and the previous varible is that the previous one is the MIN(max_nr_refs_in_memory, number of references)

Definition at line 91 of file angular_projection_matching.h.