#include <subtract_projection.h>

Inheritance diagram for ProgSubtractProjection:

Collaboration diagram for ProgSubtractProjection:

Classes
struct	Angles

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

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

void	createMask (const FileName &, Image< double > &, Image< double > &)
	Processing methods. More...

Image< double >	binarizeMask (Projection &) const

Image< double >	invertMask (const Image< double > &)

Image< double >	applyCTF (const MDRow &, Projection &)

void	processParticle (const MDRow &rowIn, int, FourierTransformer &, FourierTransformer &)

MultidimArray< std::complex< double > >	computeEstimationImage (const MultidimArray< double > &, const MultidimArray< double > &, FourierTransformer &)

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

Matrix1D< double >	checkBestModel (MultidimArray< std::complex< double > > &, const MultidimArray< std::complex< double > > &, const MultidimArray< std::complex< double > > &, const MultidimArray< std::complex< double > > &) const

	ProgSubtractProjection ()
	Empty constructor. More...

	~ProgSubtractProjection ()
	Destructor. More...

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

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

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

void	preProcess () override

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	fnVolR

FileName	fnParticles

FileName	fnImgI

FileName	fnOut

FileName	fnMaskVol

FileName	fnMask

FileName	fnProj

double	sampling

double	padFourier

double	maxResol

double	cirmaskrad

int	sigma

int	limitfreq

int	maxwiIdx

bool	nonNegative

bool	boost

bool	subtract

MultidimArray< int >	wi

Image< double >	V

Image< double >	vM

Image< double >	ivM

Image< double >	M

Image< double >	I

Image< double >	Pctf

Image< double >	iM

Image< double >	Mfinal

Image< double >	Idiff

Image< double >	cirmask

Projection	P

Projection	Pmask

Projection	PmaskVol

FourierFilter	FilterG

const MultidimArray< double > *	ctfImage = nullptr

FourierTransformer	transformerP

FourierTransformer	transformerI

MultidimArray< std::complex< double > >	IFourier

MultidimArray< std::complex< double > >	PFourier

MultidimArray< std::complex< double > >	PFourier0

MultidimArray< std::complex< double > >	PFourier1

MultidimArray< std::complex< double > >	IiMFourier

MultidimArray< std::complex< double > >	PiMFourier

FourierTransformer	transformerIiM

FourierTransformer	transformerPiM

CTFDescription	ctf

FourierFilter	FilterCTF

Image< double >	padp

Image< double >	PmaskI

Image< double >	ImgiM

MultidimArray< std::complex< double > >	ImgiMFourier

MetaDataVec	mdParticles

MDRowVec	row

Matrix1D< double >	roffset

struct Angles	part_angles

bool	disable

int	rank

FourierProjector *	projector

FourierProjector *	projectorMask

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

Subtract projections from particles

Definition at line 41 of file subtract_projection.h.

Constructor & Destructor Documentation

◆ ProgSubtractProjection()

ProgSubtractProjection::ProgSubtractProjection ( )

Empty constructor.

Definition at line 51 of file subtract_projection.cpp.

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

◆ ~ProgSubtractProjection()

ProgSubtractProjection::~ProgSubtractProjection ( )

Destructor.

Definition at line 63 of file subtract_projection.cpp.

 {
     delete projector;
     delete projectorMask;
 }

Member Function Documentation

◆ applyCTF()

Image< double > ProgSubtractProjection::applyCTF	(	const MDRow &	r,
		Projection &	proj
	)

Definition at line 188 of file subtract_projection.cpp.

                                                                                 {
     if (r.containsLabel(MDL_CTF_DEFOCUSU) || r.containsLabel(MDL_CTF_MODEL)){
         ctf.readFromMdRow(r);
         ctf.Tm = sampling;
         ctf.produceSideInfo();
         FilterCTF.FilterBand = CTF;
         FilterCTF.ctf.enable_CTFnoise = false;
         FilterCTF.ctf = ctf;
         // Padding before apply CTF
         MultidimArray <double> &mpad = padp();
         mpad.setXmippOrigin();
         MultidimArray<double> &mproj = proj();
         mproj.setXmippOrigin();
         mproj.window(mpad,STARTINGY(mproj)*(int)padFourier, STARTINGX(mproj)*(int)padFourier, FINISHINGY(mproj)*(int)padFourier, FINISHINGX(mproj)*(int)padFourier);
         FilterCTF.generateMask(mpad);
         FilterCTF.applyMaskSpace(mpad); 
         //Crop to restore original size
         mpad.window(mproj,STARTINGY(mproj), STARTINGX(mproj), FINISHINGY(mproj), FINISHINGX(mproj));
     }
     return proj;
  }

