#include <reconstruct_fourier.h>

Inheritance diagram for ProgRecFourier:

Collaboration diagram for ProgRecFourier:

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

void	defineParams ()
	Read arguments from command line. More...

void	show ()

void	run ()

void	produceSideinfo ()
	Produce side info: fill arrays with relevant transformation matrices. More...

void	finishComputations (const FileName &out_name)

void	processImages (int firstImageIndex, int lastImageIndex, bool saveFSC=false, bool reprocessFlag=false)
	Process one image. More...

void	correctWeight ()
	Method for the correction of the fourier coefficients. More...

void	forceWeightSymmetry (MultidimArray< double > &FourierWeights)
	Force the weights to be symmetrized. More...

virtual void	setIO (const FileName &fn_in, const FileName &fn_out)
	Functions of common reconstruction interface. More...

Public Member Functions inherited from ProgReconsBase
virtual	~ProgReconsBase ()

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 ()

Static Public Member Functions
static void *	processImageThread (void *threadArgs)
	Defines what a thread should do. More...

Public Attributes
FileName	fn_out

FileName	fn_sym

FileName	fn_sel

FileName	fn_doc

FileName	fn_fsc

MetaDataVec	SF

bool	do_weights

bool	useCTF

bool	phaseFlipped

double	minCTF

double	Ts

double	padding_factor_proj

double	padding_factor_vol

double	sampling_rate

double	maxResolution
	Max resolution in Angstroms. More...

int	NiterWeight
	Number of iterations for the weight. More...

int	numThreads
	Number of threads to use in parallel to process a single image. More...

pthread_t *	th_ids
	IDs for the threads. More...

ImageThreadParams *	th_args
	Contains parameters passed to each thread. More...

int	threadOpCode
	Tells the threads what to do next. More...

size_t	rowsProcessed
	Number of rows already processed on an image. More...

pthread_mutex_t	workLoadMutex
	Controls mutual exclusion on critical zones of code. More...

barrier_t	barrier
	To create a barrier synchronization for threads. More...

int *	statusArray
	A status array for each row in an image (processing, processed,etc..) More...

int	thrWidth
	How many image rows are processed at a time by a single thread. More...

int	imgSize

int	volPadSizeX

int	volPadSizeY

int	volPadSizeZ

Matrix1D< double >	Fourier_blob_table

Matrix1D< double >	blobTableSqrt

Matrix1D< double >	fourierBlobTableSqrt

double	iDeltaFourier

double	iDeltaSqrt

double	maxResolution2

struct blobtype	blob

std::vector< Matrix2D< double > >	R_repository

FourierTransformer	transformerVol

FourierTransformer	transformerImg

MultidimArray< std::complex< double > >	VoutFourier

MultidimArray< std::complex< double > >	VoutFourierTmp

MultidimArray< double >	FourierWeights

MultidimArray< double >	paddedImg

Image< double >	Vout

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

Fourier reconstruction parameters.

Definition at line 76 of file reconstruct_fourier.h.

Member Function Documentation

◆ correctWeight()

void ProgRecFourier::correctWeight ( )

Method for the correction of the fourier coefficients.

Definition at line 1056 of file reconstruct_fourier.cpp.

 {
     // If NiterWeight=0 then set the weights to one
     forceWeightSymmetry(FourierWeights);
     if (NiterWeight==0)
     {
         FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(FourierWeights)
         DIRECT_A3D_ELEM(FourierWeights, k,i,j)=1;
 
     }
     else
     {
         // Temporary save the Fourier of the volume
         MultidimArray< std::complex<double> > VoutFourierTmp;
         VoutFourierTmp=VoutFourier;
         forceWeightSymmetry(FourierWeights);
         // Prepare the VoutFourier
         FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(VoutFourier)
         {
             auto *ptrOut=(double *)&(DIRECT_A3D_ELEM(VoutFourier, k,i,j));
             if (fabs(A3D_ELEM(FourierWeights,k,i,j))>1e-3)
                 ptrOut[0] = 1.0/DIRECT_A3D_ELEM(FourierWeights, k,i,j);
         }
 
         for (int i=1;i<NiterWeight;i++)
         {
             FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(FourierWeights)
             A3D_ELEM(FourierWeights,k,i,j)=0;
             processImages(0, SF.size() - 1, !fn_fsc.empty(), true);
             forceWeightSymmetry(FourierWeights);
             FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(VoutFourier)
             {
                 auto *ptrOut=(double *)&(DIRECT_A3D_ELEM(VoutFourier, k,i,j));
                 if (fabs(A3D_ELEM(FourierWeights,k,i,j))>1e-3)
                     ptrOut[0] /= A3D_ELEM(FourierWeights,k,i,j);
             }
         }
         FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(VoutFourier)
         {
             // Put back the weights to FourierWeights from temporary variable VoutFourier
             auto *ptrOut=(double *)&(DIRECT_A3D_ELEM(VoutFourier, k,i,j));
             A3D_ELEM(FourierWeights,k,i,j) = ptrOut[0];
         }
         VoutFourier = VoutFourierTmp;
     }
 }

◆ defineParams()

void ProgRecFourier::defineParams ( )

virtual

Read arguments from command line.

Reimplemented from XmippProgram.

Definition at line 36 of file reconstruct_fourier.cpp.

 {
     //usage
     addUsageLine("Generate 3D reconstructions from projections using direct Fourier interpolation with arbitrary geometry.");
     addUsageLine("Kaisser-windows are used for interpolation in Fourier space.");
     //params
     addParamsLine("   -i <md_file>                : Metadata file with input projections");
     addParamsLine("  [-o <volume_file=\"rec_fourier.vol\">]  : Filename for output volume");
     addParamsLine("  [--iter <iterations=1>]      : Number of iterations for weight correction");
     addParamsLine("  [--sym <symfile=c1>]              : Enforce symmetry in projections");
     addParamsLine("  [--padding <proj=2.0> <vol=2.0>]  : Padding used for projections and volume");
     addParamsLine("  [--prepare_fsc <fscfile>]      : Filename root for FSC files");
     addParamsLine("  [--max_resolution <p=0.5>]     : Max resolution (Nyquist=0.5)");
     addParamsLine("  [--weight]                     : Use weights stored in the image metadata");
     addParamsLine("  [--thr <threads=1> <rows=1>]   : Number of concurrent threads and rows processed at time by a thread");
     addParamsLine("  [--blob <radius=1.9> <order=0> <alpha=15>] : Blob parameters");
     addParamsLine("                                 : radius in pixels, order of Bessel function in blob and parameter alpha");
     addParamsLine("  [--useCTF]                     : Use CTF information if present");
     addParamsLine("  [--sampling <Ts=1>]            : sampling rate of the input images in Angstroms/pixel");
     addParamsLine("                                 : It is only used when correcting for the CTF");
     addParamsLine("  [--phaseFlipped]               : Give this flag if images have been already phase flipped");
     addParamsLine("  [--minCTF <ctf=0.01>]          : Minimum value of the CTF that will be inverted");
     addParamsLine("                                 : CTF values (in absolute value) below this one will not be corrected");
     addExampleLine("For reconstruct enforcing i3 symmetry and using stored weights:", false);
     addExampleLine("   xmipp_reconstruct_fourier  -i reconstruction.sel --sym i3 --weight");
 }

◆ finishComputations()

void ProgRecFourier::finishComputations ( const FileName & out_name )

