#include <resolution_monotomo.h>

Inheritance diagram for ProgMonoTomo:

Collaboration diagram for ProgMonoTomo:

Public Member Functions
void	defineParams ()

void	readParams ()

void	produceSideInfo ()

void	amplitudeMonogenicSignal3D (MultidimArray< std::complex< double > > &myfftV, double freq, double freqH, double freqL, MultidimArray< double > &amplitude, int count, FileName fnDebug)

void	firstMonoResEstimation (MultidimArray< std::complex< double > > &myfftV, double freq, double freqH, double freqL, MultidimArray< double > &amplitude, int count, FileName fnDebug, double &mean_Signal, double &mean_noise, double &thresholdFirstEstimation)

void	median3x3x3 (MultidimArray< double > vol, MultidimArray< double > &filtered)

void	localNoise (MultidimArray< double > &noiseMap, Matrix2D< double > &noiseMatrix, int boxsize, Matrix2D< double > &thresholdMatrix)

void	postProcessingLocalResolutions (MultidimArray< double > &resolutionVol, std::vector< double > &list)

void	resolution2eval (int &count_res, double step, double &resolution, double &last_resolution, double &freq, double &freqL, int &last_fourier_idx, bool &continueIter, bool &breakIter)

void	lowestResolutionbyPercentile (MultidimArray< double > &resolutionVol, std::vector< double > &list, double &cut_value, double &resolutionThreshold)

void	run ()

void	defineParams ()

void	readParams ()

void	produceSideInfo ()

void	amplitudeMonogenicSignal3D (const std::vector< std::complex< float >> &myfftV, float freq, float freqH, float freqL, MultidimArray< float > &amplitude, int count, FileName fnDebug)

void	firstMonoResEstimation (MultidimArray< std::complex< double > > &myfftV, double freq, double freqH, double freqL, MultidimArray< double > &amplitude, int count, FileName fnDebug, double &mean_Signal, double &mean_noise, double &thresholdFirstEstimation)

void	median3x3x3 (MultidimArray< double > vol, MultidimArray< double > &filtered)

void	localNoise (MultidimArray< float > &noiseMap, Matrix2D< double > &noiseMatrix, int boxsize, Matrix2D< double > &thresholdMatrix)

void	postProcessingLocalResolutions (MultidimArray< float > &resolutionVol, std::vector< float > &list)

void	resolution2eval (int &count_res, double step, double &resolution, double &last_resolution, double &freq, double &freqL, int &last_fourier_idx, bool &continueIter, bool &breakIter)

void	smoothBorders (MultidimArray< float > &vol, MultidimArray< int > &pMask)

void	lowestResolutionbyPercentile (MultidimArray< float > &resolutionVol, std::vector< float > &list, float &cut_value, float &resolutionThreshold)

void	getFilteringResolution (size_t idx, float freq, float lastResolution, float freqL, float &resolution)

void	gaussFilter (const MultidimArray< float > &vol, const float, MultidimArray< float > &VRiesz)

void	run ()

Public Member Functions inherited from XmippProgram
const char *	getParam (const char *param, int arg=0)

const char *	getParam (const char param, const char subparam, int arg=0)

int	getIntParam (const char *param, int arg=0)

int	getIntParam (const char param, const char subparam, int arg=0)

double	getDoubleParam (const char *param, int arg=0)

double	getDoubleParam (const char param, const char subparam, int arg=0)

float	getFloatParam (const char *param, int arg=0)

float	getFloatParam (const char param, const char subparam, int arg=0)

void	getListParam (const char *param, StringVector &list)

int	getCountParam (const char *param)

bool	checkParam (const char *param)

bool	existsParam (const char *param)

void	addParamsLine (const String &line)

void	addParamsLine (const char *line)

ParamDef *	getParamDef (const char *param) const

virtual void	quit (int exit_code=0) const

virtual int	tryRun ()

void	initProgress (size_t total, size_t stepBin=60)

void	setProgress (size_t value=0)

void	endProgress ()

void	processDefaultComment (const char param, const char left)

void	setDefaultComment (const char param, const char comment)

virtual void	initComments ()

void	setProgramName (const char *name)

void	addUsageLine (const char *line, bool verbatim=false)

void	clearUsage ()

void	addExampleLine (const char *example, bool verbatim=true)

void	addSeeAlsoLine (const char *seeAlso)

void	addKeywords (const char *keywords)

const char *	name () const

virtual void	usage (int verb=0) const

virtual void	usage (const String &param, int verb=2)

int	version () const

virtual void	show () const

virtual void	read (int argc, const char **argv, bool reportErrors=true)

virtual void	read (int argc, char **argv, bool reportErrors=true)

void	read (const String &argumentsLine)

	XmippProgram ()

	XmippProgram (int argc, const char **argv)

virtual	~XmippProgram ()

Public Attributes
FileName	fnOut

FileName	fnVol

FileName	fnVol2

FileName	fnMask

FileName	fnchim

FileName	fnSpatial

FileName	fnMeanVol

FileName	fnMaskOut

FileName	fnMd

FileName	fnFilt

FileName	fnmaskWedge

double	sampling

double	minRes

double	maxRes

double	R

int	NVoxelsOriginalMask

int	Nvoxels

int	nthrs

double	freq_step

double	trimBound

double	significance

bool	noiseOnlyInHalves

bool	automaticMode

Image< int >	mask

MultidimArray< double >	iu

MultidimArray< double >	VRiesz

MultidimArray< std::complex< double > >	fftV

MultidimArray< std::complex< double > > *	fftN

FourierTransformer	transformer_inv

MultidimArray< std::complex< double > >	fftVRiesz

MultidimArray< std::complex< double > >	fftVRiesz_aux

FourierFilter	lowPassFilter

FourierFilter	FilterBand

bool	halfMapsGiven

Image< double >	Vfiltered

Image< double >	VresolutionFiltered

Matrix1D< double >	freq_fourier_z

Matrix1D< double >	freq_fourier_y

Matrix1D< double >	freq_fourier_x

Matrix2D< double >	resolutionMatrix

Matrix2D< double >	maskMatrix

size_t	xdimFT

size_t	ydimFT

size_t	zdimFT

size_t	xdim

size_t	ydim

size_t	zdim

double	resStep

std::vector< std::complex< float > >	fourierSignal

std::vector< std::complex< float > >	fourierNoise

std::vector< std::complex< float > >	fftVRiesz

std::vector< std::complex< float > >	fftVRiesz_aux

MultidimArray< float >	VRiesz

std::unique_ptr< AFT< float > >	forward_transformer

std::unique_ptr< AFT< float > >	backward_transformer

Public Attributes inherited from XmippProgram
bool	doRun

bool	runWithoutArgs

int	verbose
	Verbosity level. More...

int	debug

Additional Inherited Members
Protected Member Functions inherited from XmippProgram
void	defineCommons ()

Protected Attributes inherited from XmippProgram
int	errorCode

ProgramDef *	progDef
	Program definition and arguments parser. More...

std::map< String, CommentList >	defaultComments

int	argc
	Original command line arguments. More...

const char **	argv

Detailed Description

SSNR parameters.

Definition at line 47 of file resolution_monotomo.h.

Member Function Documentation

◆ amplitudeMonogenicSignal3D() [1/2]

void ProgMonoTomo::amplitudeMonogenicSignal3D	(	MultidimArray< std::complex< double > > &	myfftV,
		double	freq,
		double	freqH,
		double	freqL,
		MultidimArray< double > &	amplitude,
		int	count,
		FileName	fnDebug
	)

