Fourier transforms

Collaboration diagram for Fourier transforms:

Index <--> Frequency, Continuous <--> Discrete
template<typename T >
void	FFT_idx2digfreq (T &v, const Matrix1D< int > &idx, Matrix1D< double > &freq)
template<typename T >
void	digfreq2FFT_idx (T &v, const Matrix1D< double > &freq, Matrix1D< int > &idx)
void	digfreq2contfreq (const Matrix1D< double > &digfreq, Matrix1D< double > &contfreq, double pixel_size)
void	contfreq2digfreq (const Matrix1D< double > &contfreq, Matrix1D< double > &digfreq, double pixel_size)
#define	FFT_IDX2DIGFREQ(idx, size, freq)
#define	FFT_IDX2DIGFREQ_DOUBLE(idx, size, freq)
#define	FFT_IDX2DIGFREQ_FAST(idx, size, size_2, isize, freq)   freq = ( ((idx) <= (size_2)) ? (idx) : -(size) + (idx) ) * (isize);
#define	DIGFREQ2FFT_IDX(freq, size, idx)
#define	DIGFREQ2FFT_IDX_DOUBLE(freq, size, idx)

Format conversions
void	Whole2Half (const MultidimArray< std::complex< double > > &in, MultidimArray< std::complex< double > > &out)
void	Half2Whole (const MultidimArray< std::complex< double > > &in, MultidimArray< std::complex< double > > &out, size_t oridim)
void	Complex2RealImag (const MultidimArray< std::complex< double > > &in, MultidimArray< double > &real, MultidimArray< double > &imag)
void	RealImag2Complex (const MultidimArray< double > &real, const MultidimArray< double > &imag, MultidimArray< std::complex< double > > &out)

Fourier Transforms
The theoretical relationship between the Fourier transform of a discrete signal and the Fourier transform of the continuous signal is X(e^jw)=1/T*X_c(jw/T) Xmipp is not computing X(e^jw) but samples from it so that X(e^jw)=N*X_XMIPP[k] where N is the length of the signal being transformed and X_XMIPP[k] is the k-th sample. The following program illustrates how the continuous, discrete and Xmipp Fourier transform relate #include <data/matrix1d.h> #include <data/xmipp_fft.h> double discreteTransform(double w, int N1) { if (w==0) return 2*N1+1; else return sin(w*(N1+0.5))/sin(0.5*w); } double continuousTransform(double W, double T1) { if (W==0) return 2*T1; else return 2*sin(W*T1)/W; } int main() { try { MultidimArray<double> x(65); x.setXmippOrigin(); double T=0.5; double T1=6; int N1=(int)CEIL(T1/T); // Fill x with a pulse from -N1 to N1 (-T1 to T1 in continuous) FOR_ALL_ELEMENTS_IN_ARRAY1D(x) if (ABS(i)<=N1) x(i)=1; // Compute the Fourier transform MultidimArray< std::complex<double> > X; MultidimArray<double> Xmag; FourierTransform(x,X); FFT_magnitude(X,Xmag); // Compute the frequency axes MultidimArray<double> contfreq(XSIZE(X)), digfreq(XSIZE(X)); FOR_ALL_ELEMENTS_IN_ARRAY1D(X) FFT_IDX2DIGFREQ(i,XSIZE(X),digfreq(i)); digfreq*=2*PI; contfreq=digfreq/T; // Show all Fourier transforms FOR_ALL_ELEMENTS_IN_ARRAY1D(X) { if (digfreq(i)>=0) std::cout << digfreq(i) << " " << contfreq(i) << " " << XSIZE(X)*Xmag(i) << " " << ABS(discreteTransform(digfreq(i),N1)) << " " << ABS(continuousTransform(contfreq(i),T1)/T) << std::endl; } } catch (XmippError XE) { std::cout << XE << std::endl; } return 0; }
void	FourierTransform (const MultidimArray< double > &in, MultidimArray< std::complex< double > > &out)
void	InverseFourierTransform (const MultidimArray< std::complex< double > > &in, MultidimArray< double > &out)
void	FourierTransformHalf (const MultidimArray< double > &in, MultidimArray< std::complex< double > > &out)
void	InverseFourierTransformHalf (const MultidimArray< std::complex< double > > &in, MultidimArray< double > &out, int oridim)

Operations with the Fourier Transforms
void	centerFFT2 (MultidimArray< double > &v)
template<typename T >
void	CenterFFT (MultidimArray< T > &v, bool forward)
void	ShiftFFT (MultidimArray< std::complex< double > > &v, double xshift)
void	ShiftFFT (MultidimArray< std::complex< double > > &v, double xshift, double yshift)
void	ShiftFFT (MultidimArray< std::complex< double > > &v, double xshift, double yshift, double zshift)
void	CenterOriginFFT (MultidimArray< std::complex< double > > &v, bool forward)
template<typename T >
void	xmipp2PSD (const MultidimArray< T > &input, MultidimArray< T > &output, bool takeLog=true)
void	xmipp2CTF (const MultidimArray< double > &input, MultidimArray< double > &output)
#define	SWAP_ARRAY(a, b, n)   memcpy(buffer, a, n); memcpy(a, b, n); memcpy(b, buffer, n);