◆ binarizeMask()

Image< double > ProgSubtractProjection::binarizeMask ( Projection & m ) const

Definition at line 170 of file subtract_projection.cpp.

                                                                        {
     double maxMaskVol;
     double minMaskVol;
     MultidimArray<double> &mm=m();
     mm.computeDoubleMinMax(minMaskVol, maxMaskVol);
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(mm)
         DIRECT_MULTIDIM_ELEM(mm,n) = (DIRECT_MULTIDIM_ELEM(mm,n)>0.1*maxMaskVol) ? 1:0; 
     return m;
  }

◆ checkBestModel()

Matrix1D< double > ProgSubtractProjection::checkBestModel	(	MultidimArray< std::complex< double > > &	PFourierf,
		const MultidimArray< std::complex< double > > &	PFourierf0,
		const MultidimArray< std::complex< double > > &	PFourierf1,
		const MultidimArray< std::complex< double > > &	IFourierf
	)		const

Definition at line 259 of file subtract_projection.cpp.

                                                                                                                     { 
     // Compute R2 coefficient for order 0 model (R20) and order 1 model (R21)
     auto N = 2.0*(double)MULTIDIM_SIZE(PFourierf);
     double R20adj = evaluateFitting(IFourierf, PFourierf0); // adjusted R2 for an order 0 model = R2
     double R21 = evaluateFitting(IFourierf, PFourierf1); 
     double R21adj = 1.0 - (1.0 - R21) * (N - 1.0) / (N - 2.0); // adjusted R2 for an order 1 model -> p = 2
     //Decide best fitting
     Matrix1D<double> R2(2);
     if (R21adj > R20adj) { // Order 1: T(w) = b01 + b1*wi 
         PFourierf = PFourierf1;
         R2(0) = R21adj;
         R2(1) = 1;
     } 
     else { // Order 0: T(w) = b00 
         PFourierf = PFourierf0;
         R2(0) = R20adj;
         R2(1) = 0;      
     }
     return R2;
 }

◆ computeEstimationImage()

MultidimArray< std::complex< double > > ProgSubtractProjection::computeEstimationImage	(	const MultidimArray< double > &	Img,
		const MultidimArray< double > &	InvM,
		FourierTransformer &	transformerImgiM
	)

Definition at line 229 of file subtract_projection.cpp.

                                                                          {
     ImgiM().initZeros(Img);
     ImgiM().setXmippOrigin();
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(Img)
         DIRECT_MULTIDIM_ELEM(ImgiM(),n) = DIRECT_MULTIDIM_ELEM(Img,n) * DIRECT_MULTIDIM_ELEM(InvM,n);
     transformerImgiM.FourierTransform(ImgiM(),ImgiMFourier,false);
     return ImgiMFourier;
 }

◆ createMask()

void ProgSubtractProjection::createMask	(	const FileName &	fnM,
		Image< double > &	m,
		Image< double > &	im
	)

Processing methods.

Definition at line 149 of file subtract_projection.cpp.

                                                                                                  {
     if (fnM.isEmpty()) 
     {
         m().initZeros((int)XSIZE(V()),(int)YSIZE(V()),(int)ZSIZE(V()));
         im = m;
         FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(im())
             DIRECT_MULTIDIM_ELEM(im(),n) += 1; 
     }
     else 
     {
         m.read(fnM);
         m().setXmippOrigin();
         im = m;
         if (!subtract)
         {
             FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(im())
                 DIRECT_MULTIDIM_ELEM(im(),n) = (DIRECT_MULTIDIM_ELEM(m(),n)*(-1))+1;
         } 
     }
  }

◆ defineParams()

void ProgSubtractProjection::defineParams ( )

overridevirtual

Define parameters.

Reimplemented from XmippMetadataProgram.

