#include <subtomo_subtraction.h>

Inheritance diagram for ProgSubtomoSubtraction:

Collaboration diagram for ProgSubtomoSubtraction:

Classes
struct	Angles

Public Member Functions
void	POCSmask (const MultidimArray< double > &, MultidimArray< double > &)
	Processing methods. More...

void	POCSnonnegative (MultidimArray< double > &)

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

void	POCSFourierAmplitudeRadAvg (MultidimArray< std::complex< double >> &, double, const MultidimArray< double > &, int, int, int)

void	POCSMinMax (MultidimArray< double > &, double, double)

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

void	extractPhase (MultidimArray< std::complex< double >> &) const

void	computeEnergy (MultidimArray< double > &, const MultidimArray< double > &) const

void	centerFFTMagnitude (MultidimArray< double > &, MultidimArray< std::complex< double >> &, MultidimArray< double > &) const

void	radialAverage (const MultidimArray< double > &, const MultidimArray< double > &, MultidimArray< double > &)

MultidimArray< double >	computeRadQuotient (const MultidimArray< double > &, const MultidimArray< double > &, const MultidimArray< double > &, const MultidimArray< double > &)

void	createFilter (FourierFilter &, double)

Image< double >	subtraction (Image< double >, Image< double > &, const MultidimArray< double > &, const FileName &, const FileName &, FourierFilter &, double)

MultidimArray< double >	computeMagnitude (MultidimArray< double > &)

MultidimArray< double >	createMask (const Image< double > &, const FileName &, const FileName &)

void	filterMask (MultidimArray< double > &) const

MultidimArray< std::complex< double > >	computePhase (MultidimArray< double > &)

MultidimArray< double >	getSubtractionMask (const FileName &, MultidimArray< double >)

	ProgSubtomoSubtraction ()
	Empty constructor. More...

void	readParams () override
	Read argument from command line. More...

void	show () const override
	Show. More...

void	defineParams () override
	Define parameters. More...

void	readParticle (const MDRow &rowIn)
	Read and write methods. More...

void	writeParticle (MDRow &rowOut, FileName, Image< double > &)

void	preProcess () override
	MPI methods. More...

void	processImage (const FileName &fnImg, const FileName &fnImgOut, const MDRow &rowIn, MDRow &rowOut) override

void	postProcess () override

Public Member Functions inherited from XmippMetadataProgram
MetaData *	getInputMd ()

MetaDataVec &	getOutputMd ()

	XmippMetadataProgram ()
	Empty constructor. More...

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

virtual void	init ()

virtual void	setup (MetaData *md, const FileName &o="", const FileName &oroot="", bool applyGeo=false, MDLabel label=MDL_IMAGE)

virtual	~XmippMetadataProgram ()

void	setMode (WriteModeMetaData _mode)

void	setupRowOut (const FileName &fnImgIn, const MDRow &rowIn, const FileName &fnImgOut, MDRow &rowOut) const
	Prepare rowout. More...

virtual void	wait ()
	Wait for the distributor to finish. More...

virtual void	checkPoint ()
	For very long programs, it may be needed to write checkpoints. More...

virtual void	run ()
	Run over all images. More...

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)

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	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	fnVolMd

FileName	fnVolRef

FileName	fnVol2

FileName	fnOut

FileName	fnMask1

FileName	fnMask2

FileName	fnMaskSub

FileName	fnVol1F

FileName	fnVol2A

bool	computeE

size_t	iter

double	v1min

double	v1max

bool	performSubtraction

int	sigma

double	cutFreq

double	lambda

bool	radavg

bool	subtomos

MetaDataVec	mdParticles

MDRowVec	row

Matrix1D< double >	roffset

struct Angles	part_angles

Image< double >	V

Image< double >	V1

Image< double >	padv

double	pad

double	std1

size_t	n

FourierTransformer	transformer2

MultidimArray< std::complex< double > >	V2Fourier

MultidimArray< double >	V1FourierMag

MultidimArray< double >	mask

Image< double >	Vdiff

FourierFilter	filter2

int	rank

Public Attributes inherited from XmippMetadataProgram
FileName	fn_in
	Filenames of input and output Metadata. More...

FileName	fn_out

FileName	baseName

FileName	pathBaseName

FileName	oextBaseName

bool	apply_geo
	Apply geo. More...

size_t	ndimOut
	Output dimensions. More...

size_t	zdimOut

