FFTW Fourier transforms

Collaboration diagram for FFTW Fourier transforms:

Classes
class	FourierTransformer
class	CorrelationAux

Macros
#define	POWER_SPECTRUM   0
#define	AMPLITUDE_SPECTRUM   1

Functions
void	FFT_magnitude (const MultidimArray< std::complex< double > > &v, MultidimArray< double > &mag)
void	FFT_phase (const MultidimArray< std::complex< double > > &v, MultidimArray< double > &phase)
void	radial_magnitude (const MultidimArray< double > &v, MultidimArray< std::complex< double > > &V, MultidimArray< double > &radialMagnitude)
template<typename T >
void	auto_correlation_vector (const MultidimArray< T > &Img, MultidimArray< double > &R)
template<typename T >
void	correlation_vector (const MultidimArray< T > &m1, const MultidimArray< T > &m2, MultidimArray< double > &R)
void	fast_correlation_vector (const MultidimArray< std::complex< double > > &FFT1, const MultidimArray< std::complex< double > > &FFT2, MultidimArray< double > &R, FourierTransformer &transformer)
template<class T >
void	correlation_vector_no_Fourier (const MultidimArray< T > &v1, const MultidimArray< T > &v2, MultidimArray< T > &result)
void	correlation_matrix (const MultidimArray< double > &m1, const MultidimArray< double > &m2, MultidimArray< double > &R, CorrelationAux &aux, bool center=true)
void	correlation_matrix (const MultidimArray< std::complex< double > > &FFT1, const MultidimArray< double > &m2, MultidimArray< double > &R, CorrelationAux &aux, bool center=true)
void	correlation_matrix (const MultidimArray< std::complex< double > > &FFT1, const MultidimArray< std::complex< double > > &FFT2, MultidimArray< double > &R, CorrelationAux &aux, bool center=true)
template<typename T >
void	auto_correlation_matrix (const MultidimArray< T > &Img, MultidimArray< double > &R)
void	auto_correlation_matrix (const MultidimArray< double > &Img, MultidimArray< double > &R, CorrelationAux &aux)
void	convolutionFFTStack (const MultidimArray< double > &img, const MultidimArray< double > &kernel, MultidimArray< double > &result)
void	convolutionFFT (MultidimArray< double > &img, MultidimArray< double > &kernel, MultidimArray< double > &result)
void	frc_dpr (MultidimArray< double > &m1, MultidimArray< double > &m2, double sampling_rate, MultidimArray< double > &freq, MultidimArray< double > &frc, MultidimArray< double > &frc_noise, MultidimArray< double > &dpr, MultidimArray< double > &error_l2, bool skipdpr=false, bool doRfactor=false, double minFreq=-1, double maxFreq=0.5, double *rFactor=NULL)
void	scaleToSizeFourier (int Zdim, int Ydim, int Xdim, MultidimArray< double > &mdaIn, MultidimArray< double > &mdaOut, int nThreads=1)
template<typename T >
void	scaleToSizeFourier (MultidimArray< T > &mdaIn, MultidimArray< T > &mdaOut, MultidimArray< std::complex< T > > &inFourier, MultidimArray< std::complex< T > > &outFourier)
void	selfScaleToSizeFourier (int Zdim, int Ydim, int Xdim, MultidimArray< double > &mda, int nthreads=1)
void	selfScaleToSizeFourier (int Ydim, int Xdim, MultidimArray< double > &mda, int nthreads=1)
void	selfScaleToSizeFourier (int Zdim, int Ydim, int Xdim, MultidimArrayGeneric &mda, int nthreads=1)
void	selfScaleToSizeFourier (int Ydim, int Xdim, MultidimArrayGeneric &mda, int nthreads=1)
void	getSpectrum (MultidimArray< double > &Min, MultidimArray< double > &spectrum, int spectrum_type=AMPLITUDE_SPECTRUM)
void	divideBySpectrum (MultidimArray< double > &Min, MultidimArray< double > &spectrum, bool leave_origin_intact=false)
void	multiplyBySpectrum (MultidimArray< double > &Min, MultidimArray< double > &spectrum, bool leave_origin_intact=false)
void	whitenSpectrum (MultidimArray< double > &Min, MultidimArray< double > &Mout, int spectrum_type=AMPLITUDE_SPECTRUM, bool leave_origin_intact=false)
void	adaptSpectrum (MultidimArray< double > &Min, MultidimArray< double > &Mout, const MultidimArray< double > &spectrum_ref, int spectrum_type=AMPLITUDE_SPECTRUM, bool leave_origin_intact=false)
void	randomizePhases (MultidimArray< double > &Min, double wRandom)

Detailed Description

Macro Definition Documentation