Definition at line 104 of file subtract_projection.cpp.

  {
      //Usage
      addUsageLine("This program computes the subtraction between particles and a reference"); 
      addUsageLine(" volume, by computing its projections with the same angles that input particles have."); 
      addUsageLine(" Then, each particle and the correspondent projection of the reference volume are numerically");
      addUsageLine(" adjusted and subtracted using a mask which denotes the region to keep or subtract.");
      //Parameters
      XmippMetadataProgram::defineParams();
      addParamsLine("--ref <volume>\t: Reference volume to subtract");
      addParamsLine("[--mask <mask=\"\">]\t: 3D mask for region to keep, no mask implies subtraction of whole images");
      addParamsLine("[--sampling <sampling=1>]\t: Sampling rate (A/pixel)");
      addParamsLine("[--max_resolution <f=4>]\t: Maximum resolution (A)");
      addParamsLine("[--fmask_width <w=40>]\t: Extra width of final mask (A). -1 means no masking."); 
      addParamsLine("[--padding <p=2>]\t: Padding factor for Fourier projector");
      addParamsLine("[--sigma <s=2>]\t: Decay of the filter (sigma) to smooth the mask transition");
      addParamsLine("[--limit_freq <l=0>]\t: Limit frequency (= 1) or not (= 0) in adjustment process");
      addParamsLine("[--nonNegative]\t: Ignore particles with negative beta0 or R2"); 
      addParamsLine("[--boost]\t: Perform a boosting of original particles"); 
      addParamsLine("[--cirmaskrad <c=-1.0>]\t: Radius of the circular mask");
      addParamsLine("[--save <structure=\"\">]\t: Path for saving intermediate files"); 
      addParamsLine("[--subtract]\t: The mask contains the region to SUBTRACT"); 
      addExampleLine("A typical use is:",false);
      addExampleLine("xmipp_subtract_projection -i input_particles.xmd --ref input_map.mrc --mask mask_vol.mrc "
              "-o output_particles --sampling 1 --fmask_width 40 --max_resolution 4");
  }

◆ evaluateFitting()

double ProgSubtractProjection::evaluateFitting	(	const MultidimArray< std::complex< double > > &	y,
		const MultidimArray< std::complex< double > > &	yp
	)		const

Definition at line 239 of file subtract_projection.cpp.

                                                                                                     {
     double sumY = 0;
     double sumY2 = 0;
     double sumE2 = 0;
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(y) {
         double realyn = real(DIRECT_MULTIDIM_ELEM(y, n)); 
         double imagyn = imag(DIRECT_MULTIDIM_ELEM(y, n)); 
         double ereal = realyn - real(DIRECT_MULTIDIM_ELEM(yp, n)); 
         double eimag = imagyn - imag(DIRECT_MULTIDIM_ELEM(yp, n));  
         sumE2 += ereal*ereal + eimag*eimag; 
         sumY += realyn + imagyn; 
         sumY2 += realyn*realyn + imagyn*imagyn; 
     }
     auto meanY = sumY / (2.0*(double)MULTIDIM_SIZE(y));
     auto varY = sumY2 / (2.0*(double)MULTIDIM_SIZE(y)) - meanY*meanY; 
     auto R2 = 1.0 - (sumE2/(2.0*(double)MULTIDIM_SIZE(y))) / varY;
     return R2;
  }

◆ invertMask()

Image< double > ProgSubtractProjection::invertMask ( const Image< double > & m )

Definition at line 180 of file subtract_projection.cpp.

                                                                         {
     PmaskI = m;
     MultidimArray<double> &mPmaskI=PmaskI();
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(mPmaskI)
         DIRECT_MULTIDIM_ELEM(mPmaskI,n) = (DIRECT_MULTIDIM_ELEM(mPmaskI,n)*(-1))+1;
     return PmaskI;
  }

◆ postProcess()

void ProgSubtractProjection::postProcess ( )

overridevirtual

Reimplemented from XmippMetadataProgram.

Definition at line 474 of file subtract_projection.cpp.

 {
     getOutputMd().write(fn_out);
 }

◆ preProcess()

void ProgSubtractProjection::preProcess ( )

overridevirtual

Reimplemented from XmippMetadataProgram.