size_t	ydimOut

size_t	xdimOut

DataType	datatypeOut

size_t	mdInSize
	Number of input elements. More...

Public Attributes inherited from XmippProgram
bool	doRun

bool	runWithoutArgs

int	verbose
	Verbosity level. More...

int	debug

Additional Inherited Members
Protected Member Functions inherited from XmippMetadataProgram
virtual void	initComments ()

virtual bool	getImageToProcess (size_t &objId, size_t &objIndex)

void	show () const override

virtual void	startProcessing ()

virtual void	finishProcessing ()

virtual void	writeOutput ()

virtual void	showProgress ()

virtual void	defineLabelParam ()

Protected Member Functions inherited from XmippProgram
void	defineCommons ()

Protected Attributes inherited from XmippMetadataProgram
WriteModeMetaData	mode
	Metadata writing mode: OVERWRITE, APPEND. More...

FileName	oext
	Output extension and root. More...

FileName	oroot

MDLabel	image_label
	MDLabel to be used to read/write images, usually will be MDL_IMAGE. More...

bool	produces_an_output
	Indicate that a unique final output is produced. More...

bool	produces_a_metadata
	Indicate that the unique final output file is a Metadata. More...

bool	each_image_produces_an_output
	Indicate that an output is produced for each image in the input. More...

bool	allow_apply_geo

bool	decompose_stacks
	Input Metadata will treat a stack file as a set of images instead of a unique file. More...

bool	delete_output_stack
	Delete previous output stack file prior to process images. More...

bool	get_image_info
	Get the input image file dimensions to further operations. More...

bool	save_metadata_stack
	Save the associated output metadata when output file is a stack. More...

bool	track_origin
	Include the original input image filename in the output stack. More...

bool	keep_input_columns
	Keep input metadata columns. More...

bool	remove_disabled
	Remove disabled images from the input selfile. More...

bool	allow_time_bar
	Show process time bar. More...

bool	input_is_metadata
	Input is a metadata. More...

bool	single_image
	Input is a single image. More...

bool	input_is_stack
	Input is a stack. More...

bool	output_is_stack
	Output is a stack. More...

bool	create_empty_stackfile

bool	delete_mdIn

size_t	time_bar_step
	Some time bar related counters. More...

size_t	time_bar_size

size_t	time_bar_done

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

Subtomogram subtraction

Definition at line 42 of file subtomo_subtraction.h.

Constructor & Destructor Documentation

◆ ProgSubtomoSubtraction()

ProgSubtomoSubtraction::ProgSubtomoSubtraction ( )

Empty constructor.

Definition at line 37 of file subtomo_subtraction.cpp.

 {
     produces_a_metadata = true;
     each_image_produces_an_output = true;
     keep_input_columns = true;
     save_metadata_stack = true;
     remove_disabled = false;
     rank = 0;
 }

Member Function Documentation

◆ centerFFTMagnitude()

void ProgSubtomoSubtraction::centerFFTMagnitude	(	MultidimArray< double > &	VolRad,
		MultidimArray< std::complex< double >> &	VolFourierRad,
		MultidimArray< double > &	VolFourierMagRad
	)		const

Definition at line 217 of file subtomo_subtraction.cpp.

                                                       {
     FourierTransformer transformerRad;
     transformerRad.completeFourierTransform(VolRad, VolFourierRad);
     CenterFFT(VolFourierRad, true);
     FFT_magnitude(VolFourierRad, VolFourierMagRad);
     VolFourierMagRad.setXmippOrigin();
 }

◆ computeEnergy()

void ProgSubtomoSubtraction::computeEnergy	(	MultidimArray< double > &	Vdif,
		const MultidimArray< double > &	Vact
	)		const

Definition at line 209 of file subtomo_subtraction.cpp.

                                                                                                               {
     Vdif = Vdif - Vact;
     double energy = 0;
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(Vdif)
     energy += DIRECT_MULTIDIM_ELEM(Vdif, n) * DIRECT_MULTIDIM_ELEM(Vdif, n);
     energy = sqrt(energy / MULTIDIM_SIZE(Vdif));
 }

◆ computeMagnitude()

MultidimArray< double > ProgSubtomoSubtraction::computeMagnitude ( MultidimArray< double > & volume )

Definition at line 312 of file subtomo_subtraction.cpp.

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

◆ computePhase()

