Fourier projection

Collaboration diagram for Fourier projection:

Modules
	Projections (2D Image + Euler angles)

Classes
class	Projection

Functions
void	Projection::reset (int Ydim, int Xdim)
void	Projection::setAngles (double _rot, double _tilt, double _psi)
void	Projection::read (const FileName &fn, const bool only_apply_shifts=false, DataMode datamode=DATA, MDRow *row=nullptr)
void	Projection::assign (const Projection &P)
	FourierProjector::FourierProjector (double paddFactor, double maxFreq, int degree)
	FourierProjector::FourierProjector (MultidimArray< double > &V, double paddFactor, double maxFreq, int BSplinedegree)
void	FourierProjector::project (double rot, double tilt, double psi, const MultidimArray< double > *ctf=nullptr)
void	FourierProjector::updateVolume (MultidimArray< double > &V)
void	FourierProjector::produceSideInfo ()
	Prepare the Spline coefficients and projection space. More...
void	FourierProjector::produceSideInfoProjection ()
	Prepare projection space. More...

Variables
Matrix1D< double >	Projection::direction
Matrix2D< double >	Projection::euler
Matrix2D< double >	Projection::eulert
double	FourierProjector::paddingFactor
	Padding factor. More...
double	FourierProjector::maxFrequency
	Maximum Frequency for pixels. More...
double	FourierProjector::BSplineDeg
	The order of B-Spline for interpolation. More...
FourierTransformer	FourierProjector::transformer2D
MultidimArray< double > *	FourierProjector::volume
MultidimArray< double >	FourierProjector::VfourierRealCoefs
MultidimArray< double >	FourierProjector::VfourierImagCoefs
MultidimArray< std::complex< double > >	FourierProjector::projectionFourier
Image< double >	FourierProjector::projection
MultidimArray< double >	FourierProjector::phaseShiftImgB
MultidimArray< double >	FourierProjector::phaseShiftImgA
int	FourierProjector::volumeSize
int	FourierProjector::volumePaddedSize
Matrix2D< double >	FourierProjector::E

Detailed Description

Function Documentation

assign()

void Projection::assign

(

const Projection &

P

)

Another function for assignment.

Definition at line 63 of file fourier_projection.cpp.

{
    *this = P;
}

FourierProjector() [1/2]

FourierProjector::FourierProjector	(	double	paddFactor,
		double	maxFreq,
		int	degree
	)

Definition at line 68 of file fourier_projection.cpp.

{
    paddingFactor = paddFactor;
    maxFrequency = maxFreq;
    BSplineDeg = degree;
    volume = nullptr;
}

FourierProjector() [2/2]

FourierProjector::FourierProjector	(	MultidimArray< double > &	V,
		double	paddFactor,
		double	maxFreq,
		int	BSplinedegree
	)

Definition at line 76 of file fourier_projection.cpp.

{
    paddingFactor = paddFactor;
    maxFrequency = maxFreq;
    BSplineDeg = degree;
    updateVolume(V);
}

produceSideInfo()

void FourierProjector::produceSideInfo

(

)

Prepare the Spline coefficients and projection space.

Definition at line 247 of file fourier_projection.cpp.

