#include <tomo_extract_subtomograms.h>

Inheritance diagram for ProgTomoExtractSubtomograms:

Collaboration diagram for ProgTomoExtractSubtomograms:

Public Member Functions
void	defineParams ()

void	readParams ()

void	createSphere (int halfboxsize)

void	normalizeSubtomo (MultidimArray< double > &subtomo, int halfboxsize)

void	extractSubtomoFixedSize (MultidimArray< double > &subtomoExtraction)

void	writeSubtomo (int idx, int xcoor, int ycoor, int zcoor, size_t particleid)

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 ()

Static Public Member Functions
static void	upsample (const MultidimArray< std::complex< double >> &from, MultidimArray< std::complex< double >> &to)

static void	downsample (const MultidimArray< std::complex< double >> &from, MultidimArray< std::complex< double >> &to)

Public Attributes
FileName	fnOut

FileName	fnTom

FileName	fnCoor

size_t	Xdim

size_t	Ydim

size_t	Zdim

MetaDataVec	md

MetaDataVec	mdout

MDRowVec	rowout

bool	invertContrast

bool	normalize

bool	fixedBoxSize

double	scaleFactor

double	downsampleFactor

std::vector< size_t >	maskIdx

int	boxsize

int	nthrs

Image< double >	subtomoImg

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

Definition at line 44 of file tomo_extract_subtomograms.h.

Member Function Documentation

◆ createSphere()

void ProgTomoExtractSubtomograms::createSphere ( int halfboxsize )

Definition at line 62 of file tomo_extract_subtomograms.cpp.

 {
     long n=0;
 
     for (int k=0; k<boxsize; k++)
     {
         int k2 = (k-halfboxsize);
         k2 = k2*k2;
         for (int i=0; i<boxsize; i++)
         {
             int i2 = i-halfboxsize;
             int i2k2 = i2*i2 +k2 ;
             for (int j=0; j<boxsize; j++)
             {
                 int j2 = (j- halfboxsize);
                 if (sqrt(i2k2 + j2*j2)>halfboxsize)
                 {
                     maskIdx.push_back(n);
                 }
                 n++;
             }
         }
     }
 }

◆ defineParams()

void ProgTomoExtractSubtomograms::defineParams ( )

virtual

Function in which the param of each Program are defined.

Reimplemented from XmippProgram.

Definition at line 47 of file tomo_extract_subtomograms.cpp.

 {
     addUsageLine("This function takes a tomogram an extract a set of subtomogram from it. The coordinates of the subtomograms are speciffied in the metadata given by coordinates.");
     addParamsLine("  --tomogram <vol_file=\"\">               : Filename of the tomogram containing the subtomograms to be extracted");
     addParamsLine("  --coordinates <vol_file=\"\">            : Metadata (.xmd file) with the coordidanates to be extracted from the tomogram");
     addParamsLine("  --boxsize <boxsize=100>                  : Particle box size in voxels, of the particle without downsampling.");
     addParamsLine("  [--invertContrast]                       : Set this flag to invert the contrast of the extracted subtomograms");
     addParamsLine("  [--normalize]                            : This flag will set the subtomograms to have zero mean and unit standard deviation.");
     addParamsLine("  [--downsample <downsampleFactor=1.0>]    : Scale factor of the extracted subtomograms. It must be greater than 1. A downsampling 2 reduces in a factor 2 the size of the subtomos.");
     addParamsLine("  [--fixedBoxSize]                         : If selected, programs calculates the extraction box to obtain a box size of the selected size when downsampling.");
     addParamsLine("  -o <vol_file=\"\">                       : Path of the output directory. ");
     addParamsLine("  [--threads <s=4>]                        : Number of threads");
 }

◆ downsample()

void ProgTomoExtractSubtomograms::downsample	(	const MultidimArray< std::complex< double >> &	from,
		MultidimArray< std::complex< double >> &	to
	)

static

Definition at line 334 of file tomo_extract_subtomograms.cpp.

 {
     for (size_t k = 0; k < ZSIZE(to)/2; k++)
     {
         for (size_t j = 0; j < YSIZE(to)/2; j++)
         {
             for (size_t i = 0; i < XSIZE(to); i++)
             {
                 // Origin cube
                 DIRECT_A3D_ELEM(to, k, j, i) = DIRECT_A3D_ELEM(from, k, j, i);
 
                 // Oposite cube
                 DIRECT_A3D_ELEM(to,  ZSIZE(to)/2 + k, YSIZE(to)/2 + j, i) =
                 DIRECT_A3D_ELEM(from, ZSIZE(from) - ZSIZE(to)/2 + k, YSIZE(from) - YSIZE(to)/2 + j, i);
 
                 // Offset-Z cube
                 DIRECT_A3D_ELEM(to,  ZSIZE(to)/2 + k, j, i) =
                 DIRECT_A3D_ELEM(from, ZSIZE(from) - ZSIZE(to)/2 + k, j, i);
 
                 // Offset-Y cube
                 DIRECT_A3D_ELEM(to,  k, YSIZE(to)/2 + j, i) =
                 DIRECT_A3D_ELEM(from, k, YSIZE(from) - YSIZE(to)/2 + j, i);
             }
         }
     }
 }