MultidimArray< std::complex< double > > ProgSubtomoSubtraction::computePhase ( MultidimArray< double > & volume )

Definition at line 344 of file subtomo_subtraction.cpp.

                                                                                                     {
     MultidimArray<std::complex<double>> phase;
     transformer2.FourierTransform(volume, phase, true);
     extractPhase(phase);
     return phase;
 }

◆ computeRadQuotient()

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

Definition at line 263 of file subtomo_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 ProgSubtomoSubtraction::createFilter	(	FourierFilter &	filter2,
		double	cutFreq
	)

Definition at line 281 of file subtomo_subtraction.cpp.

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

◆ createMask()

MultidimArray< double > ProgSubtomoSubtraction::createMask	(	const Image< double > &	volume,
		const FileName &	fnM1,
		const FileName &	fnM2
	)

Definition at line 321 of file subtomo_subtraction.cpp.

                                                                                                                                 {
     MultidimArray<double> mask;
     if (fnM1 != "" && fnM2 != "") {
         Image<double> mask1;
         Image<double> mask2;
         mask1.read(fnM1);
         mask2.read(fnM2);
         mask = mask1() * mask2();
     } else {
         mask.resizeNoCopy(volume());
         mask.initConstant(1.0);
     }
     return mask;
 }

◆ defineParams()

void ProgSubtomoSubtraction::defineParams ( )

overridevirtual

Define parameters.

Reimplemented from XmippMetadataProgram.

Definition at line 48 of file subtomo_subtraction.cpp.

 {
     // Usage
     addUsageLine("This program modifies a volume as much as possible in order to assimilate it to another one,");
     addUsageLine("without loosing the important information in it process. Then, the subtraction of the two volumes");
     addUsageLine("can be optionally calculated. Sharpening: reference volume must be an atomic structure previously");
     addUsageLine("converted into a density map of the same specimen than in input volume 2.");
     // Parameters
     XmippMetadataProgram::defineParams();
     addParamsLine("--ref <volume>\t: Reference volume");
     addParamsLine("[--sub]\t: Perform the subtraction of the volumes. Output will be the difference");
     addParamsLine("[--sigma <s=3>]\t: Decay of the filter  (sigma) to smooth the mask transition");
     addParamsLine("[--iter <n=5>]\t: Number of iterations for the adjustment process");
     addParamsLine("[--mask1 <mask=\"\">]\t: Mask for volume 1");
     addParamsLine("[--mask2 <mask=\"\">]\t: Mask for volume 2");
     addParamsLine("[--maskSub <mask=\"\">]\t: Mask for subtraction region");
     addParamsLine("[--cutFreq <f=0>]\t: Filter both volumes with a filter which specified cutoff frequency "
             "(i.e. resolution inverse, <0.5)");
     addParamsLine("[--lambda <l=1>]\t: Relaxation factor for Fourier Amplitude POCS, i.e. 'how much modification "
             "of volume Fourier amplitudes', between 1 (full modification, recommended) and 0 (no modification)");
     addParamsLine("[--radavg]\t: Match the radially averaged Fourier amplitudes when adjusting the amplitudes instead "
             "of taking directly them from the reference volume");
     addParamsLine("[--computeEnergy]\t: Compute the energy difference between each step (energy difference gives "
             "information about the convergence of the adjustment process, while it can slightly slow the performance)");
     addParamsLine("[--saveV1 <structure=\"\"> ]\t: Save subtraction intermediate file (vol1 filtered) just when option "
             "--sub is passed, if not passed the input reference volume is not modified");
     addParamsLine("[--saveV2 <structure=\"\"> ]\t: Save subtraction intermediate file (vol2 adjusted) just when option "
             "--sub is passed, if not passed the output of the program is this file");
 }

◆ extractPhase()

void ProgSubtomoSubtraction::extractPhase ( MultidimArray< std::complex< double >> & FI ) const

Definition at line 201 of file subtomo_subtraction.cpp.

                                                                                     {
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(FI) {
         auto *ptr = (double *)&DIRECT_MULTIDIM_ELEM(FI, n);
         double phi = atan2(*(ptr + 1), *ptr);
         DIRECT_MULTIDIM_ELEM(FI, n) = std::complex<double>(cos(phi), sin(phi));
     }
 }

◆ filterMask()

void ProgSubtomoSubtraction::filterMask ( MultidimArray< double > & mask ) const

