volume_subtraction.cpp File Reference

#include "volume_subtraction.h"
#include "core/transformations.h"
#include <core/histogram.h>
#include <core/xmipp_fftw.h>
#include <core/xmipp_program.h>
#include <data/fourier_filter.h>

Include dependency graph for volume_subtraction.cpp:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions
void	POCSmask (const MultidimArray< double > &mask, MultidimArray< double > &I)
void	POCSnonnegative (MultidimArray< double > &I)
void	POCSFourierAmplitude (const MultidimArray< double > &V1FourierMag, MultidimArray< std::complex< double >> &V2Fourier, double l)
void	POCSFourierAmplitudeRadAvg (MultidimArray< std::complex< double >> &V, double l, const MultidimArray< double > &rQ, int V1size_x, int V1size_y, int V1size_z)
void	POCSMinMax (MultidimArray< double > &V, double v1m, double v1M)
void	POCSFourierPhase (const MultidimArray< std::complex< double >> &phase, MultidimArray< std::complex< double >> &FI)
void	radialAverage (const MultidimArray< double > &VolFourierMag, const MultidimArray< double > &V, MultidimArray< double > &radial_mean)
MultidimArray< double >	computeRadQuotient (const MultidimArray< double > &v1Mag, const MultidimArray< double > &vMag, const MultidimArray< double > &V1, const MultidimArray< double > &V)
void	createFilter (FourierFilter &filter2, double cutFreq)
Image< double >	subtraction (Image< double > V1, Image< double > &V, const MultidimArray< double > &mask, const FileName &fnVol1F, const FileName &fnVol2A, FourierFilter &filter2, double cutFreq)
MultidimArray< double >	computeMagnitude (MultidimArray< double > &volume)

Function Documentation

computeMagnitude()

MultidimArray<double> computeMagnitude

(

MultidimArray< double > &

volume

)

Definition at line 308 of file volume_subtraction.cpp.

                                                                      {
    FourierTransformer transformer;
    MultidimArray<std::complex<double>> fourier;
    MultidimArray<double> magnitude;
    transformer.FourierTransform(volume, fourier, false);
    FFT_magnitude(fourier, magnitude);
    return magnitude;
}

computeRadQuotient()

MultidimArray<double> computeRadQuotient	(	const MultidimArray< double > &	v1Mag,
		const MultidimArray< double > &	vMag,
		const MultidimArray< double > &	V1,
		const MultidimArray< double > &	V
	)

Definition at line 259 of file volume_subtraction.cpp.

                                        {
    // Compute the quotient of the radial mean of the volumes to use it in POCS amplitude
    MultidimArray<double> radial_meanV1;
    radialAverage(v1Mag, V1, radial_meanV1);
    MultidimArray<double> radial_meanV;
    radialAverage(vMag, V, radial_meanV);
    MultidimArray<double> radQuotient = radial_meanV1 / radial_meanV;
    FOR_ALL_ELEMENTS_IN_ARRAY1D(radQuotient)
        {
            radQuotient(i) = std::min(radQuotient(i), 1.0);
            if (radQuotient(i) != radQuotient(i)) // Check if it is NaN and change it by 0
                radQuotient(i) = 0; 
        }
    return radQuotient;
}

createFilter()

void createFilter	(	FourierFilter &	filter2,
		double	cutFreq
	)

Definition at line 277 of file volume_subtraction.cpp.

                                                          {
    filter2.FilterBand = LOWPASS;
    filter2.FilterShape = RAISED_COSINE;
    filter2.raised_w = 0.02;
    filter2.w1 = cutFreq;
}

POCSFourierAmplitude()

void POCSFourierAmplitude	(	const MultidimArray< double > &	V1FourierMag,
		MultidimArray< std::complex< double >> &	V2Fourier,
		double	l
	)

Definition at line 138 of file volume_subtraction.cpp.

                                                                {
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(V1FourierMag) {
        double mod = std::abs(DIRECT_MULTIDIM_ELEM(V2Fourier, n));
        if (mod > 1e-10) // Condition to avoid divide by zero, values smaller than
            // this threshold are considered zero
            DIRECT_MULTIDIM_ELEM(V2Fourier, n) *=
                    ((1 - l) + l * DIRECT_MULTIDIM_ELEM(V1FourierMag, n)) / mod;
    }
}

POCSFourierAmplitudeRadAvg()

void POCSFourierAmplitudeRadAvg	(	MultidimArray< std::complex< double >> &	V,
		double	l,
		const MultidimArray< double > &	rQ,
		int	V1size_x,
		int	V1size_y,
		int	V1size_z
	)

Definition at line 149 of file volume_subtraction.cpp.

                                          {
    int V1size2_x = V1size_x/2;
    double V1sizei_x = 1.0/V1size_x;
    int V1size2_y = V1size_y/2;
    double V1sizei_y = 1.0/V1size_y;
    int V1size2_z = V1size_z/2;
    double V1sizei_z = 1.0/V1size_z;
    double wx;
    double wy;
    double wz;
    for (int k=0; k<ZSIZE(V); ++k)
    {
        FFT_IDX2DIGFREQ_FAST(k,V1size_z,V1size2_z,V1sizei_z,wz)
        double wz2 = wz*wz;
        for (int i=0; i<YSIZE(V); ++i)
        {
            FFT_IDX2DIGFREQ_FAST(i,V1size_y,V1size2_y,V1sizei_y,wy)
            double wy2_wz2 = wy*wy + wz2;
            for (int j=0; j<XSIZE(V); ++j)
            {
                FFT_IDX2DIGFREQ_FAST(j,V1size_x,V1size2_x,V1sizei_x,wx)
                double w = sqrt(wx*wx + wy2_wz2);
                auto iw = std::min((int)floor(w*V1size_x), (int)XSIZE(rQ) -1);
                DIRECT_A3D_ELEM(V,k,i,j)*=(1-l)+l*DIRECT_MULTIDIM_ELEM(rQ,iw);
            }
        }
    }
}