◆ extractSubtomoFixedSize()

void ProgTomoExtractSubtomograms::extractSubtomoFixedSize ( MultidimArray< double > & subtomoExtraction )

Definition at line 144 of file tomo_extract_subtomograms.cpp.

 {
     // Downsampling
     FourierTransformer transformer1;
     FourierTransformer transformer2;
 
     MultidimArray< std::complex<double> > fftSubtomoExtraction;
     MultidimArray< std::complex<double> > fftSubtomo;
 
     fftSubtomo.initZeros(boxsize, boxsize, (boxsize/2)+1);
     transformer1.FourierTransform(subtomoExtraction, fftSubtomoExtraction, true);
 
     #ifdef DEBUG
     std::cout << "XSIZE(fftSubtomo) " << XSIZE(fftSubtomo) << std::endl;
     std::cout << "YSIZE(fftSubtomo) " << YSIZE(fftSubtomo) << std::endl;
     std::cout << "ZSIZE(fftSubtomo) " << ZSIZE(fftSubtomo) << std::endl;
 
     std::cout << "XSIZE(fftSubtomoExtraction) " << XSIZE(fftSubtomoExtraction) << std::endl;
     std::cout << "YSIZE(fftSubtomoExtraction) " << YSIZE(fftSubtomoExtraction) << std::endl;
     std::cout << "ZSIZE(fftSubtomoExtraction) " << ZSIZE(fftSubtomoExtraction) << std::endl;
     #endif
 
     if (downsampleFactor > 1)
     {
         downsample(fftSubtomoExtraction, fftSubtomo);
     }
     else  // downsampleFactor < 1
     {
         upsample(fftSubtomoExtraction, fftSubtomo);
     }
 
     subtomoExtraction.initZeros(1, boxsize, boxsize, boxsize);
 
     transformer2.inverseFourierTransform(fftSubtomo, subtomoExtraction);
 }

◆ normalizeSubtomo()

void ProgTomoExtractSubtomograms::normalizeSubtomo	(	MultidimArray< double > &	subtomo,
		int	halfboxsize
	)

Definition at line 88 of file tomo_extract_subtomograms.cpp.

 {
         MultidimArray<double> maskNormalize;
 
         double sumVal = 0;
         double sumVal2 = 0;
 
         auto counter = maskIdx.size();
         for (size_t i=0; i<maskIdx.size(); i++)
         {
                 double val = DIRECT_MULTIDIM_ELEM(subtomo, maskIdx[i]);
                 sumVal += val;
                 sumVal2 += val*val;
 
         }
 
         double mean;
         double sigma2;
         mean = sumVal/counter;
         sigma2 = sqrt(sumVal2/counter - mean*mean);
 
         FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(subtomo)
         {
             DIRECT_MULTIDIM_ELEM(subtomo, n) -= mean;
             DIRECT_MULTIDIM_ELEM(subtomo, n) /= sigma2;
         }
 }

◆ readParams()

void ProgTomoExtractSubtomograms::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 33 of file tomo_extract_subtomograms.cpp.

 {
     fnTom = getParam("--tomogram");
     fnCoor = getParam("--coordinates");
     boxsize = getIntParam("--boxsize");
     invertContrast = checkParam("--invertContrast");
     normalize = checkParam("--normalize");
     fnOut = getParam("-o");
     downsampleFactor = getDoubleParam("--downsample");
     nthrs = getIntParam("--threads");
     fixedBoxSize = checkParam("--fixedBoxSize");
 }

◆ run()

void ProgTomoExtractSubtomograms::run ( )

virtual

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

Reimplemented from XmippProgram.