AMPLITUDE_SPECTRUM

#define AMPLITUDE_SPECTRUM   1

Definition at line 669 of file xmipp_fftw.h.

POWER_SPECTRUM

#define POWER_SPECTRUM   0

Definition at line 668 of file xmipp_fftw.h.

Function Documentation

adaptSpectrum()

void adaptSpectrum	(	MultidimArray< double > &	Min,
		MultidimArray< double > &	Mout,
		const MultidimArray< double > &	spectrum_ref,
		int	spectrum_type = AMPLITUDE_SPECTRUM,
		bool	leave_origin_intact = false
	)

Adapts Min to have the same spectrum as spectrum_ref

If only_amplitudes==true, the amplitude rather than the power spectrum will be equalized

Definition at line 849 of file xmipp_fftw.cpp.

{

    Min.checkDimension(3);

    MultidimArray<double> spectrum;
    getSpectrum(Min,spectrum,spectrum_type);
    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY1D(spectrum)
    {
        dAi(spectrum, i) = (dAi(spectrum, i) > 0.) ? dAi(spectrum_ref,i)/ dAi(spectrum, i) : 1.;
    }
    Mout=Min;
    multiplyBySpectrum(Mout,spectrum,leave_origin_intact);

}

auto_correlation_matrix() [1/2]

template<typename T >

void auto_correlation_matrix	(	const MultidimArray< T > &	Img,
		MultidimArray< double > &	R
	)

Autocorrelation function of an image

Fast calcuation of the autocorrelation matrix of a given one using Fast Fourier Transform. (Using the correlation theorem)

Definition at line 594 of file xmipp_fftw.h.

{
    // Compute the Fourier Transform
    MultidimArray< std::complex< double > > FFT1;
    FourierTransformer transformer1;
    R=Img;
    transformer1.FourierTransform(R, FFT1, false);

    // Multiply FFT1 * FFT1'
    double dSize=MULTIDIM_SIZE(Img);
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(FFT1)
    {
        double *ptr=(double*)&DIRECT_MULTIDIM_ELEM(FFT1,n);
        double &realPart=*ptr;
        double &imagPart=*(ptr+1);
        realPart=dSize*(realPart*realPart+imagPart*imagPart);
        imagPart=0;
    }

    // Invert the product, in order to obtain the correlation image
    transformer1.inverseFourierTransform();

    // Center the resulting image to obtain a centered autocorrelation
    CenterFFT(R, true);
}

auto_correlation_matrix() [2/2]

void auto_correlation_matrix	(	const MultidimArray< double > &	Img,
		MultidimArray< double > &	R,
		CorrelationAux &	aux
	)

Fast autocorrelation matrix

Definition at line 957 of file xmipp_fftw.cpp.

{
    // Compute the Fourier Transform
    R=Img;
    aux.transformer1.FourierTransform(R, aux.FFT1, false);

    // Multiply FFT1 * FFT1'
    double dSize=MULTIDIM_SIZE(Img);
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(aux.FFT1)
    {
        double *ptr=(double*)&DIRECT_MULTIDIM_ELEM(aux.FFT1,n);
        double &realPart=*ptr;
        double &imagPart=*(ptr+1);
        realPart=dSize*(realPart*realPart+imagPart*imagPart);
        imagPart=0;
    }

    // Invert the product, in order to obtain the correlation image
    aux.transformer1.inverseFourierTransform();

    // Center the resulting image to obtain a centered autocorrelation
    CenterFFT(R, true);
}

auto_correlation_vector()

template<typename T >

void auto_correlation_vector	(	const MultidimArray< T > &	Img,
		MultidimArray< double > &	R
	)

Autocorrelation function of a Xmipp vector

Fast calculation of the autocorrelation vector of a given one using Fast Fourier Transform. (Using the correlation theorem)

Definition at line 452 of file xmipp_fftw.h.

{
    Img.checkDimension(1);

    // Compute the Fourier Transform
    MultidimArray< std::complex< double > > FFT1;
    FourierTransformer transformer1;
    R=Img;
    transformer1.FourierTransform(R, FFT1, false);

    // Multiply FFT1 * FFT1'
    double dSize=XSIZE(Img);
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(FFT1)
    {
        double *ptr=(double*)&DIRECT_MULTIDIM_ELEM(FFT1,n);
        double &realPart=*ptr;
        double &imagPart=*(ptr+1);
        realPart=dSize*(realPart*realPart+imagPart*imagPart);
        imagPart=0;
    }

    // Invert the product, in order to obtain the correlation image
    transformer1.inverseFourierTransform();

    // Center the resulting image to obtain a centered autocorrelation
    CenterFFT(R, true);
}

