#include <phantom_simulate_microscope.h>

Inheritance diagram for ProgSimulateMicroscope:

Collaboration diagram for ProgSimulateMicroscope:

Public Member Functions
void	readParams ()

void	defineParams ()

void	show ()

void	estimateSigma ()

void	setupFourierFilter (FourierFilter &filter, bool isBackground, double &power)

void	updateCtfs ()

void	preProcess ()

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

void	apply (MultidimArray< double > &I)

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	fn_ctf
	Filename with the CTF. More...

double	sigma
	Total noise power. More...

double	low_pass_before_CTF
	Low pass frequency before CTF. More...

bool	after_ctf_noise
	Filename with the root squared spectrum for noise after CTF. More...

double	defocus_change
	Defocus change (%) More...

double	targetSNR
	Target SNR. More...

bool	estimateSNR
	Estimate SNR. More...

bool	CTFpresent
	CTFpresent. More...

FourierFilter	ctf
	CTF. More...

FourierFilter	lowpass
	Low pass filter, if it is 0 no lowpass filter is applied. More...

FourierFilter	after_ctf
	After CTF noise root squared spectrum. More...

double	sigma_before_CTF
	Noise power before CTF. More...

double	sigma_after_CTF
	Noise power after CTF. More...

size_t	Xdim
	Input image Xdim. More...

size_t	Ydim
	Input image Ydim. More...

MetaDataVec *	pmdIn
	Particular reference to mdIn to manipulated. More...

double	downsampling

double	defocusU

double	defocusV

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

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

Parameter class for the project program

Definition at line 36 of file phantom_simulate_microscope.h.

Member Function Documentation

◆ apply()

void ProgSimulateMicroscope::apply ( MultidimArray< double > & I )

Apply to a single image. The image is modified. If the CTF is randomly selected then a new CTF is generated for each image

Definition at line 295 of file phantom_simulate_microscope.cpp.

 {
     I.setXmippOrigin();
     I.selfWindow(FIRST_XMIPP_INDEX(2*Ydim), FIRST_XMIPP_INDEX(2*Xdim),
                  LAST_XMIPP_INDEX(2*Ydim), LAST_XMIPP_INDEX(2*Xdim));
 
     // Add noise before CTF
     MultidimArray<double> noisy;
     noisy.resize(I);
     noisy.initRandom(0, sigma_before_CTF, RND_GAUSSIAN);
     if (low_pass_before_CTF < 0.5)
         lowpass.applyMaskSpace(noisy);
     I += noisy;
 
     if (CTFpresent)
     {
         // Check if the mask is a defocus changing CTF
         // In that case generate a new mask with a random defocus
         if (defocus_change != 0)
         {
             MultidimArray<double> aux;
             ctf.ctf.DeltafU = defocusU * rnd_unif(1 - defocus_change / 100, 1 + defocus_change / 100);
             ctf.ctf.DeltafV = defocusV *rnd_unif(1 - defocus_change / 100, 1 + defocus_change / 100);
             aux.initZeros(2*Ydim, 2*Xdim);
             ctf.generateMask(aux);
         }
 
         // Apply CTF
         ctf.applyMaskSpace(I);
 
         // Add noise after CTF
         noisy.initRandom(0, sigma_after_CTF, RND_GAUSSIAN);
         if (after_ctf_noise)
             after_ctf.applyMaskSpace(noisy);
         I += noisy;
     }
     I.selfWindow(FIRST_XMIPP_INDEX(Ydim), FIRST_XMIPP_INDEX(Xdim),
                  LAST_XMIPP_INDEX(Ydim), LAST_XMIPP_INDEX(Xdim));
 }

◆ defineParams()

void ProgSimulateMicroscope::defineParams ( )

virtual

Usage message. This function shows the way of introducing these parameters.

Reimplemented from XmippMetadataProgram.

Definition at line 85 of file phantom_simulate_microscope.cpp.

 {
     each_image_produces_an_output = true;
     save_metadata_stack = true;
     XmippMetadataProgram::defineParams();
     addUsageLine("Simulate the effect of the microscope on ideal projections.");
     addParamsLine("==CTF options==");
     addParamsLine(" [--ctf <CTFdescr=\"\">]       : a CTF description, if this param is not supplied it should come in metadata");
     addParamsLine(" [--after_ctf_noise]      : apply noise after CTF");
     addParamsLine(" [--defocus_change <v=0>] : change in the defocus value (percentage)");
     addParamsLine("==Noise options==");
     addParamsLine(" --noise <stddev> <w=0.5> : noise to be added, this noise is filtered at the frequency specified (<0.5).");
     addParamsLine("or --targetSNR <snr>      : the necessary noise power for a specified SNR is estimated");
     addParamsLine("or --noNoise              : do not add any noise, only simulate the CTF");
     addParamsLine(" [--downsampling <D=1>]  : Downsampling factor of the input micrograph with respect to the original");
     addParamsLine("                         : micrograph.");
     addExampleLine("Generate a set of images with the CTF applied without any noise", false);
     addExampleLine("   xmipp_phantom_simulate_microscope -i g0ta.sel --oroot g1ta --ctf untilt_ARMAavg.ctfparam");
     addExampleLine("Generate a set of images with a target SNR", false);
     addExampleLine("   xmipp_phantom_simulate_microscope -i g0ta.sel --oroot g2ta --ctf untilt_ARMAavg.ctfparam --targetSNR 0.2 --after_ctf_noise");
     addExampleLine("Generate a set of images with the CTF applied and noise before and after CTF", false);
     addExampleLine("   xmipp_phantom_simulate_microscope -i g0ta.sel --oroot g2ta --ctf untilt_ARMAavg.ctfparam --noise 4.15773 --after_ctf_noise");
 }