{
    // Zero padding
    MultidimArray<double> Vpadded;
    auto paddedDim=(int)(paddingFactor*volumeSize);
    volume->window(Vpadded,FIRST_XMIPP_INDEX(paddedDim),FIRST_XMIPP_INDEX(paddedDim),FIRST_XMIPP_INDEX(paddedDim),
                   LAST_XMIPP_INDEX(paddedDim),LAST_XMIPP_INDEX(paddedDim),LAST_XMIPP_INDEX(paddedDim));
    volume->clear();
    // Make Fourier transform, shift the volume origin to the volume center and center it
    MultidimArray< std::complex<double> > Vfourier;
    FourierTransformer transformer3D;
    transformer3D.completeFourierTransform(Vpadded,Vfourier);
    ShiftFFT(Vfourier, FIRST_XMIPP_INDEX(XSIZE(Vpadded)), FIRST_XMIPP_INDEX(YSIZE(Vpadded)), FIRST_XMIPP_INDEX(ZSIZE(Vpadded)));
    CenterFFT(Vfourier,true);
    Vfourier.setXmippOrigin();

    // Compensate for the Fourier normalization factor
    double K=(double)(XSIZE(Vpadded)*XSIZE(Vpadded)*XSIZE(Vpadded))/(double)(volumeSize*volumeSize);
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(Vfourier)
    DIRECT_MULTIDIM_ELEM(Vfourier,n)*=K;
    Vpadded.clear();
    // Compute Bspline coefficients
    if (BSplineDeg==xmipp_transformation::BSPLINE3)
    {
        MultidimArray< double > VfourierRealAux;
        MultidimArray< double > VfourierImagAux;
        Complex2RealImag(Vfourier, VfourierRealAux, VfourierImagAux);
        Vfourier.clear();
        produceSplineCoefficients(xmipp_transformation::BSPLINE3,VfourierRealCoefs,VfourierRealAux);

        // Release memory as soon as you can
        VfourierRealAux.clear();

        // Remove all those coefficients we are sure we will not use during the projections
        volumePaddedSize=XSIZE(VfourierRealCoefs);
        int idxMax=maxFrequency*XSIZE(VfourierRealCoefs)+10; // +10 is a safety guard
        idxMax=std::min(FINISHINGX(VfourierRealCoefs),idxMax);
        int idxMin=std::max(-idxMax,STARTINGX(VfourierRealCoefs));
        VfourierRealCoefs.selfWindow(idxMin,idxMin,idxMin,idxMax,idxMax,idxMax);

        produceSplineCoefficients(xmipp_transformation::BSPLINE3,VfourierImagCoefs,VfourierImagAux);
        VfourierImagAux.clear();
        VfourierImagCoefs.selfWindow(idxMin,idxMin,idxMin,idxMax,idxMax,idxMax);
    }
    else {
        Complex2RealImag(Vfourier, VfourierRealCoefs, VfourierImagCoefs);
        volumePaddedSize=XSIZE(VfourierRealCoefs);
    }

    produceSideInfoProjection();
}

produceSideInfoProjection()

void FourierProjector::produceSideInfoProjection

(

)

Prepare projection space.

Definition at line 299 of file fourier_projection.cpp.

{
    // Allocate memory for the 2D Fourier transform
    projection().initZeros(volumeSize,volumeSize);
    projection().setXmippOrigin();
    transformer2D.FourierTransform(projection(),projectionFourier,false);

    // Calculate phase shift terms
    phaseShiftImgA.initZeros(projectionFourier);
    phaseShiftImgB.initZeros(projectionFourier);
    double shift=-FIRST_XMIPP_INDEX(volumeSize);
    double xxshift = -2 * PI * shift / volumeSize;
    for (size_t i=0; i<YSIZE(projectionFourier); ++i)
    {
        double phasey=(double)(i) * xxshift;
        for (size_t j=0; j<XSIZE(projectionFourier); ++j)
        {
            // Phase shift to move the origin of the image to the corner
            double dotp = (double)(j) * xxshift + phasey;
            //sincos(dotp,&DIRECT_A2D_ELEM(phaseShiftImgB,i,j),&DIRECT_A2D_ELEM(phaseShiftImgA,i,j));
            DIRECT_A2D_ELEM(phaseShiftImgB,i,j) = sin(dotp);
            DIRECT_A2D_ELEM(phaseShiftImgA,i,j) = cos(dotp);
        }
    }
}

project()

void FourierProjector::project	(	double	rot,
		double	tilt,
		double	psi,
		const MultidimArray< double > *	ctf = nullptr
	)

This method gets the volume's Fourier and the Euler's angles as the inputs and interpolates the related projection

Definition at line 91 of file fourier_projection.cpp.

