#include <tomo_extract_particlestacks.h>

Inheritance diagram for ProgTomoExtractParticleStacks:

Collaboration diagram for ProgTomoExtractParticleStacks:

Public Member Functions
void	defineParams ()

void	createCircle (MultidimArray< double > &maskNormalize)

void	readTiltSeriesInfo (std::string &tsid)

void	getCoordinateOnTiltSeries (int xcoor, int ycoor, int zcoor, double &rot, double &tilt, double &tx, double &ty, int &x_2d, int &y_2d)

void	readParams ()

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	fnTs

FileName	fnCoor

size_t	Xdim

size_t	Ydim

size_t	Zdim

size_t	Xts

size_t	Yts

double	tilt

double	rot

double	tx

double	ty

bool	invertContrast

bool	normalize

bool	swapXY

bool	setCTF

bool	defocusPositive

std::vector< double >	tsTiltAngles

std::vector< double >	tsRotAngles

std::vector< double >	tsShiftX

std::vector< double >	tsShiftY

std::vector< double >	tsDefU

std::vector< double >	tsDefV

std::vector< double >	tsDefAng

std::vector< double >	tsDose

std::vector< MultidimArray< double > >	tsImages

double	scaleFactor

double	sampling

int	boxsize

int	nthrs

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 37 of file tomo_extract_particlestacks.h.

Member Function Documentation

◆ createCircle()

void ProgTomoExtractParticleStacks::createCircle ( MultidimArray< double > & maskNormalize )

Definition at line 65 of file tomo_extract_particlestacks.cpp.

 {
     int halfboxsize = 0.5* boxsize;
     if (normalize)
     {
         
         maskNormalize.initZeros(1, 1, boxsize, boxsize);
 
         for (int i=0; i<boxsize; i++)
         {
             int i2 = i-halfboxsize;
             int i2k2 = i2*i2;
             for (int j=0; j<boxsize; j++)
             {
                 int j2 = (j- halfboxsize);
                 if (sqrt(i2k2 + j2*j2)>halfboxsize)
                     A2D_ELEM(maskNormalize, i, j) = 1;
             }
         }
     }
 }

◆ defineParams()

void ProgTomoExtractParticleStacks::defineParams ( )

virtual

Function in which the param of each Program are defined.

Reimplemented from XmippProgram.

Definition at line 48 of file tomo_extract_particlestacks.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("  --tiltseries <xmd_file=\"\">       : Metadata (.xmd file) with the tilt series");
     addParamsLine("  --coordinates <xmd_file=\"\">      : Metadata (.xmd file) with the coordidanates to be extracted from the tomogram");
     addParamsLine("  --boxsize <boxsize=100>            : Particle box size in voxels.");
     addParamsLine("  --sampling <s=1>                   : Sampling rate in (A).");
     addParamsLine("  [--defocusPositive]                : This flag must be put if the defocus increases or decreases along the z-axis. This is requires to set the local CTF");
     addParamsLine("  [--invertContrast]                 : Put this flag if the particles to be extracted are 3D particles (subtvolumes)");
     addParamsLine("  [--swapXY]                         : Put this flag if the tomogram and the tilt series have the same dimensions but the X and Y coordinates are swaped");
     addParamsLine("  [--setCTF]                         : Put this flag if the tilt series metadata has CTF parameters. The CTF per particle will be calculated and set in the final set of particles");
     addParamsLine("  [--normalize]                      : Put this flag to normalize the background of the particle (zero mean and unit std)");
     // addParamsLine("  [--downsample <scaleFactor=0.5>]   : Scale factor of the extracted subtomograms");
     addParamsLine("  -o <mrc_file=\"\">                 : path to the output directory. ");
     addParamsLine("  [--threads <s=4>]                  : Number of threads");
 }

◆ getCoordinateOnTiltSeries()

void ProgTomoExtractParticleStacks::getCoordinateOnTiltSeries	(	int	xcoor,
		int	ycoor,
		int	zcoor,
		double &	rot,
		double &	tilt,
		double &	tx,
		double &	ty,
		int &	x_2d,
		int &	y_2d
	)

Definition at line 87 of file tomo_extract_particlestacks.cpp.

 {
     double ct = cos(tilt);
     double st = sin(tilt);
 
     double cr = cos(rot);
     double sr = sin(rot);
 
     /*
     First the picked coordinate, r, is projected on the aligned tilt series
     r' = Pr  Where r is the coordinates to be projected by the matrix P
     [x']   [ct   0   st][x]
     [y'] = [0    1    0][y]
     [z']   [0    0    0][z]
     Next is to undo the transformation This is Aligned->Unaligned
     If T is the transformation matrix unaligned->aligned, we need T^-{1}
     Let us define a rotation matrix
     R=[cos(rot) -sin(rot)]
       [ sin(rot) cos(rot)];
 
     The inverse of T is given by
     T^{-1} = [R' -R'*t;
              0  0 1]);
              
      Where t are the shifts of T
     */
 
     // Projection
     x_2d = (int) (xcoor * ct + zcoor* st);
     y_2d = (int) (ycoor);
 
     //Inverse transformation
     double x_2d_prime =   cr*x_2d  + sr*y_2d - cr*tx  - sr*ty;
     double y_2d_prime =  -sr*x_2d  + cr*y_2d + sr*tx - cr*ty;
         
     if (swapXY)
     {
         x_2d = -x_2d_prime+0.5*Xts;
         y_2d = -y_2d_prime+0.5*Yts;
     }
     else
     {
         x_2d = -x_2d_prime+0.5*Yts;
         y_2d = -y_2d_prime+0.5*Xts;
     }
 }

