ProgRecFourierAccel Class Reference

#include <reconstruct_fourier_accel.h>

Inheritance diagram for ProgRecFourierAccel:

Collaboration diagram for ProgRecFourierAccel:

Public Member Functions
	ProgRecFourierAccel ()
void	run ()
void	finishComputations (const FileName &out_name)
virtual void	setIO (const FileName &fn_in, const FileName &fn_out)
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 ()

Protected Member Functions
void	mirrorAndCropTempSpaces ()
void	forceHermitianSymmetry ()
void	processWeights ()
void	createLoadingThread ()
void	cleanLoadingThread ()
void	readParams ()
void	defineParams ()
void	show ()
void	produceSideinfo ()
void	processImages (int firstImageIndex, int lastImageIndex)
void	releaseTempSpaces ()
Protected Member Functions inherited from XmippProgram
void	defineCommons ()

Protected Attributes
LoadThreadParams	loadThread
MetaDataVec	SF
FileName	fn_out
FileName	fn_in
int	maxVolumeIndexYZ
int	maxVolumeIndexX
std::complex< float > ***	tempVolume
float ***	tempWeights
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

Additional Inherited Members
Public Attributes inherited from XmippProgram
bool	doRun
bool	runWithoutArgs
int	verbose
	Verbosity level. More...
int	debug

Detailed Description

Definition at line 99 of file reconstruct_fourier_accel.h.

Constructor & Destructor Documentation

ProgRecFourierAccel()

ProgRecFourierAccel::ProgRecFourierAccel ( )

inline

Definition at line 102 of file reconstruct_fourier_accel.h.

: tempVolume(NULL), tempWeights(NULL) {};

Member Function Documentation

cleanLoadingThread()

void ProgRecFourierAccel::cleanLoadingThread ( )

protected

Method will release all resources allocated for loading images

Definition at line 167 of file reconstruct_fourier_accel.cpp.

                                             {
    threadOpCode = EXIT_THREAD;
    // Waiting for thread to finish
    barrier_wait( &barrier );
    pthread_join(*(&loadThread.id), nullptr);
    barrier_destroy( &barrier );
}

createLoadingThread()

void ProgRecFourierAccel::createLoadingThread ( )

protected

Method will create thread used for loading images and set all necessary values/variables

Definition at line 158 of file reconstruct_fourier_accel.cpp.

                                              {
    barrier_init( &barrier, 2 ); // two barries - for main and loading thread
    loadThread.buffer1 = loadThread.buffer2 = nullptr;
    loadThread.parent = this;
    loadThread.selFile = &SF;
    pthread_create( &loadThread.id , nullptr, loadImageThread, (void *)(&loadThread) );
    threadOpCode = PRELOAD_IMAGE;
}

defineParams()

void ProgRecFourierAccel::defineParams ( )

protectedvirtual

Specify supported command line arguments

Reimplemented from XmippProgram.

Definition at line 55 of file reconstruct_fourier_accel.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("  [--sym <symfile=c1>]              : Enforce symmetry in projections");
    addParamsLine("  [--padding <proj=2.0> <vol=2.0>]  : Padding used for projections and volume");
    addParamsLine("  [--max_resolution <p=0.5>]     : Max resolution (Nyquist=0.5)");
    addParamsLine("  [--weight]                     : Use weights stored in the image metadata");
    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("  [--fast]                       : Do the blobing at the end of the computation.");
    addParamsLine("                                 : Gives slightly different results, but is faster.");
    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");
    addParamsLine("  [--bufferSize <size=25>]        : Number of projection loaded in memory (will be actually 2x as much.");
    addParamsLine("                                 : This will require up to 4*size*projSize*projSize*16B, e.g.");
    addParamsLine("                                 : 100MB for projection of 256x256 or 400MB for projection of 512x512");
    addExampleLine("For reconstruct enforcing i3 symmetry and using stored weights:", false);
    addExampleLine("   xmipp_reconstruct_fourier_accel  -i reconstruction.sel --sym i3 --weight");
}

finishComputations()

