#include <tomo_tiltseries_dose_filter.h>

Inheritance diagram for ProgTomoTSFilterDose:

Collaboration diagram for ProgTomoTSFilterDose:

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

void	show ()
	Show. More...

void	defineParams ()
	Define parameters. More...

void	run ()
	Run. More...

void	initParameters ()

void	readInputData (MetaDataVec &md)

double	doseFilter (double dose_at_end_of_frame, double critical_dose)
	Compute the dose filter, which is the signal attenuation. More...

double	criticalDose (double spatial_frequency)
	Given a spatial frequency, return the critical dose in electrons per square Angstroms. More...

double	optimalDoseGivenCriticalDose (double critical_dose)
	Given the critical dose, return an estimate of the optimal dose (at which the SNR is maximised) More...

void	applyDoseFilterToImage (int Ydim, int Xdim, const MultidimArray< std::complex< double > > &FFT1, const double dose_finish)
	Apply a dose filter to the image Fourier transform. 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
double	critical_dose_a

double	critical_dose_b

double	critical_dose_c

double	critical_dose_at_dc

double	voltage_scaling_factor

double	accVoltage

double	dose_per_image

double	pre_exposure_amount

FileName	fnTs

FileName	fnOut

int	user_supplied_first_frame

int	user_supplied_last_frame

double	sampling

double	maxFreq

bool	restore_power
	Restore noise power after filtering?', 'Renormalise the summed image after filtering. More...

std::vector< MultidimArray< std::complex< double > > *>	frameFourierVec

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 73 of file tomo_tiltseries_dose_filter.h.

Member Function Documentation

◆ applyDoseFilterToImage()

void ProgTomoTSFilterDose::applyDoseFilterToImage	(	int	Ydim,
		int	Xdim,
		const MultidimArray< std::complex< double > > &	FFT1,
		const double	dose_finish
	)

Apply a dose filter to the image Fourier transform.

Definition at line 110 of file tomo_tiltseries_dose_filter.cpp.

 {
     double ux, uy;
     double yy;
     double current_critical_dose, current_optimal_dose;
     int sizeX_2 = Xdim / 2;
     double ixsize = 1.0 / Xdim;
     int sizeY_2 = Ydim / 2;
     double iysize = 1.0 / Ydim;
 
     for (long int i = 0; i < YSIZE(FFT1); i++) //i->y,j->x xmipp convention is oposite summove
     { 
         FFT_IDX2DIGFREQ_FAST(i, Ydim, sizeY_2, iysize, uy);
         
         yy = uy * uy;
         for (long int j = 0; j < XSIZE(FFT1); j++)
         {
             FFT_IDX2DIGFREQ_FAST(j, Xdim, sizeX_2, ixsize, ux);
 
             if (i == 0 && j == 0)
             {
                 current_critical_dose = critical_dose_at_dc;
             }
             else
             {
                 current_critical_dose = criticalDose(sqrt(ux * ux + yy) / sampling);
             }
 
             current_optimal_dose = 2.51284 * current_critical_dose; //This is the optimal dose given the critical dose (see Grant, Grigorief eLife 2015)
 
             DIRECT_A2D_ELEM(FFT1, i, j ) *= doseFilter(dose_finish, current_critical_dose);
                 
 
         }
     }
 }

◆ criticalDose()

double ProgTomoTSFilterDose::criticalDose ( double spatial_frequency )

Given a spatial frequency, return the critical dose in electrons per square Angstroms.

Definition at line 73 of file tomo_tiltseries_dose_filter.cpp.

 {
     return ((critical_dose_a * pow(spatial_frequency, critical_dose_b)) +  critical_dose_c) * voltage_scaling_factor;
 }

◆ defineParams()

void ProgTomoTSFilterDose::defineParams ( )

virtual

Define parameters.

Reimplemented from XmippProgram.

Definition at line 55 of file tomo_tiltseries_dose_filter.cpp.

                                         {
     addUsageLine("AThis algorithm applies a dose filter based on Grant & Grigorief eLife 2015. The result will be a low pass filtered tilt series based on the applied dose.");
     addParamsLine("   -i <tiltseries>                  : Metadata with the input tilt series");
     addParamsLine("  [-o <fn=\"out.mrcs\">]            : Output path for the filtered tilt series.");
     addParamsLine("  [--sampling <Ts=1>]               : Sampling rate (A/pixel)");
     addParamsLine("  [--voltage <voltage=300>]         : Acceleration voltage (kV) min_value=200.0e0,max_value=300.0e0)");
     addExampleLine("Example", false);
     addExampleLine("xmipp_tomo_tiltseries_dose_filter -i tiltseries.xmd -o . --sampling 2.2 --dosePerImage 2 --voltage 300");
     addExampleLine("xmipp_tomo_tiltseries_dose_filter -i tiltseries.xmd -o . --sampling 2.2 --voltage 300");
 }

◆ doseFilter()

