#include <tomo_simulate_tilt_series.h>

Inheritance diagram for ProgTomoSimulateTiltseries:

Collaboration diagram for ProgTomoSimulateTiltseries:

Public Member Functions
void	readParams ()
	Read argument from command line. More...

void	createFiducial (MultidimArray< double > &fidImage, int boxsize)
	create fiducial More...

void	createSphere (MultidimArray< int > &mask, int boxsize)
	create a cirte to limit the projection extension More...

void	maskingRotatedSubtomo (MultidimArray< double > &subtomo, int boxsize)

void	placeSubtomoInTomo (const MultidimArray< double > &subtomo, MultidimArray< double > &tomo, const int xcoord, const int ycoord, const int zcood, const size_t halfboxsize)

void	createFiducial (MultidimArray< double > &fidImage, MultidimArray< double > &fidVol, int fidSize)

void	defineParams ()
	Define parameters. More...

void	run ()
	Run. 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)

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	fnCoords

FileName	fnVol

FileName	fnFid

double	maxTilt

double	minTilt

double	tiltStep

double	sigmaNoise

double	fidDiameter

FileName	fnTsOut

FileName	fnTomoOut

int	xdim

int	ydim

int	nfids

int	thickness

double	sampling

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

Movie alignment correlation Parameters.

Definition at line 38 of file tomo_simulate_tilt_series.h.

Member Function Documentation

◆ createFiducial() [1/2]

void ProgTomoSimulateTiltseries::createFiducial	(	MultidimArray< double > &	fidImage,
		int	boxsize
	)

create fiducial

◆ createFiducial() [2/2]

void ProgTomoSimulateTiltseries::createFiducial	(	MultidimArray< double > &	fidImage,
		MultidimArray< double > &	fidVol,
		int	fidSize
	)

Definition at line 74 of file tomo_simulate_tilt_series.cpp.

 {
     fidImage.initZeros(fidSize, fidSize);
     fidVol.initZeros(fidSize, fidSize, fidSize);
     auto halfbox = round(0.5*fidSize);
     int fidSizeVoxels;
 
     auto fidSize2 = fidSize*fidSize;
 
     // We use 5*sigmaNoise to ensure the fiducial is visible with this noise level
 
     for (size_t i= 0; i<fidSize; i++)
     {
         auto i2 = (i-halfbox)*(i-halfbox);
         for (size_t j= 0; j<fidSize; j++)
         {
             auto j2 = (j-halfbox)*(j-halfbox);
 
             if ((j2 + i2) < fidSize2)
                 A2D_ELEM(fidImage, i, j) = 5*sigmaNoise;
         }
     }
 
     for (size_t k= 0; k<fidSize; k++)
     {
         auto k2 = (k-halfbox)*(k-halfbox);
         for (size_t i= 0; i<fidSize; i++)
         {
             auto i2 = k2 + (i-halfbox)*(i-halfbox);
             for (size_t j= 0; j<fidSize; j++)
             {
                 auto j2 = i2 + (j-halfbox)*(j-halfbox);
 
                 if (i2 < fidSize2)
                     A3D_ELEM(fidVol, k, i, j) = 5*sigmaNoise;
             }
         }
     }
 
 
 }

◆ createSphere()

void ProgTomoSimulateTiltseries::createSphere	(	MultidimArray< int > &	mask,
		int	boxsize
	)

create a cirte to limit the projection extension

Definition at line 116 of file tomo_simulate_tilt_series.cpp.

 {
     int halfbox = round(0.5*boxsize);
     auto halfbox2 =halfbox*halfbox;
 
     std::cout << "halfbox = " << halfbox << std::endl;
 
     size_t idx = 0;
     long n=0;
     for (int k= 0; k<boxsize; k++)
     {
         int k2 = (k-halfbox)*(k-halfbox);
         std::cout << "k2 = " << k2 << std::endl;
         for (int i= 0; i<boxsize; i++)
         {
             int i2 = (i-halfbox)*(i-halfbox);
             int i2k2 = i2 + k2;
             for (int j= 0; j<boxsize; j++)
             {
                 int j2 = (j-halfbox)*(j-halfbox);
 
                 const bool inside = i2k2+j2 <= halfbox2;
                 A3D_ELEM(mask, k, i, j) = static_cast<int>(inside);
                 n++;
             }
         }
     }
 }

◆ defineParams()

void ProgTomoSimulateTiltseries::defineParams ( )

virtual