Definition at line 290 of file resolution_monotomo.cpp.

 {
     fftVRiesz.initZeros(myfftV);
     fftVRiesz_aux.initZeros(myfftV);
     std::complex<double> J(0,1);
 
     // Filter the input volume and add it to amplitude
     long n=0;
     double ideltal=PI/(freq-freqH);
 
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(myfftV)
     {
         double iun=DIRECT_MULTIDIM_ELEM(iu,n);
         double un=1.0/iun;
         if (freqH<=un && un<=freq)
         {
             //double H=0.5*(1+cos((un-w1)*ideltal));
             DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = DIRECT_MULTIDIM_ELEM(myfftV, n);
             DIRECT_MULTIDIM_ELEM(fftVRiesz, n) *= 0.5*(1+cos((un-freq)*ideltal));//H;
             DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n) = -J;
             DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n) *= DIRECT_MULTIDIM_ELEM(fftVRiesz, n);
             DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n) *= iun;
         } else if (un>freq)
         {
             DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = DIRECT_MULTIDIM_ELEM(myfftV, n);
             DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n) = -J;
             DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n) *= DIRECT_MULTIDIM_ELEM(fftVRiesz, n);
             DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n) *= iun;
         }
     }
 
     transformer_inv.inverseFourierTransform(fftVRiesz, VRiesz);
 
     #ifdef DEBUG
     Image<double> filteredvolume;
     filteredvolume = VRiesz;
     filteredvolume.write(formatString("Volumen_filtrado_%i.vol", count));
 
     FileName iternumber;
     iternumber = formatString("_Volume_%i.vol", count);
     Image<double> saveImg2;
     saveImg2() = VRiesz;
       if (fnDebug.c_str() != "")
       {
         saveImg2.write(fnDebug+iternumber);
       }
     saveImg2.clear(); 
     #endif
 
     amplitude.resizeNoCopy(VRiesz);
 
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(amplitude)
         DIRECT_MULTIDIM_ELEM(amplitude,n)=DIRECT_MULTIDIM_ELEM(VRiesz,n)*DIRECT_MULTIDIM_ELEM(VRiesz,n);
 
     // Calculate first component of Riesz vector
     double uz, uy, ux;
     n=0;
     for(size_t k=0; k<ZSIZE(myfftV); ++k)
     {
         for(size_t i=0; i<YSIZE(myfftV); ++i)
         {
             for(size_t j=0; j<XSIZE(myfftV); ++j)
             {
                 ux = VEC_ELEM(freq_fourier_x,j);
                 DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = ux*DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n);
                 ++n;
             }
         }
     }
     transformer_inv.inverseFourierTransform(fftVRiesz, VRiesz);
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(amplitude)
         DIRECT_MULTIDIM_ELEM(amplitude,n)+=DIRECT_MULTIDIM_ELEM(VRiesz,n)*DIRECT_MULTIDIM_ELEM(VRiesz,n);
 
     // Calculate second and third components of Riesz vector
     n=0;
     for(size_t k=0; k<ZSIZE(myfftV); ++k)
     {
         uz = VEC_ELEM(freq_fourier_z,k);
         for(size_t i=0; i<YSIZE(myfftV); ++i)
         {
             uy = VEC_ELEM(freq_fourier_y,i);
             for(size_t j=0; j<XSIZE(myfftV); ++j)
             {
                 DIRECT_MULTIDIM_ELEM(fftVRiesz, n) = uy*DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n);
                 DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n) = uz*DIRECT_MULTIDIM_ELEM(fftVRiesz_aux, n);
                 ++n;
             }
         }
     }
     transformer_inv.inverseFourierTransform(fftVRiesz, VRiesz);
 
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(amplitude)
         DIRECT_MULTIDIM_ELEM(amplitude,n)+=DIRECT_MULTIDIM_ELEM(VRiesz,n)*DIRECT_MULTIDIM_ELEM(VRiesz,n);
 
     transformer_inv.inverseFourierTransform(fftVRiesz_aux, VRiesz);
 
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(amplitude)
     {
         DIRECT_MULTIDIM_ELEM(amplitude,n)+=DIRECT_MULTIDIM_ELEM(VRiesz,n)*DIRECT_MULTIDIM_ELEM(VRiesz,n);
         DIRECT_MULTIDIM_ELEM(amplitude,n)=sqrt(DIRECT_MULTIDIM_ELEM(amplitude,n));
     }
     #ifdef DEBUG
     if (fnDebug.c_str() != "")
     {
     Image<double> saveImg;
     saveImg = amplitude;
     FileName iternumber;
     iternumber = formatString("_Amplitude_%i.vol", count);
     saveImg.write(fnDebug+iternumber);
     saveImg.clear();
     }
     #endif // DEBUG
 //
     // Low pass filter the monogenic amplitude
     transformer_inv.FourierTransform(amplitude, fftVRiesz, false);
     double raised_w = PI/(freqL-freq);
 
     n=0;
 
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(fftVRiesz)
     {
         double un=1.0/DIRECT_MULTIDIM_ELEM(iu,n);
         if (freqL>=un && un>=freq)
         {
             DIRECT_MULTIDIM_ELEM(fftVRiesz,n) *= 0.5*(1 + cos(raised_w*(un-freq)));
         }
         else
         {
             if (un>freqL)
             {
                 DIRECT_MULTIDIM_ELEM(fftVRiesz,n) = 0;
             }
         }
     }
     transformer_inv.inverseFourierTransform();
 
 
     #ifdef DEBUG
     Image<double> saveImg2;
     saveImg2 = amplitude;
     FileName fnSaveImg2;
     if (fnDebug.c_str() != "")
     {
         iternumber = formatString("_Filtered_Amplitude_%i.vol", count);
         saveImg2.write(fnDebug+iternumber);
     }
     saveImg2.clear();
     #endif // DEBUG
 }

◆ amplitudeMonogenicSignal3D() [2/2]

void ProgMonoTomo::amplitudeMonogenicSignal3D	(	const std::vector< std::complex< float >> &	myfftV,
		float	freq,
		float	freqH,
		float	freqL,
		MultidimArray< float > &	amplitude,
		int	count,
		FileName	fnDebug
	)

Definition at line 185 of file resolution_monotomo.cpp.

 {
     fftVRiesz.resize(fourierSignal.size());
     fftVRiesz_aux.resize(fourierSignal.size());
     std::complex<float> J(0,1);
 
     // Filter the input volume and add it to amplitude
     size_t n=0;
     float ideltal=PI/(freq-freqH);
     for(size_t k=0; k<zdimFT; ++k)
     {
         double uz;// = VEC_ELEM(freq_fourier_z, k);
         FFT_IDX2DIGFREQ(k, zdim, uz);
         auto uz2 = uz*uz;
         for(size_t i=0; i<ydimFT; ++i)
         {
             //const auto uy = VEC_ELEM(freq_fourier_y, i);
             double uy;
             FFT_IDX2DIGFREQ(i, ydim, uy);
             auto uy2uz2 = uz2 +uy*uy;
             for(size_t j=0; j<xdimFT; ++j)
             {
                 double ux;
                 FFT_IDX2DIGFREQ(j, xdim, ux);
                 //const auto ux = VEC_ELEM(freq_fourier_x, j);
                 const float u = std::sqrt(uy2uz2 + ux*ux);
 
                 if (freqH<=u && u<=freq)
                 {
                     fftVRiesz[n] = myfftV[n]*0.5f*(1+std::cos((u-freq)*ideltal));
                     fftVRiesz_aux[n] = -J*fftVRiesz[n]/u;
                 }
                 else
                 {
                     if (u>freq)
                     {
                         fftVRiesz[n] = myfftV[n];
                         fftVRiesz_aux[n] = -J*fftVRiesz[n]/u;
                     }
                     else
                     {
                         fftVRiesz[n] = 0.0f;
                         fftVRiesz_aux[n] = 0.0f;
                     }
                 }
                 n++;
 
             }
         }
     }
 
     backward_transformer->ifft(fftVRiesz.data(), MULTIDIM_ARRAY(VRiesz));
 
 
     amplitude = VRiesz;
 
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(VRiesz)
     {
         DIRECT_MULTIDIM_ELEM(amplitude, n) *= DIRECT_MULTIDIM_ELEM(VRiesz, n);
     }
 
 
     // Calculate first component of Riesz vector
     float uz, uy, ux;
     n=0;
     for(size_t k=0; k<zdimFT; ++k)
     {
         for(size_t i=0; i<ydimFT; ++i)
         {
             for(size_t j=0; j<xdimFT; ++j)
             {
                 //ux = VEC_ELEM(freq_fourier_x, j);
                 FFT_IDX2DIGFREQ(j, xdim, ux);
                 fftVRiesz[n] = ux*fftVRiesz_aux[n];
                 ++n;
             }
         }
     }
     backward_transformer->ifft(fftVRiesz.data(), MULTIDIM_ARRAY(VRiesz));
 
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(VRiesz)
     {
         DIRECT_MULTIDIM_ELEM(amplitude, n) += DIRECT_MULTIDIM_ELEM(VRiesz, n)*DIRECT_MULTIDIM_ELEM(VRiesz, n);
     }
 
     // Calculate second and third components of Riesz vector
     n=0;
     for(size_t k=0; k<zdimFT; ++k)
     {
         //uz = VEC_ELEM(freq_fourier_z, k);
         // = VEC_ELEM(freq_fourier_z, k);
         FFT_IDX2DIGFREQ(k, zdim, uz);
         for(size_t i=0; i<ydimFT; ++i)
         {
             //uy = VEC_ELEM(freq_fourier_y, i);
             FFT_IDX2DIGFREQ(i, ydim, uy);
             for(size_t j=0; j<xdimFT; ++j)
             {
                 fftVRiesz[n] = ux*fftVRiesz_aux[n];
                 fftVRiesz_aux[n] = uz*fftVRiesz_aux[n];
                 ++n;
             }
         }
     }
     backward_transformer->ifft(fftVRiesz.data(), MULTIDIM_ARRAY(VRiesz));
 
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(VRiesz)
     {
         DIRECT_MULTIDIM_ELEM(amplitude, n) += DIRECT_MULTIDIM_ELEM(VRiesz, n)*DIRECT_MULTIDIM_ELEM(VRiesz, n);
     }
 
 
     backward_transformer->ifft(fftVRiesz_aux.data(), MULTIDIM_ARRAY(VRiesz));
 
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(VRiesz)
     {
         DIRECT_MULTIDIM_ELEM(amplitude, n) += DIRECT_MULTIDIM_ELEM(VRiesz, n)*DIRECT_MULTIDIM_ELEM(VRiesz, n);
         DIRECT_MULTIDIM_ELEM(amplitude, n) = std::sqrt(DIRECT_MULTIDIM_ELEM(amplitude, n));
     }
 
     // Low pass filter the monogenic amplitude
     forward_transformer->fft(MULTIDIM_ARRAY(amplitude), fftVRiesz.data());
     float raised_w = PI/(freqL-freq);
 
     n = 0;
     for(size_t k=0; k<zdimFT; ++k)
     {
         //uz = VEC_ELEM(freq_fourier_z, k);
         double uz;// = VEC_ELEM(freq_fourier_z, k);
         FFT_IDX2DIGFREQ(k, zdim, uz);
         auto uz2 = uz*uz;
         for(size_t i=0; i<ydimFT; ++i)
         {
             //uy = VEC_ELEM(freq_fourier_y, i);
             double uy;
             FFT_IDX2DIGFREQ(i, ydim, uy);
             auto uy2uz2 = uz2 +uy*uy;
             for(size_t j=0; j<xdimFT; ++j)
             {
                 //ux = VEC_ELEM(freq_fourier_x, j);
                 double ux;
                 FFT_IDX2DIGFREQ(j, xdim, ux);
                 auto u = std::sqrt(uy2uz2 + ux*ux);
 
                 if (u>freqL)
                 {
                     fftVRiesz[n] = 0.0f;
                 }
                 else
                 {
                     if (freqL>=u && u>=freq)
                     {
                         fftVRiesz[n] *= 0.5f*(1 + std::cos(raised_w*(u-freq)));
                     }
                 }
                 n++;
 
             }
         }
     }
 
     backward_transformer->ifft(fftVRiesz.data(), MULTIDIM_ARRAY(amplitude));
 }