Definition at line 1103 of file reconstruct_fourier.cpp.

 {
     //#define DEBUG_VOL
 #ifdef DEBUG_VOL
     {
         VolumeXmipp save;
         save().alias( FourierWeights );
         save.write((std::string) fn_out + "Weights.vol");
 
         FourierVolume save2;
         save2().alias( VoutFourier );
         save2.write((std::string) fn_out + "FourierVol.vol");
     }
 #endif
 
     // Enforce symmetry in the Fourier values as well as the weights
     // Sjors 19aug10 enforceHermitianSymmetry first checks ndim...
     Vout().initZeros(volPadSizeZ,volPadSizeY,volPadSizeX);
     transformerVol.setReal(Vout());
     transformerVol.enforceHermitianSymmetry();
     //forceWeightSymmetry(preFourierWeights);
 
     // Tell threads what to do
     //#define DEBUG_VOL1
 #ifdef DEBUG_VOL1
 
     {
         Image<double> save;
         save().alias( FourierWeights );
         save.write((std::string) fn_out + "hermiticWeights.vol");
 
         Image< std::complex<double> > save2;
         save2().alias( VoutFourier );
         save2.write((std::string) fn_out + "hermiticFourierVol.vol");
     }
 #endif
     threadOpCode = PROCESS_WEIGHTS;
     // Awake threads
     barrier_wait( &barrier );
     // Threads are working now, wait for them to finish
     barrier_wait( &barrier );
 
     transformerVol.inverseFourierTransform();
     CenterFFT(Vout(),false);
 
     // Correct by the Fourier transform of the blob
     Vout().setXmippOrigin();
     Vout().selfWindow(FIRST_XMIPP_INDEX(imgSize),FIRST_XMIPP_INDEX(imgSize),
                       FIRST_XMIPP_INDEX(imgSize),LAST_XMIPP_INDEX(imgSize),
                       LAST_XMIPP_INDEX(imgSize),LAST_XMIPP_INDEX(imgSize));
     double pad_relation= ((double)padding_factor_proj/padding_factor_vol);
     pad_relation = (pad_relation * pad_relation * pad_relation);
 
     MultidimArray<double> &mVout=Vout();
     double ipad_relation=1.0/pad_relation;
     double meanFactor2=0;
     FOR_ALL_ELEMENTS_IN_ARRAY3D(mVout)
     {
         double radius=sqrt((double)(k*k+i*i+j*j));
         double aux=radius*iDeltaFourier;
         double factor = Fourier_blob_table(ROUND(aux));
         double factor2=(pow(Sinc(radius/(2*imgSize)),2));
         if (NiterWeight!=0)
         {
             A3D_ELEM(mVout,k,i,j) /= (ipad_relation*factor2*factor);
             meanFactor2+=factor2;
         }
         else
             A3D_ELEM(mVout,k,i,j) /= (ipad_relation*factor);
     }
     if (NiterWeight!=0)
     {
         meanFactor2/=MULTIDIM_SIZE(mVout);
         FOR_ALL_ELEMENTS_IN_ARRAY3D(mVout)
         A3D_ELEM(mVout,k,i,j) *= meanFactor2;
     }
     Vout.write(out_name);
 }

◆ forceWeightSymmetry()

void ProgRecFourier::forceWeightSymmetry ( MultidimArray< double > & FourierWeights )

Force the weights to be symmetrized.

Definition at line 1188 of file reconstruct_fourier.cpp.

 {
     int yHalf=YSIZE(FourierWeights)/2;
     if (YSIZE(FourierWeights)%2==0)
         yHalf--;
     int zHalf=ZSIZE(FourierWeights)/2;
     if (ZSIZE(FourierWeights)%2==0)
         zHalf--;
     auto zsize=(int)ZSIZE(FourierWeights);
     int zsize_1=zsize-1;
     int ysize_1=(int)YSIZE(FourierWeights)-1;
     for (int k=0; k<zsize; k++)
     {
         int ksym=intWRAP(-k,0,zsize_1);
         for (int i=1; i<=yHalf; i++)
         {
             int isym=intWRAP(-i,0,ysize_1);
             double mean=0.5*(
                             DIRECT_A3D_ELEM(FourierWeights,k,i,0)+
                             DIRECT_A3D_ELEM(FourierWeights,ksym,isym,0));
             DIRECT_A3D_ELEM(FourierWeights,k,i,0)=
                 DIRECT_A3D_ELEM(FourierWeights,ksym,isym,0)=mean;
         }
     }
     for (int k=1; k<=zHalf; k++)
     {
         int ksym=intWRAP(-k,0,zsize_1);
         double mean=0.5*(
                         DIRECT_A3D_ELEM(FourierWeights,k,0,0)+
                         DIRECT_A3D_ELEM(FourierWeights,ksym,0,0));
         DIRECT_A3D_ELEM(FourierWeights,k,0,0)=
             DIRECT_A3D_ELEM(FourierWeights,ksym,0,0)=mean;
     }
 }

◆ processImages()

void ProgRecFourier::processImages	(	int	firstImageIndex,
		int	lastImageIndex,
		bool	saveFSC = `false`,
		bool	reprocessFlag = `false`
	)

Process one image.