◆ estimateSigma()

void ProgSimulateMicroscope::estimateSigma ( )

Estimate sigma for a given SNR

Definition at line 130 of file phantom_simulate_microscope.cpp.

 {
     if (CTFpresent)
     {
         ctf.FilterBand = CTF;
         ctf.ctf.enable_CTFnoise = false;
         ctf.ctf.read(getParam("--ctf"));
         ctf.ctf.changeSamplingRate(ctf.ctf.Tm*downsampling);
         ctf.ctf.produceSideInfo();
     }
 
     size_t N_stats = pmdIn->size();
     std::cout << "N_stats=" << N_stats << std::endl;
 
     MultidimArray<double> proj_power(N_stats);
     MultidimArray<double> proj_area(N_stats);
     double power_avg=0., power_stddev, area_avg=0., area_stddev, avg, dummy;
     if (verbose!=0)
     {
         std::cerr << "Estimating noise power for target SNR=" << targetSNR << std::endl;
         init_progress_bar(N_stats);
     }
     FileName fnImg;
     Image<double> proj;
     size_t nImg=0;
     for (size_t objId : pmdIn->ids())
     {
         pmdIn->getValue(image_label, fnImg,objId);
         proj.read(fnImg);
         MultidimArray<double> mProj=proj();
 
         if (CTFpresent)
         {
             if (nImg == 0)
                 ctf.generateMask(mProj);
             ctf.applyMaskSpace(mProj);
         }
 
         // Compute projection area
         DIRECT_A1D_ELEM(proj_area,nImg) = 0;
         FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(mProj)
         if (fabs(DIRECT_MULTIDIM_ELEM(mProj,n)) > 1e-6)
             DIRECT_A1D_ELEM(proj_area,nImg)++;
 
         // Compute projection power
         mProj.computeStats(avg, DIRECT_A1D_ELEM(proj_power,nImg), dummy, dummy);
 
         if (nImg++ % 30 == 0 && verbose!=0)
             progress_bar(nImg);
     }
     progress_bar(N_stats);
     Histogram1D hist_proj, hist_area;
     compute_hist(proj_power, hist_proj, 300);
     compute_hist(proj_area, hist_area, 300);
     proj_power.computeStats(power_avg, power_stddev, dummy, dummy);
     std::cout << "# Projection power average: " << power_avg << std::endl
     << "# Projection power stddev:  " << power_stddev << std::endl
     << "# Projection percentil  2.5%: " << hist_proj.percentil(2.5) << std::endl
     << "# Projection percentil 97.5%: " << hist_proj.percentil(97.5) << std::endl;
     proj_area.computeStats(area_avg, area_stddev, dummy, dummy);
     std::cout << "# Projection area average: " << area_avg << std::endl
     << "# Projection area stddev:  " << area_stddev << std::endl
     << "# Area percentil  2.5%:    " << hist_area.percentil(2.5) << std::endl
     << "# Area percentil 97.5%:    " << hist_area.percentil(97.5) << std::endl;
 
     sigma=sqrt(power_avg*power_avg*Xdim*Ydim / (targetSNR*area_avg));
     std::cout << "Estimated sigma=" << sigma << std::endl;
     updateCtfs();
 }

◆ preProcess()

void ProgSimulateMicroscope::preProcess ( )

virtual

Produce side information.

Reimplemented from XmippMetadataProgram.

Definition at line 250 of file phantom_simulate_microscope.cpp.

 {
     randomize_random_generator();
     size_t dum, dum2;
     getImageSize(*pmdIn, Xdim, Ydim, dum, dum2);
 
     if (low_pass_before_CTF < 0.5)
     {
         lowpass.FilterBand = LOWPASS;
         lowpass.FilterShape = RAISED_COSINE;
         lowpass.w1 = low_pass_before_CTF;
     }
 
     if (estimateSNR)
         estimateSigma();
 }

◆ processImage()

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

virtual

Implements XmippMetadataProgram.