convolutionFFT()

void convolutionFFT	(	MultidimArray< double > &	img,
		MultidimArray< double > &	kernel,
		MultidimArray< double > &	result
	)

Definition at line 457 of file xmipp_fftw.cpp.

{
    FourierTransformer transformer1, transformer2;
    MultidimArray< std::complex<double> > FFT1, FFT2;
    transformer1.FourierTransform(img, FFT1, false);
    result=kernel;
    transformer2.FourierTransform(result, FFT2, false);

    // Multiply FFT1 * FFT2'
    double dSize=MULTIDIM_SIZE(img);
    double a, b, c, d; // a+bi, c+di
    double *ptrFFT2=(double*)MULTIDIM_ARRAY(FFT2);
    double *ptrFFT1=(double*)MULTIDIM_ARRAY(FFT1);
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(FFT1)
    {
        a=*ptrFFT1++;
        b=*ptrFFT1++;
        c=(*ptrFFT2)*dSize;
        d=(*(ptrFFT2+1))*dSize;
        *ptrFFT2++ = a*c-b*d;
        *ptrFFT2++ = b*c+a*d;
    }

    // Invert the product, in order to obtain the correlation image
    transformer2.inverseFourierTransform();

    // Center the resulting image to compensate for phase shift
    CenterFFT(result, false);
}

convolutionFFTStack()

void convolutionFFTStack	(	const MultidimArray< double > &	img,
		const MultidimArray< double > &	kernel,
		MultidimArray< double > &	result
	)

Definition at line 429 of file xmipp_fftw.cpp.

{
    if (&result != &img)
        result = img;

    MultidimArray<double> imgTemp;
    MultidimArray< std::complex< double> > FFTIm, FFTK;
    FourierTransformer transformer1(FFTW_BACKWARD), transformer2(FFTW_BACKWARD);

    transformer2.FourierTransform((MultidimArray<double> &)kernel, FFTK, false);

    for (size_t n = 0; n < ZSIZE(result); n++)
    {
        imgTemp.aliasSlice(result, n);
        transformer1.FourierTransform(imgTemp, FFTIm, false);

        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(FFTIm)
        DIRECT_MULTIDIM_ELEM(FFTIm,n) *= DIRECT_MULTIDIM_ELEM(FFTK,n);

        transformer1.inverseFourierTransform();

        CenterFFT(imgTemp, false);
    }

}

correlation_matrix() [1/3]

void correlation_matrix	(	const MultidimArray< double > &	m1,
		const MultidimArray< double > &	m2,
		MultidimArray< double > &	R,
		CorrelationAux &	aux,
		bool	center = true
	)

Correlation of two nD images

Fast calcuation of the correlation matrix on two matrices using Fast Fourier Transform. (Using the correlation theorem). The output matrix must be already resized

Definition at line 869 of file xmipp_fftw.cpp.

{
    aux.transformer1.FourierTransform((MultidimArray<double> &)m1, aux.FFT1, false);
    correlation_matrix(aux.FFT1,m2,R,aux,center);
}

correlation_matrix() [2/3]

void correlation_matrix	(	const MultidimArray< std::complex< double > > &	FFT1,
		const MultidimArray< double > &	m2,
		MultidimArray< double > &	R,
		CorrelationAux &	aux,
		bool	center = true
	)

Definition at line 898 of file xmipp_fftw.cpp.

{
    R=m2;
    aux.transformer2.FourierTransform(R, aux.FFT2, false);
    correlationInFourier(FF1,aux.FFT2,MULTIDIM_SIZE(R));
    aux.transformer2.inverseFourierTransform();
    if (center)
        CenterFFT(R, true);
}

correlation_matrix() [3/3]

void correlation_matrix	(	const MultidimArray< std::complex< double > > &	FFT1,
		const MultidimArray< std::complex< double > > &	FFT2,
		MultidimArray< double > &	R,
		CorrelationAux &	aux,
		bool	center = true
	)

Correlation matrix. R must already be with the right size.

Definition at line 912 of file xmipp_fftw.cpp.

{
    aux.transformer2.setReal(R);
    aux.transformer2.setFourier(FFT2);
    correlationInFourier(FFT1,aux.transformer2.fFourier,MULTIDIM_SIZE(R));
    aux.transformer2.inverseFourierTransform();
    if (center)
        CenterFFT(R, true);
}

correlation_vector()

template<typename T >

void correlation_vector	(	const MultidimArray< T > &	m1,
		const MultidimArray< T > &	m2,
		MultidimArray< double > &	R
	)

Autocorrelation function of a Xmipp vector

Fast calcuation of the correlation matrix on two matrices using Fast Fourier Transform. (Using the correlation theorem). The output matrix must be already resized