Definition at line 835 of file reconstruct_fourier.cpp.

 {
     MultidimArray< std::complex<double> > *paddedFourier;
 
     auto repaint = (int)ceil((double)SF.size()/60);
 
     bool processed;
     int imgno = 0;
     int imgIndex = firstImageIndex;
 
     // This index tells when to save work for later FSC usage
     int FSCIndex = (firstImageIndex + lastImageIndex)/2;
 
     // Index of the image that has just been processed. Used for
     // FSC purposes
     int current_index;
 
     do
     {
         threadOpCode = PRELOAD_IMAGE;
 
         for ( int nt = 0 ; nt < numThreads ; nt ++ )
         {
             if ( imgIndex <= lastImageIndex )
             {
                 th_args[nt].imageIndex = imgIndex;
                 th_args[nt].reprocessFlag = reprocessFlag;
                 imgIndex++;
             }
             else
             {
                 th_args[nt].imageIndex = -1;
             }
         }
 
         // Awaking sleeping threads
         barrier_wait( &barrier );
         // here each thread is reading a different image and compute fft
         // Threads are working now, wait for them to finish
         // processing current projection
         barrier_wait( &barrier );
 
         // each threads have read a different image and now
         // all the thread will work in a different part of a single image.
         threadOpCode = PROCESS_IMAGE;
 
         processed = false;
 
         for ( int nt = 0 ; nt < numThreads ; nt ++ )
         {
             if ( th_args[nt].read == 2 )
                 processed = true;
             else if ( th_args[nt].read == 1 )
             {
                 processed = true;
                 if (verbose) {
                     if (imgno % repaint == 0)
                         progress_bar(imgno);
                     imgno++;
                 }
 
                 double weight = th_args[nt].localweight;
                 paddedFourier = th_args[nt].localPaddedFourier;
                 current_index = th_args[nt].imageIndex;
                 Matrix2D<double> *Ainv = th_args[nt].localAInv;
 
                 //#define DEBUG22
 #ifdef DEBUG22
 
                 {
                     static int ii=0;
                     if(ii%1==0)
                     {
                         FourierImage save22;
                         //save22()=*paddedFourier;
                         save22().alias(*paddedFourier);
                         save22.write((std::string) integerToString(ii)  + "_padded_fourier.spi");
                     }
                     ii++;
                 }
 #endif
                 #undef DEBUG22
 
                 // Initialized just once
                 if ( statusArray == nullptr )
                 {
                     statusArray = (int *) malloc ( sizeof(int) * paddedFourier->ydim );
                 }
 
                 // Determine how many rows of the fourier
                 // transform are of interest for us. This is because
                 // the user can avoid to explore at certain resolutions
                 auto conserveRows=(size_t)ceil((double)paddedFourier->ydim * maxResolution * 2.0);
                 conserveRows=(size_t)ceil((double)conserveRows/2.0);
 
                 // Loop over all symmetries
                 for (size_t isym = 0; isym < R_repository.size(); isym++)
                 {
                     rowsProcessed = 0;
 
                     // Compute the coordinate axes of the symmetrized projection
                     Matrix2D<double> A_SL=R_repository[isym]*(*Ainv);
 
                     // Fill the thread arguments for each thread
                     for ( int th = 0 ; th < numThreads ; th ++ )
                     {
                         // Passing parameters to each thread
                         th_args[th].symmetry = &A_SL;
                         th_args[th].paddedFourier = paddedFourier;
                         th_args[th].weight = weight;
                         th_args[th].reprocessFlag = reprocessFlag;
                     }
 
                     // Init status array
                     for (size_t i = 0 ; i < paddedFourier->ydim ; i ++ )
                     {
                         if ( i >= conserveRows && i < (paddedFourier->ydim-conserveRows))
                         {
                             // -2 means "discarded"
                             statusArray[i] = -2;
                             rowsProcessed++;
                         }
                         else
                         {
                             statusArray[i] = 0;
                         }
                     }
 
                     // Awaking sleeping threads
                     barrier_wait( &barrier );
                     // Threads are working now, wait for them to finish
                     // processing current projection
                     barrier_wait( &barrier );
 
                     //#define DEBUG2
 #ifdef DEBUG2
 
                     {
                         static int ii=0;
                         if(ii%1==0)
                         {
                             Image<double> save;
                             save().alias( FourierWeights );
                             save.write((std::string) integerToString(ii)  + "_1_Weights.vol");
 
                             Image< std::complex<double> > save2;
                             save2().alias( VoutFourier );
                             save2.write((std::string) integerToString(ii)  + "_1_Fourier.vol");
                         }
                         ii++;
                     }
 #endif
                     #undef DEBUG2
 
                 }
 
                 if ( current_index == FSCIndex && saveFSC )
                 {
                     // Save Current Fourier, Reconstruction and Weights
                     Image<double> save;
                     save().alias( FourierWeights );
                     save.write((std::string)fn_fsc + "_1_Weights.vol");
 
                     Image< std::complex<double> > save2;
                     save2().alias( VoutFourier );
                     save2.write((std::string) fn_fsc + "_1_Fourier.vol");
 
                     finishComputations(FileName((std::string) fn_fsc + "_1_recons.vol"));
                     Vout().initZeros(volPadSizeZ, volPadSizeY, volPadSizeX);
                     transformerVol.setReal(Vout());
                     Vout().clear();
                     transformerVol.getFourierAlias(VoutFourier);
                     FourierWeights.initZeros(VoutFourier);
                     VoutFourier.initZeros();
                 }
             }
         }
     }
     while ( processed );
 
     if( saveFSC )
     {
         // Save Current Fourier, Reconstruction and Weights
         Image<double> auxVolume;
         auxVolume().alias( FourierWeights );
         auxVolume.write((std::string)fn_fsc + "_2_Weights.vol");
 
         Image< std::complex<double> > auxFourierVolume;
         auxFourierVolume().alias( VoutFourier );
         auxFourierVolume.write((std::string) fn_fsc + "_2_Fourier.vol");
 
         finishComputations(FileName((std::string) fn_fsc + "_2_recons.vol"));
 
         Vout().initZeros(volPadSizeZ, volPadSizeY, volPadSizeX);
         transformerVol.setReal(Vout());
         Vout().clear();
         transformerVol.getFourierAlias(VoutFourier);
         FourierWeights.initZeros(VoutFourier);
         VoutFourier.initZeros();
 
         auxVolume.sumWithFile(fn_fsc + "_1_Weights.vol");
         auxVolume.sumWithFile(fn_fsc + "_2_Weights.vol");
         auxFourierVolume.sumWithFile(fn_fsc + "_1_Fourier.vol");
         auxFourierVolume.sumWithFile(fn_fsc + "_2_Fourier.vol");
         remove((fn_fsc + "_1_Weights.vol").c_str());
         remove((fn_fsc + "_2_Weights.vol").c_str());
         remove((fn_fsc + "_1_Fourier.vol").c_str());
         remove((fn_fsc + "_2_Fourier.vol").c_str());
 
         /*
         //Save SUM
                                     //this is an image but not an xmipp image
                                     auxFourierVolume.write((std::string)fn_fsc + "_all_Fourier.vol",
                                             false,VDOUBLE);
                                     auxVolume.write((std::string)fn_fsc + "_all_Weight.vol",
                                             false,VDOUBLE);
         //
         */
     }
 }

◆ processImageThread()

void * ProgRecFourier::processImageThread ( void * threadArgs )

static

Defines what a thread should do.