Detailed Description

Macro Definition Documentation

DIGFREQ2FFT_IDX

#define DIGFREQ2FFT_IDX	(		freq,
			size,
			idx
	)

Value:

{ \
        (idx) = (int) (round((size) * (freq))); if ((idx) < 0) (idx) += (int) \
                    (size); }

Frequency to index (int)

Given a frequency and a size of the FFT, this macro returns the corresponding integer index

Definition at line 61 of file xmipp_fft.h.

DIGFREQ2FFT_IDX_DOUBLE

#define DIGFREQ2FFT_IDX_DOUBLE	(		freq,
			size,
			idx
	)

Value:

{ \
        (idx) = ((size) * (freq)); if ((idx) < 0) (idx) += (size); }

Frequency to index (double)

Given a frequency and a size of the FFT, this macro returns the corresponding double index

Definition at line 70 of file xmipp_fft.h.

FFT_IDX2DIGFREQ

#define FFT_IDX2DIGFREQ	(		idx,
			size,
			freq
	)

Value:

freq = (size<=1)? 0:(( (((int)idx) <= (((int)(size)) >> 1)) ? ((int)(idx)) : -((int)(size)) + ((int)(idx))) / \
           (double)(size));

Index to frequency

Given an index and a size of the FFT, this function returns the corresponding digital frequency (-1/2 to 1/2)

Definition at line 46 of file xmipp_fft.h.

FFT_IDX2DIGFREQ_DOUBLE

#define FFT_IDX2DIGFREQ_DOUBLE	(		idx,
			size,
			freq
	)

Value:

freq = (size<=1)? 0:(( (((double)idx) <= (((double)(size)) / 2.0)) ? ((double)(idx)) : -((double)(size)) + ((double)(idx))) / \
           (double)(size));

Definition at line 50 of file xmipp_fft.h.

FFT_IDX2DIGFREQ_FAST

#define FFT_IDX2DIGFREQ_FAST	(		idx,
			size,
			size_2,
			isize,
			freq
	)		freq = ( ((idx) <= (size_2)) ? (idx) : -(size) + (idx) ) * (isize);

Definition at line 54 of file xmipp_fft.h.

SWAP_ARRAY

#define SWAP_ARRAY	(		a,
			b,
			n
	)		memcpy(buffer, a, n); memcpy(a, b, n); memcpy(b, buffer, n);

Faster version of CenterFFT (now just for even images)

Definition at line 278 of file xmipp_fft.h.

Function Documentation

CenterFFT()

template<typename T >

void CenterFFT	(	MultidimArray< T > &	v,
		bool	forward
	)

CenterFFT Relation with Matlab fftshift: forward true is equals to fftshift and forward false equals to ifftshift

Definition at line 291 of file xmipp_fft.h.