void ProgRecFourierAccel::finishComputations

(

const FileName &

out_name

)

Method will take data stored in tempVolume and tempWeights (which should be cropped in X axis), calculate IFFT and store result to final destination.

Definition at line 1002 of file reconstruct_fourier_accel.cpp.

{
    if (useFast) {
        tempVolume = applyBlob(tempVolume, blob.radius, blobTableSqrt, iDeltaSqrt);
        tempWeights = applyBlob(tempWeights, blob.radius, blobTableSqrt, iDeltaSqrt);
    }

    forceHermitianSymmetry();
    processWeights();
    release(tempWeights, maxVolumeIndexYZ+1, maxVolumeIndexYZ+1);
    MultidimArray< std::complex<double> > VoutFourier;
    allocateVoutFourier(VoutFourier);
    convertToExpectedSpace(tempVolume, maxVolumeIndexYZ, VoutFourier);
    release(tempVolume, maxVolumeIndexYZ+1, maxVolumeIndexYZ+1);

    // Output volume
    Image<double> Vout;
    Vout().initZeros(paddedImgSize,paddedImgSize,paddedImgSize);
    FourierTransformer transformerVol;
    transformerVol.setThreadsNumber(2); // use just main and 'loading' thread
    transformerVol.fReal = &(Vout.data);
    transformerVol.setFourierAlias(VoutFourier);
    transformerVol.recomputePlanR2C();

    transformerVol.inverseFourierTransform();
    transformerVol.clear();
    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));
        A3D_ELEM(mVout,k,i,j) /= (ipad_relation*factor2*factor);
        meanFactor2+=factor2;
    }
    meanFactor2/=MULTIDIM_SIZE(mVout);
    FOR_ALL_ELEMENTS_IN_ARRAY3D(mVout)
    A3D_ELEM(mVout,k,i,j) *= meanFactor2;
    Vout.write(out_name);
    Vout.clear();
}

forceHermitianSymmetry()

void ProgRecFourierAccel::forceHermitianSymmetry ( )

protected

Method will enforce Hermitian symmetry, i.e will remove make sure that the values in temporal space at X=0 are complex conjugate of in respect to center of the space

Definition at line 889 of file reconstruct_fourier_accel.cpp.

                                                 {
    int x = 0;
    for (int z = 0; z <= maxVolumeIndexYZ; z++) {
        for (int y = 0; y <= maxVolumeIndexYZ/2; y++) {
            int newPos[3];
            // mirror against center of the volume, e.g. [0,0,0]->[size,size,size]. It will fit as the input space is one voxel biger
            newPos[0] = x;
            newPos[1] = maxVolumeIndexYZ - y;
            newPos[2] = maxVolumeIndexYZ - z;
            std::complex<float> tmp1 = 0.5f * (tempVolume[newPos[2]][newPos[1]][newPos[0]] + conj(tempVolume[z][y][x]));
            float tmp2 = 0.5f * (tempWeights[newPos[2]][newPos[1]][newPos[0]] + tempWeights[z][y][x]);

            tempVolume[newPos[2]][newPos[1]][newPos[0]] = tmp1;
            tempVolume[z][y][x] = conj(tmp1);
            tempWeights[newPos[2]][newPos[1]][newPos[0]] = tempWeights[z][y][x] = tmp2;
        }
    }
}

mirrorAndCropTempSpaces()

void ProgRecFourierAccel::mirrorAndCropTempSpaces ( )

protected

Method will take temp spaces (containing complex conjugate values in the 'right X side'), transfer them to 'left X side' and remove the 'right X side'. As a result, the X dimension of the temp spaces will be half of the original.

Definition at line 853 of file reconstruct_fourier_accel.cpp.

                                                  {
    maxVolumeIndexX = maxVolumeIndexYZ/2; // just half of the space is necessary, the rest is complex conjugate
    mirrorAndCrop(tempWeights, &identity<float>);
    mirrorAndCrop(tempVolume, &conjugate);
}

processImages()

void ProgRecFourierAccel::processImages ( int  firstImageIndex, int  lastImageIndex  )