Definition at line 336 of file subtomo_subtraction.cpp.

                                                                         {
     FourierFilter Filter;
     Filter.FilterShape = REALGAUSSIAN;
     Filter.FilterBand = LOWPASS;
     Filter.w1 = sigma;
     Filter.applyMaskSpace(mask);
 }

◆ getSubtractionMask()

MultidimArray< double > ProgSubtomoSubtraction::getSubtractionMask	(	const FileName &	fnMSub,
		MultidimArray< double >	mask
	)

Definition at line 351 of file subtomo_subtraction.cpp.

                                                                                                                   {
     if (fnMSub.isEmpty()){
         filterMask(mask);
         return mask;
     }
     else {
         Image<double> masksub;
         masksub.read(fnMSub);
         return masksub();
     }
 }

◆ POCSFourierAmplitude()

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

Definition at line 142 of file subtomo_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 ProgSubtomoSubtraction::POCSFourierAmplitudeRadAvg	(	MultidimArray< std::complex< double >> &	V,
		double	l,
		const MultidimArray< double > &	rQ,
		int	V1size_x,
		int	V1size_y,
		int	V1size_z
	)

Definition at line 153 of file subtomo_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 ProgSubtomoSubtraction::POCSFourierPhase	(	const MultidimArray< std::complex< double >> &	phase,
		MultidimArray< std::complex< double >> &	FI
	)

Definition at line 193 of file subtomo_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 ProgSubtomoSubtraction::POCSmask	(	const MultidimArray< double > &	mask,
		MultidimArray< double > &	I
	)

Processing methods.

Definition at line 132 of file subtomo_subtraction.cpp.

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

◆ POCSMinMax()

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

Definition at line 183 of file subtomo_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 ProgSubtomoSubtraction::POCSnonnegative ( MultidimArray< double > & I )

Definition at line 137 of file subtomo_subtraction.cpp.

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

◆ postProcess()

void ProgSubtomoSubtraction::postProcess ( )

overridevirtual

Reimplemented from XmippMetadataProgram.

Definition at line 491 of file subtomo_subtraction.cpp.

 {
     getOutputMd().write(fn_out);
 }

◆ preProcess()

void ProgSubtomoSubtraction::preProcess ( )

overridevirtual

MPI methods.

Reimplemented from XmippMetadataProgram.

Definition at line 363 of file subtomo_subtraction.cpp.

                                         {
     show();
  }

◆ processImage()

void ProgSubtomoSubtraction::processImage	(	const FileName &	fnImg,
		const FileName &	fnImgOut,
		const MDRow &	rowIn,
		MDRow &	rowOut
	)

overridevirtual

Implements XmippMetadataProgram.