{
    bool firstTime=true;                        // First time executing inner loops.

    // Check dimension is between 1 and 3 inclusive.
    if ( v.getDim() > 0 && v.getDim() <= 3)
    {
        // Shift in the X direction
        size_t l=0;

        // Shift in the X direction
        if (((l = XSIZE(v)) > 1) && (YSIZE(v) == 1) && (ZSIZE(v) == 1))
        {
            size_t firstHalfSize=0;
            size_t secondHalfSize=0;
            MultidimArray< T > aux;

            if (forward)
            {
                firstHalfSize  = (l + 1) / 2;
                secondHalfSize = l - firstHalfSize;
                aux.resizeNoCopy( firstHalfSize);

                for (size_t k = 0; k < ZSIZE(v); k++)
                {
                    for (size_t i = 0; i < YSIZE(v); i++)
                    {
                        memcpy( &dAi(aux, 0), &dAkij(v, k, i, 0), sizeof(T)*firstHalfSize);
                        memcpy( &dAkij(v, k, i, 0), &dAkij(v, k, i, firstHalfSize), sizeof(T)*secondHalfSize);
                        memcpy( &dAkij(v, k, i, secondHalfSize), &dAi(aux, 0), sizeof(T)*firstHalfSize);
                    }
                }
            }
            else
            {
                secondHalfSize = (l + 1) / 2;
                firstHalfSize  = l - secondHalfSize;
                aux.resizeNoCopy( secondHalfSize);

                for (size_t k = 0; k < ZSIZE(v); k++)
                {
                    for (size_t i = 0; i < YSIZE(v); i++)
                    {
                        memcpy( &dAi(aux, 0), &dAkij(v, k, i, firstHalfSize), sizeof(T)*secondHalfSize);
                        memcpy( &dAkij(v, k, i, secondHalfSize), &dAkij(v, k, i, 0), sizeof(T)*firstHalfSize);
                        memcpy( &dAkij(v, k, i, 0), &dAi(aux, 0), sizeof(T)*secondHalfSize);
                    }
                }
            }
        }
        else
        {
            // 3D
            MultidimArray< T > aux;
            size_t  l=0;
            long int shift;

            int     i=0;                            // Loop counter.
            int     halfRows=0;                     // Half rows of the matrix.
            int     rowSize=0;                      // Size in bytes of the row.

            // Size in bytes of first and second half row.
            int     firstHalfRowSize=0, secondHalfRowSize=0;

            // # elements in first and second half row.
            int     nElemsFirstHalf_X=0, nElemsSecondHalf_X=0;

            MultidimArray< T >  tempVector;         // Temporary vector.

            bool    isOdd=false;                    // Odd # rows in matrix.
            MultidimArray< T >  savedRow;           // Temporary vector in odd matrix.

            // Execute FFT in 2 dimensions before 3D.
            if (forward)
            {
                for (size_t k=0; k<ZSIZE(v); k++)
                {
                    // Pointers to upper and lower half matrix.
                    T       *upperHalfPtr, *lowerHalfSrcPtr, *lowerHalfDestPtr;

                    // This is computed first time only.
                    if (firstTime)
                    {
                        firstTime = false;

                        // Compute # elements in each half vector in X-dimension.
                        nElemsFirstHalf_X = (XSIZE(v) + 1) / 2;
                        nElemsSecondHalf_X = XSIZE(v) - nElemsFirstHalf_X;

                        // Compute X dimension size in bytes.
                        rowSize = XSIZE(v)*sizeof(T);

                        // Allocate temporary vector.
                        firstHalfRowSize = nElemsFirstHalf_X*sizeof(T);
                        secondHalfRowSize = rowSize - firstHalfRowSize;
                        tempVector.resizeNoCopy( XSIZE(v));

                        // Compute # iterations.
                        halfRows = YSIZE(v) / 2;

                        // If even # rows then save row located in the middle of the matrix.
                        if ((YSIZE(v) % 2) != 0)
                        {
                            isOdd = true;
                            savedRow.resizeNoCopy( XSIZE(v));
                        }
                    }

                    // Initialize pointers to upper and lower matrix rows.
                    upperHalfPtr     = &dAkij(v, k, 0, 0);
                    lowerHalfSrcPtr  = upperHalfPtr + (halfRows*XSIZE(v));
                    lowerHalfDestPtr = lowerHalfSrcPtr;

                    // If even # rows then save row located in the middle of the matrix.
                    if (isOdd)
                    {
                        memcpy( &dAi(savedRow,0), lowerHalfSrcPtr, rowSize);

                        // Source and destination rows are not the same in odd matrix.
                        lowerHalfSrcPtr = lowerHalfSrcPtr + XSIZE(v);
                    }

                    // Half of the rows must be exchanged.
                    for (i=0; i<halfRows ;i++)
                    {
                        memcpy( &dAi(tempVector,0), upperHalfPtr, rowSize);
                        memcpy( upperHalfPtr, lowerHalfSrcPtr + nElemsFirstHalf_X, secondHalfRowSize);
                        memcpy( upperHalfPtr + nElemsSecondHalf_X, lowerHalfSrcPtr, firstHalfRowSize);
                        memcpy( lowerHalfDestPtr, &dAi(tempVector, nElemsFirstHalf_X), secondHalfRowSize);
                        memcpy( lowerHalfDestPtr + nElemsSecondHalf_X, &dAi(tempVector,0), firstHalfRowSize);

                        upperHalfPtr     = upperHalfPtr + XSIZE(v);
                        lowerHalfSrcPtr  = lowerHalfSrcPtr + XSIZE(v);
                        lowerHalfDestPtr = lowerHalfDestPtr + XSIZE(v);
                    }

                    // If # rows is odd then restore last row.
                    if (isOdd)
                    {
                        memcpy( lowerHalfDestPtr + nElemsSecondHalf_X, &dAi(savedRow,0), firstHalfRowSize);
                        memcpy( lowerHalfDestPtr, &dAi(savedRow, nElemsFirstHalf_X), secondHalfRowSize);
                    }
                }
            }
            else
            {
                for (size_t k = 0; k<ZSIZE(v); k++)
                {
                    // Pointers to upper and lower half matrix.
                    T       *upperHalfSrcPtr, *upperHalfDestPtr, *lowerHalfPtr;

                    // This is computed first time only.
                    if (firstTime)
                    {
                        firstTime = false;

                        // Compute # elements in each half vector in X-dimension.
                        nElemsSecondHalf_X = (XSIZE(v) + 1) / 2;
                        nElemsFirstHalf_X  = XSIZE(v) - nElemsSecondHalf_X;

                        // Allocate temporary vector.
                        rowSize = XSIZE(v)*sizeof(T);
                        firstHalfRowSize = nElemsFirstHalf_X*sizeof(T);
                        secondHalfRowSize = rowSize - firstHalfRowSize;
                        tempVector.resizeNoCopy(rowSize);

                        // Compute # iterations.
                        halfRows = YSIZE(v) / 2;

                        // If odd # rows then save last row.
                        if ((YSIZE(v) % 2) != 0)
                        {
                            isOdd = true;
                            savedRow.resizeNoCopy(XSIZE(v));
                        }
                    }

                    // Initialize pointers to quadrants.
                    upperHalfSrcPtr  = &dAkij(v, k, halfRows-1, 0);
                    lowerHalfPtr     = &dAkij(v, k, YSIZE(v)-1, 0);
                    upperHalfDestPtr = upperHalfSrcPtr;

                    // If odd # rows then save last row.
                    if (isOdd)
                    {
                        // Source and destination rows are not the same in odd matrix.
                        upperHalfDestPtr = upperHalfDestPtr + XSIZE(v);

                        memcpy( &dAi(savedRow,0), upperHalfDestPtr, rowSize);
                    }

                    // Half of the rows must be exchanged.
                    for (i=0; i<halfRows ;i++)
                    {
                        memcpy( &dAi(tempVector, 0), upperHalfSrcPtr, rowSize);
                        memcpy( upperHalfDestPtr, lowerHalfPtr + nElemsFirstHalf_X, secondHalfRowSize);
                        memcpy( upperHalfDestPtr + nElemsSecondHalf_X, lowerHalfPtr, firstHalfRowSize);
                        memcpy( lowerHalfPtr, &dAi(tempVector, nElemsFirstHalf_X), secondHalfRowSize);
                        memcpy( lowerHalfPtr + nElemsSecondHalf_X, &dAi(tempVector,0), firstHalfRowSize);

                        upperHalfSrcPtr  = upperHalfSrcPtr - XSIZE(v);
                        lowerHalfPtr     = lowerHalfPtr - XSIZE(v);
                        upperHalfDestPtr = upperHalfDestPtr - XSIZE(v);
                    }

                    // If # rows is odd then restore row at medium position.
                    if (isOdd)
                    {
                        memcpy( upperHalfDestPtr + nElemsSecondHalf_X, &dAi(savedRow, 0), firstHalfRowSize);
                        memcpy( upperHalfDestPtr, &dAi(savedRow, nElemsFirstHalf_X), secondHalfRowSize);
                    }
                }
            }

            // Shift in the Z direction
            if ((l = ZSIZE(v)) > 1)
            {
                aux.resizeNoCopy(l);
                shift = (long int)(l / 2);

                if (!forward)
                    shift = -shift;
                size_t lmax=(l/4)*4;
                for (size_t i = 0; i < YSIZE(v); i++)
                {
                    for (size_t j = 0; j < XSIZE(v); j++)
                    {
                        // Shift the input in an auxiliary vector
                        for (size_t k = 0; k < l; k++)
                        {
                            size_t kp = k + shift;
                            if (-shift > (long int)k)
                                kp += l;
                            else if (kp >= l)
                                kp -= l;

                            dAi(aux,kp) = dAkij(v, k, i, j);
                        }

                        // Copy the vector
                        const T* ptrAux=&dAi(aux,0);
                        for (size_t k = 0; k < lmax; k+=4,ptrAux+=4)
                        {
                            dAkij(v, k,   i, j) = *ptrAux;
                            dAkij(v, k+1, i, j) = *(ptrAux+1);
                            dAkij(v, k+2, i, j) = *(ptrAux+2);
                            dAkij(v, k+3, i, j) = *(ptrAux+3);
                        }
                        for (size_t k = lmax; k < l; ++k, ++ptrAux)
                        {
                            dAkij(v, k, i, j) = *ptrAux;
                        }
                    }
                }
            }
        }
    }
    else
    {
        REPORT_ERROR(ERR_MULTIDIM_DIM,"CenterFFT ERROR: Dimension should be 1, 2 or 3");
    }
}