Definition at line 180 of file tomo_extract_subtomograms.cpp.

 {
     std::cout << "Starting ... "<< std::endl;
 
     md.read(fnCoor);
 
     Image<double> tomImg;
     auto &tom = tomImg();
     tomImg.read(fnTom);
 
     size_t Xtom;
     size_t Ytom;
     size_t Ztom;
 
     Xtom = XSIZE(tom);
     Ytom = YSIZE(tom);
     Ztom = ZSIZE(tom);
 
     auto &subtomo = subtomoImg();
 
     MultidimArray<double> subtomoExtraction;
 
     int xcoor;
     int ycoor;
     int zcoor;
     int xinit;
     int yinit;
     int zinit;
 
     size_t idx=1;
 
     int halfboxsize = floor(0.5*boxsize);
 
     double invertSign = 1.0;
 
     if (invertContrast)
     {
         invertSign = -1;
     }
 
     if (normalize)
     {
         createSphere(halfboxsize);
     }
 
     double dsFactorTolerance = 0.01;
     double dsFactorDiff = abs(downsampleFactor - 1);
 
     size_t boxSizeExtraction;
     if (fixedBoxSize && dsFactorDiff > dsFactorTolerance)
     {
         #ifdef DEBUG
         std::cout << "Entering fixed box size mode" << std::endl;
         #endif
 
 
         boxSizeExtraction = boxsize * downsampleFactor;
         halfboxsize = floor(0.5*boxSizeExtraction);
     }
 
     for (const auto& row : md)
     {
         row.getValue(MDL_XCOOR, xcoor);
         row.getValue(MDL_YCOOR, ycoor);
         row.getValue(MDL_ZCOOR, zcoor);
 
         int xlim = xcoor + halfboxsize;
         int ylim = ycoor + halfboxsize;
         int zlim = zcoor + halfboxsize;
 
         xinit = xcoor - halfboxsize;
         yinit = ycoor - halfboxsize;
         zinit = zcoor - halfboxsize;
 
         if ((xlim>Xtom) || (ylim>Ytom) || (zlim>Ztom) || (xinit<0) || (yinit<0) || (zinit<0))
             continue;
 
         if (fixedBoxSize && dsFactorDiff > dsFactorTolerance)
         {
             subtomo.initZeros(1, boxSizeExtraction, boxSizeExtraction, boxSizeExtraction);
         }
         else
         {
             subtomo.initZeros(1, boxsize, boxsize, boxsize);
         }
 
         // Contrast inversion
         for (int k=zinit; k<zlim; k++)
         {
             int kk = k - zcoor;
             for (int i=yinit; i<ylim; i++)
             {
                 int ii = i-ycoor;
                 for (int j=xinit; j<xlim; j++)
                 {
                     A3D_ELEM(subtomo, kk+halfboxsize, ii+halfboxsize, j+halfboxsize-xcoor) = invertSign * A3D_ELEM(tom, k, i, j);
                 }
             }
         }
 
         if (fixedBoxSize && dsFactorDiff > dsFactorTolerance)
         {
             extractSubtomoFixedSize(subtomo);
         }
 
         // Normalization
         if (normalize)
         {
             normalizeSubtomo(subtomo, halfboxsize);
         }
 
         size_t particleid = -1;
         if (row.containsLabel(MDL_PARTICLE_ID))
         {
             row.getValue(MDL_PARTICLE_ID, particleid);
         }
 
         writeSubtomo(idx, xcoor, ycoor, zcoor, particleid);
 
         ++idx;
     }
 
 
     std::cout << "Subtomo substraction finished succesfully!!" << std::endl;
 }

◆ upsample()

void ProgTomoExtractSubtomograms::upsample	(	const MultidimArray< std::complex< double >> &	from,
		MultidimArray< std::complex< double >> &	to
	)

static

Definition at line 307 of file tomo_extract_subtomograms.cpp.

 {
     for (size_t k = 0; k < ZSIZE(from)/2; k++)
     {
         for (size_t j = 0; j < YSIZE(from)/2; j++)
         {
             for (size_t i = 0; i < XSIZE(from); i++)
             {
                 // Origin cube
                 DIRECT_A3D_ELEM(to, k, j, i) = DIRECT_A3D_ELEM(from, k, j, i);
 
                 // Oposite cube
                 DIRECT_A3D_ELEM(to,  ZSIZE(to) - ZSIZE(from)/2 + k, YSIZE(to) - YSIZE(from)/2 + j, i) =
                 DIRECT_A3D_ELEM(from, ZSIZE(from)/2 + k, YSIZE(from)/2 + j, i);
 
                 // Offset-Z cube
                 DIRECT_A3D_ELEM(to,  ZSIZE(to) - ZSIZE(from)/2 + k, j, i) =
                 DIRECT_A3D_ELEM(from, ZSIZE(from)/2 + k, j, i);
 
                 // Offset-Y cube
                 DIRECT_A3D_ELEM(to,  k, YSIZE(to) - YSIZE(from)/2 + j, i) =
                 DIRECT_A3D_ELEM(from, k, YSIZE(from)/2 + j, i);
             }
         }
     }
 }

◆ writeSubtomo()

void ProgTomoExtractSubtomograms::writeSubtomo	(	int	idx,
		int	xcoor,
		int	ycoor,
		int	zcoor,
		size_t	particleid
	)

Definition at line 116 of file tomo_extract_subtomograms.cpp.

 {
     FileName fn;
     fn = fnCoor.getBaseName() + formatString("-%i.mrc", idx);
 
     #ifdef DEBUG
     std::cout << fn << std::endl;
     std::cout << fnOut << std::endl;
     std::cout << fnOut+fn << std::endl;
     #endif
 
     subtomoImg.write(fnOut+"/"+fn);
 
     if (particleid != -1)
     {
         rowout.setValue(MDL_PARTICLE_ID, particleid);
     }
 
     rowout.setValue(MDL_XCOOR, xcoor);
     rowout.setValue(MDL_YCOOR, ycoor);
     rowout.setValue(MDL_ZCOOR, zcoor);
     rowout.setValue(MDL_IMAGE, fn);
     mdout.addRow(rowout);
 
     fn = fnCoor.getBaseName() + "_extracted.xmd";
     mdout.write(fnOut+"/"+fn);
 }