Definition at line 281 of file subtract_projection.cpp.

                                          {
     // Read input volume, mask and particles metadata
     show();
     V.read(fnVolR);
     V().setXmippOrigin();
 
     // Create 2D circular mask to avoid edge artifacts after wrapping
     size_t Xdim;
     size_t Ydim;
     size_t Zdim;
     size_t Ndim;
     V.getDimensions(Xdim, Ydim, Zdim, Ndim);
     cirmask().initZeros((int)Ydim, (int)Xdim);
     cirmask().setXmippOrigin();
     if (cirmaskrad == -1.0)
         cirmaskrad = (double)XSIZE(V())/2;
     RaisedCosineMask(cirmask(), cirmaskrad*0.8, cirmaskrad*0.9);
     cirmask.write(formatString("%s/cirmask.mrc", fnProj.c_str()));
     
     // Create mock image of same size as particles (and referencce volume) to get
     I().initZeros((int)Ydim, (int)Xdim);
     I().initConstant(1);
     transformerI.FourierTransform(I(), IFourier, false);
 
     // Initialize Gaussian LPF to smooth mask
     FilterG.FilterShape=REALGAUSSIAN;
     FilterG.FilterBand=LOWPASS;
     FilterG.w1=sigma;
 
     // Construct frequencies image
     wi.initZeros(IFourier);
     Matrix1D<double> w(2);  
     double cutFreq = sampling/maxResol;
     for (int i=0; i<YSIZE(wi); i++) {
         FFT_IDX2DIGFREQ(i,YSIZE(IFourier),YY(w)) 
         for (int j=0; j<XSIZE(wi); j++)  {
             FFT_IDX2DIGFREQ(j,XSIZE(IFourier),XX(w))
             DIRECT_A2D_ELEM(wi,i,j) = (int)round((sqrt(YY(w)*YY(w) + XX(w)*XX(w))) * (int)XSIZE(IFourier)); // indexes
         }
     }
     if (limitfreq == 0)
         maxwiIdx = (int)XSIZE(wi); 
     else
         DIGFREQ2FFT_IDX(cutFreq, (int)YSIZE(IFourier), maxwiIdx)
 
     if (rank==0)
     {
         // Read or create mask keep and compute inverse of mask keep (mask subtract)
         createMask(fnMask, vM, ivM);
         FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(V())
             DIRECT_MULTIDIM_ELEM(V(),n) = DIRECT_MULTIDIM_ELEM(V(),n)*DIRECT_MULTIDIM_ELEM(ivM(),n); 
         // Initialize Fourier projectors
         std::cout << "-------Initializing projectors-------" << std::endl;
         projector = new FourierProjector(V(), padFourier, cutFreq, xmipp_transformation::BSPLINE3);
         projectorMask = new FourierProjector(vM(), padFourier, cutFreq, xmipp_transformation::BSPLINE3);
         std::cout << "-------Projectors initialized-------" << std::endl;
     }
     else
     {
         projector = new FourierProjector(padFourier,cutFreq,xmipp_transformation::BSPLINE3);
         projectorMask = new FourierProjector(padFourier,cutFreq,xmipp_transformation::BSPLINE3);
     }
  }

◆ processImage()

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

overridevirtual

Implements XmippMetadataProgram.