centerFFT2()

void centerFFT2

(

MultidimArray< double > &

v

)

Center FFT for 2D arrays. The function is optimized for the particular case of 2D.

Definition at line 276 of file xmipp_fft.cpp.

{
    if (v.getDim() == 2)
    {
        //Just separe the even and odd dimensions case
        if (XSIZE(v) % 2 == 0 && YSIZE(v) % 2 == 0)
        {
            int xsize = XSIZE(v);
            int xhalf = xsize / 2;
            int yhalf = YSIZE(v) / 2;
            double * posA = MULTIDIM_ARRAY(v);
            double * posB = posA + xhalf;
            double * posC = posA + xsize * yhalf;
            double * posD = posC + xhalf;
            double * buffer=new double[xhalf];
            size_t bytes = xhalf * sizeof(double);

            for (int i = 0; i < yhalf; ++i,
                 posA += xsize, posB += xsize, posC += xsize, posD += xsize)
            {
                SWAP_ARRAY(posA, posD, bytes);
                SWAP_ARRAY(posB, posC, bytes);
            }
            delete []buffer;
        }
        else
        {
            //todo: implementation for the odd case needed
            CenterFFT(v, true);
        }
    }
    else
        std::cerr <<"bad dim: " << v.getDim() << std::endl;
}