Definition at line 370 of file subtomo_subtraction.cpp.

                                                                                                                             {
     Image<double> V1;
     V1.read(fnVolRef);
     V1().setXmippOrigin();
     mask = createMask(V1, fnMask1, fnMask2);
     POCSmask(mask, V1());
     POCSnonnegative(V1());
     V1().computeDoubleMinMax(v1min, v1max);
     std1 = V1().computeStddev();
     createFilter(filter2, cutFreq);
     
     readParticle(rowIn);
     MultidimArray<double> &mv = V();
     mv.setXmippOrigin();
 
     // Window subtomo (padding) before apply alignment
     MultidimArray <double> &mpad = padv();
     mpad.setXmippOrigin();
     pad = XSIZE(mv);
     mv.window(mpad, STARTINGZ(mv)-(int)pad/2, STARTINGY(mv)-(int)pad/2, STARTINGX(mv)-(int)pad/2, FINISHINGZ(mv)+(int)pad/2, FINISHINGZ(mv)+(int)pad/2, FINISHINGZ(mv)+(int)pad/2);
     // Read alignment
     rowIn.getValueOrDefault(MDL_ANGLE_ROT, part_angles.rot, 0);
     rowIn.getValueOrDefault(MDL_ANGLE_TILT, part_angles.tilt, 0);
     rowIn.getValueOrDefault(MDL_ANGLE_PSI, part_angles.psi, 0);
     roffset.initZeros(3);
     rowIn.getValueOrDefault(MDL_SHIFT_X, roffset(0), 0);
     rowIn.getValueOrDefault(MDL_SHIFT_Y, roffset(1), 0);
     rowIn.getValueOrDefault(MDL_SHIFT_Z, roffset(2), 0);
     // Apply alignment
     Image<double> Vaux;
     MultidimArray<double> &mvaux = Vaux();
     mvaux.setXmippOrigin();
     Vaux = padv;
     mvaux.initZeros();
     Euler_rotate(mpad, part_angles.rot, part_angles.tilt, part_angles.psi, mvaux);
     selfTranslate(xmipp_transformation::LINEAR, mvaux, roffset, xmipp_transformation::WRAP);
     // Crop to restore original size
     mvaux.window(mv, STARTINGZ(mv), STARTINGY(mv), STARTINGX(mv), FINISHINGZ(mv), FINISHINGY(mv), FINISHINGX(mv));
     // Preprocessing
     POCSmask(mask, mv);
     POCSnonnegative(mv);
     Vdiff = V;
     auto V2FourierPhase = computePhase(mv);
     auto V1FourierMag = computeMagnitude(V1());
     auto V2FourierMag = computeMagnitude(mv);
     auto radQuotient = computeRadQuotient(V1FourierMag, V2FourierMag, V1(), mv);
 
     for (n = 0; n < iter; ++n) 
     {
         transformer2.FourierTransform(mv, V2Fourier, false);
         if (radavg) {
             auto V1size_x = (int)XSIZE(V1());
             auto V1size_y = (int)YSIZE(V1());
             auto V1size_z = (int)ZSIZE(V1());
             POCSFourierAmplitudeRadAvg(V2Fourier, lambda, radQuotient, V1size_x, V1size_y, V1size_z);
         }
         else {
             POCSFourierAmplitude(V1FourierMag, V2Fourier, lambda);
         }
         transformer2.inverseFourierTransform();
         if (computeE) {
             computeEnergy(Vdiff(), mv);
             Vdiff = V;
         }
         POCSMinMax(mv, v1min, v1max);
         if (computeE) {
             computeEnergy(Vdiff(), mv);
             Vdiff = V;
         }
         POCSmask(mask, mv);
         if (computeE) {
             computeEnergy(Vdiff(), mv);
             Vdiff = V;
         }
         transformer2.FourierTransform();
         POCSFourierPhase(V2FourierPhase, V2Fourier);
         transformer2.inverseFourierTransform();
         if (computeE) {
             computeEnergy(Vdiff(), mv);
             Vdiff = V;
         }
         POCSnonnegative(mv);
         if (computeE) {
             computeEnergy(Vdiff(), mv);
             Vdiff = V;
         }
         double std2 = mv.computeStddev();
         mv *= std1 / std2;
         if (computeE) {
             computeEnergy(Vdiff(), mv);
             Vdiff = V;
         }
         if (cutFreq != 0) {
             filter2.generateMask(mv);
             filter2.do_generate_3dmask = true;
             filter2.applyMaskSpace(mv);
             if (computeE) {
                 computeEnergy(Vdiff(), mv);
                 Vdiff = V;
             }
         }
     }
 
     if (performSubtraction) {
         auto masksub = getSubtractionMask(fnMaskSub, mask);
         V1.read(fnVolRef);
         V = subtraction(V1, V, masksub, fnVol1F, fnVol2A, filter2, cutFreq);
     }
 
     // Recover original alignment
     Image<double> Vf;
     MultidimArray<double> &mvf = Vf();
     mvf.setXmippOrigin();
     Vf = V;
     Matrix2D<double> A;
     A.initIdentity(4);
     geo2TransformationMatrix(rowIn,A);
     applyGeometry(xmipp_transformation::LINEAR, mvf, mv, A, xmipp_transformation::IS_INV, xmipp_transformation::DONT_WRAP);
     writeParticle(rowOut, fnImgOut, Vf); 
 }

◆ radialAverage()

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

Definition at line 227 of file subtomo_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;
 }

◆ readParams()

void ProgSubtomoSubtraction::readParams ( )

overridevirtual

Read argument from command line.

Reimplemented from XmippMetadataProgram.

Definition at line 79 of file subtomo_subtraction.cpp.

 {
     XmippMetadataProgram::readParams();
     fnVolRef = getParam("--ref");
     performSubtraction = checkParam("--sub");
     iter = getIntParam("--iter");
     sigma = getIntParam("--sigma");
     fnMask1 = getParam("--mask1");
     fnMask2 = getParam("--mask2");
     fnMaskSub = getParam("--maskSub");
     cutFreq = getDoubleParam("--cutFreq");
     lambda = getDoubleParam("--lambda");
     fnVol1F = getParam("--saveV1");
     if (fnVol1F.isEmpty())
         fnVol1F = "volume1_filtered.mrc";
     fnVol2A = getParam("--saveV2");
     if (fnVol2A.isEmpty())
         fnVol2A = "volume2_adjusted.mrc";
     radavg = checkParam("--radavg");
     computeE = checkParam("--computeEnergy");
 }