◆ readParams()

void ProgTomoExtractParticleStacks::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 31 of file tomo_extract_particlestacks.cpp.

 {
     fnTs = getParam("--tiltseries");
     fnCoor = getParam("--coordinates");
     boxsize = getIntParam("--boxsize");
     sampling = getDoubleParam("--sampling");
     defocusPositive = checkParam("--defocusPositive");
     invertContrast = checkParam("--invertContrast");
     // scaleFactor = getDoubleParam("--downsample");
     normalize = checkParam("--normalize");
     swapXY = checkParam("--swapXY");
     setCTF = checkParam("--setCTF");
     fnOut = getParam("-o");
     nthrs = getIntParam("--threads");
 }

◆ readTiltSeriesInfo()

void ProgTomoExtractParticleStacks::readTiltSeriesInfo ( std::string & tsid )

Definition at line 134 of file tomo_extract_particlestacks.cpp.

 {
     //TODO: Ensure there is only a single TSId
     MetaDataVec mdts;
     mdts.read(fnTs);
 
     Image<double> tiltImg;
     auto &ptrtiltImg = tiltImg();
 
     // The tilt series is stored as a stack of images;
     MultidimArray<double> tsImg;
     std::vector<FileName> tsNames(0);
 
     FileName fnImg;
 
     size_t Nimages = 0;
     double defU=0, defV=0, defAng=0, dose = 0;
 
     for (const auto& row : mdts)
     {
         row.getValue(MDL_IMAGE, fnImg);
         row.getValue(MDL_ANGLE_TILT, tilt);
         row.getValue(MDL_ANGLE_ROT, rot);
 
         if (setCTF)
         {
             row.getValue(MDL_CTF_DEFOCUSU, defU);
             row.getValue(MDL_CTF_DEFOCUSV, defV);
             row.getValue(MDL_CTF_DEFOCUS_ANGLE, defAng);
             row.getValue(MDL_CTF_DEFOCUS_ANGLE, dose);
             tsDefU.push_back(defU);
             tsDefV.push_back(defV);
             tsDefAng.push_back(defAng);
 
             tsDose.push_back(dose);
         }
 
 
         row.getValue(MDL_SHIFT_X, tx);
         row.getValue(MDL_SHIFT_Y, ty);
         row.getValue(MDL_TSID, tsid);
 
         tiltImg.read(fnImg);
 
         tsImages.push_back(ptrtiltImg);
 
         tsTiltAngles.push_back(tilt);
         tsRotAngles.push_back(rot);
         tsShiftX.push_back(tx);
         tsShiftY.push_back(ty);
         tsNames.push_back(fnImg);
         Nimages +=1;
     }
 
     Xts = XSIZE(ptrtiltImg);
     Yts = YSIZE(ptrtiltImg);
 }

◆ run()

void ProgTomoExtractParticleStacks::run ( )

virtual

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

Reimplemented from XmippProgram.