Definition at line 345 of file subtract_projection.cpp.

  { 
     // Initialize aux variable
     disable = false;
     // Project volume and process projections 
     const auto sizeI = (int)XSIZE(I());
     processParticle(rowIn, sizeI, transformerP, transformerI);
     // Build projected and final masks
     if (fnMask.isEmpty()) { // If there is no provided mask
         M().initZeros(P());
         // inverse mask (iM) is all 1s
         iM = invertMask(M);
     }
     else { // If a mask has been provided
         projectVolume(*projectorMask, Pmask, sizeI, sizeI, part_angles.rot, part_angles.tilt, part_angles.psi, ctfImage);   
         // Apply binarization, shift and gaussian filter to the projected mask
         M = binarizeMask(Pmask);
         selfTranslate(xmipp_transformation::LINEAR, M(), roffset, xmipp_transformation::DONT_WRAP);
         FilterG.applyMaskSpace(M());
         if (subtract) // If the mask contains the part to SUBTRACT: iM = input mask
             iM = M;
         else // If the mask contains the part to KEEP: iM = INVERSE of original mask
             iM = invertMask(M);
     }
     
     // Compute estimation images: IiM = I*iM and PiM = P*iM 
     IiMFourier = computeEstimationImage(I(), iM(), transformerIiM);
     PiMFourier = computeEstimationImage(Pctf(), iM(), transformerPiM);  
 
     // Estimate transformation with model of order 0: T(w) = beta00 and model of order 1: T(w) = beta01 + beta1*w
     MultidimArray<double> num0;
     num0.initZeros(maxwiIdx+1); 
     MultidimArray<double> den0;
     den0.initZeros(maxwiIdx+1);
     Matrix2D<double> A1;
     A1.initZeros(2,2);
     Matrix1D<double> b1;
     b1.initZeros(2);
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(PiMFourier) {
         int win = DIRECT_MULTIDIM_ELEM(wi, n);
         if (win < maxwiIdx) 
         {
             double realPiMFourier = real(DIRECT_MULTIDIM_ELEM(PiMFourier,n));
             double imagPiMFourier = imag(DIRECT_MULTIDIM_ELEM(PiMFourier,n));
             DIRECT_MULTIDIM_ELEM(num0,win) += real(DIRECT_MULTIDIM_ELEM(IiMFourier,n)) * realPiMFourier
                                             + imag(DIRECT_MULTIDIM_ELEM(IiMFourier,n)) * imagPiMFourier;
             DIRECT_MULTIDIM_ELEM(den0,win) += realPiMFourier*realPiMFourier + imagPiMFourier*imagPiMFourier;
             A1(0,0) += realPiMFourier*realPiMFourier + imagPiMFourier*imagPiMFourier;
             A1(0,1) += win*(realPiMFourier + imagPiMFourier);
             A1(1,1) += 2*win;
             b1(0) += real(DIRECT_MULTIDIM_ELEM(IiMFourier,n)) * realPiMFourier + imag(DIRECT_MULTIDIM_ELEM(IiMFourier,n)) * imagPiMFourier;
             b1(1) += win*(real(DIRECT_MULTIDIM_ELEM(IiMFourier,n))+imag(DIRECT_MULTIDIM_ELEM(IiMFourier,n)));
         }
     }
     A1(1,0) = A1(0,1);
 
     // Compute beta00 from order 0 model
     double beta00 = num0.sum()/den0.sum();
     if (nonNegative && beta00 < 0) 
     {
         disable = true;
     }
     // Apply adjustment order 0: PFourier0 = T(w) * PFourier = beta00 * PFourier
     PFourier0 = PFourier;
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(PFourier0) 
         DIRECT_MULTIDIM_ELEM(PFourier0,n) *= beta00; 
     PFourier0(0,0) = IiMFourier(0,0); 
 
     // Compute beta01 and beta1 from order 1 model
     PseudoInverseHelper h;
     h.A = A1;
     h.b = b1;
     Matrix1D<double> betas1;
     solveLinearSystem(h,betas1); 
     double beta01 = betas1(0);
     double beta1 = betas1(1);
 
     // Apply adjustment order 1: PFourier1 = T(w) * PFourier = (beta01 + beta1*w) * PFourier
     PFourier1 = PFourier;
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(PFourier1)
         DIRECT_MULTIDIM_ELEM(PFourier1,n) *= (beta01+beta1*DIRECT_MULTIDIM_ELEM(wi,n)); 
     PFourier1(0,0) = IiMFourier(0,0); 
 
     // Check best model
     Matrix1D<double> R2adj = checkBestModel(PFourier, PFourier0, PFourier1, IFourier);
     double beta0save;
     double beta1save;       
     if (R2adj(1) == 0)
     {
         beta0save = beta00;
         beta1save = 0;
     }
     else
     {
         beta0save = beta01;
         beta1save = beta1;
     }
 
     // Create empty new image for output particle
     MultidimArray<double> &mIdiff=Idiff();
     mIdiff.initZeros(I());
     mIdiff.setXmippOrigin();
 
     if (boost) // Boosting of original particles
     {
         if (R2adj(1) == 0)
         {
             FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(IFourier) 
                 DIRECT_MULTIDIM_ELEM(IFourier,n) /= beta00; 
         } 
         else if (R2adj(1) == 1)
         {
             FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(IFourier)
                 DIRECT_MULTIDIM_ELEM(IFourier,n) /= (beta01+beta1*DIRECT_MULTIDIM_ELEM(wi,n)); 
         }
         transformerI.inverseFourierTransform(IFourier, Idiff());
     } 
     else  // Subtraction
     {
         // Recover adjusted projection (P) in real space
         transformerP.inverseFourierTransform(PFourier, P());
         mIdiff.initZeros(I());
         mIdiff.setXmippOrigin();
         FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(mIdiff)
             DIRECT_MULTIDIM_ELEM(mIdiff,n) = DIRECT_MULTIDIM_ELEM(I(),n)-DIRECT_MULTIDIM_ELEM(P(),n);
     }
     writeParticle(rowOut, fnImgOut, Idiff, R2adj(0), beta0save, beta1save); 
 }

◆ processParticle()