Definition at line 289 of file reconstruct_fourier.cpp.

 {
 
     auto * threadParams = (ImageThreadParams *) threadArgs;
     ProgRecFourier * parent = threadParams->parent;
     barrier_t * barrier = &(parent->barrier);
 
     int minSeparation;
 
     if ( (int)ceil(parent->blob.radius) > parent->thrWidth )
         minSeparation = (int)ceil(parent->blob.radius);
     else
         minSeparation = parent->thrWidth;
 
     minSeparation+=1;
 
     Matrix2D<double>  localA(3, 3), localAinv(3, 3);
     MultidimArray< std::complex<double> > localPaddedFourier;
     MultidimArray<double> localPaddedImg;
     FourierTransformer localTransformerImg;
 
     std::vector<size_t> objId;
 
     threadParams->selFile->findObjects(objId);
     ApplyGeoParams params;
     params.only_apply_shifts = true;
     MultidimArray<int> zWrapped(3*parent->volPadSizeZ),yWrapped(3*parent->volPadSizeY),xWrapped(3*parent->volPadSizeX),
     zNegWrapped, yNegWrapped, xNegWrapped;
     zWrapped.initConstant(-1);
     yWrapped.initConstant(-1);
     xWrapped.initConstant(-1);
     zWrapped.setXmippOrigin();
     yWrapped.setXmippOrigin();
     xWrapped.setXmippOrigin();
     zNegWrapped=zWrapped;
     yNegWrapped=yWrapped;
     xNegWrapped=xWrapped;
 
     MultidimArray<double> x2precalculated(XSIZE(xWrapped)), y2precalculated(XSIZE(yWrapped)), z2precalculated(XSIZE(zWrapped));
     x2precalculated.initConstant(-1);
     y2precalculated.initConstant(-1);
     z2precalculated.initConstant(-1);
     x2precalculated.setXmippOrigin();
     y2precalculated.setXmippOrigin();
     z2precalculated.setXmippOrigin();
 
     bool hasCTF=(threadParams->selFile->containsLabel(MDL_CTF_MODEL) || threadParams->selFile->containsLabel(MDL_CTF_DEFOCUSU)) &&
                 parent->useCTF;
     if (hasCTF)
     {
         threadParams->ctf.enable_CTF=true;
         threadParams->ctf.enable_CTFnoise=false;
     }
     do
     {
         barrier_wait( barrier );
 
         switch ( parent->threadOpCode )
         {
         case PRELOAD_IMAGE:
             {
 
                 threadParams->read = 0;
 
                 if ( threadParams->imageIndex >= 0 )
                 {
                     // Read input image
                     double rot, tilt, psi, weight;
                     Projection proj;
 
                     //Read projection from selfile, read also angles and shifts if present
                     //but only apply shifts
 
                     proj.readApplyGeo(*(threadParams->selFile), objId[threadParams->imageIndex], params);
                     rot  = proj.rot();
                     tilt = proj.tilt();
                     psi  = proj.psi();
                     weight = proj.weight();
                     if (hasCTF)
                     {
                         threadParams->ctf.readFromMetadataRow(*(threadParams->selFile),objId[threadParams->imageIndex]);
                         // threadParams->ctf.Tm=threadParams->parent->Ts;
                         threadParams->ctf.produceSideInfo();
                     }
 
                     threadParams->weight = 1.;
 
                     if(parent->do_weights)
                         threadParams->weight = weight;
                     else if (!parent->do_weights)
                     {
                         weight=1.0;
                     }
                     else if (weight==0.0)
                     {
                         threadParams->read = 2;
                         break;
                     }
 
                     // Copy the projection to the center of the padded image
                     // and compute its Fourier transform
                     proj().setXmippOrigin();
                     auto localPaddedImgSize=(size_t)(parent->imgSize*parent->padding_factor_proj);
                     if (threadParams->reprocessFlag)
                         localPaddedFourier.initZeros(localPaddedImgSize,localPaddedImgSize/2+1);
                     else
                     {
                         localPaddedImg.initZeros(localPaddedImgSize,localPaddedImgSize);
                         localPaddedImg.setXmippOrigin();
                         const MultidimArray<double> &mProj=proj();
                         FOR_ALL_ELEMENTS_IN_ARRAY2D(mProj)
                         A2D_ELEM(localPaddedImg,i,j)=A2D_ELEM(mProj,i,j);
                         // COSS A2D_ELEM(localPaddedImg,i,j)=weight*A2D_ELEM(mProj,i,j);
                         CenterFFT(localPaddedImg,true);
 
                         // Fourier transformer for the images
                         localTransformerImg.setReal(localPaddedImg);
                         localTransformerImg.FourierTransform();
                         localTransformerImg.getFourierAlias(localPaddedFourier);
                     }
 
                     // Compute the coordinate axes associated to this image
                     Euler_angles2matrix(rot, tilt, psi, localA);
                     localAinv=localA.transpose();
 
                     threadParams->localweight = weight;
                     threadParams->localAInv = &localAinv;
                     threadParams->localPaddedFourier = &localPaddedFourier;
                     //#define DEBUG22
 #ifdef DEBUG22
 
                     {//CORRECTO
 
                         if(threadParams->myThreadID%1==0)
                         {
                             proj.write((std::string) integerToString(threadParams->myThreadID)  + "_" +\
                                        integerToString(threadParams->imageIndex) + "proj.spi");
 
                             ImageXmipp save44;
                             save44()=localPaddedImg;
                             save44.write((std::string) integerToString(threadParams->myThreadID)  + "_" +\
                                          integerToString(threadParams->imageIndex) + "local_padded_img.spi");
 
                             FourierImage save33;
                             save33()=localPaddedFourier;
                             save33.write((std::string) integerToString(threadParams->myThreadID)  + "_" +\
                                          integerToString(threadParams->imageIndex) + "local_padded_fourier.spi");
                             FourierImage save22;
                             //save22()=*paddedFourier;
                             save22().alias(*(threadParams->localPaddedFourier));
                             save22.write((std::string) integerToString(threadParams->myThreadID)  + "_" +\
                                          integerToString(threadParams->imageIndex) + "_padded_fourier.spi");
                         }
 
                     }
 #endif
                     #undef DEBUG22
 
                     threadParams->read = 1;
                 }
                 break;
             }
         case EXIT_THREAD:
             return nullptr;
         case PROCESS_WEIGHTS:
             {
 
                 // Get a first approximation of the reconstruction
                 double corr2D_3D=pow(parent->padding_factor_proj,2.)/
                                  (parent->imgSize* pow(parent->padding_factor_vol,3.));
                 // Divide by Zdim because of the
                 // the extra dimension added
                 // and padding differences
                 MultidimArray<double> &mFourierWeights=parent->FourierWeights;
                 for (int k=threadParams->myThreadID; k<=FINISHINGZ(mFourierWeights); k+=parent->numThreads)
                     for (int i=STARTINGY(mFourierWeights); i<=FINISHINGY(mFourierWeights); i++)
                         for (int j=STARTINGX(mFourierWeights); j<=FINISHINGX(mFourierWeights); j++)
                         {
                             if (parent->NiterWeight==0)
                                 A3D_ELEM(parent->VoutFourier,k,i,j)*=corr2D_3D;
                             else
                             {
                                 double weight_kij=A3D_ELEM(mFourierWeights,k,i,j);
                                 if (1.0/weight_kij>ACCURACY)
                                     A3D_ELEM(parent->VoutFourier,k,i,j)*=corr2D_3D*A3D_ELEM(mFourierWeights,k,i,j);
                                 else
                                     A3D_ELEM(parent->VoutFourier,k,i,j)=0;
                             }
                         }
                 break;
             }
         case PROCESS_IMAGE:
             {
                 MultidimArray< std::complex<double> > *paddedFourier = threadParams->paddedFourier;
                 if (threadParams->weight==0.0)
                     break;
                 bool reprocessFlag = threadParams->reprocessFlag;
                 int * statusArray = parent->statusArray;
 
                 int minAssignedRow;
                 int maxAssignedRow;
                 bool breakCase;
                 bool assigned;
 
                 // Get the inverse of the sampling rate
                 // double iTs=parent->padding_factor_proj/parent->Ts;
                 double iTs=1.0/parent->Ts; // The padding factor is not considered here, but later when the indexes
                 //                         // are converted to digital frequencies
                 do
                 {
                     minAssignedRow = -1;
                     maxAssignedRow = -1;
                     breakCase = false;
                     assigned = false;
 
                     do
                     {
                         pthread_mutex_lock( &(parent->workLoadMutex) );
 
                         if ( parent->rowsProcessed == YSIZE(*paddedFourier) )
                         {
                             pthread_mutex_unlock( &(parent->workLoadMutex) );
                             breakCase = true;
                             break;
                         }
 
                         for (size_t w = 0 ; w < YSIZE(*paddedFourier) ; w++ )
                         {
                             if ( statusArray[w]==0 )
                             {
                                 assigned = true;
                                 minAssignedRow = w;
                                 maxAssignedRow = w+minSeparation-1;
 
                                 if ( maxAssignedRow > (int)(YSIZE(*paddedFourier)-1) )
                                     maxAssignedRow = YSIZE(*paddedFourier)-1;
 
                                 for ( int in = (minAssignedRow - minSeparation) ; in < (int)minAssignedRow ; in ++ )
                                 {
                                     if ( ( in >= 0 ) && ( in < (int)YSIZE(*paddedFourier) ))
                                     {
                                         if ( statusArray[in] > -1 )
                                             statusArray[in]++;
                                     }
                                 }
 
                                 for ( int in = minAssignedRow ; in <= (int)maxAssignedRow ; in ++ )
                                 {
                                     if ( ( in >= 0 ) && ( in < (int)YSIZE(*paddedFourier) ))
                                     {
                                         if ( statusArray[in] == 0 )
                                         {
                                             statusArray[in] = -1;
                                             parent->rowsProcessed++;
                                         }
                                     }
                                 }
 
                                 for ( int in = maxAssignedRow+1 ; in <= (maxAssignedRow+minSeparation) ; in ++ )
                                 {
                                     if ( ( in >= 0 ) && ( in < (int)YSIZE(*paddedFourier) ))
                                     {
                                         if ( statusArray[in] > -1 )
                                             statusArray[in]++;
                                     }
                                 }
 
                                 break;
                             }
                         }
 
                         pthread_mutex_unlock( &(parent->workLoadMutex) );
                     }
                     while ( !assigned );
 
                     if ( breakCase == true )
                     {
                         break;
                     }
 
                     Matrix2D<double> * A_SL = threadParams->symmetry;
 
                     // Loop over all Fourier coefficients in the padded image
                     Matrix1D<double> freq(3), gcurrent(3), real_position(3), contFreq(3);
                     Matrix1D<int> corner1(3), corner2(3);
 
                     // Some alias and calculations moved from heavy loops
                     double wCTF=1, wModulator=1.0;
                     double blobRadiusSquared = parent->blob.radius * parent->blob.radius;
                     double iDeltaSqrt = parent->iDeltaSqrt;
                     Matrix1D<double> & blobTableSqrt = parent->blobTableSqrt;
                     int xsize_1 = XSIZE(parent->VoutFourier) - 1;
                     int zsize_1 = ZSIZE(parent->VoutFourier) - 1;
                     MultidimArray< std::complex<double> > &VoutFourier=parent->VoutFourier;
                     MultidimArray<double> &fourierWeights = parent->FourierWeights;
                     //                    MultidimArray<double> &prefourierWeights = parent->preFourierWeights;
                     // Get i value for the thread
                     for (int i = minAssignedRow; i <= maxAssignedRow ; i ++ )
                     {
                         // Discarded rows can be between minAssignedRow and maxAssignedRow, check
                         if ( statusArray[i] == -1 )
                             for (int j=STARTINGX(*paddedFourier); j<=FINISHINGX(*paddedFourier); j++)
                             {
                                 // Compute the frequency of this coefficient in the
                                 // universal coordinate system
                                 FFT_IDX2DIGFREQ(j,XSIZE(parent->paddedImg),XX(freq));
                                 FFT_IDX2DIGFREQ(i,YSIZE(parent->paddedImg),YY(freq));
                                 ZZ(freq)=0;
                                 if (XX(freq)*XX(freq)+YY(freq)*YY(freq)>parent->maxResolution2)
                                     continue;
                                 wModulator=1.0;
                                 if (hasCTF && !reprocessFlag)
                                 {
                                     XX(contFreq)=XX(freq)*iTs;
                                     YY(contFreq)=YY(freq)*iTs;
                                     threadParams->ctf.precomputeValues(XX(contFreq),YY(contFreq));
                                     //wCTF=threadParams->ctf.getValueAt();
                                     wCTF=threadParams->ctf.getValuePureNoKAt();
                                     //wCTF=threadParams->ctf.getValuePureWithoutDampingAt();
 
                                     if (std::isnan(wCTF))
                                     {
                                         if (i==0 && j==0)
                                             wModulator=wCTF=1.0;
                                         else
                                             wModulator=wCTF=0.0;
                                     }
                                     if (fabs(wCTF)<parent->minCTF)
                                     {
                                         wModulator=fabs(wCTF);
                                         wCTF=SGN(wCTF);
                                     }
                                     else
                                         wCTF=1.0/wCTF;
                                     if (parent->phaseFlipped)
                                         wCTF=fabs(wCTF);
                                 }
 
                                 SPEED_UP_temps012;
                                 M3x3_BY_V3x1(freq,*A_SL,freq);
 
                                 // Look for the corresponding index in the volume Fourier transform
                                 DIGFREQ2FFT_IDX_DOUBLE(XX(freq),parent->volPadSizeX,XX(real_position));
                                 DIGFREQ2FFT_IDX_DOUBLE(YY(freq),parent->volPadSizeY,YY(real_position));
                                 DIGFREQ2FFT_IDX_DOUBLE(ZZ(freq),parent->volPadSizeZ,ZZ(real_position));
 
                                 // Put a box around that coefficient
                                 XX(corner1)=CEIL (XX(real_position)-parent->blob.radius);
                                 YY(corner1)=CEIL (YY(real_position)-parent->blob.radius);
                                 ZZ(corner1)=CEIL (ZZ(real_position)-parent->blob.radius);
                                 XX(corner2)=FLOOR(XX(real_position)+parent->blob.radius);
                                 YY(corner2)=FLOOR(YY(real_position)+parent->blob.radius);
                                 ZZ(corner2)=FLOOR(ZZ(real_position)+parent->blob.radius);
 
 #ifdef DEBUG
 
                                 std::cout << "Idx Img=(0," << i << "," << j << ") -> Freq Img=("
                                 << freq.transpose() << ") ->\n    Idx Vol=("
                                 << real_position.transpose() << ")\n"
                                 << "   Corner1=" << corner1.transpose() << std::endl
                                 << "   Corner2=" << corner2.transpose() << std::endl;
 #endif
                                 // Loop within the box
                                 auto *ptrIn=(double *)&(A2D_ELEM(*paddedFourier, i,j));
 
                                 // Some precalculations
                                 for (int intz = ZZ(corner1); intz <= ZZ(corner2); ++intz)
                                 {
                                     double z = intz - ZZ(real_position);
                                     A1D_ELEM(z2precalculated,intz)=z*z;
                                     if (A1D_ELEM(zWrapped,intz)<0)
                                     {
                                         int iz, izneg;
                                         fastIntWRAP(iz, intz, 0, zsize_1);
                                         A1D_ELEM(zWrapped,intz)=iz;
                                         int miz=-iz;
                                         fastIntWRAP(izneg, miz,0,zsize_1);
                                         A1D_ELEM(zNegWrapped,intz)=izneg;
                                     }
                                 }
                                 for (int inty = YY(corner1); inty <= YY(corner2); ++inty)
                                 {
                                     double y = inty - YY(real_position);
                                     A1D_ELEM(y2precalculated,inty)=y*y;
                                     if (A1D_ELEM(yWrapped,inty)<0)
                                     {
                                         int iy, iyneg;
                                         fastIntWRAP(iy, inty, 0, zsize_1);
                                         A1D_ELEM(yWrapped,inty)=iy;
                                         int miy=-iy;
                                         fastIntWRAP(iyneg, miy,0,zsize_1);
                                         A1D_ELEM(yNegWrapped,inty)=iyneg;
                                     }
                                 }
                                 for (int intx = XX(corner1); intx <= XX(corner2); ++intx)
                                 {
                                     double x = intx - XX(real_position);
                                     A1D_ELEM(x2precalculated,intx)=x*x;
                                     if (A1D_ELEM(xWrapped,intx)<0)
                                     {
                                         int ix, ixneg;
                                         fastIntWRAP(ix, intx, 0, zsize_1);
                                         A1D_ELEM(xWrapped,intx)=ix;
                                         int mix=-ix;
                                         fastIntWRAP(ixneg, mix,0,zsize_1);
                                         A1D_ELEM(xNegWrapped,intx)=ixneg;
                                     }
                                 }
 
                                 // Actually compute
                                 for (int intz = ZZ(corner1); intz <= ZZ(corner2); ++intz)
                                 {
                                     double z2 = A1D_ELEM(z2precalculated,intz);
                                     int iz=A1D_ELEM(zWrapped,intz);
                                     int izneg=A1D_ELEM(zNegWrapped,intz);
 
                                     for (int inty = YY(corner1); inty <= YY(corner2); ++inty)
                                     {
                                         double y2z2 = A1D_ELEM(y2precalculated,inty) + z2;
                                         if (y2z2 > blobRadiusSquared)
                                             continue;
                                         int iy=A1D_ELEM(yWrapped,inty);
                                         int iyneg=A1D_ELEM(yNegWrapped,inty);
 
                                         int size1=YXSIZE(VoutFourier)*izneg+(iyneg*XSIZE(VoutFourier));
                                         int size2=YXSIZE(VoutFourier)*iz+(iy*XSIZE(VoutFourier));
                                         int fixSize=0;
 
                                         for (int intx = XX(corner1); intx <= XX(corner2); ++intx)
                                         {
                                             // Compute distance to the center of the blob
                                             // Compute blob value at that distance
                                             double d2 = A1D_ELEM(x2precalculated,intx) + y2z2;
 
                                             if (d2 > blobRadiusSquared)
                                                 continue;
                                             auto aux = (int)(d2 * iDeltaSqrt + 0.5);//Same as ROUND but avoid comparison
                                             double w = VEC_ELEM(blobTableSqrt, aux)*threadParams->weight *wModulator;
 
                                             // Look for the location of this logical index
                                             // in the physical layout
 #ifdef DEBUG
 
                                             std::cout << "   gcurrent=" << gcurrent.transpose()
                                             << " d=" << d << std::endl;
                                             std::cout << "   1: intx=" << intx
                                             << " inty=" << inty
                                             << " intz=" << intz << std::endl;
 #endif
 
                                             int ix=A1D_ELEM(xWrapped,intx);
 #ifdef DEBUG
 
                                             std::cout << "   2: ix=" << ix << " iy=" << iy
                                             << " iz=" << iz << std::endl;
 #endif
 
                                             bool conjugate=false;
                                             int izp, iyp, ixp;
                                             if (ix > xsize_1)
                                             {
                                                 izp = izneg;
                                                 iyp = iyneg;
                                                 ixp = A1D_ELEM(xNegWrapped,intx);
                                                 conjugate=true;
                                                 fixSize = size1;
                                             }
                                             else
                                             {
                                                 izp=iz;
                                                 iyp=iy;
                                                 ixp=ix;
                                                 fixSize = size2;
                                             }
 #ifdef DEBUG
                                             std::cout << "   3: ix=" << ix << " iy=" << iy
                                             << " iz=" << iz << " conj="
                                             << conjugate << std::endl;
 #endif
 
                                             // Add the weighted coefficient
                                             if (reprocessFlag)
                                             {
                                                 // Use VoutFourier as temporary to save the memory
                                                 auto *ptrOut=(double *)&(DIRECT_A3D_ELEM(VoutFourier, izp,iyp,ixp));
                                                 DIRECT_A3D_ELEM(fourierWeights, izp,iyp,ixp) += (w * ptrOut[0]);
                                             }
                                             else
                                             {
                                                 double wEffective=w*wCTF;
                                                 size_t memIdx=fixSize + ixp;//YXSIZE(VoutFourier)*(izp)+((iyp)*XSIZE(VoutFourier))+(ixp);
                                                 auto *ptrOut=(double *)&(DIRECT_A1D_ELEM(VoutFourier, memIdx));
                                                 ptrOut[0] += wEffective * ptrIn[0];
                                                 DIRECT_A1D_ELEM(fourierWeights, memIdx) += w;
 
                                                 if (conjugate)
                                                     ptrOut[1]-=wEffective*ptrIn[1];
                                                 else
                                                     ptrOut[1]+=wEffective*ptrIn[1];
                                             }
                                         }
                                     }
                                 }
                             }
                     }
 
                     pthread_mutex_lock( &(parent->workLoadMutex) );
 
                     for ( int w = (minAssignedRow - minSeparation) ; w < minAssignedRow ; w ++ )
                     {
                         if ( ( w >= 0 ) && ( w < (int)YSIZE(*paddedFourier) ))
                         {
                             if ( statusArray[w] > 0 )
                             {
                                 statusArray[w]--;
                             }
                         }
                     }
 
                     for ( int w = maxAssignedRow+1 ; w <= (maxAssignedRow+minSeparation) ; w ++ )
                     {
                         if ( ( w >= 0 ) && ( w < (int)YSIZE(*paddedFourier) ))
                         {
                             if ( statusArray[w] > 0 )
                             {
                                 statusArray[w]--;
                             }
                         }
                     }
 
                     pthread_mutex_unlock( &(parent->workLoadMutex) );
 
                 }
                 while (!breakCase);
                 break;
             }
         default:
             break;
         }
 
         barrier_wait( barrier );
     }
     while ( 1 );
 }