Definition at line 488 of file xmipp_fftw.h.

{
    m1.checkDimension(2);
    m2.checkDimension(2);

    // Compute the Fourier Transforms
    MultidimArray< std::complex< double > > FFT1, FFT2;
    FourierTransformer transformer1, transformer2;
    R=m1;
    transformer1.FourierTransform(R, FFT1, false);
    transformer2.FourierTransform((MultidimArray<T> &)m2, FFT2, false);

    // Multiply FFT1 * FFT2'
    double dSize=XSIZE(m1);
    FOR_ALL_ELEMENTS_IN_ARRAY1D(FFT1)
    FFT1(i) *= dSize * conj(FFT2(i));

    // Invert the product, in order to obtain the correlation image
    transformer1.inverseFourierTransform();

    // Center the resulting image to obtain a centered autocorrelation
    CenterFFT(R, true);
}

correlation_vector_no_Fourier()

template<class T >

void correlation_vector_no_Fourier	(	const MultidimArray< T > &	v1,
		const MultidimArray< T > &	v2,
		MultidimArray< T > &	result
	)

Compute the correlation vector without using Fourier.

The results are the same as the previous ones but this function is threadsafe while the previous one is not.

It is assumed that the two vectors v1, and v2 are of the same size.

Definition at line 534 of file xmipp_fftw.h.

{
    v1.checkDimension(1);
    v2.checkDimension(1);

    result.initZeros(v1);
    result.setXmippOrigin();
    int N=XSIZE(v1)-1;
    FOR_ALL_ELEMENTS_IN_ARRAY1D(result)
    for (int k=0; k<XSIZE(v1); ++k)
        A1D_ELEM(result,i)+=DIRECT_A1D_ELEM(v1,intWRAP(k+i,0,N))*
                            DIRECT_A1D_ELEM(v2,k);
    STARTINGX(result)=0;
}

divideBySpectrum()

void divideBySpectrum	(	MultidimArray< double > &	Min,
		MultidimArray< double > &	spectrum,
		bool	leave_origin_intact = false
	)

Divide the input map in Fourier-space by the spectrum provided.

If leave_origin_intact==true, the origin pixel will remain untouched

Definition at line 788 of file xmipp_fftw.cpp.

{

    Min.checkDimension(3);

    MultidimArray<double> div_spec(spectrum);
    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY1D(spectrum)
    {
        if (ABS(dAi(spectrum,i)) > 0.)
            dAi(div_spec,i) = 1./dAi(spectrum,i);
        else
            dAi(div_spec,i) = 1.;
    }
    multiplyBySpectrum(Min,div_spec,leave_origin_intact);
}

fast_correlation_vector()

void fast_correlation_vector	(	const MultidimArray< std::complex< double > > &	FFT1,
		const MultidimArray< std::complex< double > > &	FFT2,
		MultidimArray< double > &	R,
		FourierTransformer &	transformer
	)

Autocorrelation function of a Xmipp vector

Same as correlation_vector, but the Fourier transforms have already been computed and the transformer is reused. The transformer internal variables must have already been resized.

Definition at line 926 of file xmipp_fftw.cpp.

{
    transformer.setFourier(FFT1);

    // Multiply FFT1 * FFT2'
    double a, b, c, d; // a+bi, c+di
    double *ptrFFT2=(double*)MULTIDIM_ARRAY(FFT2);
    double *ptrFFT1=(double*)&(A1D_ELEM(transformer.fFourier,0));
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(FFT1)
    {
        a=*ptrFFT1;
        b=*(ptrFFT1+1);
        c=(*ptrFFT2++);
        d=(*ptrFFT2++)*(-1);
        *ptrFFT1++ = a*c-b*d;
        *ptrFFT1++ = b*c+a*d;
    }

    // Invert the product, in order to obtain the correlation image
    transformer.inverseFourierTransform();

    // Center the resulting image to obtain a centered autocorrelation
    R=*transformer.fReal;
    CenterFFT(R, true);
    R.setXmippOrigin();
}

FFT_magnitude()

void FFT_magnitude	(	const MultidimArray< std::complex< double > > &	v,
		MultidimArray< double > &	mag
	)

FFT Magnitude 1D

Definition at line 393 of file xmipp_fftw.cpp.

{
    mag.resizeNoCopy(v);
    double * ptrv=(double *)MULTIDIM_ARRAY(v);
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(v)
    {
        double re=*ptrv;
        double im=*(ptrv+1);
        DIRECT_MULTIDIM_ELEM(mag, n) = sqrt(re*re+im*im);
        ptrv+=2;
    }
}

FFT_phase()