◆ defineParams() [1/2]

void ProgMonoTomo::defineParams ( )

virtual

Function in which the param of each Program are defined.

Reimplemented from XmippProgram.

Definition at line 55 of file resolution_monotomo.cpp.

 {
     addUsageLine("This function determines the local resolution of a tomogram. It makes use of two reconstructions, odd and even. The difference between them"
             "gives a noise reconstruction. Thus, by computing the local amplitude of the signal at different frequencies and establishing a comparison with"
             "the noise, the local resolution is computed");
     addParamsLine("  --vol <vol_file=\"\">              : Half volume 1");
     addParamsLine("  --vol2 <vol_file=\"\">             : Half volume 2");
     addParamsLine("  [--mask <vol_file=\"\">]           : Mask defining the signal. ");
     addParamsLine("  -o <output=\"MGresolution.vol\">   : Local resolution volume (in Angstroms)");
     addParamsLine("  [--maskWedge <vol_file=\"\">]  : Mask containing the missing wedge in Fourier space");
     addParamsLine("  [--filteredMap <vol_file=\"\">]                : Local resolution volume filtered considering the missing wedge (in Angstroms)");
     addParamsLine("  --meanVol <vol_file=\"\">          : Mean volume of half1 and half2 (only it is necessary the two haves are used)");
     addParamsLine("  [--sampling_rate <s=1>]            : Sampling rate (A/px)");
     addParamsLine("  [--step <s=0.25>]                  : The resolution is computed at a number of frequencies between minimum and");
     addParamsLine("                                     : maximum resolution px/A. This parameter determines that number");
     addParamsLine("  [--minRes <s=30>]                  : Minimum resolution (A)");
     addParamsLine("  [--maxRes <s=1>]                   : Maximum resolution (A)");
     addParamsLine("  [--trimmed <s=0.5>]                : Trimming percentile");
     addParamsLine("  [--significance <s=0.95>]          : The level of confidence for the hypothesis test.");
     addParamsLine("  [--threads <s=4>]                  : Number of threads");
 }

◆ defineParams() [2/2]

void ProgMonoTomo::defineParams ( )

virtual

Function in which the param of each Program are defined.

Reimplemented from XmippProgram.

◆ firstMonoResEstimation() [1/2]

void ProgMonoTomo::firstMonoResEstimation	(	MultidimArray< std::complex< double > > &	myfftV,
		double	freq,
		double	freqH,
		double	freqL,
		MultidimArray< double > &	amplitude,
		int	count,
		FileName	fnDebug,
		double &	mean_Signal,
		double &	mean_noise,
		double &	thresholdFirstEstimation
	)

◆ firstMonoResEstimation() [2/2]

void ProgMonoTomo::firstMonoResEstimation	(	MultidimArray< std::complex< double > > &	myfftV,
		double	freq,
		double	freqH,
		double	freqL,
		MultidimArray< double > &	amplitude,
		int	count,
		FileName	fnDebug,
		double &	mean_Signal,
		double &	mean_noise,
		double &	thresholdFirstEstimation
	)

◆ gaussFilter()

void ProgMonoTomo::gaussFilter	(	const MultidimArray< float > &	vol,
		const float	sigma,
		MultidimArray< float > &	VRiesz
	)

Definition at line 567 of file resolution_monotomo.cpp.

 {
     float isigma2 = (sigma*sigma);
 
     forward_transformer->fft(MULTIDIM_ARRAY(vol), fftVRiesz.data());
     size_t n=0;
     for(size_t k=0; k<zdimFT; ++k)
     {
         //const auto uz = VEC_ELEM(freq_fourier_z, k);
         double uz;
         FFT_IDX2DIGFREQ(k, zdim, uz);
         auto uz2 = uz*uz;
 
         for(size_t i=0; i<ydimFT; ++i)
         {
             //const auto uy = VEC_ELEM(freq_fourier_y, i);
             double uy;
             FFT_IDX2DIGFREQ(i, ydim, uy);
             double uy2uz2 = uz2 +uy*uy;
             for(size_t j=0; j<xdimFT; ++j)
             {
                 //const auto ux = VEC_ELEM(freq_fourier_x, j);
                 double ux;
                 FFT_IDX2DIGFREQ(j, xdim, ux);
                 const float u2 = (float) (uy2uz2 + ux*ux);
 
                 fftVRiesz[n] *= std::exp(-PI*PI*u2*isigma2);
 
                 n++;
 
             }
         }
     }
 
     backward_transformer->ifft(fftVRiesz.data(), MULTIDIM_ARRAY(VRiesz));
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(VRiesz)
     {
         DIRECT_MULTIDIM_ELEM(VRiesz, n) /= xdim*ydim*zdim;
     }
 }

◆ getFilteringResolution()

void ProgMonoTomo::getFilteringResolution	(	size_t	idx,
		float	freq,
		float	lastResolution,
		float	freqL,
		float &	resolution
	)

◆ localNoise() [1/2]

void ProgMonoTomo::localNoise	(	MultidimArray< float > &	noiseMap,
		Matrix2D< double > &	noiseMatrix,
		int	boxsize,
		Matrix2D< double > &	thresholdMatrix
	)