◆ produceSideinfo()

void ProgRecFourier::produceSideinfo ( )

Produce side info: fill arrays with relevant transformation matrices.

Definition at line 184 of file reconstruct_fourier.cpp.

 {
     // Translate the maximum resolution to digital frequency
     // maxResolution=sampling_rate/maxResolution;
     maxResolution2=maxResolution*maxResolution;
 
     // Read the input images
     SF.read(fn_sel);
     SF.removeDisabled();
 
     // Ask for memory for the output volume and its Fourier transform
     size_t objId = SF.firstRowId();
     FileName fnImg;
     SF.getValue(MDL_IMAGE,fnImg,objId);
     Image<double> I;
     I.read(fnImg, HEADER);
     int Ydim=YSIZE(I());
     int Xdim=XSIZE(I());
     if (Ydim!=Xdim)
         REPORT_ERROR(ERR_MULTIDIM_SIZE,"This algorithm only works for squared images");
     imgSize=Xdim;
     volPadSizeX = volPadSizeY = volPadSizeZ=(int)(Xdim*padding_factor_vol);
     Vout().initZeros(volPadSizeZ,volPadSizeY,volPadSizeX);
 
     //use threads for volume inverse fourier transform, plan is created in setReal()
     transformerVol.setThreadsNumber(numThreads);
     transformerVol.setReal(Vout());
 
     Vout().clear(); // Free the memory so that it is available for FourierWeights
     transformerVol.getFourierAlias(VoutFourier);
     VoutFourier.initZeros();
     FourierWeights.initZeros(VoutFourier);
 
     // Ask for memory for the padded images
     auto paddedImgSize=(size_t)(Xdim*padding_factor_proj);
     paddedImg.resize(paddedImgSize,paddedImgSize);
     paddedImg.setXmippOrigin();
     transformerImg.setReal(paddedImg);
 
     // Build a table of blob values
     blobTableSqrt.resize(BLOB_TABLE_SIZE_SQRT);
     fourierBlobTableSqrt.resize(BLOB_TABLE_SIZE_SQRT);
     Fourier_blob_table.resize(BLOB_TABLE_SIZE_SQRT);
 
     struct blobtype blobFourier,blobnormalized;
     blobFourier=blob;
     //Sjors 18aug10 blobFourier.radius/=(padding_factor_proj*Xdim);
     blobFourier.radius/=(padding_factor_vol*Xdim);
     blobnormalized=blob;
     blobnormalized.radius/=((double)padding_factor_proj/padding_factor_vol);
     double deltaSqrt     = (blob.radius*blob.radius) /(BLOB_TABLE_SIZE_SQRT-1);
     double deltaFourier  = (sqrt(3.)*Xdim/2.)/(BLOB_TABLE_SIZE_SQRT-1);
 
     // The interpolation kernel must integrate to 1
     double iw0 = 1.0 / blob_Fourier_val(0.0, blobnormalized);
     //Sjors 18aug10 double padXdim3 = padding_factor_proj * Xdim;
     double padXdim3 = padding_factor_vol * Xdim;
     padXdim3 = padXdim3 * padXdim3 * padXdim3;
     double blobTableSize = blob.radius*sqrt(1./ (BLOB_TABLE_SIZE_SQRT-1));
     //***
     //Following commented line seems to be the right thing but I do not understand it
     //double fourierBlobTableSize = (sqrt(3.)*Xdim*Xdim/2.)*blobFourier.radius *sqrt(1./ (BLOB_TABLE_SIZE_SQRT-1));
     FOR_ALL_ELEMENTS_IN_MATRIX1D(blobTableSqrt)
 
     {
         //use a r*r sample instead of r
         //DIRECT_VEC_ELEM(blob_table,i)         = blob_val(delta*i, blob)  *iw0;
         VEC_ELEM(blobTableSqrt,i)    = blob_val(blobTableSize*sqrt((double)i), blob)  *iw0;
         //***
         //DIRECT_VEC_ELEM(fourierBlobTableSqrt,i) =
         //     blob_Fourier_val(fourierBlobTableSize*sqrt(i), blobFourier)*padXdim3  *iw0;
         VEC_ELEM(Fourier_blob_table,i) =
             blob_Fourier_val(deltaFourier*i, blobFourier)*padXdim3  *iw0;
         //#define DEBUG
 #ifdef DEBUG
 
         std::cout << VEC_ELEM(Fourier_blob_table,i)
         << " " << VEC_ELEM(fourierBlobTableSqrt,i)
         << std::endl;
 #endif
   #undef DEBUG
 
     }
     //iDelta        = 1/delta;
     iDeltaSqrt    = 1/deltaSqrt;
     iDeltaFourier = 1/deltaFourier;
 
     // Get symmetries
     Matrix2D<double>  Identity(3,3);
     Identity.initIdentity();
     R_repository.push_back(Identity);
     if (fn_sym != "")
     {
         SymList SL;
         SL.readSymmetryFile(fn_sym);
         for (int isym = 0; isym < SL.trueSymsNo(); isym++)
         {
             Matrix2D<double>  L(4, 4), R(4, 4);
             SL.getMatrices(isym, L, R);
             R.resize(3, 3);
             R_repository.push_back(R);
         }
     }
 }