CenterOriginFFT()

void CenterOriginFFT	(	MultidimArray< std::complex< double > > &	v,
		bool	forward
	)

Place the origin of the FFT at the center of the vector and back

Changes the real and the fourier space origin

Definition at line 386 of file xmipp_fft.cpp.

{
    if ( v.getDim() == 1 )
    {
        // 1D
        double xshift = -(double)(int)(XSIZE(v) / 2);
        if (forward)
        {
            ShiftFFT(v, xshift);
            CenterFFT(v, forward);
        }
        else
        {
            CenterFFT(v, forward);
            ShiftFFT(v, -xshift);
        }
    }
    else if ( v.getDim() == 2)
    {
        // 2D

        double xshift = -(double)(int)(XSIZE(v) / 2);
        double yshift = -(double)(int)(YSIZE(v) / 2);
        if (forward)
        {
            ShiftFFT(v, xshift, yshift);
            CenterFFT(v, forward);
        }
        else
        {
            CenterFFT(v, forward);
            ShiftFFT(v, -xshift, -yshift);
        }
    }
    else if ( v.getDim() == 3)
    {
        // 3D
        double xshift = -(double)(int)(XSIZE(v) / 2);
        double yshift = -(double)(int)(YSIZE(v) / 2);
        double zshift = -(double)(int)(ZSIZE(v) / 2);
        if (forward)
        {
            ShiftFFT(v, xshift, yshift, zshift);
            CenterFFT(v, forward);
        }
        else
        {
            CenterFFT(v, forward);
            ShiftFFT(v, -xshift, -yshift, -zshift);
        }
    }
    else
        REPORT_ERROR(ERR_MULTIDIM_DIM,"CenterOriginFFT ERROR: only valis for 1D or 2D or 3D");
}

Complex2RealImag()

void Complex2RealImag	(	const MultidimArray< std::complex< double > > &	in,
		MultidimArray< double > &	real,
		MultidimArray< double > &	imag
	)

Conversion from complex -> real,imag

Convert complex -> real,imag Fourier transforms 2D. –

Definition at line 103 of file xmipp_fft.cpp.

{
    real.resizeNoCopy(in);
    imag.resizeNoCopy(in);
    Complex2RealImag(MULTIDIM_ARRAY(in), MULTIDIM_ARRAY(real),
                     MULTIDIM_ARRAY(imag), MULTIDIM_SIZE(in));
}

contfreq2digfreq()

void contfreq2digfreq ( const Matrix1D< double > &  contfreq, Matrix1D< double > &  digfreq, double  pixel_size  )

inline

Continuous to Digital frequency

The pixel size must be given in Amstrongs. The digital frequency is between [-1/2,1/2]

Definition at line 139 of file xmipp_fft.h.

{
    digfreq.resizeNoCopy(contfreq);
    FOR_ALL_ELEMENTS_IN_MATRIX1D(contfreq)
    VEC_ELEM(digfreq,i) = VEC_ELEM(contfreq,i) * pixel_size;
}

digfreq2contfreq()

void digfreq2contfreq ( const Matrix1D< double > &  digfreq, Matrix1D< double > &  contfreq, double  pixel_size  )

inline

Digital to Continuous frequency

The pixel size must be given in Amstrongs. The digital frequency is between [-1/2,1/2]

Definition at line 125 of file xmipp_fft.h.

{
    contfreq.resizeNoCopy(digfreq);
    FOR_ALL_ELEMENTS_IN_MATRIX1D(digfreq)
    VEC_ELEM(contfreq,i) = VEC_ELEM(digfreq,i) / pixel_size;
}

digfreq2FFT_idx()

template<typename T >

void digfreq2FFT_idx	(	T &	v,
		const Matrix1D< double > &	freq,
		Matrix1D< int > &	idx
	)

Frequency to index

This function can be used with vectors of any size (1,2,3). The Digital spectrum is lim:confirm bd ited between -1/2 and 1/2. If the vector has got more than 3 coordinates, then an exception is thrown

Definition at line 106 of file xmipp_fft.h.

{
    if (VEC_XSIZE(freq) < 1 || VEC_XSIZE(freq) > 3)
        REPORT_ERROR(ERR_MATRIX_SIZE, "digfreq2FFT_idx: freq is not of the correct size");

    idx.resizeNoCopy(VEC_XSIZE(freq));

    int size[3];
    v.getSize(size);

    FOR_ALL_ELEMENTS_IN_MATRIX1D(idx)
    DIGFREQ2FFT_IDX(VEC_ELEM(freq,i), size[i], VEC_ELEM(idx,i));
}

FFT_idx2digfreq()

template<typename T >

void FFT_idx2digfreq	(	T &	v,
		const Matrix1D< int > &	idx,
		Matrix1D< double > &	freq
	)

Index to frequency

This function can be used with vectors of any size (1,2,3). The Digital spectrum is limited between -1/2 and 1/2. If the vector has got more than 3 coordinates, then an exception is thrown

Definition at line 80 of file xmipp_fft.h.