Definition at line 351 of file resolution_monotomo.cpp.

 {
 //  std::cout << "Analyzing local noise" << std::endl;
 
     int xdim = XSIZE(noiseMap);
     int ydim = YSIZE(noiseMap);
 
     int Nx = xdim/boxsize;
     int Ny = ydim/boxsize;
 
 
 
     // For the spline regression
     int lX=std::min(8,Nx-2), lY=std::min(8,Ny-2);
     WeightedLeastSquaresHelper helper;
     helper.A.initZeros(Nx*Ny,lX*lY);
     helper.b.initZeros(Nx*Ny);
     helper.w.initZeros(Nx*Ny);
     helper.w.initConstant(1);
     double hX = xdim / (double)(lX-3);
     double hY = ydim / (double)(lY-3);
 
     if ( (xdim<boxsize) || (ydim<boxsize) )
         std::cout << "Error: The tomogram in x-direction or y-direction is too small" << std::endl;
 
     std::vector<double> noiseVector(1);
     std::vector<double> x,y,t;
 
     int xLimit, yLimit, xStart, yStart;
 
     long counter;
     int idxBox=0;
 
     for (int X_boxIdx=0; X_boxIdx<Nx; ++X_boxIdx)
     {
         if (X_boxIdx==Nx-1)
         {
             xStart = STARTINGX(noiseMap) + X_boxIdx*boxsize;
             xLimit = FINISHINGX(noiseMap);
         }
         else
         {
             xStart = STARTINGX(noiseMap) + X_boxIdx*boxsize;
             xLimit = STARTINGX(noiseMap) + (X_boxIdx+1)*boxsize;
         }
 
         for (int Y_boxIdx=0; Y_boxIdx<Ny; ++Y_boxIdx)
         {
             if (Y_boxIdx==Ny-1)
             {
                 yStart = STARTINGY(noiseMap) + Y_boxIdx*boxsize;
                 yLimit =  FINISHINGY(noiseMap);
             }
             else
             {
                 yStart = STARTINGY(noiseMap) + Y_boxIdx*boxsize;
                 yLimit = STARTINGY(noiseMap) + (Y_boxIdx+1)*boxsize;
             }
 
 
 
             counter = 0;
             for (int i = yStart; i<yLimit; i++)
             {
                 for (int j = xStart; j<xLimit; j++)
                 {
                     for (int k = STARTINGZ(noiseMap); k<FINISHINGZ(noiseMap); k++)
                     {
                         if (counter%257 == 0) //we take one voxel each 257 (prime number) points to reduce noise data
                             noiseVector.push_back( A3D_ELEM(noiseMap, k, i, j) );
                         ++counter;
                     }
                 }
             }
 
             std::sort(noiseVector.begin(),noiseVector.end());
             noiseMatrix.initZeros(Ny, Nx);
 
             MAT_ELEM(noiseMatrix, Y_boxIdx, X_boxIdx) = noiseVector[size_t(noiseVector.size()*significance)];
 
             double tileCenterY=0.5*(yLimit+yStart)-STARTINGY(noiseMap); // Translated to physical coordinates
             double tileCenterX=0.5*(xLimit+xStart)-STARTINGX(noiseMap);
             // Construction of the spline equation system
             long idxSpline=0;
             for(int controlIdxY = -1; controlIdxY < (lY - 1); ++controlIdxY)
             {
                 double tmpY = Bspline03((tileCenterY / hY) - controlIdxY);
                 VEC_ELEM(helper.b,idxBox)=MAT_ELEM(noiseMatrix, Y_boxIdx, X_boxIdx);
                 if (tmpY == 0.0)
                 {
                     idxSpline+=lX;
                     continue;
                 }
 
                 for(int controlIdxX = -1; controlIdxX < (lX - 1); ++controlIdxX)
                 {
                     double tmpX = Bspline03((tileCenterX / hX) - controlIdxX);
                     MAT_ELEM(helper.A,idxBox,idxSpline) = tmpY * tmpX;
                     idxSpline+=1;
                 }
             }
             x.push_back(tileCenterX);
             y.push_back(tileCenterY);
             t.push_back(MAT_ELEM(noiseMatrix, Y_boxIdx, X_boxIdx));
             noiseVector.clear();
             idxBox+=1;
         }
     }
 
 
     // Spline coefficients
     Matrix1D<double> cij;
     weightedLeastSquares(helper, cij);
 
     thresholdMatrix.initZeros(ydim, xdim);
 
     for (int i=0; i<ydim; ++i)
     {
         for (int j=0; j<xdim; ++j)
         {
             long idxSpline=0;
 
             for(int controlIdxY = -1; controlIdxY < (lY - 1); ++controlIdxY)
             {
                 double tmpY = Bspline03((i / hY) - controlIdxY);
 
                 if (tmpY == 0.0)
                 {
                     idxSpline+=lX;
                     continue;
                 }
 
                 double xContrib=0.0;
                 for(int controlIdxX = -1; controlIdxX < (lX - 1); ++controlIdxX)
                 {
                     double tmpX = Bspline03((j / hX) - controlIdxX);
                     xContrib+=VEC_ELEM(cij,idxSpline) * tmpX;// *tmpY;
                     idxSpline+=1;
                 }
                 MAT_ELEM(thresholdMatrix,i,j)+=xContrib*tmpY;
             }
         }
     }
 }

◆ localNoise() [2/2]

void ProgMonoTomo::localNoise	(	MultidimArray< double > &	noiseMap,
		Matrix2D< double > &	noiseMatrix,
		int	boxsize,
		Matrix2D< double > &	thresholdMatrix
	)

Definition at line 442 of file resolution_monotomo.cpp.

 {
 //  std::cout << "Analyzing local noise" << std::endl;
 
     int xdim = XSIZE(noiseMap);
     int ydim = YSIZE(noiseMap);
 
     int Nx = xdim/boxsize;
     int Ny = ydim/boxsize;
 
     noiseMatrix.initZeros(Ny, Nx);
 
     // For the spline regression
     int lX=std::min(8,Nx-2), lY=std::min(8,Ny-2);
     WeightedLeastSquaresHelper helper;
     helper.A.initZeros(Nx*Ny,lX*lY);
     helper.b.initZeros(Nx*Ny);
     helper.w.initZeros(Nx*Ny);
     helper.w.initConstant(1);
     double hX = xdim / (double)(lX-3);
     double hY = ydim / (double)(lY-3);
 
     if ( (xdim<boxsize) || (ydim<boxsize) )
         std::cout << "Error: The tomogram in x-direction or y-direction is too small" << std::endl;
 
     std::vector<double> noiseVector(1);
     std::vector<double> x,y,t;
 
     int xLimit, yLimit, xStart, yStart;
 
     long counter;
     int idxBox=0;
 
     for (int X_boxIdx=0; X_boxIdx<Nx; ++X_boxIdx)
     {
         if (X_boxIdx==Nx-1)
         {
             xStart = STARTINGX(noiseMap) + X_boxIdx*boxsize;
             xLimit = FINISHINGX(noiseMap);
         }
         else
         {
             xStart = STARTINGX(noiseMap) + X_boxIdx*boxsize;
             xLimit = STARTINGX(noiseMap) + (X_boxIdx+1)*boxsize;
         }
 
         for (int Y_boxIdx=0; Y_boxIdx<Ny; ++Y_boxIdx)
         {
             if (Y_boxIdx==Ny-1)
             {
                 yStart = STARTINGY(noiseMap) + Y_boxIdx*boxsize;
                 yLimit =  FINISHINGY(noiseMap);
             }
             else
             {
                 yStart = STARTINGY(noiseMap) + Y_boxIdx*boxsize;
                 yLimit = STARTINGY(noiseMap) + (Y_boxIdx+1)*boxsize;
             }
 
 
 
             counter = 0;
             for (int i = yStart; i<yLimit; i++)
             {
                 for (int j = xStart; j<xLimit; j++)
                 {
                     for (int k = STARTINGZ(noiseMap); k<FINISHINGZ(noiseMap); k++)
                     {
                         if (counter%257 == 0) //we take one voxel each 257 (prime number) points to reduce noise data
                             noiseVector.push_back( A3D_ELEM(noiseMap, k, i, j) );
                         ++counter;
                     }
                 }
             }
 
             std::sort(noiseVector.begin(),noiseVector.end());
             MAT_ELEM(noiseMatrix, Y_boxIdx, X_boxIdx) = noiseVector[size_t(noiseVector.size()*significance)];
 
             double tileCenterY=0.5*(yLimit+yStart)-STARTINGY(noiseMap); // Translated to physical coordinates
             double tileCenterX=0.5*(xLimit+xStart)-STARTINGX(noiseMap);
             // Construction of the spline equation system
             long idxSpline=0;
             for(int controlIdxY = -1; controlIdxY < (lY - 1); ++controlIdxY)
             {
                 double tmpY = Bspline03((tileCenterY / hY) - controlIdxY);
                 VEC_ELEM(helper.b,idxBox)=MAT_ELEM(noiseMatrix, Y_boxIdx, X_boxIdx);
                 if (tmpY == 0.0)
                 {
                     idxSpline+=lX;
                     continue;
                 }
 
                 for(int controlIdxX = -1; controlIdxX < (lX - 1); ++controlIdxX)
                 {
                     double tmpX = Bspline03((tileCenterX / hX) - controlIdxX);
                     MAT_ELEM(helper.A,idxBox,idxSpline) = tmpY * tmpX;
                     idxSpline+=1;
                 }
             }
             x.push_back(tileCenterX);
             y.push_back(tileCenterY);
             t.push_back(MAT_ELEM(noiseMatrix, Y_boxIdx, X_boxIdx));
             noiseVector.clear();
             idxBox+=1;
         }
     }
 
 
     // Spline coefficients
     Matrix1D<double> cij;
     weightedLeastSquares(helper, cij);
 
     thresholdMatrix.initZeros(ydim, xdim);
 
     for (int i=0; i<ydim; ++i)
     {
         for (int j=0; j<xdim; ++j)
         {
             long idxSpline=0;
 
             for(int controlIdxY = -1; controlIdxY < (lY - 1); ++controlIdxY)
             {
                 double tmpY = Bspline03((i / hY) - controlIdxY);
 
                 if (tmpY == 0.0)
                 {
                     idxSpline+=lX;
                     continue;
                 }
 
                 double xContrib=0.0;
                 for(int controlIdxX = -1; controlIdxX < (lX - 1); ++controlIdxX)
                 {
                     double tmpX = Bspline03((j / hX) - controlIdxX);
                     xContrib+=VEC_ELEM(cij,idxSpline) * tmpX;// *tmpY;
                     idxSpline+=1;
                 }
                 MAT_ELEM(thresholdMatrix,i,j)+=xContrib*tmpY;
             }
         }
     }
 }

◆ lowestResolutionbyPercentile() [1/2]

void ProgMonoTomo::lowestResolutionbyPercentile	(	MultidimArray< double > &	resolutionVol,
		std::vector< double > &	list,
		double &	cut_value,
		double &	resolutionThreshold
	)