Define parameters.

Reimplemented from XmippProgram.

Definition at line 54 of file tomo_simulate_tilt_series.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("  --coordinates <xmd_file=\"\">           : Metadata (.xmd file) with the tilt series");
     addParamsLine("  --vol <xmd_file=\"\">                   : Volume (.mrc file) with the coordidanates to be extracted from the tomogram");
     addParamsLine("[  --xdim <xdim=1000>]                    : Particle box size in voxels.");
     addParamsLine("[  --ydim <ydim=1000> ]                   : Particle box size in voxels.");
     addParamsLine("[  --minTilt <minTilt=-60.0>]             : Minimum tilt angle");
     addParamsLine("[  --maxTilt <maxTilt=60.0>]              : Maximum tilt angle");
     addParamsLine("[  --tiltStep <tiltStep=60.0>]            : Step of the angular sampling. Distance between two consequtive tilt angles");
     addParamsLine("[  --thickness <thickness=300>]           : Tomogram thickness in pixel");
     addParamsLine("[  --sampling <s=1>]                      : Sampling rate in (A).");
     addParamsLine("[  --fiducialCoordinates <xmd_file=\"\">] : Metadata (.xmd file) with fiducial coordinates in the tomogram");
     addParamsLine("[  --fiducialDiameter <xmd_file=\"\">]    : Fiducial diameter in (A)");
     addParamsLine("[  --sigmaNoise <sigmaNoise=-1>]          : Standard deviation of the noise.");
     addParamsLine("  --tiltseries  <mrc_file=\"\">           : path to the output directory. ");
     addParamsLine("  --tomogram  <mrc_file=\"\">             : path to the output directory. ");
 }

◆ maskingRotatedSubtomo()

void ProgTomoSimulateTiltseries::maskingRotatedSubtomo	(	MultidimArray< double > &	subtomo,
		int	boxsize
	)

Definition at line 145 of file tomo_simulate_tilt_series.cpp.

 {
     //subtomo.resetOrigin();
     auto halfbox = round(0.5*boxsize);
     auto halfbox2 =halfbox*halfbox;
 
     std::cout << " boxsize = " << boxsize << "    halfboxsize = " << halfbox << std::endl;
 
     long n=0;
     for (size_t k= 0; k<boxsize; k++)
     {
         auto k2 = (k-halfbox)*(k-halfbox);
         for (size_t j= 0; j<boxsize; j++)
         {
             auto j2 = k2+(j-halfbox)*(j-halfbox);
             for (size_t i= 0; i<boxsize; i++)
             {
                 auto i2 = (i-halfbox)*(i-halfbox);
                 if (i2+j2>halfbox2)
                 {
                     DIRECT_MULTIDIM_ELEM(subtomo, n) = 0.0;
 
                 }
                 n++;
             }
         }
     }
 }

◆ placeSubtomoInTomo()

void ProgTomoSimulateTiltseries::placeSubtomoInTomo	(	const MultidimArray< double > &	subtomo,
		MultidimArray< double > &	tomo,
		const int	xcoord,
		const int	ycoord,
		const int	zcood,
		const size_t	halfboxsize
	)

Definition at line 175 of file tomo_simulate_tilt_series.cpp.

 {
     auto initzcoord = thickness/2 + zcoord-boxsize/2;
     auto initxcoord = xdim/2 + xcoord-boxsize/2;
     auto initycoord = ydim/2 + ycoord-boxsize/2;
 
     for (int k = initzcoord; k<(initzcoord+boxsize); k++)
     {
         for (int i = initycoord; i<(initycoord+boxsize); i++)
         {
             for (int j = initxcoord; j<(initxcoord+boxsize); j++)
             {
                 DIRECT_A3D_ELEM(tomo, k, i, j) =  -DIRECT_A3D_ELEM(subtomo, k-initzcoord, i-initycoord, j-initxcoord);
             }
         }
     }
 }

◆ readParams()

void ProgTomoSimulateTiltseries::readParams ( )

virtual

Read argument from command line.

Reimplemented from XmippProgram.

