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

#include "stir/scatter/ScatterSimulation.h"

Inheritance diagram for stir::ScatterSimulation:

[legend]

Public Member Functions
	ScatterSimulation ()
	Default constructor.

virtual Succeeded	process_data ()

virtual std::string	method_info () const =0
	gives method information

virtual void	ask_parameters ()
	prompts the user to enter parameter values manually

void	downsample_density_image_for_scatter_points (float _zoom_xy, float _zoom_z, int _size_xy=-1, int _size_z=-1)
	This function is a less powerfull tool than directly zooming the image. However it will check that the downsampling is done in manner compatible with the ScatterSimulation.

Succeeded	downsample_scanner (int new_num_rings=-1, int new_num_dets=-1)
	Downsample the scanner keeping the total axial length the same. More...

Succeeded	downsample_images_to_scanner_size ()
	Downsamples activity and attenuation images to voxel sizes appropriate for the (downsampled) scanner. More...

float	detection_efficiency (const float energy) const
	gamma-energy-part of the detection efficiency More...

virtual Succeeded	set_up ()

virtual void	write_log (const double simulation_time, const float total_scatter)
	Output the log of the process.

void	set_use_cache (const bool)
	Enable/disable caching of line integrals.

bool	get_use_cache () const
	Return if line integrals are cached or not.

check functions
bool	has_template_proj_data_info () const

bool	has_exam_info () const
	Returns true if template_exam_info_sptr has been set.

get functions
shared_ptr< ProjData >	get_output_proj_data_sptr () const

int	get_num_scatter_points () const

shared_ptr< const ProjDataInfo >	get_template_proj_data_info_sptr () const
	Get the template ProjDataInfo.

shared_ptr< const ExamInfo >	get_exam_info_sptr () const
	Get the ExamInfo.

const DiscretisedDensity< 3, float > &	get_activity_image () const

const DiscretisedDensity< 3, float > &	get_attenuation_image () const

const DiscretisedDensity< 3, float > &	get_attenuation_image_for_scatter_points () const

shared_ptr< const DiscretisedDensity< 3, float > >	get_density_image_for_scatter_points_sptr () const

set functions
void	set_template_proj_data_info (const std::string &)

void	set_template_proj_data_info (const ProjDataInfo &)

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

void	set_activity_image (const std::string &filename)

void	set_exam_info (const ExamInfo &)

void	set_exam_info_sptr (const shared_ptr< const ExamInfo >)

void	set_output_proj_data_sptr (shared_ptr< ProjData >)

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

void	set_density_image (const std::string &)

void	set_output_proj_data (const std::string &)
	This function depends on the ProjDataInfo of the scanner. You first have to set that.

void	set_output_proj_data_sptr (const shared_ptr< const ExamInfo >, const shared_ptr< const ProjDataInfo >, const std::string &)

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

void	set_image_downsample_factors (float factor_xy=1.f, float factor_z=1.f, int _size_zoom_xy=-1, int _size_zoom_z=-1)

void	set_density_image_for_scatter_points (const std::string &)
	set_density_image_for_scatter_points

void	set_attenuation_threshold (const float)
	set the attenuation threshold

void	set_randomly_place_scatter_points (const bool)
	The scattering point in the voxel will be chosen randomly, instead of choosing the centre. More...

void	set_cache_enabled (const bool)


void	set_downsample_scanner_bool (const bool arg)
	Get and set methods for the downsample_scanner_bool.

bool	get_downsample_scanner_bool () const


int	get_num_downsample_scanner_rings () const
	Get and set methods for downsample_scanner_rings.

void	set_num_downsample_scanner_rings (const int arg)


int	get_num_downsample_scanner_dets () const
	Get and set methods for downsample_scanner_dets.

void	set_num_downsample_scanner_dets (const int arg)

Public Member Functions inherited from stir::RegisteredObjectBase
virtual std::string	get_registered_name () const =0
	Returns the name of the type of the object. More...

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
Compton scatter cross sections
static float	dif_Compton_cross_section (const float cos_theta, float energy)

static float	total_Compton_cross_section (float energy)

static float	photon_energy_after_Compton_scatter (const float cos_theta, const float energy)

static float	photon_energy_after_Compton_scatter_511keV (const float cos_theta)

static float	total_Compton_cross_section_relative_to_511keV (const float energy)