◆ readParams()

void ProgRecFourier::readParams ( )

virtual

Read arguments from command line.

Reimplemented from XmippProgram.

Definition at line 64 of file reconstruct_fourier.cpp.

 {
     fn_sel = getParam("-i");
     fn_out = getParam("-o");
     fn_sym = getParam("--sym");
     if(checkParam("--prepare_fsc"))
         fn_fsc = getParam("--prepare_fsc");
     do_weights = checkParam("--weight");
     padding_factor_proj = getDoubleParam("--padding", 0);
     padding_factor_vol = getDoubleParam("--padding", 1);
     blob.radius   = getDoubleParam("--blob", 0);
     blob.order    = getIntParam("--blob", 1);
     blob.alpha    = getDoubleParam("--blob", 2);
     maxResolution = getDoubleParam("--max_resolution");
     numThreads = getIntParam("--thr");
     thrWidth = getIntParam("--thr", 1);
     NiterWeight = getIntParam("--iter");
     useCTF = checkParam("--useCTF");
     phaseFlipped = checkParam("--phaseFlipped");
     minCTF = getDoubleParam("--minCTF");
     if (useCTF)
         Ts=getDoubleParam("--sampling");
 }

◆ run()

void ProgRecFourier::run ( )

virtual

Run.

Reimplemented from XmippProgram.