{
    double freqy;
    double freqx;
    std::complex< double > f;
    Euler_angles2matrix(rot,tilt,psi,E);
    projectionFourier.initZeros();
    double maxFreq2=maxFrequency*maxFrequency;
    auto Xdim=(int)XSIZE(VfourierRealCoefs);
    auto Ydim=(int)YSIZE(VfourierRealCoefs);
    auto Zdim=(int)ZSIZE(VfourierRealCoefs);

    for (size_t i=0; i<YSIZE(projectionFourier); ++i)
    {
        FFT_IDX2DIGFREQ(i,volumeSize,freqy);
        double freqy2=freqy*freqy;

        double freqYvol_X=MAT_ELEM(E,1,0)*freqy;
        double freqYvol_Y=MAT_ELEM(E,1,1)*freqy;
        double freqYvol_Z=MAT_ELEM(E,1,2)*freqy;
        for (size_t j=0; j<XSIZE(projectionFourier); ++j)
        {
            // The frequency of pairs (i,j) in 2D
            FFT_IDX2DIGFREQ(j,volumeSize,freqx);

            // Do not consider pixels with high frequency
            if ((freqy2+freqx*freqx)>maxFreq2)
                continue;

            // Compute corresponding frequency in the volume
            double freqvol_X=freqYvol_X+MAT_ELEM(E,0,0)*freqx;
            double freqvol_Y=freqYvol_Y+MAT_ELEM(E,0,1)*freqx;
            double freqvol_Z=freqYvol_Z+MAT_ELEM(E,0,2)*freqx;

            double c;
            double d;
            if (BSplineDeg==xmipp_transformation::NEAREST)
            {
                // 0 order interpolation
                // Compute corresponding index in the volume
                auto kVolume=(int)round(freqvol_Z*volumePaddedSize);
                auto iVolume=(int)round(freqvol_Y*volumePaddedSize);
                auto jVolume=(int)round(freqvol_X*volumePaddedSize);
                c = A3D_ELEM(VfourierRealCoefs,kVolume,iVolume,jVolume);
                d = A3D_ELEM(VfourierImagCoefs,kVolume,iVolume,jVolume);
            }
            else if (BSplineDeg==xmipp_transformation::LINEAR)
            {
                // B-spline linear interpolation
                double kVolume=freqvol_Z*volumePaddedSize;
                double iVolume=freqvol_Y*volumePaddedSize;
                double jVolume=freqvol_X*volumePaddedSize;
                c=VfourierRealCoefs.interpolatedElement3D(jVolume,iVolume,kVolume);
                d=VfourierImagCoefs.interpolatedElement3D(jVolume,iVolume,kVolume);
            }
            else
            {
                // B-spline cubic interpolation
                double kVolume=freqvol_Z*volumePaddedSize;
                double iVolume=freqvol_Y*volumePaddedSize;
                double jVolume=freqvol_X*volumePaddedSize;

                // Commented for speed-up, the corresponding code is below
                // c=VfourierRealCoefs.interpolatedElementBSpline3D(jVolume,iVolume,kVolume);
                // d=VfourierImagCoefs.interpolatedElementBSpline3D(jVolume,iVolume,kVolume);

                // The code below is a replicate for speed reasons of interpolatedElementBSpline3D
                double z=kVolume;
                double y=iVolume;
                double x=jVolume;

                // Logical to physical
                z -= STARTINGZ(VfourierRealCoefs);
                y -= STARTINGY(VfourierRealCoefs);
                x -= STARTINGX(VfourierRealCoefs);

                auto l1 = (int)ceil(x - 2);
                int l2 = l1 + 3;

                auto m1 = (int)ceil(y - 2);
                int m2 = m1 + 3;

                auto n1 = (int)ceil(z - 2);
                int n2 = n1 + 3;

                c = d = 0.0;
                double aux;
                for (int nn = n1; nn <= n2; nn++)
                {
                    int equivalent_nn=nn;
                    if      (nn<0)
                        equivalent_nn=-nn-1;
                    else if (nn>=Zdim)
                        equivalent_nn=2*Zdim-nn-1;
                    double yxsumRe = 0.0;
                    double yxsumIm = 0.0;
                    for (int m = m1; m <= m2; m++)
                    {
                        int equivalent_m=m;
                        if      (m<0)
                            equivalent_m=-m-1;
                        else if (m>=Ydim)
                            equivalent_m=2*Ydim-m-1;
                        double xsumRe = 0.0;
                        double xsumIm = 0.0;
                        for (int l = l1; l <= l2; l++)
                        {
                            double xminusl = x - (double) l;
                            int equivalent_l=l;
                            if      (l<0)
                                equivalent_l=-l-1;
                            else if (l>=Xdim)
                                equivalent_l=2*Xdim-l-1;
                            auto CoeffRe = (double) DIRECT_A3D_ELEM(VfourierRealCoefs,equivalent_nn,equivalent_m,equivalent_l);
                            auto CoeffIm = (double) DIRECT_A3D_ELEM(VfourierImagCoefs,equivalent_nn,equivalent_m,equivalent_l);
                            BSPLINE03(aux,xminusl);
                            xsumRe += CoeffRe * aux;
                            xsumIm += CoeffIm * aux;
                        }

                        double yminusm = y - (double) m;
                        BSPLINE03(aux,yminusm);
                        yxsumRe += xsumRe * aux;
                        yxsumIm += xsumIm * aux;
                    }

                    double zminusn = z - (double) nn;
                    BSPLINE03(aux,zminusn);
                    c += yxsumRe * aux;
                    d += yxsumIm * aux;
                }
            }
            // Phase shift to move the origin of the image to the corner
            double a=DIRECT_A2D_ELEM(phaseShiftImgA,i,j);
            double b=DIRECT_A2D_ELEM(phaseShiftImgB,i,j);
            if (ctf!=nullptr)
            {
                double ctfij=DIRECT_A2D_ELEM(*ctf,i,j);
                a*=ctfij;
                b*=ctfij;
            }

            // Multiply Fourier coefficient in volume times phase shift
            double ac = a * c;
            double bd = b * d;
            double ab_cd = (a + b) * (c + d);

            // And store the multiplication
            auto *ptrI_ij=(double *)&DIRECT_A2D_ELEM(projectionFourier,i,j);
            *ptrI_ij = ac - bd;
            *(ptrI_ij+1) = ab_cd - ac - bd;
        }
    }
    transformer2D.inverseFourierTransform();
}