protected

Method will initialize temporal storage (if necessary), processes images in given interval (including boundary indexes) and store results in temporal storage (

See also

tempVolume and tempWeights)

Definition at line 969 of file reconstruct_fourier_accel.cpp.

{
    int loops = ceil((lastImageIndex-firstImageIndex+1)/(float)bufferSize);

    // the +1 is to prevent outOfBound reading when mirroring the result (later)
    if (nullptr == tempVolume) {
        allocate(tempVolume, maxVolumeIndexYZ+1, maxVolumeIndexYZ+1, maxVolumeIndexYZ+1);
    }
    if (nullptr == tempWeights) {
        allocate(tempWeights, maxVolumeIndexYZ+1, maxVolumeIndexYZ+1, maxVolumeIndexYZ+1);
    }

    int startLoadIndex = firstImageIndex;

    loadImages(startLoadIndex, std::min(lastImageIndex+1, startLoadIndex+bufferSize));
    barrier_wait( &barrier );
    for(int i = 0; i < loops; i++) {
        swapLoadBuffers();
        startLoadIndex += bufferSize;
        loadImages(startLoadIndex, std::min(lastImageIndex+1, startLoadIndex+bufferSize));
        processBuffer(loadThread.buffer2);
        barrier_wait( &barrier );
    }
    delete[] loadThread.buffer1;
    delete[] loadThread.buffer2;
    loadThread.buffer1 = loadThread.buffer2 = nullptr;
}

processWeights()

void ProgRecFourierAccel::processWeights ( )

protected

Method will in effect do the point-wise division of tempVolume and tempWeights (i.e. correct Fourier coefficients by proper weight)

Definition at line 908 of file reconstruct_fourier_accel.cpp.

                                         {
    // Get a first approximation of the reconstruction
    float corr2D_3D=pow(padding_factor_proj,2.)/
                     (imgSize* pow(padding_factor_vol,3.));
    for (int z = 0; z <= maxVolumeIndexYZ; z++) {
        for (int y = 0; y <= maxVolumeIndexYZ; y++) {
            for (int x = 0; x <= maxVolumeIndexX; x++) {
                float weight = tempWeights[z][y][x];

                if (weight > ACCURACY)
                    tempVolume[z][y][x] *= corr2D_3D / weight;
                else
                    tempVolume[z][y][x] = 0;
            }
        }
    }
}

produceSideinfo()

void ProgRecFourierAccel::produceSideinfo ( )

protected

Method will fill other help structures and variables

Definition at line 175 of file reconstruct_fourier_accel.cpp.

{
    // Translate the maximum resolution to digital frequency
    // maxResolution=sampling_rate/maxResolution;
    maxResolutionSqr=maxResolution*maxResolution;

    // Read the input images
    SF.read(fn_in);
    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;
    paddedImgSize = Xdim*padding_factor_vol;
    auto conserveRows = (size_t) ceil((double) paddedImgSize * maxResolution * 2.0);
    conserveRows = (size_t) ceil((double) conserveRows / 2.0);
    maxVolumeIndexX = maxVolumeIndexYZ = 2 * conserveRows;

    // Build a table of blob values
    Fourier_blob_table.resize(BLOB_TABLE_SIZE_SQRT);

    struct blobtype blobFourier,blobnormalized;
    blobFourier=blob;
    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);
    double padXdim3 = padding_factor_vol * Xdim;
    padXdim3 = padXdim3 * padXdim3 * padXdim3;
    double blobTableSize = blob.radius*sqrt(1./ (BLOB_TABLE_SIZE_SQRT-1));
    for (int i = 0; i < BLOB_TABLE_SIZE_SQRT; i++)
    {
        //use a r*r sample instead of r
        //DIRECT_VEC_ELEM(blob_table,i)         = blob_val(delta*i, blob)  *iw0;
        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.symsNo(); 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 ProgRecFourierAccel::readParams ( )

protectedvirtual

Read arguments from command line

Reimplemented from XmippProgram.

Definition at line 85 of file reconstruct_fourier_accel.cpp.

{
    fn_in = getParam("-i");
    fn_out = getParam("-o");
    fn_sym = getParam("--sym");
    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);
    useFast       = checkParam("--fast");
    maxResolution = getDoubleParam("--max_resolution");
    useCTF = checkParam("--useCTF");
    isPhaseFlipped = checkParam("--phaseFlipped");
    minCTF = getDoubleParam("--minCTF");
    if (useCTF)
        iTs = 1 / getDoubleParam("--sampling");
    bufferSize = getIntParam("--bufferSize");
}