double ProgTomoTSFilterDose::doseFilter	(	double	dose_at_end_of_frame,
		double	critical_dose
	)

Compute the dose filter, which is the signal attenuation.

Definition at line 67 of file tomo_tiltseries_dose_filter.cpp.

 {
     return exp((-0.5 * dose_at_end_of_frame) / critical_dose);
 }

◆ initParameters()

void ProgTomoTSFilterDose::initParameters ( )

Set up the critical curve function (summovie)

Definition at line 78 of file tomo_tiltseries_dose_filter.cpp.

 {
     critical_dose_a = 0.24499;
     critical_dose_b = -1.6649;
     critical_dose_c = 2.8141;
     // init with a very large number
     critical_dose_at_dc = std::numeric_limits<double>::max() * 0.001;
 
     // Defining the scaling factor given the acceleration voltage. The methods considers only two voltage 200keV and 300keV
     // Any other voltage would result in an error
     if (accVoltage < 301 && accVoltage > 299.)
     {
             accVoltage = 300.0;
             voltage_scaling_factor = 1.0;
     }
     else
     {
         if (accVoltage < 201.0 && accVoltage > 199.0)
         {
             accVoltage = 200.0;
             voltage_scaling_factor = 0.8;
         }
         else
         {
             REPORT_ERROR(ERR_ARG_INCORRECT, "Bad acceleration voltage (must be 200 or 300 kV");
         }
             
     }
 
 }

◆ optimalDoseGivenCriticalDose()

double ProgTomoTSFilterDose::optimalDoseGivenCriticalDose ( double critical_dose )

Given the critical dose, return an estimate of the optimal dose (at which the SNR is maximised)

◆ readInputData()

void ProgTomoTSFilterDose::readInputData ( MetaDataVec & md )

Definition at line 148 of file tomo_tiltseries_dose_filter.cpp.

 {
     size_t Xdim, Ydim, Zdim, Ndim;
 
     //if input is an stack create a metadata.
     if (fnTs.isMetaData())
         mdts.read(fnTs);
     else
     {
         REPORT_ERROR(ERR_ARG_INCORRECT, "The input must be a metadata with the image filenames and the accumulated dose");
     }
 }

◆ readParams()

void ProgTomoTSFilterDose::readParams ( )

virtual

Read argument from command line.

Reimplemented from XmippProgram.

Definition at line 35 of file tomo_tiltseries_dose_filter.cpp.

                                       {
     fnTs = getParam("-i");
     fnOut = getParam("-o");
     sampling = getDoubleParam("--sampling");
     accVoltage = getDoubleParam("--voltage");
 }

◆ run()

void ProgTomoTSFilterDose::run ( )

virtual

Run.

Reimplemented from XmippProgram.

Definition at line 162 of file tomo_tiltseries_dose_filter.cpp.

 {
     
     show();
 
     initParameters();
     
     MetaDataVec mdts;
 
     readInputData(mdts);
 
     std::cout << " Starting filtering..."<< std::endl;
 
     
     FourierTransformer transformer;
     FileName fnti;
     FileName fnOutFrame;
     
     Image<double> tiltImage;
     MultidimArray<double> &ptrtiltImage = tiltImage();
 
     int mode = WRITE_OVERWRITE;
     double doseTi;
 
     FileName fnk;
     std::vector<FileName> ti_fn;
     std::vector<double> ti_dose;
 
     for (size_t objId : mdts.ids())
     {
         mdts.getValue(MDL_IMAGE, fnti, objId);
 
         mdts.getValue(MDL_DOSE, doseTi, objId);
 
         ti_fn.push_back(fnti);
         ti_dose.push_back(doseTi);
     }
 
     FileName fnTsOut;
 
     size_t n = 0;
     MultidimArray<std::complex<double> > FFT1;
     
     for (size_t objId : mdts.ids())
     {
         mdts.getValue(MDL_IMAGE, fnti, objId);
         mdts.getValue(MDL_DOSE, doseTi, objId);
         
         tiltImage.read(fnti);
 
         // Now do the Fourier transform and filter
         transformer.FourierTransform(ptrtiltImage, FFT1, false);
 
         applyDoseFilterToImage((int) YSIZE(ptrtiltImage), (int) XSIZE(ptrtiltImage), FFT1, doseTi);
 
         transformer.inverseFourierTransform();
 
         fnTsOut = fnOut;
         tiltImage.write(fnTsOut, n+FIRST_IMAGE, true, mode);
         mode = WRITE_APPEND;
 
         ++n;
 
     }
 }

◆ show()

void ProgTomoTSFilterDose::show ( )

Show.

Definition at line 43 of file tomo_tiltseries_dose_filter.cpp.

 {
     if (!verbose)
         return;
     std::cout << "Input movie:           " << fnTs
             << std::endl << "Output movie:          "
             << fnOut << std::endl
             << "Sampling:              " << sampling << std::endl
             << "Acceleration_Voltage (kV):   " << accVoltage << std::endl;
 }