void FFT_phase	(	const MultidimArray< std::complex< double > > &	v,
		MultidimArray< double > &	phase
	)

FFT Phase 1D

Definition at line 408 of file xmipp_fftw.cpp.

{
    phase.resizeNoCopy(v);
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(v)
    DIRECT_MULTIDIM_ELEM(phase, n) = atan2(DIRECT_MULTIDIM_ELEM(v,n).imag(), DIRECT_MULTIDIM_ELEM(v, n).real());
}

frc_dpr()

void frc_dpr	(	MultidimArray< double > &	m1,
		MultidimArray< double > &	m2,
		double	sampling_rate,
		MultidimArray< double > &	freq,
		MultidimArray< double > &	frc,
		MultidimArray< double > &	frc_noise,
		MultidimArray< double > &	dpr,
		MultidimArray< double > &	error_l2,
		bool	skipdpr = false,
		bool	doRfactor = false,
		double	minFreq = -1,
		double	maxFreq = 0.5,
		double *	rFactor = NULL
	)

Fourier-Ring-Correlation between two multidimArrays using FFT

Definition at line 491 of file xmipp_fftw.cpp.

{
    if (!m1.sameShape(m2))
        REPORT_ERROR(ERR_MULTIDIM_SIZE,"MultidimArrays have different shapes!");

    int m1sizeX = XSIZE(m1), m1sizeY = YSIZE(m1), m1sizeZ = ZSIZE(m1);

    MultidimArray< std::complex< double > > FT1;
    FourierTransformer transformer1(FFTW_BACKWARD);
    transformer1.FourierTransform(m1, FT1, false);

//#define DEBUG
#ifdef DEBUG
     std::cerr << "FT1" << FT1 << std::endl; 
#endif
#undef DEBUG


    m1.clear(); // Free memory

    MultidimArray< std::complex< double > > FT2;
    FourierTransformer transformer2(FFTW_BACKWARD);
    transformer2.FourierTransform(m2, FT2, false);
    m2.clear(); // Free memory

    MultidimArray< int > radial_count(m1sizeX/2+1);
    MultidimArray<double> num, den1, den2, den_dpr;
    Matrix1D<double> f(3);
    //to calculate r-factor
    double rFactorNumerator=0;
    double rFactorDenominator=0;


    num.initZeros(radial_count);
    den1.initZeros(radial_count);
    den2.initZeros(radial_count);

    //dpr calculation takes for ever in large volumes
    //since atan2 is called many times
    //until atan2 is changed by a table let us make dpr an option
    if (dodpr)
    {
        dpr.initZeros(radial_count);
        den_dpr.initZeros(radial_count);
    }
    freq.initZeros(radial_count);
    frc.initZeros(radial_count);
    frc_noise.initZeros(radial_count);
    error_l2.initZeros(radial_count);
#ifdef SAVE_REAL_PART

    std::vector<double> *realPart=
        new std::vector<double>[XSIZE(radial_count)];
#endif

    int sizeZ_2 = m1sizeZ/2;
    if (sizeZ_2==0)
        sizeZ_2=1;
    double isizeZ = 1.0/m1sizeZ;
    int sizeY_2 = m1sizeY/2;
    double iysize = 1.0/m1sizeY;
    int sizeX_2 = m1sizeX/2;
    double ixsize = 1.0/m1sizeX;
    double R;

    int ZdimFT1=(int)ZSIZE(FT1);
    int YdimFT1=(int)YSIZE(FT1);
    int XdimFT1=(int)XSIZE(FT1);

    double maxFreq_2 =0.;
    maxFreq_2 = maxFreq * maxFreq;
    for (int k=0; k<ZdimFT1; k++)
    {
        FFT_IDX2DIGFREQ_FAST(k,m1sizeZ,sizeZ_2,isizeZ,ZZ(f));
        double fz2=ZZ(f)*ZZ(f);
        for (int i=0; i<YdimFT1; i++)
        {
            FFT_IDX2DIGFREQ_FAST(i,YSIZE(m1),sizeY_2, iysize, YY(f));
            double fz2_fy2=fz2 + YY(f)*YY(f);
            for (int j=0; j<XdimFT1; j++)
            {
                FFT_IDX2DIGFREQ_FAST(j,m1sizeX, sizeX_2, ixsize, XX(f));

                double R2 =fz2_fy2 + XX(f)*XX(f);
                if (R2>maxFreq_2)
                    continue;

                R = sqrt(R2);
                int idx = (int)round(R * m1sizeX);
                std::complex<double> &z1 = dAkij(FT1, k, i, j);
                std::complex<double> &z2 = dAkij(FT2, k, i, j);
                double absz1 = abs(z1);
                double absz2 = abs(z2);
                dAi(num,idx) += real(conj(z1) * z2);
                dAi(den1,idx) += absz1*absz1;
                dAi(den2,idx) += absz2*absz2;
                dAi(error_l2,idx) += abs(z1-z2);
                //for calculating r-factor
                if ( R > minFreq && R < maxFreq)
                {
                    rFactorNumerator += fabs(absz1 - absz2);
                    rFactorDenominator += absz1;
                }

                if (dodpr) //this takes to long for a huge volume
                {
                    double phaseDiff=atan2(z1.imag(),z1.real()) - atan2(z2.imag(),z2.real());
                    phaseDiff = RAD2DEG(phaseDiff);
                    phaseDiff = realWRAP(phaseDiff,-180, 180);
                    dAi(dpr,idx) += ((absz1+absz2)*phaseDiff*phaseDiff);
                    dAi(den_dpr,idx) += (absz1+absz2);
#ifdef SAVE_REAL_PART

                    realPart[idx].push_back(z1.real());
#endif

                }
                dAi(radial_count,idx)++;
            }
        }
    }
    //to calculate r-factor
    if (doRfactor)
    {
        //std::cout << "rFactorNumerator = " << rFactorNumerator << std::endl;
        //std::cout << "rFactorDenominator = " << rFactorDenominator << std::endl;
        *rFactor = rFactorNumerator / rFactorDenominator;
        
        //std::cout << "R-factor = " << rFactorValue << std::endl;
        //*rFactor = rFactorValue;
        //std::cout << "R-rFactor = " << *rFactor << std::endl;
    }

    FOR_ALL_ELEMENTS_IN_ARRAY1D(freq)
    {
        dAi(freq,i) = (double) i / (m1sizeX * sampling_rate);
        dAi(frc,i) = dAi(num,i)/sqrt(dAi(den1,i)*dAi(den2,i));
        dAi(frc_noise,i) = 2 / sqrt((double) dAi(radial_count,i));
        dAi(error_l2,i) /= dAi(radial_count,i);

        if (dodpr)
            dAi(dpr,i) = sqrt(dAi(dpr,i) / dAi(den_dpr,i));
#ifdef SAVE_REAL_PART

        std::ofstream fhOut;
        fhOut.open(((std::string)"PPP_RealPart_"+integerToString(i)+".txt").
                   c_str());
        for (int j=0; j<realPart[i].size(); j++)
            fhOut << realPart[i][j] << std::endl;
        fhOut.close();
#endif

    }
}