void ProgSubtractProjection::processParticle	(	const MDRow &	rowIn,
		int	sizeImg,
		FourierTransformer &	transformerPf,
		FourierTransformer &	transformerIf
	)

Definition at line 210 of file subtract_projection.cpp.

                                                                                                                                                        {
     readParticle(rowprocess);
     rowprocess.getValueOrDefault(MDL_ANGLE_ROT, part_angles.rot, 0);
     rowprocess.getValueOrDefault(MDL_ANGLE_TILT, part_angles.tilt, 0);
     rowprocess.getValueOrDefault(MDL_ANGLE_PSI, part_angles.psi, 0);
     roffset.initZeros(2);
     rowprocess.getValueOrDefault(MDL_SHIFT_X, roffset(0), 0);
     rowprocess.getValueOrDefault(MDL_SHIFT_Y, roffset(1), 0);
     roffset *= -1;
     projectVolume(*projector, P, sizeImg, sizeImg, part_angles.rot, part_angles.tilt, part_angles.psi, ctfImage);
     selfTranslate(xmipp_transformation::LINEAR, P(), roffset, xmipp_transformation::WRAP);
     Pctf = applyCTF(rowprocess, P);
     MultidimArray<double> &mPctf = Pctf();
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(mPctf)
         DIRECT_MULTIDIM_ELEM(mPctf,n) = DIRECT_MULTIDIM_ELEM(mPctf,n) * DIRECT_MULTIDIM_ELEM(cirmask(),n);
     transformerPf.FourierTransform(Pctf(), PFourier, false);
     transformerIf.FourierTransform(I(), IFourier, false);
 }

◆ readParams()

void ProgSubtractProjection::readParams ( )

overridevirtual

Read argument from command line.

Reimplemented from XmippMetadataProgram.

Definition at line 70 of file subtract_projection.cpp.

  {
     XmippMetadataProgram::readParams();
     fnVolR = getParam("--ref");
     fnMask=getParam("--mask");
     sigma=getIntParam("--sigma");
     sampling = getDoubleParam("--sampling");
     padFourier = getDoubleParam("--padding");
     maxResol = getDoubleParam("--max_resolution");
     limitfreq = getIntParam("--limit_freq");
     cirmaskrad = getDoubleParam("--cirmaskrad");
     fnProj = getParam("--save"); 
     nonNegative = checkParam("--nonNegative");
     boost = checkParam("--boost");
     subtract = checkParam("--subtract");
  }

◆ readParticle()

void ProgSubtractProjection::readParticle ( const MDRow & rowIn )

Read and write methods.

Definition at line 131 of file subtract_projection.cpp.

                                                          {
     r.getValueOrDefault(MDL_IMAGE, fnImgI, "no_filename");
     I.read(fnImgI);
     I().setXmippOrigin();
  }

◆ show()

void ProgSubtractProjection::show ( ) const

overridevirtual

Show.

Reimplemented from XmippProgram.

Definition at line 88 of file subtract_projection.cpp.

                                         {
     if (!verbose)
         return;
     std::cout
     << "Input particles:\t" << fnParticles << std::endl
     << "Reference volume:\t" << fnVolR << std::endl
     << "Mask of the region to keep:\t" << fnMask << std::endl
     << "Sigma of low pass filter:\t" << sigma << std::endl
     << "Sampling rate:\t" << sampling << std::endl
     << "Padding factor:\t" << padFourier << std::endl
     << "Max. Resolution:\t" << maxResol << std::endl
     << "Limit frequency:\t" << limitfreq << std::endl
     << "Output particles:\t" << fnOut << std::endl;
  }

◆ writeParticle()

void ProgSubtractProjection::writeParticle	(	MDRow &	rowOut,
		FileName	fnImgOut,
		Image< double > &	img,
		double	R2a,
		double	b0save,
		double	b1save
	)

Definition at line 137 of file subtract_projection.cpp.

                                                                                                                                           {
     img.write(fnImgOut);
     rowOut.setValue(MDL_IMAGE, fnImgOut);
     rowOut.setValue(MDL_SUBTRACTION_R2, R2a); 
     rowOut.setValue(MDL_SUBTRACTION_BETA0, b0save); 
     rowOut.setValue(MDL_SUBTRACTION_BETA1, b1save); 
     if (nonNegative && (disable || R2a < 0)) 
     {
         rowOut.setValue(MDL_ENABLED, -1);
     }
  }

Member Data Documentation

◆ boost

bool ProgSubtractProjection::boost

