Estimate the scatter probability using a model-based approach. More...

#include "stir/scatter/ScatterEstimation.h"

Inheritance diagram for stir::ScatterEstimation:

Public Member Functions
	ScatterEstimation ()
	Default constructor (calls set_defaults())

	ScatterEstimation (const std::string &parameter_filename)
	Overloaded constructor with parameter file and initialisation.

virtual Succeeded	process_data ()
	Full process_data which performs set_up() before beginning.

shared_ptr< ProjData >	get_output () const
	Get current scatter estimate.

shared_ptr< ProjData >	make_2D_projdata_sptr (const shared_ptr< ProjData > in_3d_sptr)
	make projdata 2D shared pointer

shared_ptr< ProjData >	make_2D_projdata_sptr (const shared_ptr< ProjData > in_3d_sptr, string template_filename)

virtual Succeeded	set_up ()
	set_up More...

void	set_input_proj_data_sptr (const shared_ptr< ProjData >)
	Set the input projdata.

void	set_input_data (const shared_ptr< ProjData > &data)
	Set the input projdata. More...

shared_ptr< const ProjData >	get_input_data () const

void	set_reconstruction_method_sptr (const shared_ptr< Reconstruction< DiscretisedDensity< 3, float >>>)
	Set the reconstruction method for the scatter estimation.

void	set_attenuation_image_sptr (const shared_ptr< const DiscretisedDensity< 3, float >>)
	Set the full resolution attenuation image.

void	set_attenuation_correction_proj_data_sptr (const shared_ptr< ProjData >)
	set projection data that contains the attenuation correction factors

void	set_normalisation_sptr (const shared_ptr< BinNormalisation >)
	set normalisation object (excluding attenuation)

void	set_background_proj_data_sptr (const shared_ptr< ProjData >)

void	set_initial_activity_image_sptr (const shared_ptr< const DiscretisedDensity< 3, float >>)

void	set_mask_image_sptr (const shared_ptr< const DiscretisedDensity< 3, float >>)

void	set_mask_proj_data_sptr (const shared_ptr< ProjData > arg)
	set mask for tail-fitting More...

void	set_recompute_mask_image (bool arg)

void	set_recompute_mask_projdata (bool arg)

void	set_scatter_simulation_method_sptr (const shared_ptr< ScatterSimulation >)

void	set_num_iterations (int)

void	set_output_scatter_estimate_prefix (const std::string &)

void	set_export_scatter_estimates_of_each_iteration (bool)

void	set_max_scale_value (float value)

void	set_min_scale_value (float value)

void	set_mask_projdata_filename (std::string name)

void	set_mask_image_filename (std::string name)

void	set_output_additive_estimate_prefix (std::string name)

void	set_run_debug_mode (bool debug)

void	set_restart_reconstruction_every_scatter_iteration (bool setting)

bool	get_restart_reconstruction_every_scatter_iteration () const

int	get_iterations_num () const
	Set the zoom factor in the XY plane for the downsampling of the activity and attenuation image. More...

int	get_num_iterations () const
	Get the number of iterations for the scatter estimation.

shared_ptr< const DiscretisedDensity< 3, float > >	get_estimated_activity_image_sptr () const
	Get the (low resolution) estimate of the activity image.

virtual bool	already_setup () const
	allows checking if we have called set_up()

Public Member Functions inherited from stir::ParsingObject
	ParsingObject (const ParsingObject &)

ParsingObject &	operator= (const ParsingObject &)

void	ask_parameters ()

virtual std::string	parameter_info ()

bool	parse (std::istream &f)

bool	parse (const char *const filename)

Static Public Member Functions
static void	upsample_and_fit_scatter_estimate (ProjData &scaled_scatter_proj_data, const ProjData &emission_proj_data, const ProjData &scatter_proj_data, BinNormalisation &scatter_normalisation, const ProjData &weights_proj_data, const float min_scale_factor, const float max_scale_factor, const unsigned half_filter_width, BSpline::BSplineType spline_type=BSpline::BSplineType::linear, const bool remove_interleaving=true)
	upsample coarse scatter estimate and fit it to tails of the emission data More...