Definition at line 35 of file tomo_simulate_tilt_series.cpp.

 {
     fnCoords = getParam("--coordinates");
     fnVol = getParam("--vol");
     xdim = getIntParam("--xdim");
     ydim = getIntParam("--ydim");
     minTilt = getDoubleParam("--minTilt");
     maxTilt = getDoubleParam("--maxTilt");
     tiltStep = getDoubleParam("--tiltStep");
     thickness = getIntParam("--thickness");
     sampling = getDoubleParam("--sampling");
     fnFid = getParam("--fiducialCoordinates");
     fidDiameter = getDoubleParam("--fiducialDiameter");
     sigmaNoise = getDoubleParam("--sigmaNoise");
     fnTsOut = getParam("--tiltseries");
     fnTomoOut = getParam("--tomogram");
 }

◆ run()

void ProgTomoSimulateTiltseries::run ( )

virtual

Run.

Reimplemented from XmippProgram.

Definition at line 195 of file tomo_simulate_tilt_series.cpp.

 {
     std::cout << "Starting ... " << std::endl;
 
     // Reading the input map
     Image<double> proteinImg;
     auto &ptrVol = proteinImg();
     proteinImg.read(fnVol);
 
     MultidimArray<double> rotatedVol;   // This map is the rotated volume
     rotatedVol.resizeNoCopy(ptrVol);
 
     // Defining input coordinates and subtomogram orientation
     MetaDataVec mdCoords;
     int xcoor, ycoor, zcoor;
     double rot =0.0, tilt=0.0, psi=0, sx=0, sy=0;
     Matrix2D<double> eulerMat_proj, eulerMat_VolRotation, eulerMat;
     Projection imgPrj;
 
     MultidimArray<double> &ptrProj = imgPrj();
 
     int BSplinedegree = 3;
     MultidimArray<double> fidImage;
     MultidimArray<double> fidVol;
 
     if (fnFid != "")
     {
         int fidBoxsize;
         fidBoxsize = round(fidDiameter/sampling);
         createFiducial(fidImage, fidVol, fidBoxsize);
     }
 
 
     // Estimating the number of projections adn ouput tilt series
     std::vector<double> tiltAngles;
     Image<double> tiltseriesImg;
     MultidimArray<double> &tiltseries = tiltseriesImg();
     int numberOfProjections;
     numberOfProjections = (maxTilt - minTilt)/tiltStep + 1;
 
     size_t boxsize, halfboxsize;
     boxsize = XSIZE(ptrVol);
     halfboxsize = round(0.5*boxsize);
     tiltseries.initZeros(numberOfProjections, 1, ydim, xdim+2*boxsize);
 
     for (size_t idx = 0; idx<numberOfProjections; idx++)
     {
         tiltAngles.push_back(minTilt+idx*tiltStep);
     }
 
     MultidimArray<int> mask;
     mask.resizeNoCopy(ptrVol);
 
     // Reading the coordinates and generating tilt series, tomograms and coordinates.
     MetaDataVec mdTiltSeries;
     mdCoords.read(fnCoords);
     MDRowVec row, rowOut;
 
     // This is the output tomogram
     Image<double> tomogram;
     auto &ptrTomo = tomogram();
     ptrTomo.initZeros(1, thickness, ydim, xdim);
 
     std::cout << "Starting rows " << std::endl;
     int addCoords = 0;
 
     for (const auto& row : mdCoords)
     {
         row.getValue(MDL_XCOOR, xcoor);
         row.getValue(MDL_YCOOR, ycoor);
         row.getValue(MDL_ZCOOR, zcoor);
 
         double theta, phi, xi;
 
 
 
         if (row.containsLabel(MDL_ANGLE_ROT) && row.containsLabel(MDL_ANGLE_TILT) && row.containsLabel(MDL_ANGLE_PSI))
         {
 //          row.getValue(MDL_ANGLE_ROT, theta);
 //          row.getValue(MDL_ANGLE_TILT, phi);
 //          row.getValue(MDL_ANGLE_PSI, xi);
 
             row.getValue(MDL_ANGLE_PSI, theta);
             row.getValue(MDL_ANGLE_TILT, phi);
             row.getValue(MDL_ANGLE_ROT, xi);
 
 //          theta *= -1;
 //          phi  *= -1;
 //          xi  *= -1;
         }
         else
         {
             double u = (double) rand()/RAND_MAX;
             double v = (double) rand()/RAND_MAX;
             double w = (double) rand()/RAND_MAX;
 
             theta = 360.0*u;
             phi = acos(2*v - 1.0)*180.0/PI;
             xi = 360*w;
 
             //eulerMat_proj.initIdentity(4);
         }
 
         Euler_angles2matrix(theta, phi, xi, eulerMat_VolRotation, true);
 
         applyGeometry(xmipp_transformation::BSPLINE3, rotatedVol, ptrVol, eulerMat_VolRotation, xmipp_transformation::IS_NOT_INV, true, 0.);
         maskingRotatedSubtomo(rotatedVol, boxsize);
 
         rotatedVol.setXmippOrigin();
         //CenterFFT(rotatedVol,true);
         placeSubtomoInTomo(rotatedVol, ptrTomo, xcoor, ycoor, zcoor, boxsize);
 
 
         auto projector = FourierProjector(rotatedVol, 2.0, 0.5, BSplinedegree);
 
         for (size_t idx = 0; idx<tiltAngles.size(); idx++)
         {
             double tiltProj = tiltAngles[idx];
             if (addCoords<1)
             {
                 rowOut.setValue(MDL_ANGLE_TILT, tiltProj);
                 FileName fnTi;
                 fnTi = formatString("%i@", idx+1);
                 rowOut.setValue(MDL_IMAGE, fnTi + fnTsOut);
                 mdTiltSeries.addRow(rowOut);
             }
 
             projectVolume(projector, imgPrj, boxsize, boxsize, 0.0, tiltProj, 0.0);
 
             Matrix1D<double> shifts(2);
             tiltProj *= PI/180.0;
 
             double ct = cos(tiltProj);
             double st = sin(tiltProj);
 
 
             // Projection
             auto x_2d = (int) (xcoor * ct + zcoor* st);
             auto y_2d = (int) (ycoor);
 
             y_2d = y_2d+0.5*ydim;
             x_2d = x_2d+0.5*xdim;
 
             int xlim = x_2d + halfboxsize;
             int ylim = y_2d + halfboxsize;
 
             int xinit = x_2d - halfboxsize;
             int yinit = y_2d - halfboxsize;
 
             if ((xlim>xdim) || (ylim>ydim) || (xinit<0) || (yinit<0))
             {
                 //std::cout << "Particle out of the tilt image" << std::endl;
                 continue;
             }
 
             for (int i=0; i<boxsize; i++)
             {
                 if (yinit+i<ydim)
                 {
                     for (int j=0; j<boxsize; j++)
                     {
                         if (xinit+j<xdim)
                         {
                             NZYX_ELEM(tiltseries, idx, 0, yinit+i, xinit+j) +=  DIRECT_A2D_ELEM(ptrProj, i, j);
                         }
                     }
                 }
             }
         }
         addCoords++;
     }
 
     MultidimArray<double> tiltseriesOut;
     tiltseriesOut.initZeros(numberOfProjections, 1, ydim, xdim);
 
     if (sigmaNoise>0)
     {
         std::random_device rd{};
         std::mt19937 gen{rd()};
 
         // values near the mean are the most likely
         // standard deviation affects the dispersion of generated values from the mean
         std::normal_distribution dts{0.0, sigmaNoise};
         std::normal_distribution dtom{0.0, sigmaNoise/boxsize};
 
         // draw a sample from the normal distribution and round it to an integer
         auto randomNumberTs = [&dts, &gen]{ return dts(gen); };
         auto randomNumberTom = [&dtom, &gen]{ return dtom(gen); };
 
 
         for (size_t idx = 0; idx<tiltAngles.size(); idx++)
         {
             for (int i=0; i<ydim; i++)
             {
                 for (int j=0; j<xdim; j++)
                 {
                     NZYX_ELEM(tiltseriesOut, idx, 0, i, j) =  -NZYX_ELEM(tiltseries, idx, 0, i, j)+randomNumberTs();
                 }
             }
         }
 
         FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(ptrTomo)
         {
             DIRECT_MULTIDIM_ELEM(ptrTomo, n) += randomNumberTom();
         }
     }
     else
     {
 
         for (size_t idx = 0; idx<tiltAngles.size(); idx++)
         {
             for (int i=0; i<ydim; i++)
             {
                 for (int j=0; j<xdim; j++)
                 {
                     NZYX_ELEM(tiltseriesOut, idx, 0, i, j) =  -NZYX_ELEM(tiltseries, idx, 0, i, j);
                 }
             }
         }
     }
 
     tiltseriesImg() = tiltseriesOut;
     tiltseriesImg.write(fnTsOut);
     FileName fnOutXmd;
     fnOutXmd = fnTsOut.removeLastExtension() + ".xmd";
     tomogram.write(fnTomoOut);
     std::cout << fnOutXmd << std::endl;
     mdTiltSeries.write(fnOutXmd);
 }