◆ readParticle()

void ProgSubtomoSubtraction::readParticle ( const MDRow & rowIn )

Read and write methods.

Definition at line 119 of file subtomo_subtraction.cpp.

                                                         {
     r.getValueOrDefault(MDL_IMAGE, fnVol2, "no_filename");
     V.read(fnVol2);
     V().setXmippOrigin();
  }

◆ show()

void ProgSubtomoSubtraction::show ( ) const

overridevirtual

Show.

Reimplemented from XmippProgram.

Definition at line 102 of file subtomo_subtraction.cpp.

                                         {
     if (!verbose)
         return;
     std::cout 
     << "Input volume(s):\t" << fnVolMd << std::endl
     << "Reference volume:\t" << fnVolRef << std::endl
     << "Input mask 1:\t" << fnMask1 << std::endl
     << "Input mask 2:\t" << fnMask2 << std::endl
     << "Input mask subtract:\t" << fnMaskSub << std::endl
     << "Sigma:\t" << sigma << std::endl
     << "Iterations:\t" << iter << std::endl
     << "Cutoff frequency:\t" << cutFreq << std::endl
     << "Relaxation factor:\t" << lambda << std::endl
     << "Match radial averages:\t" << radavg << std::endl
     << "Output:\t" << fnOut << std::endl;
 }

◆ subtraction()

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

Definition at line 288 of file subtomo_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;
 }

◆ writeParticle()

void ProgSubtomoSubtraction::writeParticle	(	MDRow &	rowOut,
		FileName	fnImgOut,
		Image< double > &	img
	)

Definition at line 125 of file subtomo_subtraction.cpp.

                                                                                                 {
     img.write(fnImgOut);
     rowOut.setValue(MDL_IMAGE, fnImgOut);
  }

Member Data Documentation

◆ computeE

bool ProgSubtomoSubtraction::computeE

Definition at line 56 of file subtomo_subtraction.h.

◆ cutFreq

double ProgSubtomoSubtraction::cutFreq

Definition at line 62 of file subtomo_subtraction.h.

◆ filter2

FourierFilter ProgSubtomoSubtraction::filter2

Definition at line 91 of file subtomo_subtraction.h.

◆ fnMask1

FileName ProgSubtomoSubtraction::fnMask1

Definition at line 50 of file subtomo_subtraction.h.

◆ fnMask2

FileName ProgSubtomoSubtraction::fnMask2

Definition at line 51 of file subtomo_subtraction.h.

◆ fnMaskSub

FileName ProgSubtomoSubtraction::fnMaskSub

Definition at line 52 of file subtomo_subtraction.h.

◆ fnOut

FileName ProgSubtomoSubtraction::fnOut

Definition at line 49 of file subtomo_subtraction.h.

◆ fnVol1F

FileName ProgSubtomoSubtraction::fnVol1F

Definition at line 53 of file subtomo_subtraction.h.

◆ fnVol2

FileName ProgSubtomoSubtraction::fnVol2

Definition at line 48 of file subtomo_subtraction.h.

◆ fnVol2A

FileName ProgSubtomoSubtraction::fnVol2A

Definition at line 54 of file subtomo_subtraction.h.

◆ fnVolMd

FileName ProgSubtomoSubtraction::fnVolMd

Definition at line 46 of file subtomo_subtraction.h.

◆ fnVolRef

FileName ProgSubtomoSubtraction::fnVolRef

Definition at line 47 of file subtomo_subtraction.h.

◆ iter

size_t ProgSubtomoSubtraction::iter

Definition at line 57 of file subtomo_subtraction.h.

◆ lambda

double ProgSubtomoSubtraction::lambda

Definition at line 63 of file subtomo_subtraction.h.

◆ mask

MultidimArray<double> ProgSubtomoSubtraction::mask

Definition at line 89 of file subtomo_subtraction.h.

◆ mdParticles

MetaDataVec ProgSubtomoSubtraction::mdParticles