Definition at line 194 of file tomo_extract_particlestacks.cpp.

 {
     std::cout << "Starting ... "<< std::endl;
     double signValue = 1;
     if (invertContrast)
     {
         signValue = -1;
     }
 
     MetaDataVec mdcoords, mdparticlestack;
     mdcoords.read(fnCoor);
     MetaDataVec mdts;
 
     int xcoor, ycoor, zcoor, xinit, yinit;
     size_t idx=1;
 
     mdts.read(fnTs);
 
     std::string tsid;
     readTiltSeriesInfo(tsid);
 
     double tilt, rot, tx, ty;
 
     size_t Nimages = 0;
 
     int halfboxsize = floor(0.5*boxsize);
     Image<double> finalStack;
     auto &particlestack = finalStack();
 
     MultidimArray<double> maskNormalize;
     if (normalize)
     {
         createCircle(maskNormalize);
     }
 
     size_t elem = 0;
     double signDef = -1.0;
 
     FileName fnXmd;
     fnXmd = tsid + formatString(".xmd");
 
     for (const auto& row : mdcoords)
     {
         row.getValue(MDL_XCOOR, xcoor);
         row.getValue(MDL_YCOOR, ycoor);
         row.getValue(MDL_ZCOOR, zcoor);
         
         std::vector<MultidimArray<double>> imgVec;
         MultidimArray<double> singleImage;
         singleImage.initZeros(1,1,boxsize, boxsize);
 
         FileName fnMrc;
         fnMrc = tsid + formatString("-%i.mrcs", elem);
 
         size_t imgNumber = 0;
         for (size_t idx = 0; idx<tsImages.size(); idx++)
         {
             auto tsImg = tsImages[idx];
             tilt = tsTiltAngles[idx]*PI/180;
             rot = tsRotAngles[idx]*PI/180;
             tx = tsShiftX[idx];
             ty = tsShiftY[idx];
             int x_2d, y_2d;
             getCoordinateOnTiltSeries(xcoor, ycoor, zcoor, rot, tilt, tx, ty, x_2d, y_2d);
 
             int xlim = x_2d + halfboxsize;
             int ylim = y_2d + halfboxsize;
 
             xinit = x_2d - halfboxsize;
             yinit = y_2d - halfboxsize;
 
             if ((xlim>Xts) || (ylim>Yts) || (xinit<0) || (yinit<0))
             {
                 continue;
             }
 
             for (int i=yinit; i<ylim; i++)
             {
                 int ii = i-y_2d;
                 for (int j=xinit; j<xlim; j++)
                 {
                     A2D_ELEM(singleImage, ii+halfboxsize, j+halfboxsize-x_2d) = signValue*A2D_ELEM(tsImages[idx], i, j);
                 }
             }
 
             if (normalize)
             {
                 double sumVal = 0, sumVal2 = 0, counter = 0;
 
                 long n = 0;
                 for (int i=0; i<boxsize; i++)
                 {
                     for (int j=0; j<boxsize; j++)
                     {
                         if (DIRECT_MULTIDIM_ELEM(maskNormalize, n)>0)
                         {
                             double val = A2D_ELEM(singleImage, i, j);
                             sumVal += val;
                             sumVal2 += val*val;
                             counter = counter + 1;
                         }
                         n++;
 
                     }
                 }
 
                 double mean, sigma2;
                 mean = sumVal/counter;
                 sigma2 = sqrt(sumVal2/counter - mean*mean);
 
                 for (int i=0; i<boxsize; i++)
                 {
                     for (int j=0; j<boxsize; j++)
                     {
                         A2D_ELEM(singleImage, i, j) -= mean;
                         A2D_ELEM(singleImage, i, j) /= sigma2;
                     }
                 }
             }
 
             imgVec.push_back(singleImage);
             MDRowVec rowParticleStack;
             rowParticleStack.setValue(MDL_TSID, tsid);
             FileName idxstr;
             idxstr = formatString("%i@",imgNumber+1);
             rowParticleStack.setValue(MDL_IMAGE, idxstr+fnMrc);
             rowParticleStack.setValue(MDL_ANGLE_TILT, tsTiltAngles[idx]);
             rowParticleStack.setValue(MDL_ANGLE_ROT, tsRotAngles[idx]);
             rowParticleStack.setValue(MDL_SHIFT_X, tsShiftX[idx]);
             rowParticleStack.setValue(MDL_SHIFT_Y, tsShiftY[idx]);
             rowParticleStack.setValue(MDL_DOSE, tsDose[idx]);
 
             double defU=0, defV=0, defAng=0;
             if (setCTF)
             {
 //              if (defocusPositive)
 //              {
 //                  signDef = 1.0;
 //              }
 
                 double Df = (xcoor * cos(tilt) + zcoor* sin(tilt))*sampling*sin(tilt);
 
                 defU = tsDefU[idx]  + signDef*Df;
                 defV = tsDefV[idx]  + signDef*Df;
             }
             rowParticleStack.setValue(MDL_CTF_DEFOCUSU, defU);
             rowParticleStack.setValue(MDL_CTF_DEFOCUSV, defV);
             rowParticleStack.setValue(MDL_CTF_DEFOCUS_ANGLE, tsRotAngles[idx]);
             rowParticleStack.setValue(MDL_XCOOR, x_2d);
             rowParticleStack.setValue(MDL_YCOOR, y_2d);
 
             mdparticlestack.addRow(rowParticleStack);
             imgNumber++;
         }
         mdparticlestack.write(fnOut+"/"+fnXmd);
         Nimages = imgVec.size();
         particlestack.initZeros(Nimages, 1, boxsize, boxsize);
         for (int k=0; k<Nimages; k++)
         {
             singleImage = imgVec[k];
             for (int i=0; i<boxsize; i++)
             {
                 for (int j=0; j<boxsize; j++)
                 {
                     NZYX_ELEM(particlestack, k, 0, i, j) = A2D_ELEM(singleImage, i, j);
                 }
             }
         }
 
 
 
         finalStack.write(fnOut+"/"+fnMrc);
 
         elem += 1;
 
     }
 
 
 }