{
    if (VEC_XSIZE(idx) < 1 || VEC_XSIZE(idx) > 3)
        REPORT_ERROR(ERR_MATRIX_SIZE, "FFT_idx2digfreq: Index is not of the correct size");

    freq.resizeNoCopy(VEC_XSIZE(idx));

    switch (VEC_XSIZE(idx))
    {
    case 3:
        FFT_IDX2DIGFREQ(VEC_ELEM(idx,2), ZSIZE(v), VEC_ELEM(freq,2));
    case 2:
        FFT_IDX2DIGFREQ(VEC_ELEM(idx,1), YSIZE(v), VEC_ELEM(freq,1));
    case 1:
        FFT_IDX2DIGFREQ(VEC_ELEM(idx,0), XSIZE(v), VEC_ELEM(freq,0));
    }
}

FourierTransform()

void FourierTransform	(	const MultidimArray< double > &	in,
		MultidimArray< std::complex< double > > &	out
	)

Direct Fourier Transform

Direct Fourier Transform nD ----------------------------------------—

Definition at line 124 of file xmipp_fft.cpp.

{
    if ( in.getDim() == 1 )
    {
        // 1D
        int N = XSIZE(in);
        MultidimArray<double> re(in), tmp(N), im(N), cas(N);
        out.resizeNoCopy(N);

        GetCaS(MULTIDIM_ARRAY(cas), N);
        DftRealToRealImaginary(MULTIDIM_ARRAY(re), MULTIDIM_ARRAY(im),
                               MULTIDIM_ARRAY(tmp), MULTIDIM_ARRAY(cas), N);
        RealImag2Complex(MULTIDIM_ARRAY(re), MULTIDIM_ARRAY(im),
                         MULTIDIM_ARRAY(out), N);
    }
    else
    {
        // 2D and 3D
        int Status;
        MultidimArray<double> re(in), im;
        im.resizeNoCopy(in);
        out.resizeNoCopy(in);
        VolumeDftRealToRealImaginary(MULTIDIM_ARRAY(re),
                                     MULTIDIM_ARRAY(im), XSIZE(in), YSIZE(in), ZSIZE(in), &Status);
        RealImag2Complex(re,im,out);
    }

}

FourierTransformHalf()

void FourierTransformHalf	(	const MultidimArray< double > &	in,
		MultidimArray< std::complex< double > > &	out
	)

Direct Fourier Transform, output half of (centro-symmetric) transform

Direct Fourier Transform 1D / 2D, output half of (centro-symmetric) transform -—

Definition at line 187 of file xmipp_fft.cpp.

{

    MultidimArray<std::complex<double> > aux;
    FourierTransform(in, aux);
    Whole2Half(aux, out);
}

Half2Whole()

void Half2Whole	(	const MultidimArray< std::complex< double > > &	in,
		MultidimArray< std::complex< double > > &	out,
		size_t	oridim
	)

Conversion from half -> whole

Convert half -> whole of (centro-symmetric) Fourier transforms 2D. –

Definition at line 66 of file xmipp_fft.cpp.

{
    if (in.getDim() == 1)
    {
        // 1D
        out.resizeNoCopy(oridim);
        for (size_t j = 0; j < XSIZE(in); j++)
            DIRECT_A1D_ELEM(out,j) = DIRECT_A1D_ELEM(in,j);
        for (size_t j = XSIZE(in); j < oridim; j++)
            DIRECT_A1D_ELEM(out,j) = conj(DIRECT_A1D_ELEM(in,oridim - j));
    }
    else if (in.getDim() == 2)
    {
        // 2D
        out.resizeNoCopy(oridim, XSIZE(in));

        // Old part
        for (size_t i = 0; i < YSIZE(in); i++)
            for (size_t j = 0; j < XSIZE(in); j++)
                dAij(out, i, j) = dAij(in, i, j);

        // Complete first column of old part
        for (size_t j = YSIZE(in); j < XSIZE(in); j++)
            dAij(out, 0, j) = conj(dAij(in, 0, XSIZE(in) - j));

        // New part
        for (size_t i = YSIZE(in); i < oridim; i++)
        {
            dAij(out, i, 0) = conj(dAij(in, oridim - i, 0));
            for (size_t j = 1; j < XSIZE(in); j++)
                dAij(out, i, j) = conj(dAij(in, oridim - i, XSIZE(in) - j));
        }
    }
}

InverseFourierTransform()

void InverseFourierTransform	(	const MultidimArray< std::complex< double > > &	in,
		MultidimArray< double > &	out
	)

Inverse Fourier Transform

Inverse Fourier Transform nD. --------------------------------------—

Definition at line 155 of file xmipp_fft.cpp.

