xmipp_fft.h

Go to the documentation of this file.

/***************************************************************************
*
* Authors:     Carlos Oscar S. Sorzano (coss@cnb.csic.es)
*
* Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
* 02111-1307  USA
*
*  All comments concerning this program package may be sent to the
*  e-mail address 'xmipp@cnb.csic.es'
***************************************************************************/

#ifndef CORE_FFT_H
#define CORE_FFT_H

#include <complex>

#include "multidim_array.h"
#include "xmipp_funcs.h"

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

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

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

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

#define DIGFREQ2FFT_IDX_DOUBLE(freq, size, idx) { \
        (idx) = ((size) * (freq)); if ((idx) < 0) (idx) += (size); }

template <typename T>
void FFT_idx2digfreq(T& v, const Matrix1D< int >& idx, Matrix1D< double >& freq)
{
    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));
    }
}

template <typename T>
void digfreq2FFT_idx(T& v, const Matrix1D< double >& freq, Matrix1D< int >& idx)
{
    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));
}

inline void digfreq2contfreq(const Matrix1D< double >& digfreq,
                             Matrix1D< double >& contfreq,
                             double pixel_size)
{
    contfreq.resizeNoCopy(digfreq);
    FOR_ALL_ELEMENTS_IN_MATRIX1D(digfreq)
    VEC_ELEM(contfreq,i) = VEC_ELEM(digfreq,i) / pixel_size;
}

inline void contfreq2digfreq(const Matrix1D< double >& contfreq,
                             Matrix1D< double >& digfreq,
                             double pixel_size)
{
    digfreq.resizeNoCopy(contfreq);
    FOR_ALL_ELEMENTS_IN_MATRIX1D(contfreq)
    VEC_ELEM(digfreq,i) = VEC_ELEM(contfreq,i) * pixel_size;
}

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

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


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

void centerFFT2(MultidimArray<double> &v);


template <typename T>
void CenterFFT(MultidimArray< T >& v, bool forward)
{
    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");
    }
}

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

#endif