Static Public Member Functions inherited from stir::RegisteredObject< ScatterSimulation >
static ScatterSimulation *	read_registered_object (std::istream *in, const std::string &registered_name)
	Construct a new object (of a type derived from Root, its actual type determined by the registered_name parameter) by parsing the istream. More...

static ScatterSimulation *	ask_type_and_parameters ()
	ask the user for the type, and then calls read_registered_object(0, type) More...

static void	list_registered_names (std::ostream &stream)
	List all possible registered names to the stream. More...

Protected Types
enum	image_type { act_image_type, att_image_type }

Protected Types inherited from stir::RegisteredObject< ScatterSimulation >
typedef ScatterSimulation (	RootFactory) (std::istream *)
	The type of a root factory is a function, taking an istream* as argument, and returning a Root*.

typedef FactoryRegistry< std::string, RootFactory, interfile_less >	RegistryType
	The type of the registry.

Protected Member Functions
virtual double	process_data_for_view_segment_num (const ViewSegmentNumbers &vs_num)
	computes scatter for one viewgram More...

float	compute_emis_to_det_points_solid_angle_factor (const CartesianCoordinate3D< float > &emis_point, const CartesianCoordinate3D< float > &detector_coord)

void	set_defaults () override
	Set defaults before parsing.

void	initialise_keymap () override
	Initialise all keywords.

bool	post_processing () override

void	sample_scatter_points ()
	find scatter points More...

virtual void	remove_cache_for_integrals_over_attenuation ()
	remove cached attenuation integrals More...

virtual void	remove_cache_for_integrals_over_activity ()
	reset cached activity integrals More...

virtual double	scatter_estimate (const Bin &bin)=0
	virtual function that computes the scatter for one (downsampled) bin

void	initialise_cache_for_scattpoint_det_integrals_over_attenuation ()
	set-up cache for attenuation integrals More...

void	initialise_cache_for_scattpoint_det_integrals_over_activity ()
	set-up cache for activity integrals 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
std::vector< ScatterPoint >	scatt_points_vector

float	scatter_volume

std::string	template_proj_data_filename

shared_ptr< ProjDataInfo >	proj_data_info_sptr

shared_ptr< ExamInfo >	template_exam_info_sptr

std::string	density_image_filename

std::string	density_image_for_scatter_points_filename

std::string	density_image_for_scatter_points_output_filename

shared_ptr< const DiscretisedDensity< 3, float > >	density_image_sptr

shared_ptr< const DiscretisedDensity< 3, float > >	activity_image_sptr
	Pointer to hold the current activity estimation.

std::string	output_proj_data_filename
	Output proj_data fileanme prefix.

shared_ptr< ProjData >	output_proj_data_sptr
	Shared ptr to hold the simulated data.

float	attenuation_threshold
	threshold below which a voxel in the attenuation image will not be considered as a candidate scatter point

bool	randomly_place_scatter_points
	boolean to see if we need to move the scatter point randomly within in its voxel

bool	use_cache
	boolean to see if we need to cache the integrals More...

std::string	activity_image_filename
	Filename for the initial activity estimate.

float	zoom_xy
	Zoom factor on plane XY. Defaults on 1.f.

float	zoom_z
	Zoom factor on Z axis. Defaults on 1.f.

int	zoom_size_xy
	Zoomed image size on plane XY. Defaults on -1.

int	zoom_size_z
	Zoomed image size on Z axis. Defaults on -1.

int	downsample_scanner_rings
	Number of rings of downsampled scanner.

int	downsample_scanner_dets
	Number of detectors per ring of downsampled scanner.

bool	downsample_scanner_bool

bool	_already_set_up

Protected Attributes inherited from stir::ParsingObject
KeyParser	parser

detection related functions
CartesianCoordinate3D< float >	shift_detector_coordinates_to_origin

std::vector< CartesianCoordinate3D< float > >	detection_points_vector

virtual void	find_detectors (unsigned &det_num_A, unsigned &det_num_B, const Bin &bin) const

unsigned	find_in_detection_points_vector (const CartesianCoordinate3D< float > &coord) const

double	detection_efficiency_no_scatter (const unsigned det_num_A, const unsigned det_num_B) const
	average detection efficiency of unscattered counts

static float	max_cos_angle (const float low, const float approx, const float resolution_at_511keV)
	maximum angle to consider above which detection after Compton scatter is considered too small