Definition at line 623 of file resolution_monotomo.cpp.

 {
     double last_resolution_2 = list[(list.size()-1)];
 
     double lowest_res;
     lowest_res = list[0];
 
     // Count number of voxels with resolution
 
     double rVol;
     std::vector<double> resolVec(0);
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(resolutionVol)
     {
         rVol = DIRECT_MULTIDIM_ELEM(resolutionVol, n);
         if (rVol>=(last_resolution_2-0.001))//&& (DIRECT_MULTIDIM_ELEM(resolutionVol, n)<lowest_res) ) //the value 0.001 is a tolerance
         {
             resolVec.push_back(rVol);
         }
 
     }
     size_t N;
     N = resolVec.size();
 
     std::sort(resolVec.begin(), resolVec.end());
 
     resolutionThreshold = resolVec[(int)((0.95)*N)];
 
     std::cout << "resolutionThreshold = " << resolutionThreshold <<  std::endl;
 }

◆ lowestResolutionbyPercentile() [2/2]

void ProgMonoTomo::lowestResolutionbyPercentile	(	MultidimArray< float > &	resolutionVol,
		std::vector< float > &	list,
		float &	cut_value,
		float &	resolutionThreshold
	)

Definition at line 535 of file resolution_monotomo.cpp.

 {
     double last_resolution_2 = list[(list.size()-1)];
 
     double lowest_res;
     lowest_res = list[0];
 
     // Count number of voxels with resolution
 
     double rVol;
     std::vector<double> resolVec(0);
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(resolutionVol)
     {
         rVol = DIRECT_MULTIDIM_ELEM(resolutionVol, n);
         if (rVol>=(last_resolution_2-0.001))//&& (DIRECT_MULTIDIM_ELEM(resolutionVol, n)<lowest_res) ) //the value 0.001 is a tolerance
         {
             resolVec.push_back(rVol);
         }
 
     }
     size_t N;
     N = resolVec.size();
 
     std::sort(resolVec.begin(), resolVec.end());
 
     resolutionThreshold = resolVec[(int)((0.95)*N)];
 
     std::cout << "resolutionThreshold = " << resolutionThreshold <<  std::endl;
 }

◆ median3x3x3() [1/2]

void ProgMonoTomo::median3x3x3	(	MultidimArray< double >	vol,
		MultidimArray< double > &	filtered
	)

◆ median3x3x3() [2/2]

void ProgMonoTomo::median3x3x3	(	MultidimArray< double >	vol,
		MultidimArray< double > &	filtered
	)

◆ postProcessingLocalResolutions() [1/2]

void ProgMonoTomo::postProcessingLocalResolutions	(	MultidimArray< double > &	resolutionVol,
		std::vector< double > &	list
	)

Definition at line 587 of file resolution_monotomo.cpp.

 {
     MultidimArray<double> resolutionVol_aux = resolutionVol;
     double init_res, last_res;
 
     init_res = list[0];
     last_res = list[(list.size()-1)];
     
     double last_resolution_2 = list[last_res];
 
     double lowest_res;
     lowest_res = list[1]; //Example resolutions between 10-300, list(0)=300, list(1)=290, it is used list(1) due to background
     //is at 300 and the smoothing cast values of 299 and they must be removed.
 
     // Count number of voxels with resolution
     std::vector<double> resolVec(0);
     double rVol;
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(resolutionVol)
     {
         rVol = DIRECT_MULTIDIM_ELEM(resolutionVol, n);
         if ( (rVol>=(last_resolution_2-0.001)) && (rVol<=lowest_res) ) //the value 0.001 is a tolerance
         {
             resolVec.push_back(rVol);
         }
     }
 
     size_t N;
     N = resolVec.size();
     std::sort(resolVec.begin(), resolVec.end());
 
     std::cout << "median Resolution = " << resolVec[(int)(0.5*N)] << std::endl;
 }

◆ postProcessingLocalResolutions() [2/2]

void ProgMonoTomo::postProcessingLocalResolutions	(	MultidimArray< float > &	resolutionVol,
		std::vector< float > &	list
	)

Definition at line 499 of file resolution_monotomo.cpp.

 {
     MultidimArray<float> resolutionVol_aux = resolutionVol;
     float init_res, last_res;
 
     init_res = list[0];
     last_res = list[(list.size()-1)];
     
     auto last_resolution_2 = list[last_res];
 
     double lowest_res;
     lowest_res = list[1]; //Example resolutions between 10-300, list(0)=300, list(1)=290, it is used list(1) due to background
     //is at 300 and the smoothing cast values of 299 and they must be removed.
 
     // Count number of voxels with resolution
     std::vector<float> resolVec(0);
     float rVol;
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(resolutionVol)
     {
         rVol = DIRECT_MULTIDIM_ELEM(resolutionVol, n);
         if ( (rVol>=(last_resolution_2-0.001)) && (rVol<=lowest_res) ) //the value 0.001 is a tolerance
         {
             resolVec.push_back(rVol);
         }
     }
 
     size_t N;
     N = resolVec.size();
     std::sort(resolVec.begin(), resolVec.end());
 
     std::cout << "median Resolution = " << resolVec[(int)(0.5*N)] << std::endl;
 }

◆ produceSideInfo() [1/2]

void ProgMonoTomo::produceSideInfo ( )

Definition at line 78 of file resolution_monotomo.cpp.

 {
     std::cout << "Starting..." << std::endl;
     std::cout << "           " << std::endl;
     std::cout << "IMPORTANT: If the angular step of the tilt series is higher than 3 degrees"<< std::endl;
     std::cout << "           then, the tomogram is not properly for MonoTomo. Despite this is not "<< std::endl;
     std::cout << "           optimal, MonoTomo will try to compute the local resolution." << std::endl;
     std::cout << "           " << std::endl;
 
     Image<double> V;
     Image<double> V1, V2;
     V1.read(fnVol);
     V2.read(fnVol2);
     V()=0.5*(V1()+V2());
     V.write(fnMeanVol);
 
     V().setXmippOrigin();
 
     transformer_inv.setThreadsNumber(nthrs);
 
     FourierTransformer transformer;
     MultidimArray<double> &inputVol = V();
     VRiesz.resizeNoCopy(inputVol);
 
     #ifdef TEST_FRINGES
 
     double modulus, xx, yy, zz;
 
     long nnn=0;
     for(size_t k=0; k<ZSIZE(inputVol); ++k)
     {
         zz = (k-ZSIZE(inputVol)*0.5)*(k-ZSIZE(inputVol)*0.5);
         for(size_t i=0; i<YSIZE(inputVol); ++i)
         {
             yy = (i-YSIZE(inputVol)*0.5)*(i-YSIZE(inputVol)*0.5);
             for(size_t j=0; j<XSIZE(inputVol); ++j)
             {
                 xx = (j-XSIZE(inputVol)*0.5)*(j-XSIZE(inputVol)*0.5);
                 DIRECT_MULTIDIM_ELEM(inputVol,nnn) = cos(0.1*sqrt(xx+yy+zz));
                 ++nnn;
             }
         }
     }
 
     Image<double> saveiu;
     saveiu = inputVol;
     saveiu.write("franjas.vol");
     exit(0);
     #endif
 
 
     transformer.FourierTransform(inputVol, fftV);
     iu.initZeros(fftV);
 
     // Calculate u and first component of Riesz vector
     double uz, uy, ux, uz2, u2, uz2y2;
     long n=0;
     for(size_t k=0; k<ZSIZE(fftV); ++k)
     {
         FFT_IDX2DIGFREQ(k,ZSIZE(inputVol),uz);
         uz2=uz*uz;
 
         for(size_t i=0; i<YSIZE(fftV); ++i)
         {
             FFT_IDX2DIGFREQ(i,YSIZE(inputVol),uy);
             uz2y2=uz2+uy*uy;
 
             for(size_t j=0; j<XSIZE(fftV); ++j)
             {
                 FFT_IDX2DIGFREQ(j,XSIZE(inputVol),ux);
                 u2=uz2y2+ux*ux;
                 if ((k != 0) || (i != 0) || (j != 0))
                     DIRECT_MULTIDIM_ELEM(iu,n) = 1.0/sqrt(u2);
                 else
                     DIRECT_MULTIDIM_ELEM(iu,n) = 1e38;
                 ++n;
             }
         }
     }
     #ifdef DEBUG
     Image<double> saveiu;
     saveiu = 1/iu;
     saveiu.write("iu.vol");
     #endif
 
     // Prepare low pass filter
     lowPassFilter.FilterShape = RAISED_COSINE;
     lowPassFilter.raised_w = 0.01;
     lowPassFilter.do_generate_3dmask = false;
     lowPassFilter.FilterBand = LOWPASS;
 
     // Prepare mask
     MultidimArray<int> &pMask=mask();
 
     if (fnMask != "")
     {
         mask.read(fnMask);
         mask().setXmippOrigin();
     }
     else
     {
         size_t Zdim, Ydim, Xdim, Ndim;
         inputVol.getDimensions(Xdim, Ydim, Zdim, Ndim);
         mask().resizeNoCopy(Ndim, Zdim, Ydim, Xdim);
         mask().initConstant(1);
     }
 
     NVoxelsOriginalMask = 0;
 
     FOR_ALL_ELEMENTS_IN_ARRAY3D(pMask)
     {
         if (A3D_ELEM(pMask, k, i, j) == 1)
             ++NVoxelsOriginalMask;
 //      if (i*i+j*j+k*k > R*R)
 //          A3D_ELEM(pMask, k, i, j) = -1;
     }
 
     #ifdef DEBUG_MASK
     mask.write("mask.vol");
     #endif
 
 
     V1.read(fnVol);
     V2.read(fnVol2);
 
     V1()-=V2();
     V1()/=2;
 
     fftN=new MultidimArray< std::complex<double> >;
     FourierTransformer transformer2;
 
     #ifdef DEBUG
       V1.write("diff_volume.vol");
     #endif
 
     int N_smoothing = 10;
 
     int siz_z = ZSIZE(inputVol)*0.5;
     int siz_y = YSIZE(inputVol)*0.5;
     int siz_x = XSIZE(inputVol)*0.5;
 
 
     int limit_distance_x = (siz_x-N_smoothing);
     int limit_distance_y = (siz_y-N_smoothing);
     int limit_distance_z = (siz_z-N_smoothing);
 
     n=0;
     for(int k=0; k<ZSIZE(inputVol); ++k)
     {
         uz = (k - siz_z);
         for(int i=0; i<YSIZE(inputVol); ++i)
         {
             uy = (i - siz_y);
             for(int j=0; j<XSIZE(inputVol); ++j)
             {
                 ux = (j - siz_x);
 
                 if (abs(ux)>=limit_distance_x)
                 {
                     DIRECT_MULTIDIM_ELEM(V1(), n) *= 0.5*(1+cos(PI*(limit_distance_x - abs(ux))/N_smoothing));
                     DIRECT_MULTIDIM_ELEM(pMask, n) = 0;
                 }
                 if (abs(uy)>=limit_distance_y)
                 {
                     DIRECT_MULTIDIM_ELEM(V1(), n) *= 0.5*(1+cos(PI*(limit_distance_y - abs(uy))/N_smoothing));
                     DIRECT_MULTIDIM_ELEM(pMask, n) = 0;
                 }
                 if (abs(uz)>=limit_distance_z)
                 {
                     DIRECT_MULTIDIM_ELEM(V1(), n) *= 0.5*(1+cos(PI*(limit_distance_z - abs(uz))/N_smoothing));
                     DIRECT_MULTIDIM_ELEM(pMask, n) = 0;
                 }
                 ++n;
             }
         }
     }
 
 
     transformer2.FourierTransform(V1(), *fftN);
 
     V.clear();
 
     double u;
 
     freq_fourier_z.initZeros(ZSIZE(fftV));
     freq_fourier_x.initZeros(XSIZE(fftV));
     freq_fourier_y.initZeros(YSIZE(fftV));
 
     VEC_ELEM(freq_fourier_z,0) = 1e-38;
     for(size_t k=0; k<ZSIZE(fftV); ++k)
     {
         FFT_IDX2DIGFREQ(k,ZSIZE(pMask), u);
         VEC_ELEM(freq_fourier_z,k) = u;
     }
 
     VEC_ELEM(freq_fourier_y,0) = 1e-38;
     for(size_t k=0; k<YSIZE(fftV); ++k)
     {
         FFT_IDX2DIGFREQ(k,YSIZE(pMask), u);
         VEC_ELEM(freq_fourier_y,k) = u;
     }
 
     VEC_ELEM(freq_fourier_x,0) = 1e-38;
     for(size_t k=0; k<XSIZE(fftV); ++k)
     {
         FFT_IDX2DIGFREQ(k,XSIZE(pMask), u);
         VEC_ELEM(freq_fourier_x,k) = u;
     }
 
 }