Protected Member Functions
void	set_defaults () override
	All recomputes_** will default true.

void	initialise_keymap () override
	Initialise all keywords.

bool	post_processing () override
	This will be called at the end of the parsing. More...

Protected Member Functions inherited from stir::ParsingObject
virtual void	set_key_values ()
	This will be called before parsing or parameter_info is called. More...

Protected Attributes
bool	recompute_mask_image
	Recompute or load the mask image.

bool	recompute_mask_projdata
	If set the mask projdata will be recomputed.

bool	recompute_atten_projdata
	If set to 1 the attenuation coefficients are going to be recalculated.

bool	restart_reconstruction_every_scatter_iteration
	If set to true, the activity image will be reset to 1 in each iteration of the scatter estimation. Therefore, more reconstruction subiterations will be required for convergence.

shared_ptr< Reconstruction< DiscretisedDensity< 3, float > > >	reconstruction_template_sptr
	This is the reconstruction object which is going to be used for the scatter estimation and the calculation of the initial activity image (if recompute set). It can be defined in the same parameters file as the scatter parameters or in an external defined in the reconstruction_template_par_filename.

shared_ptr< DiscretisedDensity< 3, float > >	current_activity_image_sptr
	The current activity estimate.

shared_ptr< const DiscretisedDensity< 3, float > >	atten_image_sptr
	Image with attenuation values.

shared_ptr< BinNormalisation >	norm_3d_sptr
	normalisation components in 3D (without atten)

shared_ptr< ProjData >	mask_projdata_sptr
	Mask proj_data.

shared_ptr< ProjData >	input_projdata_sptr
	The full 3D projdata are used for the calculation of the 2D and later for the upsampling back to 3D.

shared_ptr< ProjData >	input_projdata_2d_sptr
	The 2D projdata are used for the scatter estimation.

shared_ptr< ProjData >	add_projdata_2d_sptr
	Additive projection data after SSRB – Randoms.

shared_ptr< ProjData >	data_to_fit_projdata_sptr
	Prompts - randoms.

shared_ptr< ProjData >	add_projdata_sptr

shared_ptr< ProjData >	back_projdata_2d_sptr
	(Additive + Scatter Estimate) * Mult in 2D

shared_ptr< ProjData >	back_projdata_sptr
	Initially this points to the un-normalised randoms.

std::string	mask_image_filename
	Filename of mask image.

std::string	mask_projdata_filename
	Filename of mask's projdata.

std::string	back_projdata_filename
	Filename of background projdata.

std::string	tail_mask_par_filename
	Optional parameter file for the tail fitting. More...

std::string	input_projdata_filename
	Filename of the measured emission 3D data.

std::string	atten_image_filename
	This is the image file name with the anatomic information.

std::string	recon_template_par_filename
	The filename for the parameters file of the reconstruction method.

std::string	atten_coeff_filename
	The file name for the attenuation coefficients. If they are to be recalculated they will be stored here, if set.

MaskingParameters	masking_parameters

int	num_scatter_iterations

std::string	output_scatter_estimate_prefix
	Output file name prefix.

std::string	output_additive_estimate_prefix

bool	_already_setup
	variable to check if we have called set_up()

Protected Attributes inherited from stir::ParsingObject
KeyParser	parser

Detailed Description

Estimate the scatter probability using a model-based approach.

Author: Nikos Efthimiou; Kris Thielemans

This is an implementation of Watson's SSS iterative algorithm to estimate scatter. Simulation is normally at low resolution, ideally reconstruction as well (to save time). The code also does tail-fitting.

A particular feature is that the first iterations happen in 2D (after SSRB of all the data). Currently only the last simulation is in 3D. This could be changed later.

Todo:

The code should throw an error if 2D input data are loaded. It should just deactivate the final upsampling.

Currently FBP reconstruction is not working and just throws an error.

This code needs far more documentation.