POCSFourierPhase()

void POCSFourierPhase	(	const MultidimArray< std::complex< double >> &	phase,
		MultidimArray< std::complex< double >> &	FI
	)

Definition at line 189 of file volume_subtraction.cpp.

                                               {
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(phase)
                            DIRECT_MULTIDIM_ELEM(FI, n) =
                                    std::abs(DIRECT_MULTIDIM_ELEM(FI, n)) * DIRECT_MULTIDIM_ELEM(phase, n);
}

POCSmask()

void POCSmask	(	const MultidimArray< double > &	mask,
		MultidimArray< double > &	I
	)

Definition at line 128 of file volume_subtraction.cpp.

                                                                           {
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(I)
        DIRECT_MULTIDIM_ELEM(I, n) *= DIRECT_MULTIDIM_ELEM(mask, n);
}

POCSMinMax()

void POCSMinMax	(	MultidimArray< double > &	V,
		double	v1m,
		double	v1M
	)

Definition at line 179 of file volume_subtraction.cpp.

                                                                  {
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(V) {
        double val = DIRECT_MULTIDIM_ELEM(V, n);
        if (val < v1m)
            DIRECT_MULTIDIM_ELEM(V, n) = v1m;
        else if (val > v1M)
            DIRECT_MULTIDIM_ELEM(V, n) = v1M;
    }
}

POCSnonnegative()

void POCSnonnegative

(

MultidimArray< double > &

I

)

Definition at line 133 of file volume_subtraction.cpp.

                                               {
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(I)
        DIRECT_MULTIDIM_ELEM(I, n) = std::max(0.0, DIRECT_MULTIDIM_ELEM(I, n));
}

radialAverage()

void radialAverage	(	const MultidimArray< double > &	VolFourierMag,
		const MultidimArray< double > &	V,
		MultidimArray< double > &	radial_mean
	)

Definition at line 223 of file volume_subtraction.cpp.

                                                                            {
    MultidimArray<double> radial_count;
    int Vsize2_x = int(XSIZE(V))/2;
    double Vsizei_x = 1.0/int(XSIZE(V));
    int Vsize2_y = int(YSIZE(V))/2;
    double Vsizei_y = 1.0/int(YSIZE(V));
    int Vsize2_z = int(ZSIZE(V))/2;
    double Vsizei_z = 1.0/int(ZSIZE(V));
    double wx;
    double wy;
    double wz;
    auto maxrad = int(floor(sqrt(Vsize2_x*Vsize2_x + Vsize2_y*Vsize2_y + Vsize2_z*Vsize2_z)));
    radial_count.initZeros(maxrad);
    radial_mean.initZeros(maxrad);
    for (int k=0; k<Vsize2_z; ++k)
        {
            FFT_IDX2DIGFREQ_FAST(k,ZSIZE(V),Vsize2_z,Vsizei_z,wz)
            double wz2 = wz*wz;
            for (int i=0; i<Vsize2_y; ++i)
            {
                FFT_IDX2DIGFREQ_FAST(i,YSIZE(V),Vsize2_y,Vsizei_y,wy)
                double wy2_wz2 = wy*wy + wz2;
                for (int j=0; j<Vsize2_x; ++j)
                {
                    FFT_IDX2DIGFREQ_FAST(j,XSIZE(V),Vsize2_x,Vsizei_x,wx)
                    double w = sqrt(wx*wx + wy2_wz2);
                    auto iw = (int)round(w*int(XSIZE(V)));
                    DIRECT_A1D_ELEM(radial_mean,iw)+=DIRECT_A3D_ELEM(VolFourierMag,k,i,j);
                    DIRECT_A1D_ELEM(radial_count,iw)+=1.0;
                }
            }
        }
    radial_mean/= radial_count;
}

subtraction()

Image<double> subtraction	(	Image< double >	V1,
		Image< double > &	V,
		const MultidimArray< double > &	mask,
		const FileName &	fnVol1F,
		const FileName &	fnVol2A,
		FourierFilter &	filter2,
		double	cutFreq
	)

Definition at line 284 of file volume_subtraction.cpp.

                                                                         {
    Image<double> V1Filtered;
    V1Filtered() = V1();
    if (cutFreq != 0){
        filter2.applyMaskSpace(V1Filtered());
    }
    if (!fnVol1F.isEmpty()) {
        V1Filtered.write(fnVol1F);
    }
    if (!fnVol2A.isEmpty()) {
        V.write(fnVol2A);
    }
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(V1())
                            DIRECT_MULTIDIM_ELEM(V1(), n) =
                                    DIRECT_MULTIDIM_ELEM(V1(), n) * (1 - DIRECT_MULTIDIM_ELEM(mask, n)) +
                                    (DIRECT_MULTIDIM_ELEM(V1Filtered(), n) -
                                            std::min(DIRECT_MULTIDIM_ELEM(V(), n),
                                                    DIRECT_MULTIDIM_ELEM(V1Filtered(), n))) *
                                                    DIRECT_MULTIDIM_ELEM(mask, n);
    return V1;
}