◆ produceSideInfo() [2/2]

void ProgMonoTomo::produceSideInfo ( )

◆ readParams() [1/2]

void ProgMonoTomo::readParams ( )

virtual

Function in which each program will read parameters that it need. If some error occurs the usage will be printed out.

Reimplemented from XmippProgram.

Definition at line 36 of file resolution_monotomo.cpp.

 {
     fnVol = getParam("--vol");
     fnVol2 = getParam("--vol2");
     fnMeanVol = getParam("--meanVol");
     fnOut = getParam("-o");
     fnFilt = getParam("--filteredMap");
     fnMask = getParam("--mask");
     sampling = getDoubleParam("--sampling_rate");
     fnmaskWedge = getParam("--maskWedge");
     minRes = getDoubleParam("--minRes");
     maxRes = getDoubleParam("--maxRes");
     freq_step = getDoubleParam("--step");
     trimBound = getDoubleParam("--trimmed");
     significance = getDoubleParam("--significance");
     nthrs = getIntParam("--threads");
 }

◆ readParams() [2/2]

void ProgMonoTomo::readParams ( )

virtual

Function in which each program will read parameters that it need. If some error occurs the usage will be printed out.

Reimplemented from XmippProgram.

◆ resolution2eval() [1/2]

void ProgMonoTomo::resolution2eval	(	int &	count_res,
		double	step,
		double &	resolution,
		double &	last_resolution,
		double &	freq,
		double &	freqL,
		int &	last_fourier_idx,
		bool &	continueIter,
		bool &	breakIter
	)

Definition at line 655 of file resolution_monotomo.cpp.

 {
     resolution = maxRes - count_res*step;
     freq = sampling/resolution;
     ++count_res;
 
     double Nyquist = 2*sampling;
     double aux_frequency;
     int fourier_idx;
 
     DIGFREQ2FFT_IDX(freq, ZSIZE(VRiesz), fourier_idx);
 
     FFT_IDX2DIGFREQ(fourier_idx, ZSIZE(VRiesz), aux_frequency);
 
     freq = aux_frequency;
 
     if (fourier_idx == last_fourier_idx)
     {
         continueIter = true;
         return;
     }
 
     last_fourier_idx = fourier_idx;
     resolution = sampling/aux_frequency;
 
 
     if (count_res == 0)
         last_resolution = resolution;
 
     if (resolution<Nyquist)
     {
         breakIter = true;
         return;
     }
 
     freqL = sampling/(resolution + step);
 
     int fourier_idx_2;
 
     DIGFREQ2FFT_IDX(freqL, ZSIZE(VRiesz), fourier_idx_2);
 
     if (fourier_idx_2 == fourier_idx)
     {
         if (fourier_idx > 0){
             FFT_IDX2DIGFREQ(fourier_idx - 1, ZSIZE(VRiesz), freqL);
         }
         else{
             freqL = sampling/(resolution + step);
         }
     }
 
 }

◆ resolution2eval() [2/2]

void ProgMonoTomo::resolution2eval	(	int &	count_res,
		double	step,
		double &	resolution,
		double &	last_resolution,
		double &	freq,
		double &	freqL,
		int &	last_fourier_idx,
		bool &	continueIter,
		bool &	breakIter
	)

◆ run() [1/2]

void ProgMonoTomo::run ( )

virtual

This function will be start running the program. it also should be implemented by derived classes.

Reimplemented from XmippProgram.