Deprecated:: The current parsing code still allows the bin normalisation type keyword to specify a "chained" normalisation object with first the norm and second the atten. As there is no way to guarantee that this will be the case, we have now replaced this with explicit setting of the norm factors only (keyword normalisation type) and the attenuation factors. The naming of the keywords is confusing however. (Ensuring backwards compatibility was not so easy, so the code might look confusing.)

Member Function Documentation

◆ upsample_and_fit_scatter_estimate()

void stir::ScatterEstimation::upsample_and_fit_scatter_estimate	(	ProjData &	scaled_scatter_proj_data,
		const ProjData &	emission_proj_data,
		const ProjData &	scatter_proj_data,
		BinNormalisation &	scatter_normalisation,
		const ProjData &	weights_proj_data,
		const float	min_scale_factor,
		const float	max_scale_factor,
		const unsigned	half_filter_width,
		BSpline::BSplineType	spline_type = `BSpline::BSplineType::linear`,
		const bool	remove_interleaving = `true`
	)

static

upsample coarse scatter estimate and fit it to tails of the emission data

Current procedure:

interpolate segment 0 of scatter_proj_data to size of segment 0 of emission_proj_data
inverseSSRB to create oblique segments
undo normalisation (as measured data is not normalised)
find scale factors with get_scale_factors_per_sinogram()
apply thresholds
filter scale-factors in axial direction (independently for every segment)
apply scale factors using scale_sinograms()

References stir::VectorWithOffset< T >::begin(), stir::Array< num_dimensions, elemT >::begin_all(), stir::VectorWithOffset< T >::end(), stir::Array< num_dimensions, elemT >::end_all(), stir::error(), stir::Array< num_dimensions, elemT >::fill(), stir::ProjData::fill(), stir::VectorWithOffset< T >::fill(), stir::ExamData::get_exam_info_sptr(), stir::ProjDataInfo::get_max_axial_pos_num(), stir::ProjDataInfo::get_max_segment_num(), stir::ProjDataInfo::get_min_axial_pos_num(), stir::ProjDataInfo::get_min_segment_num(), stir::ProjData::get_proj_data_info_sptr(), stir::get_scale_factors_per_sinogram(), stir::Array< num_dimensions, elemT >::grow(), stir::info(), stir::interpolate_projdata(), stir::inverse_SSRB(), stir::BinNormalisation::is_trivial(), stir::VectorWithOffset< T >::resize(), stir::scale_sinograms(), stir::BinNormalisation::set_up(), stir::threshold_lower(), stir::threshold_upper(), stir::BinNormalisation::undo(), and stir::warning().

Referenced by process_data().

◆ set_up()

Succeeded stir::ScatterEstimation::set_up ( )

virtual

set_up

Returns

This function will take care most of the initialisation needed:

Procedure:
1. Check if debug mode and activate all export flags.
2. Load the input_projdata_3d_sptr and perform SSRB
3. Initialise (partially for the moment) the reconstruction method:
4. Load Normalisation data and perform SSRB
5. Load the background data (randoms) and do normalisation (to get the additive data)

References _already_setup, add_projdata_2d_sptr, stir::FilePath::append(), atten_image_sptr, back_projdata_2d_sptr, back_projdata_sptr, current_activity_image_sptr, data_to_fit_projdata_sptr, stir::error(), stir::FilePath::get_current_working_directory(), stir::info(), input_projdata_2d_sptr, input_projdata_sptr, make_2D_projdata_sptr(), mask_image_filename, mask_projdata_sptr, masking_parameters, distributed::min_threshold, stir::ProjData::read_from_file(), recompute_mask_image, recompute_mask_projdata, reconstruction_template_sptr, and stir::warning().

◆ set_input_data()

void stir::ScatterEstimation::set_input_data ( const shared_ptr< ProjData > & data )

Set the input projdata.

Using same name as Reconstruction

References stir::error(), input_projdata_sptr, and set_input_proj_data_sptr().

◆ set_mask_proj_data_sptr()