releaseTempSpaces()

void ProgRecFourierAccel::releaseTempSpaces ( )

protected

Method will release temporal spaces for weights and Fourier coefs.

Definition at line 997 of file reconstruct_fourier_accel.cpp.

                                            {
    release(tempWeights, maxVolumeIndexYZ+1, maxVolumeIndexYZ+1);
    release(tempVolume, maxVolumeIndexYZ+1, maxVolumeIndexYZ+1);
}

run()

void ProgRecFourierAccel::run ( )

virtual

Run the image processing. Method will load data, process them and store result to final destination.

Reimplemented from XmippProgram.

Definition at line 139 of file reconstruct_fourier_accel.cpp.

{
    show();
    produceSideinfo();
    // Process all images in the selfile
    if (verbose) {
        init_progress_bar(SF.size());
    }
    // Create loading thread stuff
    createLoadingThread();
    //Computing interpolated volume
    processImages(0, SF.size() - 1);
    // remove complex conjugate of the intermediate result
    mirrorAndCropTempSpaces();
    //Saving the volume
    finishComputations(fn_out);
    cleanLoadingThread();
}

setIO()

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

virtual

Functions of common reconstruction interface

Implements ProgReconsBase.

Definition at line 1057 of file reconstruct_fourier_accel.cpp.

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

show()

void ProgRecFourierAccel::show ( )

protected

Show basic info to standard output

Definition at line 107 of file reconstruct_fourier_accel.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_in << 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 (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: " << 1 / iTs << std::endl
            << "Phase flipped: " << isPhaseFlipped << 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 max_resolution          : "  << maxResolution
        << "\n -----------------------------------------------------------------" << std::endl;
    }
}

Member Data Documentation

fn_in

FileName ProgRecFourierAccel::fn_in

protected

Input file name

Definition at line 132 of file reconstruct_fourier_accel.h.

fn_out

FileName ProgRecFourierAccel::fn_out

protected

Output file name

Definition at line 129 of file reconstruct_fourier_accel.h.

loadThread

LoadThreadParams ProgRecFourierAccel::loadThread

protected

Thread used for loading input data

Definition at line 123 of file reconstruct_fourier_accel.h.

maxVolumeIndexX

int ProgRecFourierAccel::maxVolumeIndexX

protected

maximal index in the temporal volumes, X axis

Definition at line 138 of file reconstruct_fourier_accel.h.

maxVolumeIndexYZ

int ProgRecFourierAccel::maxVolumeIndexYZ

protected

maximal index in the temporal volumes, Y and Z axis

Definition at line 135 of file reconstruct_fourier_accel.h.

SF

MetaDataVec ProgRecFourierAccel::SF

protected

SelFile containing all projections

Definition at line 126 of file reconstruct_fourier_accel.h.

tempVolume

std::complex<float>*** ProgRecFourierAccel::tempVolume

protected

3D volume holding the cropped (without high frequencies) Fourier space. Lowest frequencies are in the center, i.e. Fourier space creates a sphere within a cube.

Definition at line 146 of file reconstruct_fourier_accel.h.

tempWeights

float*** ProgRecFourierAccel::tempWeights

protected

3D volume holding the weights of the Fourier coefficients stored in tempVolume.

Definition at line 152 of file reconstruct_fourier_accel.h.

---

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

• xmipp/libraries/reconstruction/reconstruct_fourier_accel.h
• xmipp/libraries/reconstruction/reconstruct_fourier_accel.cpp