Definition at line 124 of file reconstruct_fourier.cpp.

 {
     show();
     produceSideinfo();
     // Process all images in the selfile
     if (verbose)
     {
         if (NiterWeight!=0)
             init_progress_bar(NiterWeight*SF.size());
         else
             init_progress_bar(SF.size());
     }
     // Create threads stuff
     barrier_init( &barrier, numThreads+1 );
     pthread_mutex_init( &workLoadMutex, nullptr);
     statusArray = nullptr;
     th_ids = (pthread_t *)malloc( numThreads * sizeof( pthread_t));
     th_args = (ImageThreadParams *) malloc ( numThreads * sizeof( ImageThreadParams ) );
 
     // Create threads
     for ( int nt = 0 ; nt < numThreads ; nt ++ )
     {
         // Passing parameters to each thread
         th_args[nt].parent = this;
         th_args[nt].myThreadID = nt;
         th_args[nt].selFile = new MetaDataVec(SF);
         pthread_create( (th_ids+nt) , nullptr, processImageThread, (void *)(th_args+nt) );
     }
 
     //Computing interpolated volume
     processImages(0, SF.size() - 1, !fn_fsc.empty(), false);
 
     // Correcting the weights
     correctWeight();
 
     //Saving the volume
     finishComputations(fn_out);
 
     threadOpCode = EXIT_THREAD;
 
     // Waiting for threads to finish
     barrier_wait( &barrier );
     for ( int nt = 0 ; nt < numThreads ; nt ++ )
         pthread_join(*(th_ids+nt), nullptr);
     barrier_destroy( &barrier );
 
     // Deallocate resources.
     if ( statusArray != nullptr)
     {
         free(statusArray);
     }
     for (int nt=0; nt<numThreads ;nt++)
     {
         delete(th_args[nt].selFile);
     }
     free(th_ids);
     free(th_args);
 }

◆ setIO()

void ProgRecFourier::setIO	(	const FileName &	fn_in,
		const FileName &	fn_out
	)

virtual

Functions of common reconstruction interface.

Implements ProgReconsBase.

Definition at line 1182 of file reconstruct_fourier.cpp.

 {
     this->fn_sel = fn_in;
     this->fn_out = fn_out;
 }

◆ show()

void ProgRecFourier::show ( )

Show.

Definition at line 89 of file reconstruct_fourier.cpp.

 {
     if (verbose > 0)
     {
         std::cout << " =====================================================================" << std::endl;
         std::cout << " Direct 3D reconstruction method using Kaiser windows as interpolators" << std::endl;
         std::cout << " =====================================================================" << std::endl;
         std::cout << " Input selfile             : "  << fn_sel << std::endl;
         std::cout << " padding_factor_proj       : "  << padding_factor_proj << std::endl;
         std::cout << " padding_factor_vol        : "  << padding_factor_vol << std::endl;
         std::cout << " Output volume             : "  << fn_out << std::endl;
         if (fn_sym != "")
             std::cout << " Symmetry file for projections : "  << fn_sym << std::endl;
         if (fn_fsc != "")
             std::cout << " File root for FSC files: " << fn_fsc << std::endl;
         if (do_weights)
             std::cout << " Use weights stored in the image headers or doc file" << std::endl;
         else
             std::cout << " Do NOT use weights" << std::endl;
         if (useCTF)
             std::cout << "Using CTF information" << std::endl
             << "Sampling rate: " << Ts << std::endl
             << "Phase flipped: " << phaseFlipped << std::endl
             << "Minimum CTF: " << minCTF << std::endl;
         std::cout << "\n Interpolation Function"
         << "\n   blrad                 : "  << blob.radius
         << "\n   blord                 : "  << blob.order
         << "\n   blalpha               : "  << blob.alpha
         //<< "\n sampling_rate           : "  << sampling_rate
         << "\n max_resolution          : "  << maxResolution
         << "\n -----------------------------------------------------------------" << std::endl;
     }
 }

Member Data Documentation

◆ barrier

barrier_t ProgRecFourier::barrier

To create a barrier synchronization for threads.

Definition at line 139 of file reconstruct_fourier.h.

◆ blob

struct blobtype ProgRecFourier::blob

Definition at line 170 of file reconstruct_fourier.h.

◆ blobTableSqrt

Matrix1D<double> ProgRecFourier::blobTableSqrt

Definition at line 160 of file reconstruct_fourier.h.

◆ do_weights

bool ProgRecFourier::do_weights