Definition at line 60 of file subtract_projection.h.

◆ cirmask

Image<double> ProgSubtractProjection::cirmask

Definition at line 75 of file subtract_projection.h.

◆ cirmaskrad

double ProgSubtractProjection::cirmaskrad

Definition at line 55 of file subtract_projection.h.

◆ ctf

CTFDescription ProgSubtractProjection::ctf

Definition at line 97 of file subtract_projection.h.

◆ ctfImage

const MultidimArray<double>* ProgSubtractProjection::ctfImage = nullptr

Definition at line 82 of file subtract_projection.h.

◆ disable

bool ProgSubtractProjection::disable

Definition at line 116 of file subtract_projection.h.

◆ FilterCTF

FourierFilter ProgSubtractProjection::FilterCTF

Definition at line 98 of file subtract_projection.h.

◆ FilterG

FourierFilter ProgSubtractProjection::FilterG

Definition at line 80 of file subtract_projection.h.

◆ fnImgI

FileName ProgSubtractProjection::fnImgI

Definition at line 47 of file subtract_projection.h.

◆ fnMask

FileName ProgSubtractProjection::fnMask

Definition at line 50 of file subtract_projection.h.

◆ fnMaskVol

FileName ProgSubtractProjection::fnMaskVol

Definition at line 49 of file subtract_projection.h.

◆ fnOut

FileName ProgSubtractProjection::fnOut

Definition at line 48 of file subtract_projection.h.

◆ fnParticles

FileName ProgSubtractProjection::fnParticles

Definition at line 46 of file subtract_projection.h.

◆ fnProj

FileName ProgSubtractProjection::fnProj

Definition at line 51 of file subtract_projection.h.

◆ fnVolR

FileName ProgSubtractProjection::fnVolR

Definition at line 45 of file subtract_projection.h.

◆ I

Image<double> ProgSubtractProjection::I

Definition at line 70 of file subtract_projection.h.

◆ Idiff

Image<double> ProgSubtractProjection::Idiff

Definition at line 74 of file subtract_projection.h.

◆ IFourier

MultidimArray< std::complex<double> > ProgSubtractProjection::IFourier

Definition at line 86 of file subtract_projection.h.

◆ IiMFourier

MultidimArray< std::complex<double> > ProgSubtractProjection::IiMFourier

Definition at line 90 of file subtract_projection.h.

◆ iM

Image<double> ProgSubtractProjection::iM

Definition at line 72 of file subtract_projection.h.

◆ ImgiM

Image<double> ProgSubtractProjection::ImgiM

Definition at line 101 of file subtract_projection.h.

◆ ImgiMFourier

MultidimArray< std::complex<double> > ProgSubtractProjection::ImgiMFourier

Definition at line 102 of file subtract_projection.h.

◆ ivM

Image<double> ProgSubtractProjection::ivM

Definition at line 67 of file subtract_projection.h.

◆ limitfreq

int ProgSubtractProjection::limitfreq

Definition at line 57 of file subtract_projection.h.

◆ M

Image<double> ProgSubtractProjection::M

Definition at line 69 of file subtract_projection.h.

◆ maxResol

double ProgSubtractProjection::maxResol

Definition at line 54 of file subtract_projection.h.

◆ maxwiIdx

int ProgSubtractProjection::maxwiIdx

Definition at line 58 of file subtract_projection.h.

◆ mdParticles

MetaDataVec ProgSubtractProjection::mdParticles

Definition at line 105 of file subtract_projection.h.

◆ Mfinal

Image<double> ProgSubtractProjection::Mfinal

Definition at line 73 of file subtract_projection.h.

◆ nonNegative

bool ProgSubtractProjection::nonNegative

Definition at line 59 of file subtract_projection.h.

◆ P

Projection ProgSubtractProjection::P

Definition at line 77 of file subtract_projection.h.

◆ padFourier

double ProgSubtractProjection::padFourier

Definition at line 53 of file subtract_projection.h.

◆ padp

Image<double> ProgSubtractProjection::padp

Definition at line 99 of file subtract_projection.h.

◆ part_angles

struct Angles ProgSubtractProjection::part_angles

Definition at line 114 of file subtract_projection.h.

◆ Pctf

Image<double> ProgSubtractProjection::Pctf

Definition at line 71 of file subtract_projection.h.

◆ PFourier

MultidimArray< std::complex<double> > ProgSubtractProjection::PFourier

Definition at line 87 of file subtract_projection.h.