read()

void Projection::read	(	const FileName &	fn,
		const bool	only_apply_shifts = false,
		DataMode	datamode = DATA,
		MDRow *	row = nullptr
	)

Read a Projection from file.

When a projection is read, the Euler matrices and perpendicular direction is computed and stored in the Projection structure.

Definition at line 50 of file fourier_projection.cpp.

{
    Image<double>::read(fn, datamode);
    if (row != nullptr)
        applyGeo(*row, only_apply_shifts);
    Euler_angles2matrix(rot(), tilt(), psi(), euler);
    eulert = euler.transpose();
    euler.getRow(2, direction);
    direction.selfTranspose();
}

reset()

void Projection::reset	(	int	Ydim,
		int	Xdim
	)

Init_zeros and move origin to center.

This function initialises the projection plane with 0's, and then moves the logical origin of the image to the physical center of the image (using the Xmipp conception of image center).

Definition at line 33 of file fourier_projection.cpp.

{
    data.initZeros(Ydim, Xdim);
    data.setXmippOrigin();
}

setAngles()

void Projection::setAngles	(	double	_rot,
		double	_tilt,
		double	_psi
	)

Set Euler angles for this projection.

The Euler angles are stored in the Xmipp header, then the pass matrices (Universe <—> Projection coordinate system) are computed, and the vector perpendicular to this projection plane is also calculated.

Definition at line 40 of file fourier_projection.cpp.

{
    setEulerAngles(_rot, _tilt, _psi);
    Euler_angles2matrix(_rot, _tilt, _psi, euler);
    eulert = euler.transpose();
    euler.getRow(2, direction);
    direction.selfTranspose();
}

updateVolume()

void FourierProjector::updateVolume

(

MultidimArray< double > &

V

)

Update volume

Definition at line 84 of file fourier_projection.cpp.

{
    volume = &V;
    volumeSize=XSIZE(*volume);
    produceSideInfo();
}

Variable Documentation

BSplineDeg

double FourierProjector::BSplineDeg

The order of B-Spline for interpolation.

Definition at line 119 of file fourier_projection.h.

direction

Matrix1D< double > Projection::direction

Vector perpendicular to the projection plane. It is calculated as a function of rot and tilt.

Definition at line 62 of file fourier_projection.h.

E

Matrix2D<double> FourierProjector::E

Definition at line 149 of file fourier_projection.h.

euler

Matrix2D< double > Projection::euler

Matrix used to pass from the Universal coordinate system to the projection coordinate system.

Rp = euler * Ru

Definition at line 71 of file fourier_projection.h.

eulert

Matrix2D< double > Projection::eulert

Just the opposite.

Ru = eulert * Rp.

Definition at line 79 of file fourier_projection.h.

maxFrequency

double FourierProjector::maxFrequency

Maximum Frequency for pixels.

Definition at line 117 of file fourier_projection.h.

paddingFactor

double FourierProjector::paddingFactor

Padding factor.

Definition at line 115 of file fourier_projection.h.

phaseShiftImgA

MultidimArray<double> FourierProjector::phaseShiftImgA

Definition at line 140 of file fourier_projection.h.

phaseShiftImgB

MultidimArray<double> FourierProjector::phaseShiftImgB

Definition at line 139 of file fourier_projection.h.

projection

Image<double> FourierProjector::projection

Definition at line 136 of file fourier_projection.h.

projectionFourier

MultidimArray< std::complex<double> > FourierProjector::projectionFourier

Definition at line 133 of file fourier_projection.h.

transformer2D

FourierTransformer FourierProjector::transformer2D

Definition at line 123 of file fourier_projection.h.

VfourierImagCoefs

MultidimArray< double > FourierProjector::VfourierImagCoefs

Definition at line 130 of file fourier_projection.h.

VfourierRealCoefs

MultidimArray< double > FourierProjector::VfourierRealCoefs

Definition at line 129 of file fourier_projection.h.

volume

MultidimArray<double>* FourierProjector::volume

Definition at line 126 of file fourier_projection.h.

volumePaddedSize

int FourierProjector::volumePaddedSize

Definition at line 146 of file fourier_projection.h.

volumeSize

int FourierProjector::volumeSize

Definition at line 143 of file fourier_projection.h.