Definition at line 713 of file resolution_monotomo.cpp.

 {
     produceSideInfo();
 
     Image<double> outputResolution;
 
     outputResolution().resizeNoCopy(VRiesz);
     outputResolution().initConstant(maxRes);
 
     MultidimArray<int> &pMask = mask();
     MultidimArray<double> &pOutputResolution = outputResolution();
     MultidimArray<double> &pVfiltered = Vfiltered();
     MultidimArray<double> &pVresolutionFiltered = VresolutionFiltered();
     MultidimArray<double> amplitudeMS, amplitudeMN;
 
     double criticalZ=icdf_gauss(significance);
     double criticalW=-1;
     double resolution, resolution_2, last_resolution = 10000;  //A huge value for achieving
                                                 //last_resolution < resolution
     double freq, freqH, freqL;
     double max_meanS = -1e38;
     double cut_value = 0.025;
     int boxsize = 50;
 
 
     double R_ = freq_step;
 
     if (R_<0.25)
         R_=0.25;
 
     double Nyquist = 2*sampling;
     if (minRes<2*sampling)
         minRes = Nyquist;
 
     bool doNextIteration=true;
 
     bool lefttrimming = false;
     int last_fourier_idx = -1;
 
     int count_res = 0;
     FileName fnDebug;
 
     int iter=0;
     std::vector<double> list;
 
     std::cout << "Analyzing frequencies" << std::endl;
     std::cout << "                     " << std::endl;
     std::vector<double> noiseValues;
 
     int xdim = XSIZE(pOutputResolution);
     int ydim = YSIZE(pOutputResolution);
 
     do
     {
         bool continueIter = false;
         bool breakIter = false;
 
         resolution2eval(count_res, R_,
                         resolution, last_resolution,
                         freq, freqH,
                         last_fourier_idx, continueIter, breakIter);
 
         if (continueIter)
             continue;
 
         if (breakIter)
             break;
 
         std::cout << "resolution = " << resolution << std::endl;
 
 
         list.push_back(resolution);
 
         if (iter <2)
             resolution_2 = list[0];
         else
             resolution_2 = list[iter - 2];
 
         fnDebug = "Signal";
 
         freqL = freq + 0.01;
 
         amplitudeMonogenicSignal3D(fftV, freq, freqH, freqL, amplitudeMS, iter, fnDebug);
         fnDebug = "Noise";
         amplitudeMonogenicSignal3D(*fftN, freq, freqH, freqL, amplitudeMN, iter, fnDebug);
 
         Matrix2D<double> noiseMatrix;
 
         Matrix2D<double> thresholdMatrix;
         localNoise(amplitudeMN, noiseMatrix, boxsize, thresholdMatrix);
 
 
         double sumS=0, sumS2=0, sumN=0, sumN2=0, NN = 0, NS = 0;
         noiseValues.clear();
 
 
         FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(amplitudeMS)
         {
             if (DIRECT_MULTIDIM_ELEM(pMask, n)>=1)
             {
                 double amplitudeValue=DIRECT_MULTIDIM_ELEM(amplitudeMS, n);
                 double amplitudeValueN=DIRECT_MULTIDIM_ELEM(amplitudeMN, n);
                 sumS  += amplitudeValue;
                 noiseValues.push_back(amplitudeValueN);
                 sumN  += amplitudeValueN;
                 ++NS;
                 ++NN;
             }
         }
     
         #ifdef DEBUG
         std::cout << "NS" << NS << std::endl;
         std::cout << "NVoxelsOriginalMask" << NVoxelsOriginalMask << std::endl;
         std::cout << "NS/NVoxelsOriginalMask = " << NS/NVoxelsOriginalMask << std::endl;
         #endif
         
 
             if (NS == 0)
             {
                 std::cout << "There are no points to compute inside the mask" << std::endl;
                 std::cout << "If the number of computed frequencies is low, perhaps the provided"
                         "mask is not enough tight to the volume, in that case please try another mask" << std::endl;
                 break;
             }
 
             double meanS=sumS/NS;
 
             #ifdef DEBUG
               std::cout << "Iteration = " << iter << ",   Resolution= " << resolution <<
                       ",   Signal = " << meanS << ",   Noise = " << meanN << ",  Threshold = "
                       << thresholdNoise <<std::endl;
             #endif
 
 
             if (meanS>max_meanS)
                 max_meanS = meanS;
 
             if (meanS<0.001*max_meanS)
             {
                 std::cout << "Search of resolutions stopped due to too low signal" << std::endl;
                 break;
             }
 
             FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(amplitudeMS)
             {
                 if (DIRECT_A3D_ELEM(pMask, k,i,j)>=1)
                 {
                     if ( DIRECT_A3D_ELEM(amplitudeMS, k,i,j)>MAT_ELEM(thresholdMatrix, i, j) )
                     {
 
                         DIRECT_A3D_ELEM(pMask,  k,i,j) = 1;
                         DIRECT_A3D_ELEM(pOutputResolution, k,i,j) = resolution;
                     }
                     else{
                         DIRECT_A3D_ELEM(pMask,  k,i,j) += 1;
                         if (DIRECT_A3D_ELEM(pMask,  k,i,j) >2)
                         {
                             DIRECT_A3D_ELEM(pMask,  k,i,j) = -1;
                             DIRECT_A3D_ELEM(pOutputResolution,  k,i,j) = resolution_2;
                         }
                     }
                 }
             }
 
 
 
             #ifdef DEBUG_MASK
             FileName fnmask_debug;
             fnmask_debug = formatString("maske_%i.vol", iter);
             mask.write(fnmask_debug);
             #endif
 
 
             if (doNextIteration)
             {
                 if (resolution <= (minRes-0.001))
                     doNextIteration = false;
             }
 
 //      }
         iter++;
         last_resolution = resolution;
     } while (doNextIteration);
 
     Image<double> outputResolutionImage2;
     outputResolutionImage2() = pOutputResolution;
     outputResolutionImage2.write("resultado.vol");
 
 
     amplitudeMN.clear();
     amplitudeMS.clear();
 
     //Convolution with a real gaussian to get a smooth map
     MultidimArray<double> FilteredResolution = pOutputResolution;
     double sigma = 25.0;
 
     double resolutionThreshold;
 
     sigma = 3;
 
     realGaussianFilter(FilteredResolution, sigma);
 
 
     lowestResolutionbyPercentile(FilteredResolution, list, cut_value, resolutionThreshold);
 
 
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(FilteredResolution)
     {
         if ( (DIRECT_MULTIDIM_ELEM(FilteredResolution, n)<resolutionThreshold) && (DIRECT_MULTIDIM_ELEM(FilteredResolution, n)>DIRECT_MULTIDIM_ELEM(pOutputResolution, n)) )
             DIRECT_MULTIDIM_ELEM(FilteredResolution, n) = DIRECT_MULTIDIM_ELEM(pOutputResolution, n);
         if ( DIRECT_MULTIDIM_ELEM(FilteredResolution, n)<Nyquist)
             DIRECT_MULTIDIM_ELEM(FilteredResolution, n) = Nyquist;
     }
 
     realGaussianFilter(FilteredResolution, sigma);
 
 
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(FilteredResolution)
     {
         if ((DIRECT_MULTIDIM_ELEM(FilteredResolution, n)<resolutionThreshold)  && (DIRECT_MULTIDIM_ELEM(FilteredResolution, n)>DIRECT_MULTIDIM_ELEM(pOutputResolution, n)))
             DIRECT_MULTIDIM_ELEM(FilteredResolution, n) = DIRECT_MULTIDIM_ELEM(pOutputResolution, n);
         if ( DIRECT_MULTIDIM_ELEM(FilteredResolution, n)<Nyquist)
             DIRECT_MULTIDIM_ELEM(FilteredResolution, n) = Nyquist;
     }
     Image<double> outputResolutionImage;
     MultidimArray<double> resolutionFiltered, resolutionChimera;
 
     postProcessingLocalResolutions(FilteredResolution, list);
     outputResolutionImage() = FilteredResolution;
     outputResolutionImage.write(fnOut);
 }

◆ run() [2/2]

void ProgMonoTomo::run ( )

virtual

This function will be start running the program. it also should be implemented by derived classes.

Reimplemented from XmippProgram.

◆ smoothBorders()

void ProgMonoTomo::smoothBorders	(	MultidimArray< float > &	vol,
		MultidimArray< int > &	pMask
	)

Definition at line 139 of file resolution_monotomo.cpp.

 {
     int N_smoothing = 10;
 
     int siz_z = zdim*0.5;
     int siz_y = ydim*0.5;
     int siz_x = xdim*0.5;
 
 
     int limit_distance_x = (siz_x-N_smoothing);
     int limit_distance_y = (siz_y-N_smoothing);
     int limit_distance_z = (siz_z-N_smoothing);
 
     long n=0;
     for(int k=0; k<zdim; ++k)
     {
         auto uz = (k - siz_z);
         for(int i=0; i<ydim; ++i)
         {
             auto uy = (i - siz_y);
             for(int j=0; j<xdim; ++j)
             {
                 auto ux = (j - siz_x);
 
                 if (abs(ux)>=limit_distance_x)
                 {
                     DIRECT_MULTIDIM_ELEM(vol, n) *= 0.5*(1+std::cos(PI*(limit_distance_x - std::abs(ux))/N_smoothing));
                     DIRECT_MULTIDIM_ELEM(pMask, n) = 0;
                 }
                 if (abs(uy)>=limit_distance_y)
                 {
                     DIRECT_MULTIDIM_ELEM(vol, n) *= 0.5*(1+std::cos(PI*(limit_distance_y - std::abs(uy))/N_smoothing));
                     DIRECT_MULTIDIM_ELEM(pMask, n) = 0;
                 }
                 if (abs(uz)>=limit_distance_z)
                 {
                     DIRECT_MULTIDIM_ELEM(vol, n) *= 0.5*(1+std::cos(PI*(limit_distance_z - std::abs(uz))/N_smoothing));
                     DIRECT_MULTIDIM_ELEM(pMask, n) = 0;
                 }
                 ++n;
             }
         }
     }
 }

Member Data Documentation

◆ automaticMode

bool ProgMonoTomo::automaticMode

Definition at line 64 of file resolution_monotomo.h.

◆ backward_transformer

std::unique_ptr<AFT<float> > ProgMonoTomo::backward_transformer

Definition at line 113 of file resolution_monotomo.h.

◆ fftN

MultidimArray< std::complex<double> > * ProgMonoTomo::fftN

Definition at line 105 of file resolution_monotomo.h.

◆ fftV

MultidimArray< std::complex<double> > ProgMonoTomo::fftV

Definition at line 105 of file resolution_monotomo.h.

◆ fftVRiesz [1/2]

MultidimArray< std::complex<double> > ProgMonoTomo::fftVRiesz

Definition at line 107 of file resolution_monotomo.h.

◆ fftVRiesz [2/2]

std::vector<std::complex<float> > ProgMonoTomo::fftVRiesz

Definition at line 110 of file resolution_monotomo.h.

◆ fftVRiesz_aux [1/2]

MultidimArray< std::complex<double> > ProgMonoTomo::fftVRiesz_aux