getSpectrum()

void getSpectrum	(	MultidimArray< double > &	Min,
		MultidimArray< double > &	spectrum,
		int	spectrum_type = AMPLITUDE_SPECTRUM
	)

Get the amplitude or power spectrum of the map in Fourier space.

i.e. the radial average of the (squared) amplitudes of all Fourier components

Definition at line 752 of file xmipp_fftw.cpp.

{
    MultidimArray<std::complex<double> > Faux;
    int xsize = XSIZE(Min);
    Matrix1D<double> f(3);
    MultidimArray<double> count(xsize);
    FourierTransformer transformer;

    spectrum.initZeros(xsize);
    count.initZeros();
    transformer.FourierTransform(Min, Faux, false);
    if (ZSIZE(Faux)==1)
        ZZ(f)=0;
    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(Faux)
    {
        FFT_IDX2DIGFREQ(j,xsize,XX(f));
        FFT_IDX2DIGFREQ(i,YSIZE(Faux),YY(f));
        if (ZSIZE(Faux)>1)
            FFT_IDX2DIGFREQ(k,ZSIZE(Faux),ZZ(f));
        double R=f.module();
        //if (R>0.5) continue;
        int idx = round(R*xsize);
        double F=abs(dAkij(Faux, k, i, j));
        if (spectrum_type == AMPLITUDE_SPECTRUM)
            A1D_ELEM(spectrum,idx) += F;
        else
            A1D_ELEM(spectrum,idx) += F*F;
        A1D_ELEM(count,idx) += 1.;
    }
    for (int i = 0; i < xsize; i++)
        if (A1D_ELEM(count,i) > 0.)
            A1D_ELEM(spectrum,i) /= A1D_ELEM(count,i);
}

multiplyBySpectrum()

void multiplyBySpectrum	(	MultidimArray< double > &	Min,
		MultidimArray< double > &	spectrum,
		bool	leave_origin_intact = false
	)

Multiply the input map in Fourier-space by the spectrum provided.

If leave_origin_intact==true, the origin pixel will remain untouched

Definition at line 806 of file xmipp_fftw.cpp.