Definition at line 68 of file subtomo_subtraction.h.

◆ n

size_t ProgSubtomoSubtraction::n

Definition at line 85 of file subtomo_subtraction.h.

◆ pad

double ProgSubtomoSubtraction::pad

Definition at line 83 of file subtomo_subtraction.h.

◆ padv

Image<double> ProgSubtomoSubtraction::padv

Definition at line 82 of file subtomo_subtraction.h.

◆ part_angles

struct Angles ProgSubtomoSubtraction::part_angles

Definition at line 77 of file subtomo_subtraction.h.

◆ performSubtraction

bool ProgSubtomoSubtraction::performSubtraction

Definition at line 60 of file subtomo_subtraction.h.

◆ radavg

bool ProgSubtomoSubtraction::radavg

Definition at line 64 of file subtomo_subtraction.h.

◆ rank

int ProgSubtomoSubtraction::rank

Definition at line 92 of file subtomo_subtraction.h.

◆ roffset

Matrix1D<double> ProgSubtomoSubtraction::roffset

Definition at line 70 of file subtomo_subtraction.h.

◆ row

MDRowVec ProgSubtomoSubtraction::row

Definition at line 69 of file subtomo_subtraction.h.

◆ sigma

int ProgSubtomoSubtraction::sigma

Definition at line 61 of file subtomo_subtraction.h.

◆ std1

double ProgSubtomoSubtraction::std1

Definition at line 84 of file subtomo_subtraction.h.

◆ subtomos

bool ProgSubtomoSubtraction::subtomos

Definition at line 65 of file subtomo_subtraction.h.

◆ transformer2

FourierTransformer ProgSubtomoSubtraction::transformer2

Definition at line 86 of file subtomo_subtraction.h.

◆ V

Image<double> ProgSubtomoSubtraction::V

Definition at line 80 of file subtomo_subtraction.h.

◆ V1

Image<double> ProgSubtomoSubtraction::V1

Definition at line 81 of file subtomo_subtraction.h.

◆ V1FourierMag

MultidimArray<double> ProgSubtomoSubtraction::V1FourierMag

Definition at line 88 of file subtomo_subtraction.h.

◆ v1max

double ProgSubtomoSubtraction::v1max

Definition at line 59 of file subtomo_subtraction.h.

◆ v1min

double ProgSubtomoSubtraction::v1min

Definition at line 58 of file subtomo_subtraction.h.

◆ V2Fourier

MultidimArray<std::complex<double> > ProgSubtomoSubtraction::V2Fourier

Definition at line 87 of file subtomo_subtraction.h.

◆ Vdiff

Image<double> ProgSubtomoSubtraction::Vdiff

Definition at line 90 of file subtomo_subtraction.h.

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

xmipp/libraries/reconstruction/subtomo_subtraction.h
xmipp/libraries/reconstruction/subtomo_subtraction.cpp

Classes

Public Member Functions

Public Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ ProgSubtomoSubtraction()

Member Function Documentation

◆ centerFFTMagnitude()

◆ computeEnergy()

◆ computeMagnitude()

◆ computePhase()

◆ computeRadQuotient()

◆ createFilter()

◆ createMask()

◆ defineParams()

◆ extractPhase()

◆ filterMask()

◆ getSubtractionMask()

◆ POCSFourierAmplitude()

◆ POCSFourierAmplitudeRadAvg()

◆ POCSFourierPhase()

◆ POCSmask()

◆ POCSMinMax()

◆ POCSnonnegative()

◆ postProcess()

◆ preProcess()

◆ processImage()

◆ radialAverage()

◆ readParams()

◆ readParticle()

◆ show()

◆ subtraction()

◆ writeParticle()

Member Data Documentation

◆ computeE

◆ cutFreq

◆ filter2

◆ fnMask1

◆ fnMask2

◆ fnMaskSub

◆ fnOut

◆ fnVol1F

◆ fnVol2

◆ fnVol2A

◆ fnVolMd

◆ fnVolRef

◆ iter

◆ lambda

◆ mask

◆ mdParticles

◆ n

◆ pad

◆ padv

◆ part_angles

◆ performSubtraction

◆ radavg

◆ rank

◆ roffset

◆ row

◆ sigma

◆ std1

◆ subtomos

◆ transformer2

◆ V

◆ V1

◆ V1FourierMag

◆ v1max

◆ v1min

◆ V2Fourier

◆ Vdiff