Definition at line 107 of file resolution_monotomo.h.

◆ fftVRiesz_aux [2/2]

std::vector<std::complex<float> > ProgMonoTomo::fftVRiesz_aux

Definition at line 110 of file resolution_monotomo.h.

◆ FilterBand

FourierFilter ProgMonoTomo::FilterBand

Definition at line 108 of file resolution_monotomo.h.

◆ fnchim

FileName ProgMonoTomo::fnchim

Definition at line 51 of file resolution_monotomo.h.

◆ fnFilt

FileName ProgMonoTomo::fnFilt

Definition at line 51 of file resolution_monotomo.h.

◆ fnMask

FileName ProgMonoTomo::fnMask

Definition at line 51 of file resolution_monotomo.h.

◆ fnMaskOut

FileName ProgMonoTomo::fnMaskOut

Definition at line 51 of file resolution_monotomo.h.

◆ fnmaskWedge

FileName ProgMonoTomo::fnmaskWedge

Definition at line 51 of file resolution_monotomo.h.

◆ fnMd

FileName ProgMonoTomo::fnMd

Definition at line 51 of file resolution_monotomo.h.

◆ fnMeanVol

FileName ProgMonoTomo::fnMeanVol

Definition at line 51 of file resolution_monotomo.h.

◆ fnOut

FileName ProgMonoTomo::fnOut

Filenames

Definition at line 51 of file resolution_monotomo.h.

◆ fnSpatial

FileName ProgMonoTomo::fnSpatial

Definition at line 51 of file resolution_monotomo.h.

◆ fnVol

FileName ProgMonoTomo::fnVol

Definition at line 51 of file resolution_monotomo.h.

◆ fnVol2

FileName ProgMonoTomo::fnVol2

Definition at line 51 of file resolution_monotomo.h.

◆ forward_transformer

std::unique_ptr<AFT<float> > ProgMonoTomo::forward_transformer

Definition at line 113 of file resolution_monotomo.h.

◆ fourierNoise

std::vector<std::complex<float> > ProgMonoTomo::fourierNoise

Definition at line 65 of file resolution_monotomo.h.

◆ fourierSignal

std::vector<std::complex<float> > ProgMonoTomo::fourierSignal

Definition at line 65 of file resolution_monotomo.h.

◆ freq_fourier_x

Matrix1D<double> ProgMonoTomo::freq_fourier_x

Definition at line 111 of file resolution_monotomo.h.

◆ freq_fourier_y

Matrix1D<double> ProgMonoTomo::freq_fourier_y

Definition at line 111 of file resolution_monotomo.h.

◆ freq_fourier_z

Matrix1D<double> ProgMonoTomo::freq_fourier_z

Definition at line 111 of file resolution_monotomo.h.

◆ freq_step

double ProgMonoTomo::freq_step

Step in digital frequency

Definition at line 61 of file resolution_monotomo.h.

◆ halfMapsGiven

bool ProgMonoTomo::halfMapsGiven

Definition at line 109 of file resolution_monotomo.h.

◆ iu

MultidimArray<double> ProgMonoTomo::iu

Definition at line 104 of file resolution_monotomo.h.

◆ lowPassFilter

FourierFilter ProgMonoTomo::lowPassFilter

Definition at line 108 of file resolution_monotomo.h.

◆ mask

Image< int > ProgMonoTomo::mask

Definition at line 103 of file resolution_monotomo.h.

◆ maskMatrix

Matrix2D< double > ProgMonoTomo::maskMatrix

Definition at line 112 of file resolution_monotomo.h.

◆ maxRes

double ProgMonoTomo::maxRes

Definition at line 55 of file resolution_monotomo.h.

◆ minRes

double ProgMonoTomo::minRes

Definition at line 55 of file resolution_monotomo.h.

◆ noiseOnlyInHalves

bool ProgMonoTomo::noiseOnlyInHalves

The search for resolutions is linear or inverse

Definition at line 64 of file resolution_monotomo.h.

◆ nthrs

int ProgMonoTomo::nthrs

Definition at line 58 of file resolution_monotomo.h.

◆ Nvoxels

int ProgMonoTomo::Nvoxels

Definition at line 58 of file resolution_monotomo.h.

◆ NVoxelsOriginalMask

int ProgMonoTomo::NVoxelsOriginalMask

Is the volume previously masked?

Definition at line 58 of file resolution_monotomo.h.

◆ R

double ProgMonoTomo::R

Definition at line 55 of file resolution_monotomo.h.

◆ resolutionMatrix

Matrix2D< double > ProgMonoTomo::resolutionMatrix

Definition at line 112 of file resolution_monotomo.h.

◆ resStep

double ProgMonoTomo::resStep

Step in digital frequency

Definition at line 60 of file resolution_monotomo.h.

◆ sampling

double ProgMonoTomo::sampling

sampling rate, minimum resolution, and maximum resolution

Definition at line 55 of file resolution_monotomo.h.

◆ significance

double ProgMonoTomo::significance

Definition at line 61 of file resolution_monotomo.h.

◆ transformer_inv

FourierTransformer ProgMonoTomo::transformer_inv

Definition at line 106 of file resolution_monotomo.h.

◆ trimBound

double ProgMonoTomo::trimBound

Definition at line 61 of file resolution_monotomo.h.

◆ Vfiltered

Image<double> ProgMonoTomo::Vfiltered

Definition at line 110 of file resolution_monotomo.h.

◆ VresolutionFiltered

Image<double> ProgMonoTomo::VresolutionFiltered

Definition at line 110 of file resolution_monotomo.h.

◆ VRiesz [1/2]

MultidimArray<double> ProgMonoTomo::VRiesz

Definition at line 104 of file resolution_monotomo.h.

◆ VRiesz [2/2]

MultidimArray<float> ProgMonoTomo::VRiesz

Definition at line 112 of file resolution_monotomo.h.

◆ xdim

size_t ProgMonoTomo::xdim

Definition at line 57 of file resolution_monotomo.h.

◆ xdimFT

size_t ProgMonoTomo::xdimFT

Definition at line 57 of file resolution_monotomo.h.

◆ ydim

size_t ProgMonoTomo::ydim

Definition at line 57 of file resolution_monotomo.h.

◆ ydimFT

size_t ProgMonoTomo::ydimFT

Definition at line 57 of file resolution_monotomo.h.

◆ zdim

size_t ProgMonoTomo::zdim

Definition at line 57 of file resolution_monotomo.h.

◆ zdimFT

size_t ProgMonoTomo::zdimFT

Definition at line 57 of file resolution_monotomo.h.

The documentation for this class was generated from the following files:

xmipp/legacy/libraries/tomo/resolution_monotomo.h
xmipp/legacy/libraries/tomo/resolution_monotomo.cpp

Public Member Functions

Public Attributes

Additional Inherited Members

Detailed Description

Member Function Documentation

◆ amplitudeMonogenicSignal3D() [1/2]

◆ amplitudeMonogenicSignal3D() [2/2]

◆ defineParams() [1/2]

◆ defineParams() [2/2]

◆ firstMonoResEstimation() [1/2]

◆ firstMonoResEstimation() [2/2]

◆ gaussFilter()

◆ getFilteringResolution()

◆ localNoise() [1/2]

◆ localNoise() [2/2]

◆ lowestResolutionbyPercentile() [1/2]

◆ lowestResolutionbyPercentile() [2/2]

◆ median3x3x3() [1/2]

◆ median3x3x3() [2/2]

◆ postProcessingLocalResolutions() [1/2]

◆ postProcessingLocalResolutions() [2/2]

◆ produceSideInfo() [1/2]

◆ produceSideInfo() [2/2]

◆ readParams() [1/2]

◆ readParams() [2/2]

◆ resolution2eval() [1/2]

◆ resolution2eval() [2/2]

◆ run() [1/2]

◆ run() [2/2]

◆ smoothBorders()

Member Data Documentation

◆ automaticMode

◆ backward_transformer

◆ fftN

◆ fftV

◆ fftVRiesz [1/2]

◆ fftVRiesz [2/2]

◆ fftVRiesz_aux [1/2]

◆ fftVRiesz_aux [2/2]

◆ FilterBand

◆ fnchim

◆ fnFilt

◆ fnMask

◆ fnMaskOut

◆ fnmaskWedge

◆ fnMd

◆ fnMeanVol

◆ fnOut

◆ fnSpatial

◆ fnVol

◆ fnVol2

◆ forward_transformer

◆ fourierNoise

◆ fourierSignal

◆ freq_fourier_x

◆ freq_fourier_y

◆ freq_fourier_z

◆ freq_step

◆ halfMapsGiven

◆ iu

◆ lowPassFilter

◆ mask

◆ maskMatrix

◆ maxRes

◆ minRes

◆ noiseOnlyInHalves

◆ nthrs

◆ Nvoxels

◆ NVoxelsOriginalMask

◆ R

◆ resolutionMatrix

◆ resStep

◆ sampling

◆ significance

◆ transformer_inv

◆ trimBound

◆ Vfiltered

◆ VresolutionFiltered

◆ VRiesz [1/2]

◆ VRiesz [2/2]