{
    Min.checkDimension(3);

    MultidimArray<std::complex<double> > Faux;
    Matrix1D<double> f(3);
    MultidimArray<double> lspectrum;
    FourierTransformer transformer;
    double dim3 = XSIZE(Min)*YSIZE(Min)*ZSIZE(Min);

    transformer.FourierTransform(Min, Faux, false);
    lspectrum=spectrum;
    if (leave_origin_intact)
        lspectrum(0)=1.;
    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(Faux)
    {
        FFT_IDX2DIGFREQ(j,XSIZE(Min), XX(f));
        FFT_IDX2DIGFREQ(i,YSIZE(Faux),YY(f));
        FFT_IDX2DIGFREQ(k,ZSIZE(Faux),ZZ(f));
        double R=f.module();
        //if (R > 0.5) continue;
        int idx=ROUND(R*XSIZE(Min));
        dAkij(Faux, k, i, j) *=  lspectrum(idx) * dim3;
    }
    transformer.inverseFourierTransform();
}

radial_magnitude()

void radial_magnitude	(	const MultidimArray< double > &	v,
		MultidimArray< std::complex< double > > &	V,
		MultidimArray< double > &	radialMagnitude
	)

FFT Radial average. This function computes the Fourier transform (in a relatively inefficient way), the amplitude and radial average of the input image. This is primarily meant for debugging. It also returns the Fourier transform just in case it might be useful outside.

Definition at line 416 of file xmipp_fftw.cpp.

{
    FourierTransform(v,V);
    MultidimArray<double> Vmag;
    FFT_magnitude(V,Vmag);
    CenterFFT(Vmag,true);
    Vmag.setXmippOrigin();
    MultidimArray<int> radialCount;
    Matrix1D<int> center(3);
    radialAverage(Vmag,center,radialMagnitude,radialCount);
}

randomizePhases()

void randomizePhases	(	MultidimArray< double > &	Min,
		double	wRandom
	)

Randomize phases beyond a certain frequency

Definition at line 981 of file xmipp_fftw.cpp.

{
    FourierTransformer transformer;
    MultidimArray< std::complex<double> > F;
    transformer.FourierTransform(Min,F,false);

    Matrix1D<double> f(3);
    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(F)
    {
        FFT_IDX2DIGFREQ(j,XSIZE(Min), XX(f));
        FFT_IDX2DIGFREQ(i,YSIZE(F),YY(f));
        FFT_IDX2DIGFREQ(k,ZSIZE(F),ZZ(f));
        double w=f.module();
        if (w > wRandom)
        {
            double alpha=rnd_unif(0,2*PI);
            double c,s;
            //sincos(alpha,&s,&c);
            s = sin(alpha);
            c = cos(alpha);

            DIRECT_A3D_ELEM(F,k,i,j)*=std::complex<double>(s,c);
        }
    }
    transformer.inverseFourierTransform();
}

scaleToSizeFourier() [1/2]

void scaleToSizeFourier	(	int	Zdim,
		int	Ydim,
		int	Xdim,
		MultidimArray< double > &	mdaIn,
		MultidimArray< double > &	mdaOut,
		int	nThreads = 1
	)

Scale matrix using Fourier transform

Ydim and Xdim define the output size, mda is the MultidimArray to scale

Definition at line 658 of file xmipp_fftw.cpp.

{
    //Mmem = *this
    //memory for fourier transform output
    MultidimArray<std::complex<double> > MmemFourier;
    // Perform the Fourier transform
    FourierTransformer transformerM;
    transformerM.setThreadsNumber(nThreads);
    transformerM.FourierTransform(mdaIn, MmemFourier, false);

    // Create space for the downsampled image and its Fourier transform
    mdaOut.resizeNoCopy(Zdim, Ydim, Xdim);
    MultidimArray<std::complex<double> > MpmemFourier;
    FourierTransformer transformerMp;
    transformerMp.setReal(mdaOut);
    transformerMp.getFourierAlias(MpmemFourier);

    scaleToSizeFourier(mdaIn, mdaOut, MmemFourier, MpmemFourier);

    // Transform data
    transformerMp.inverseFourierTransform();
}

scaleToSizeFourier() [2/2]

template<typename T >

void scaleToSizeFourier	(	MultidimArray< T > &	mdaIn,
		MultidimArray< T > &	mdaOut,
		MultidimArray< std::complex< T > > &	inFourier,
		MultidimArray< std::complex< T > > &	outFourier
	)

Scale Fourier transform mdaIn and mdaOut define the sizes of the respective image, inFourier is transform to be scaled