static float	energy_lower_limit (const float low, const float approx, const float resolution_at_511keV)
	mimumum energy to consider above which detection after Compton scatter is considered too small

integrating functions
float	exp_integral_over_attenuation_image_between_scattpoint_det (const CartesianCoordinate3D< float > &scatter_point, const CartesianCoordinate3D< float > &detector_coord)

float	integral_over_activity_image_between_scattpoint_det (const CartesianCoordinate3D< float > &scatter_point, const CartesianCoordinate3D< float > &detector_coord)

float	cached_integral_over_activity_image_between_scattpoint_det (const unsigned scatter_point_num, const unsigned det_num)

float	cached_exp_integral_over_attenuation_image_between_scattpoint_det (const unsigned scatter_point_num, const unsigned det_num)

static float	integral_between_2_points (const DiscretisedDensity< 3, float > &density, const CartesianCoordinate3D< float > &point1, const CartesianCoordinate3D< float > &point2)

Additional Inherited Members
Static Protected Member Functions inherited from stir::RegisteredObject< ScatterSimulation >
static RegistryType &	registry ()
	Static function returning the registry. More...

Detailed Description

Simulate the scatter probability using a model-based approach.

This base class intends to computes a Compton scatter estimate using an analytical approximation of an integral. It takes as input an emission image and an attenuation image. This is effectively an implementation of the simulation part of the algorithms of Watson and Ollinger. The base class takes care of downsampling and upsampling, book-keeping, caching and looping over bins.

One non-standard feature is that you can specify a different attenuation image to find the scatter points and one to compute the integrals over the attenuation image. The idea is that maybe you want to compute the integrals on a finer grid then you sample the attenuation image. This is probably not very useful though.

Todo:: detector coordinates are derived from ProjDataInfo, but areas and orientations are determined by using a cylindrical scanner.

Todo:: variables/function named density really should use attenuation. This is currently only done for a few variables, but parsing keywords are correct.

References

This implementation is described in the following

C. Tsoumpas, P. Aguiar, K. S. Nikita, D. Ros, K. Thielemans, Evaluation of the Single Scatter Simulation Algorithm Implemented in the STIR Library, Proc. of IEEE Medical Imaging Conf. 2004, vol. 6 pp. 3361 - 3365.

Please refer to the above if you use this implementation.

See also these papers for extra evaluation

N. Dikaios , T. J. Spinks , K. Nikita , K. Thielemans, Evaluation of the Scatter Simulation Algorithm Implemented in STIR, proc. 5th ESBME, Patras, Greece.
P. Aguiar, Ch. Tsoumpas, C. Crespo, J. Pavia, C. Falcon, A. Cot, K. Thielemans and D. Ros, Assessment of scattered photons in the quantification of the small animal PET studies, Eur J Nucl Med Mol I 33:S315-S315 Sep 2006, Proc. EANM 2006, Athens, Greece.
I Polycarpou, P K Marsden and C Tsoumpas, A comparative investigation of scatter correction in 3D PET, J Phys: Conference Series (317), conference 1, 2011

Member Function Documentation

◆ get_density_image_for_scatter_points_sptr()

shared_ptr< const DiscretisedDensity< 3, float > > stir::ScatterSimulation::get_density_image_for_scatter_points_sptr ( ) const

Deprecated:

◆ set_exam_info()

void stir::ScatterSimulation::set_exam_info ( const ExamInfo & arg )

Since July 2016, the information for the energy window and energy resolution are stored in ExamInfo.

References stir::ExamInfo::create_shared_clone(), and template_exam_info_sptr.

Referenced by get_exam_info_sptr().

◆ set_randomly_place_scatter_points()

void stir::ScatterSimulation::set_randomly_place_scatter_points ( const bool arg )

The scattering point in the voxel will be chosen randomly, instead of choosing the centre.

This was first recommended by Watson. It is recommended to leave this on, as otherwise discretisation errors are more obvious.

Note that the random generator is seeded via date/time, so re-running the scatter simulation will give a slightly different result if this boolean is on.

References randomly_place_scatter_points, and use_cache.

◆ downsample_scanner()

Succeeded stir::ScatterSimulation::downsample_scanner	(	int	new_num_rings = `-1`,
		int	new_num_dets = `-1`
	)