Definition at line 267 of file phantom_simulate_microscope.cpp.

 {
     static Image<double> img;
     static FileName last_ctf;
     static bool firstImage=true;
     last_ctf = fn_ctf;
     img.read(fnImg);
 
     rowOut=rowIn;
 
     rowIn.getValue(MDL_CTF_MODEL, fn_ctf);
     rowOut.setValue(MDL_IMAGE, fnImgOut);
     rowOut.setValue(MDL_CTF_MODEL, fn_ctf);
     if (fn_ctf != last_ctf || firstImage)
     {
         updateCtfs();
         firstImage=false;
     }
 
     if (ZSIZE(img())!=1)
         REPORT_ERROR(ERR_MULTIDIM_DIM,"This process is not intended for volumes");
 
     apply(img());
 
     img.write(fnImgOut);
 }

◆ readParams()

void ProgSimulateMicroscope::readParams ( )

virtual

Read from a command line. An exception might be thrown by any of the internal conversions, this would mean that there is an error in the command line and you might show a usage message.

Reimplemented from XmippMetadataProgram.

Definition at line 33 of file phantom_simulate_microscope.cpp.

 {
     XmippMetadataProgram::readParams();
 
     fn_ctf = "";//initialize empty, force recalculation of first time
     pmdIn = dynamic_cast<MetaDataVec*>(getInputMd());
     CTFpresent=true;
     if (checkParam("--ctf"))
     {
         //Fill the input metadata with the value of 'fn_ctf'
         fn_ctf=getParam("--ctf");
         MDConstGenerator generator(fn_ctf);
         generator.label = MDL_CTF_MODEL;
         generator.fill(*pmdIn);
     }
     else
     {
         if (pmdIn->containsLabel(MDL_CTF_MODEL))
         {
             //sort the images according to the ctf to avoid the recaculation of it
             //beeten images of the same ctf group
             MetaDataDb md(*pmdIn);
             pmdIn->sort(md, MDL_CTF_MODEL);
         }
         else
             CTFpresent=false;
     }
     after_ctf_noise = checkParam("--after_ctf_noise");
     defocus_change = getDoubleParam("--defocus_change");
     if (checkParam("--noise"))
     {
         sigma = getDoubleParam("--noise");
         low_pass_before_CTF = getDoubleParam("--noise",1);
         estimateSNR=false;
     }
     else if (checkParam("--noNoise"))
     {
         estimateSNR=0;
         sigma=0;
         low_pass_before_CTF=0.5;
     }
     else
     {
         targetSNR = getDoubleParam("--targetSNR");
         low_pass_before_CTF=0.5;
         estimateSNR=true;
     }
     downsampling = getDoubleParam("--downsampling");
 
 }

◆ setupFourierFilter()

void ProgSimulateMicroscope::setupFourierFilter	(	FourierFilter &	filter,
		bool	isBackground,
		double &	power
	)

Prepare fourier filter

Definition at line 200 of file phantom_simulate_microscope.cpp.

 {
     static int dXdim = 2 * Xdim, dYdim = 2 * Ydim;
     static MultidimArray<double> aux;
 
     filter.FilterBand = CTF;
     filter.ctf.read(fn_ctf);
     defocusU = ctf.ctf.DeltafU;
     defocusV = ctf.ctf.DeltafV;
     filter.ctf.enable_CTF = !isBackground;
     filter.ctf.enable_CTFnoise = isBackground;
     filter.ctf.produceSideInfo();
     aux.resizeNoCopy(dYdim, dXdim);
     aux.setXmippOrigin();
     filter.do_generate_3dmask=true;
     filter.generateMask(aux);
     power = filter.maskPower();
 }

◆ show()

void ProgSimulateMicroscope::show ( )

Show parameters.

Definition at line 110 of file phantom_simulate_microscope.cpp.

 {
     if (!verbose)
         return;
     XmippMetadataProgram::show();
     std::cout
     << "Noise: " << sigma << std::endl
     << "Low pass freq: " << low_pass_before_CTF << std::endl
     << "After CTF noise: " << after_ctf_noise << std::endl
     << "Defocus change: " << defocus_change << std::endl
     ;
     if (estimateSNR)
         std::cout
         << "Target SNR: " << targetSNR << std::endl;
     if (CTFpresent)
         std::cout << "CTF file: " << fn_ctf << std::endl;
     else
         std::cout << "No CTF provided\n";
 }

◆ updateCtfs()

void ProgSimulateMicroscope::updateCtfs ( )

Update ctfs

Definition at line 220 of file phantom_simulate_microscope.cpp.

 {
     double before_power = 0, after_power = 0;
     sigma_after_CTF = sigma_before_CTF = 0;
 
     if (CTFpresent)
     {
         setupFourierFilter(ctf, false, before_power);
         if (after_ctf_noise)
             setupFourierFilter(after_ctf, true, after_power);
     }
 
     // Compute noise balance
     if ((after_power != 0 || before_power != 0) && sigma!=0)
     {
         double p = after_power / (after_power + before_power);
         double K = 1 / sqrt(p * after_power + (1 - p) * before_power);
         sigma_after_CTF = sqrt(p) * K * sigma;
         sigma_before_CTF = sqrt(1 - p) * K * sigma;
     }
     else if (sigma != 0)
     {
         sigma_before_CTF = sigma;
         sigma_after_CTF = 0;
     }
     else
         sigma_before_CTF=sigma_after_CTF=0;
 }