◆ PFourier0

MultidimArray< std::complex<double> > ProgSubtractProjection::PFourier0

Definition at line 88 of file subtract_projection.h.

◆ PFourier1

MultidimArray< std::complex<double> > ProgSubtractProjection::PFourier1

Definition at line 89 of file subtract_projection.h.

◆ PiMFourier

MultidimArray< std::complex<double> > ProgSubtractProjection::PiMFourier

Definition at line 91 of file subtract_projection.h.

◆ Pmask

Projection ProgSubtractProjection::Pmask

Definition at line 78 of file subtract_projection.h.

◆ PmaskI

Image<double> ProgSubtractProjection::PmaskI

Definition at line 100 of file subtract_projection.h.

◆ PmaskVol

Projection ProgSubtractProjection::PmaskVol

Definition at line 79 of file subtract_projection.h.

◆ projector

FourierProjector* ProgSubtractProjection::projector

Definition at line 133 of file subtract_projection.h.

◆ projectorMask

FourierProjector* ProgSubtractProjection::projectorMask

Definition at line 134 of file subtract_projection.h.

◆ rank

int ProgSubtractProjection::rank

Definition at line 132 of file subtract_projection.h.

◆ roffset

Matrix1D<double> ProgSubtractProjection::roffset

Definition at line 107 of file subtract_projection.h.

◆ row

MDRowVec ProgSubtractProjection::row

Definition at line 106 of file subtract_projection.h.

◆ sampling

double ProgSubtractProjection::sampling

Definition at line 52 of file subtract_projection.h.

◆ sigma

int ProgSubtractProjection::sigma

Definition at line 56 of file subtract_projection.h.

◆ subtract

bool ProgSubtractProjection::subtract

Definition at line 61 of file subtract_projection.h.

◆ transformerI

FourierTransformer ProgSubtractProjection::transformerI

Definition at line 84 of file subtract_projection.h.

◆ transformerIiM

FourierTransformer ProgSubtractProjection::transformerIiM

Definition at line 94 of file subtract_projection.h.

◆ transformerP

FourierTransformer ProgSubtractProjection::transformerP

Definition at line 83 of file subtract_projection.h.

◆ transformerPiM

FourierTransformer ProgSubtractProjection::transformerPiM

Definition at line 95 of file subtract_projection.h.

◆ V

Image<double> ProgSubtractProjection::V

Definition at line 65 of file subtract_projection.h.

◆ vM

Image<double> ProgSubtractProjection::vM

Definition at line 66 of file subtract_projection.h.

◆ wi

MultidimArray<int> ProgSubtractProjection::wi

Definition at line 62 of file subtract_projection.h.

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

xmipp/libraries/reconstruction/subtract_projection.h
xmipp/libraries/reconstruction/subtract_projection.cpp

Classes

Public Member Functions

Public Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ ProgSubtractProjection()

◆ ~ProgSubtractProjection()

Member Function Documentation

◆ applyCTF()

◆ binarizeMask()

◆ checkBestModel()

◆ computeEstimationImage()

◆ createMask()

◆ defineParams()

◆ evaluateFitting()

◆ invertMask()

◆ postProcess()

◆ preProcess()

◆ processImage()

◆ processParticle()

◆ readParams()

◆ readParticle()

◆ show()

◆ writeParticle()

Member Data Documentation

◆ boost

◆ cirmask

◆ cirmaskrad

◆ ctf

◆ ctfImage

◆ disable

◆ FilterCTF

◆ FilterG

◆ fnImgI

◆ fnMask

◆ fnMaskVol

◆ fnOut

◆ fnParticles

◆ fnProj

◆ fnVolR

◆ I

◆ Idiff

◆ IFourier

◆ IiMFourier

◆ iM

◆ ImgiM

◆ ImgiMFourier

◆ ivM

◆ limitfreq

◆ M

◆ maxResol

◆ maxwiIdx

◆ mdParticles

◆ Mfinal

◆ nonNegative

◆ P

◆ padFourier

◆ padp

◆ part_angles

◆ Pctf

◆ PFourier

◆ PFourier0

◆ PFourier1

◆ PiMFourier

◆ Pmask

◆ PmaskI

◆ PmaskVol

◆ projector

◆ projectorMask

◆ rank

◆ roffset

◆ row

◆ sampling

◆ sigma

◆ subtract

◆ transformerI

◆ transformerIiM

◆ transformerP

◆ transformerPiM