void stir::ScatterEstimation::set_mask_proj_data_sptr ( const shared_ptr< ProjData > arg )

inline

set mask for tail-fitting

arg will not be modified

◆ get_iterations_num()

int stir::ScatterEstimation::get_iterations_num ( ) const

Set the zoom factor in the XY plane for the downsampling of the activity and attenuation image.

Set the zoom factor in the Z axis for the downsampling of the activity and attenuation image. Get the number of iterations for the scatter estimation

Deprecated:: Use get_num_iterations()

References stir::error(), stir::FilePath::get_filename(), stir::FilePath::get_path(), input_projdata_sptr, norm_3d_sptr, num_scatter_iterations, stir::FilePath::prepend_directory_name(), and stir::warning().

◆ post_processing()

bool stir::ScatterEstimation::post_processing ( )

overrideprotectedvirtual

This will be called at the end of the parsing.

Returns: false if everything OK, true if not

Reimplemented from stir::ParsingObject.

References stir::KeyParser::add_parsing_key(), stir::KeyParser::add_start_key(), stir::KeyParser::add_stop_key(), atten_coeff_filename, atten_image_filename, atten_image_sptr, back_projdata_filename, back_projdata_sptr, stir::error(), stir::MaskingParameters::filter_filename, stir::MaskingParameters::filter_sptr, stir::info(), input_projdata_filename, input_projdata_sptr, mask_image_filename, mask_projdata_filename, mask_projdata_sptr, masking_parameters, output_scatter_estimate_prefix, stir::KeyParser::parse(), stir::read_from_file(), stir::ProjData::read_from_file(), recompute_mask_image, recompute_mask_projdata, recon_template_par_filename, reconstruction_template_sptr, set_attenuation_correction_proj_data_sptr(), and stir::warning().

Member Data Documentation

◆ tail_mask_par_filename

std::string stir::ScatterEstimation::tail_mask_par_filename

protected

Optional parameter file for the tail fitting.

If not provided, sensible defaults are used

Referenced by initialise_keymap(), and process_data().

◆ masking_parameters

MaskingParameters stir::ScatterEstimation::masking_parameters

protected

the set of parameters to obtain a mask from the attenuation image

The attenuation image will be thresholded to find a plausible mask for where there can be emission data. This mask will be then forward projected to find the tails in the projection data.

This is a simple strategy that can fail due to motion etc, so the attenuation image is first blurred, and the default threshold is low.

Note that there is currently no attempt to eliminate the bed from the attenuation image first. Tails are therefore going to be too small, which could create trouble.

By default, a Gaussian filter of FWHM (15,20,20) will be applied before thresholding with a value 0.003 cm^-1

Referenced by initialise_keymap(), post_processing(), process_data(), set_defaults(), and set_up().

◆ num_scatter_iterations

int stir::ScatterEstimation::num_scatter_iterations

protected

The number of iterations the scatter estimation will perform. Default = 5.

Referenced by get_iterations_num(), get_num_iterations(), initialise_keymap(), process_data(), and set_defaults().

The documentation for this class was generated from the following files:

/home/sirfuser/devel/STIRdistrib/STIR/src/include/stir/scatter/ScatterEstimation.h
/home/sirfuser/devel/STIRdistrib/STIR/src/include/stir/scatter/ScatterEstimation.inl
/home/sirfuser/devel/STIRdistrib/STIR/src/scatter_buildblock/ScatterEstimation.cxx
/home/sirfuser/devel/STIRdistrib/STIR/src/scatter_buildblock/upsample_and_fit_scatter_estimate.cxx

Public Member Functions

Static Public Member Functions

Protected Member Functions

Protected Attributes

Detailed Description

Member Function Documentation

◆ upsample_and_fit_scatter_estimate()

◆ set_up()

◆ set_input_data()

◆ set_mask_proj_data_sptr()

◆ get_iterations_num()

◆ post_processing()

Member Data Documentation

◆ tail_mask_par_filename

◆ masking_parameters

◆ num_scatter_iterations