{
    if ( in.getDim() == 1 )
    {
        // 1D
        int N = XSIZE(in);
        MultidimArray<double> tmp(N), im(N), cas(N);
        out.resizeNoCopy(N);

        GetCaS(MULTIDIM_ARRAY(cas), N);
        Complex2RealImag(MULTIDIM_ARRAY(in), MULTIDIM_ARRAY(out),
                         MULTIDIM_ARRAY(im), N);
        InvDftRealImaginaryToReal(MULTIDIM_ARRAY(out), MULTIDIM_ARRAY(im),
                                  MULTIDIM_ARRAY(tmp), MULTIDIM_ARRAY(cas), N);
    }
    else
    {
        // 2D and 3D
        int Status;
        MultidimArray<double> im;
        out.resizeNoCopy(in);
        im.resizeNoCopy(in);
        Complex2RealImag(in, out, im);
        VolumeInvDftRealImaginaryToReal(MULTIDIM_ARRAY(out),
                                        MULTIDIM_ARRAY(im),
                                        XSIZE(in), YSIZE(in), ZSIZE(in), &Status);
    }
}

InverseFourierTransformHalf()

void InverseFourierTransformHalf	(	const MultidimArray< std::complex< double > > &	in,
		MultidimArray< double > &	out,
		int	oridim
	)

Inverse Fourier Transform 1D, input half of (centro-symmetric) transform

Inverse Fourier Transform 1D / 2D, input half of (centro-symmetric) transform -—

Definition at line 197 of file xmipp_fft.cpp.

{

    MultidimArray< std::complex<double> > aux;
    Half2Whole(in, aux, oridim);
    InverseFourierTransform(aux, out);
    out.setXmippOrigin();
}

RealImag2Complex()

void RealImag2Complex	(	const MultidimArray< double > &	real,
		const MultidimArray< double > &	imag,
		MultidimArray< std::complex< double > > &	out
	)

Conversion from real,imag -> complex

Convert real,imag -> complex Fourier transforms 3D. –

Definition at line 114 of file xmipp_fft.cpp.

{
    out.resizeNoCopy(real);
    RealImag2Complex(MULTIDIM_ARRAY(real), MULTIDIM_ARRAY(imag),
                     MULTIDIM_ARRAY(out), MULTIDIM_SIZE(real));
}

ShiftFFT() [1/3]

void ShiftFFT	(	MultidimArray< std::complex< double > > &	v,
		double	xshift
	)

FFT shift 1D

Calculates the Fourier Transform of the shifted real-space vector by phase shifts in Fourier space

Definition at line 311 of file xmipp_fft.cpp.

{
    v.checkDimension(1);
    double dotp, a, b, c, d, ac, bd, ab_cd;
    double xxshift = xshift / (double)XSIZE(v);
    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY1D(v)
    {
        dotp = -2 * PI * ((double)(i) * xxshift);
        a = cos(dotp);
        b = sin(dotp);
        c = DIRECT_A1D_ELEM(v,i).real();
        d = DIRECT_A1D_ELEM(v,i).imag();
        ac = a * c;
        bd = b * d;
        ab_cd = (a + b) * (c + d); // (ab_cd-ac-bd = ad+bc : but needs 4 multiplications)
        DIRECT_A1D_ELEM(v,i) = std::complex<double>(ac - bd, ab_cd - ac - bd);
    }
}

ShiftFFT() [2/3]

void ShiftFFT	(	MultidimArray< std::complex< double > > &	v,
		double	xshift,
		double	yshift
	)

FFT shift 2D

Calculates the Fourier Transform of the shifted real-space vector by phase shifts in Fourier space

Definition at line 331 of file xmipp_fft.cpp.

{
    v.checkDimension(2);
    double dotp, a, b, c, d, ac, bd, ab_cd;
    double xxshift = xshift / (double)XSIZE(v);
    double yyshift = yshift / (double)YSIZE(v);
    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(v)
    {
        dotp = -2 * PI * ((double)(j) * xxshift + (double)(i) * yyshift);
        a = cos(dotp);
        b = sin(dotp);
        c = DIRECT_A2D_ELEM(v,i,j).real();
        d = DIRECT_A2D_ELEM(v,i,j).imag();
        ac = a * c;
        bd = b * d;
        ab_cd = (a + b) * (c + d);
        DIRECT_A2D_ELEM(v,i,j) = std::complex<double>(ac - bd, ab_cd - ac - bd);
    }
}

ShiftFFT() [3/3]

void ShiftFFT	(	MultidimArray< std::complex< double > > &	v,
		double	xshift,
		double	yshift,
		double	zshift
	)

FFT shift 3D

Calculates the Fourier Transform of the shifted real-space vector by phase shifts in Fourier space

Definition at line 352 of file xmipp_fft.cpp.