Flag whether to use the weights in the image metadata

Definition at line 86 of file reconstruct_fourier.h.

◆ fn_doc

FileName ProgRecFourier::fn_doc

Definition at line 80 of file reconstruct_fourier.h.

◆ fn_fsc

FileName ProgRecFourier::fn_fsc

Definition at line 80 of file reconstruct_fourier.h.

◆ fn_out

FileName ProgRecFourier::fn_out

Filenames

Definition at line 80 of file reconstruct_fourier.h.

◆ fn_sel

FileName ProgRecFourier::fn_sel

Definition at line 80 of file reconstruct_fourier.h.

◆ fn_sym

FileName ProgRecFourier::fn_sym

Definition at line 80 of file reconstruct_fourier.h.

◆ Fourier_blob_table

Matrix1D<double> ProgRecFourier::Fourier_blob_table

Definition at line 157 of file reconstruct_fourier.h.

◆ fourierBlobTableSqrt

Matrix1D<double> ProgRecFourier::fourierBlobTableSqrt

Definition at line 160 of file reconstruct_fourier.h.

◆ FourierWeights

MultidimArray<double> ProgRecFourier::FourierWeights

Definition at line 185 of file reconstruct_fourier.h.

◆ iDeltaFourier

double ProgRecFourier::iDeltaFourier

Definition at line 164 of file reconstruct_fourier.h.

◆ iDeltaSqrt

double ProgRecFourier::iDeltaSqrt

Definition at line 164 of file reconstruct_fourier.h.

◆ imgSize

int ProgRecFourier::imgSize

Definition at line 149 of file reconstruct_fourier.h.

◆ maxResolution

double ProgRecFourier::maxResolution

Max resolution in Angstroms.

Definition at line 112 of file reconstruct_fourier.h.

◆ maxResolution2

double ProgRecFourier::maxResolution2

Definition at line 167 of file reconstruct_fourier.h.

◆ minCTF

double ProgRecFourier::minCTF

Minimum CTF value to invert

Definition at line 97 of file reconstruct_fourier.h.

◆ NiterWeight

int ProgRecFourier::NiterWeight

Number of iterations for the weight.

Definition at line 115 of file reconstruct_fourier.h.

◆ numThreads

int ProgRecFourier::numThreads

Number of threads to use in parallel to process a single image.

Definition at line 118 of file reconstruct_fourier.h.

◆ paddedImg

MultidimArray<double> ProgRecFourier::paddedImg

Definition at line 188 of file reconstruct_fourier.h.

◆ padding_factor_proj

double ProgRecFourier::padding_factor_proj

Projection padding Factor

Definition at line 103 of file reconstruct_fourier.h.

◆ padding_factor_vol

double ProgRecFourier::padding_factor_vol

Volume padding Factor

Definition at line 106 of file reconstruct_fourier.h.

◆ phaseFlipped

bool ProgRecFourier::phaseFlipped

Phase flipped. True if the images have been phase flipped before entering.

Definition at line 94 of file reconstruct_fourier.h.

◆ R_repository

std::vector<Matrix2D<double> > ProgRecFourier::R_repository

Definition at line 173 of file reconstruct_fourier.h.

◆ rowsProcessed

size_t ProgRecFourier::rowsProcessed

Number of rows already processed on an image.

Definition at line 130 of file reconstruct_fourier.h.

◆ sampling_rate

double ProgRecFourier::sampling_rate

Sampling rate in Angstroms/pixel

Definition at line 109 of file reconstruct_fourier.h.

◆ SF

MetaDataVec ProgRecFourier::SF

SelFile containing all projections

Definition at line 83 of file reconstruct_fourier.h.

◆ statusArray

int* ProgRecFourier::statusArray

A status array for each row in an image (processing, processed,etc..)

Definition at line 142 of file reconstruct_fourier.h.

◆ th_args

ImageThreadParams* ProgRecFourier::th_args

Contains parameters passed to each thread.

Definition at line 124 of file reconstruct_fourier.h.

◆ th_ids

pthread_t* ProgRecFourier::th_ids

IDs for the threads.

Definition at line 121 of file reconstruct_fourier.h.

◆ threadOpCode

int ProgRecFourier::threadOpCode

Tells the threads what to do next.

Definition at line 127 of file reconstruct_fourier.h.

◆ thrWidth

int ProgRecFourier::thrWidth

How many image rows are processed at a time by a single thread.

Definition at line 145 of file reconstruct_fourier.h.

◆ transformerImg

FourierTransformer ProgRecFourier::transformerImg

Definition at line 179 of file reconstruct_fourier.h.

◆ transformerVol

FourierTransformer ProgRecFourier::transformerVol

Definition at line 176 of file reconstruct_fourier.h.

◆ Ts

double ProgRecFourier::Ts

Sampling rate

Definition at line 100 of file reconstruct_fourier.h.

◆ useCTF

bool ProgRecFourier::useCTF

Use CTF

Definition at line 89 of file reconstruct_fourier.h.

◆ volPadSizeX

int ProgRecFourier::volPadSizeX

Definition at line 152 of file reconstruct_fourier.h.

◆ volPadSizeY

int ProgRecFourier::volPadSizeY

Definition at line 153 of file reconstruct_fourier.h.

◆ volPadSizeZ

int ProgRecFourier::volPadSizeZ

Definition at line 154 of file reconstruct_fourier.h.

◆ Vout

Image<double> ProgRecFourier::Vout

Definition at line 191 of file reconstruct_fourier.h.

◆ VoutFourier

MultidimArray< std::complex<double> > ProgRecFourier::VoutFourier

Definition at line 182 of file reconstruct_fourier.h.

◆ VoutFourierTmp

MultidimArray< std::complex<double> > ProgRecFourier::VoutFourierTmp

Definition at line 182 of file reconstruct_fourier.h.

◆ workLoadMutex

pthread_mutex_t ProgRecFourier::workLoadMutex

Controls mutual exclusion on critical zones of code.

Definition at line 136 of file reconstruct_fourier.h.

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

xmipp/libraries/reconstruction/reconstruct_fourier.h
xmipp/libraries/reconstruction/reconstruct_fourier.cpp

Public Member Functions

Static Public Member Functions

Public Attributes

Additional Inherited Members

Detailed Description

Member Function Documentation

◆ correctWeight()

◆ defineParams()

◆ finishComputations()

◆ forceWeightSymmetry()

◆ processImages()

◆ processImageThread()

◆ produceSideinfo()

◆ readParams()

◆ run()

◆ setIO()

◆ show()

Member Data Documentation

◆ barrier

◆ blob

◆ blobTableSqrt

◆ do_weights

◆ fn_doc

◆ fn_fsc

◆ fn_out

◆ fn_sel

◆ fn_sym

◆ Fourier_blob_table

◆ fourierBlobTableSqrt

◆ FourierWeights

◆ iDeltaFourier

◆ iDeltaSqrt

◆ imgSize

◆ maxResolution

◆ maxResolution2

◆ minCTF

◆ NiterWeight

◆ numThreads

◆ paddedImg

◆ padding_factor_proj

◆ padding_factor_vol

◆ phaseFlipped

◆ R_repository

◆ rowsProcessed

◆ sampling_rate

◆ SF

◆ statusArray

◆ th_args

◆ th_ids

◆ threadOpCode

◆ thrWidth

◆ transformerImg

◆ transformerVol

◆ Ts

◆ useCTF

◆ volPadSizeX

◆ volPadSizeY

◆ volPadSizeZ

◆ Vout

◆ VoutFourier

◆ VoutFourierTmp

◆ workLoadMutex