Downsample the scanner keeping the total axial length the same.

If new_num_rings<=0, use rings of approximately 2 cm thickness. If new_num_dets <=0, use the default set (currently set in set_defaults())

References _PI, downsample_scanner_dets, downsample_scanner_rings, stir::Scanner::get_axial_length(), stir::Scanner::get_num_detectors_per_ring(), stir::info(), output_proj_data_filename, stir::ProjDataInfo::ProjDataInfoCTI(), stir::round(), and set_output_proj_data().

Referenced by post_processing().

◆ downsample_images_to_scanner_size()

Succeeded stir::ScatterSimulation::downsample_images_to_scanner_size ( )

Downsamples activity and attenuation images to voxel sizes appropriate for the (downsampled) scanner.

This step is not necessary but could result in a speed-up in computing the line integrals. It also avoids problems with too high resolution images compared to the downsampled scanner.

Another way to resolve that is to smooth the images before the scatter simulation. This is currently not implemented in this class.

Warning: This function should be called after having set all data.

References activity_image_sptr, stir::VoxelsOnCartesianGrid< elemT >::get_empty_copy(), remove_cache_for_integrals_over_activity(), remove_cache_for_integrals_over_attenuation(), and stir::zoom_image().

◆ detection_efficiency()

float stir::ScatterSimulation::detection_efficiency ( const float energy ) const

gamma-energy-part of the detection efficiency

Formula used is based on a Gaussian pdf for the efficiency, with

$FWHM_E = \alpha \sqrt(E)$

where the proportionality constant $\alpha$ is determined by the energy FWHM of the Scanner at the reference energy.

This pdf is then integrated from the lower to the upper energy window limits.

See also: Scanner::get_energy_resolution; Scanner::get_reference_energy

References stir::error().

◆ process_data_for_view_segment_num()

double stir::ScatterSimulation::process_data_for_view_segment_num ( const ViewSegmentNumbers & vs_num )

protectedvirtual

computes scatter for one viewgram

Returns: total scatter estimated for this viewgram

References stir::Bin::axial_pos_num(), stir::error(), output_proj_data_sptr, scatter_estimate(), stir::SegmentIndices::segment_num(), stir::Bin::tangential_pos_num(), and stir::ViewgramIndices::view_num().

Referenced by set_use_cache().

◆ post_processing()

bool stir::ScatterSimulation::post_processing ( )

overrideprotectedvirtual

Warning: post_processing will set everything that has a file name in the par file. The corresponding set functions should be used either for files that are not stored in the drive.

Reimplemented from stir::ParsingObject.

Reimplemented in stir::SingleScatterSimulation.

References activity_image_filename, activity_image_sptr, downsample_density_image_for_scatter_points(), downsample_scanner(), stir::error(), stir::ProjDataInfoGenericNoArcCorr::find_cartesian_coordinates_of_detection(), stir::VectorWithOffset< Array< num_dimensions - 1, elemT > >::get_max_index(), stir::VectorWithOffset< Array< num_dimensions - 1, elemT > >::get_min_index(), stir::Scanner::get_num_rings(), stir::Scanner::get_ring_spacing(), stir::VoxelsOnCartesianGrid< elemT >::get_voxel_size(), initialise_cache_for_scattpoint_det_integrals_over_activity(), initialise_cache_for_scattpoint_det_integrals_over_attenuation(), output_proj_data_filename, stir::read_from_file(), remove_cache_for_integrals_over_activity(), remove_cache_for_integrals_over_attenuation(), sample_scatter_points(), set_density_image_for_scatter_points(), set_output_proj_data(), template_exam_info_sptr, zoom_size_xy, zoom_size_z, zoom_xy, and zoom_z.

Referenced by stir::SingleScatterSimulation::post_processing().

◆ sample_scatter_points()

void stir::ScatterSimulation::sample_scatter_points ( )

protected

find scatter points

This function sets scatt_points_vector and scatter_volume. It will also remove any cached integrals as they would be incorrect otherwise.

References stir::convert_int_to_float(), stir::error(), stir::VectorWithOffset< Array< num_dimensions - 1, elemT > >::get_max_index(), stir::VectorWithOffset< Array< num_dimensions - 1, elemT > >::get_min_index(), stir::DiscretisedDensity< num_dimensions, elemT >::get_origin(), and stir::VoxelsOnCartesianGrid< elemT >::get_voxel_size().