Definition at line 685 of file xmipp_fftw.cpp.

                                                                                             {
    size_t xsize = std::min(XSIZE(inFourier),XSIZE(outFourier))*sizeof(std::complex<T>);
    size_t yhalf = std::min(std::min(YSIZE(mdaIn)/2,YSIZE(mdaOut)/2),YSIZE(mdaIn)-1);
    size_t zhalf = std::min(std::min(ZSIZE(mdaIn)/2,ZSIZE(mdaOut)/2),ZSIZE(mdaIn)-1);

    size_t kp0=0;
    size_t kpF=zhalf;
    size_t ip0=0;
    size_t ipF=yhalf;
    size_t km0=ZSIZE(inFourier)>=ZSIZE(outFourier)?(kpF+1):(ZSIZE(outFourier)-(ZSIZE(inFourier)-(zhalf+1)));
    size_t kmF=ZSIZE(outFourier)-1;
    size_t im0=YSIZE(inFourier)>=YSIZE(outFourier)?(ipF+1):(YSIZE(outFourier)-(YSIZE(inFourier)-(yhalf+1)));
    size_t imF=YSIZE(outFourier)-1;

    //Init with zero
    outFourier.initZeros();

    for (size_t k = kp0; k<=kpF; ++k)
    {
        for (size_t i=ip0; i<=ipF; ++i)
            memcpy(&dAkij(outFourier,k,i,0),&dAkij(inFourier,k,i,0),xsize);
        for (size_t i=im0; i<=imF; ++i) {
            size_t ip = i + YSIZE(inFourier)-YSIZE(outFourier) ;
            memcpy(&dAkij(outFourier,k,i,0),&dAkij(inFourier,k,ip,0),xsize);
        }
    }
    for (size_t k = km0; k<=kmF; ++k) {
        size_t kp = k + ZSIZE(inFourier)-ZSIZE(outFourier) ;
        for (size_t i=ip0; i<=ipF; ++i)
            memcpy(&dAkij(outFourier,k,i,0),&dAkij(inFourier,kp,i,0),xsize);
        for (size_t i=im0; i<=imF; ++i) {
            size_t ip = i + YSIZE(inFourier)-YSIZE(outFourier) ;
            memcpy(&dAkij(outFourier,k,i,0),&dAkij(inFourier,kp,ip,0),xsize);
        }
    }
}

selfScaleToSizeFourier() [1/4]

void selfScaleToSizeFourier	(	int	Zdim,
		int	Ydim,
		int	Xdim,
		MultidimArray< double > &	mda,
		int	nthreads = 1
	)

Definition at line 723 of file xmipp_fftw.cpp.

{
  MultidimArray<double> aux;
  scaleToSizeFourier(Zdim, Ydim, Xdim, mda, aux, nThreads);
  mda=aux;
}

selfScaleToSizeFourier() [2/4]

void selfScaleToSizeFourier	(	int	Ydim,
		int	Xdim,
		MultidimArray< double > &	mda,
		int	nthreads = 1
	)

Definition at line 731 of file xmipp_fftw.cpp.

{
    selfScaleToSizeFourier(1, Ydim, Xdim, Mpmem, nThreads);
}

selfScaleToSizeFourier() [3/4]

void selfScaleToSizeFourier	(	int	Zdim,
		int	Ydim,
		int	Xdim,
		MultidimArrayGeneric &	mda,
		int	nthreads = 1
	)

MultidimArrayGeneric version

Definition at line 736 of file xmipp_fftw.cpp.

{
    MultidimArray<double> aux;
    Mpmem.getImage(aux);
    selfScaleToSizeFourier(Zdim, Ydim, Xdim, aux, nThreads);
    Mpmem.setImage(aux);
}

selfScaleToSizeFourier() [4/4]

void selfScaleToSizeFourier	(	int	Ydim,
		int	Xdim,
		MultidimArrayGeneric &	mda,
		int	nthreads = 1
	)

Definition at line 744 of file xmipp_fftw.cpp.

{
    MultidimArray<double> aux;
    Mpmem.getImage(aux);
    selfScaleToSizeFourier(1, Ydim, Xdim, aux, nThreads);
    Mpmem.setImage(aux);
}

whitenSpectrum()

void whitenSpectrum	(	MultidimArray< double > &	Min,
		MultidimArray< double > &	Mout,
		int	spectrum_type = AMPLITUDE_SPECTRUM,
		bool	leave_origin_intact = false
	)

Perform a whitening of the amplitude/power spectrum of a 3D map

If leave_origin_intact==true, the origin pixel will remain untouched

Definition at line 835 of file xmipp_fftw.cpp.

{
    Min.checkDimension(3);

    MultidimArray<double> spectrum;
    getSpectrum(Min,spectrum,spectrum_type);
    Mout=Min;
    divideBySpectrum(Mout,spectrum,leave_origin_intact);

}