{
    v.checkDimension(3);
    double dotp, a, b, c, d, ac, bd, ab_cd;
    double xxshift = -2 * PI * xshift / (double)XSIZE(v);
    double yyshift = -2 * PI * yshift / (double)YSIZE(v);
    double zzshift = -2 * PI * zshift / (double)ZSIZE(v);
    for (size_t k=0; k<ZSIZE(v); ++k)
    {
        double zdot=(double)(k) * zzshift;
        for (size_t i=0; i<YSIZE(v); ++i)
        {
            double zydot=zdot+(double)(i) * yyshift;
            for (size_t j=0; j<XSIZE(v); ++j)
            {
                double *ptrv_kij=(double *)&DIRECT_A3D_ELEM(v,k,i,j);
                dotp = (double)(j) * xxshift + zydot;
                //sincos(dotp,&b,&a);
                b = sin(dotp);
                a = cos(dotp);
                c = *ptrv_kij;
                d = *(ptrv_kij+1);
                ac = a * c;
                bd = b * d;
                ab_cd = (a + b) * (c + d);
                *ptrv_kij = ac - bd;
                *(ptrv_kij+1) = ab_cd - ac - bd;
            }
        }
    }
}

Whole2Half()

void Whole2Half	(	const MultidimArray< std::complex< double > > &	in,
		MultidimArray< std::complex< double > > &	out
	)

Conversion from whole -> half

Convert whole -> half of (centro-symmetric) Fourier transforms 1D. –

Definition at line 34 of file xmipp_fft.cpp.

{
    if (in.getDim() == 1)
    {
        // 1D
        int ldim = (int)(XSIZE(in) / 2) + 1;
        out.resize(ldim);
        for (int j = 0; j < ldim; j++)
            out(j) = in(j);
    }
    else if (in.getDim() == 2)
    {
        // 2D
        // This assumes squared images...
        int ldim = (int)(YSIZE(in) / 2) + 1;

        out.initZeros(ldim, XSIZE(in));
        // Fill first column only half
        for (int j = 0; j < ldim; j++)
            dAij(out, 0, j) = dAij(in, 0, j);
        // Fill rest
        for (int i = 1; i < ldim; i++)
            for (size_t j = 0; j < XSIZE(in); j++)
                dAij(out, i, j) = dAij(in, i, j);
    }
    else
        REPORT_ERROR(ERR_MULTIDIM_DIM,"ERROR: Whole2Half only implemented for 1D and 2D multidimArrays");

}

xmipp2CTF()

void xmipp2CTF	(	const MultidimArray< double > &	input,
		MultidimArray< double > &	output
	)

Xmipp image -> Xmipp CTF. The log10 is taken, outliers rejected and the image is reorganized.

Definition at line 472 of file xmipp_fft.cpp.

{
    output = input;
    CenterFFT(output, true);

    // Prepare PSD part
    double min_val = output(0, XSIZE(output) - 1);
    double max_val = min_val;
    bool first = true;
    int Xdim = XSIZE(output);
    int Ydim = YSIZE(output);
    FOR_ALL_ELEMENTS_IN_ARRAY2D(output)
    {
        if ((i < Ydim / 2 && j >= Xdim / 2) || (i >= Ydim / 2 && j < Xdim / 2))
        {
            if (output(i, j) > XMIPP_EQUAL_ACCURACY &&
                (output(i, j) < min_val || first))
                min_val = output(i, j);
            if (output(i, j) > XMIPP_EQUAL_ACCURACY &&
                (output(i, j) > max_val || first))
            {
                max_val = output(i, j);
                first = false;
            }
        }
    }
    MultidimArray<double> left(YSIZE(output), XSIZE(output));
    min_val = 10 * log10(min_val);
    FOR_ALL_ELEMENTS_IN_ARRAY2D(output)
    {
        if ((i < Ydim / 2 && j >= Xdim / 2) || (i >= Ydim / 2 && j < Xdim / 2))
        {
            if (output(i, j) > XMIPP_EQUAL_ACCURACY)
                left(i, j) = 10 * log10(output(i, j));
            else
                left(i, j) = min_val;
        }
    }
    reject_outliers(left);

    // Join both parts
    FOR_ALL_ELEMENTS_IN_ARRAY2D(output)
    if ((i < Ydim / 2 && j >= Xdim / 2) || (i >= Ydim / 2 && j < Xdim / 2))
        output(i, j) = left(i, j);
    else
        output(i, j) = ABS(output(i, j));
}

xmipp2PSD()

template<typename T >

void xmipp2PSD	(	const MultidimArray< T > &	input,
		MultidimArray< T > &	output,
		bool	takeLog = true
	)

Xmipp image -> Xmipp PSD. The log10 is taken, outliers rejected and the image is reorganized.

Definition at line 443 of file xmipp_fft.cpp.

{
    output = input;
    CenterFFT(output, true);
    if (takeLog) {
        auto min_val = std::numeric_limits<T>::max();
        FOR_ALL_ELEMENTS_IN_ARRAY2D(output)
        {
            auto pixval=A2D_ELEM(output,i,j);
            if (pixval > 0 && pixval < min_val)
                min_val = pixval;
        }
        min_val = 10 * log10(min_val);
        FOR_ALL_ELEMENTS_IN_ARRAY2D(output)
        {
            auto pixval=A2D_ELEM(output,i,j);
            if (pixval > 0)
                A2D_ELEM(output,i,j) = 10 * log10(pixval);
            else
                A2D_ELEM(output,i,j) = min_val;
        }
    }
    reject_outliers(output);
}