Referenced by downsample_density_image_for_scatter_points(), and post_processing().

◆ remove_cache_for_integrals_over_attenuation()

void stir::ScatterSimulation::remove_cache_for_integrals_over_attenuation ( )

protectedvirtual

remove cached attenuation integrals

should be used before recalculating scatter for a new attenuation image or when changing the sampling of the detector etc

Referenced by downsample_density_image_for_scatter_points(), downsample_images_to_scanner_size(), get_exam_info_sptr(), post_processing(), set_defaults(), and set_use_cache().

◆ remove_cache_for_integrals_over_activity()

void stir::ScatterSimulation::remove_cache_for_integrals_over_activity ( )

protectedvirtual

reset cached activity integrals

should be used before recalculating scatter for a new activity image or when changing the sampling of the detector etc

Referenced by downsample_images_to_scanner_size(), get_exam_info_sptr(), post_processing(), set_defaults(), and set_use_cache().

◆ initialise_cache_for_scattpoint_det_integrals_over_attenuation()

void stir::ScatterSimulation::initialise_cache_for_scattpoint_det_integrals_over_attenuation ( )

protected

set-up cache for attenuation integrals

Warning: This will not remove existing cached data (if the sizes match). If you need this, call remove_cache_for_scattpoint_det_integrals_over_attenuation() first.

References stir::VectorWithOffset< T >::resize().

Referenced by post_processing().

◆ initialise_cache_for_scattpoint_det_integrals_over_activity()

void stir::ScatterSimulation::initialise_cache_for_scattpoint_det_integrals_over_activity ( )

protected

set-up cache for activity integrals

Warning: This will not remove existing cached data (if the sizes match). If you need this, call remove_cache_for_scattpoint_det_integrals_over_activity() first.

References stir::VectorWithOffset< T >::resize().

Referenced by post_processing().

Member Data Documentation

◆ template_exam_info_sptr

shared_ptr<ExamInfo> stir::ScatterSimulation::template_exam_info_sptr

protected

Exam info extracted from the scanner template

Referenced by get_exam_info_sptr(), post_processing(), set_exam_info(), and set_output_proj_data().

◆ use_cache

bool stir::ScatterSimulation::use_cache

protected

boolean to see if we need to cache the integrals

By default, we cache the integrals over the emission and attenuation image. If you run out of memory, you can switch this off, but performance will suffer dramatically.

Referenced by ask_parameters(), get_use_cache(), initialise_keymap(), set_defaults(), set_randomly_place_scatter_points(), and set_use_cache().

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

/home/sirfuser/devel/STIRdistrib/STIR/src/include/stir/scatter/ScatterSimulation.h
/home/sirfuser/devel/STIRdistrib/STIR/src/include/stir/scatter/ScatterSimulation.inl
/home/sirfuser/devel/STIRdistrib/STIR/src/scatter_buildblock/cached_single_scatter_integrals.cxx
/home/sirfuser/devel/STIRdistrib/STIR/src/scatter_buildblock/sample_scatter_points.cxx
/home/sirfuser/devel/STIRdistrib/STIR/src/scatter_buildblock/scatter_detection_modelling.cxx
/home/sirfuser/devel/STIRdistrib/STIR/src/scatter_buildblock/ScatterSimulation.cxx
/home/sirfuser/devel/STIRdistrib/STIR/src/scatter_buildblock/single_scatter_integrals.cxx

Public Member Functions

Static Public Member Functions

Protected Types

Protected Member Functions

Protected Attributes

detection related functions

integrating functions

Additional Inherited Members

Detailed Description

Member Function Documentation

◆ get_density_image_for_scatter_points_sptr()

◆ set_exam_info()

◆ set_randomly_place_scatter_points()

◆ downsample_scanner()

◆ downsample_images_to_scanner_size()

◆ detection_efficiency()

◆ process_data_for_view_segment_num()

◆ post_processing()

◆ sample_scatter_points()

◆ remove_cache_for_integrals_over_attenuation()

◆ remove_cache_for_integrals_over_activity()

◆ initialise_cache_for_scattpoint_det_integrals_over_attenuation()

◆ initialise_cache_for_scattpoint_det_integrals_over_activity()

Member Data Documentation

◆ template_exam_info_sptr

◆ use_cache