STIR 6.4.0
stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData< TargetT > Class Template Reference

An objective function class appropriate for PET list mode data. More...

#include "stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeData.h"

Inheritance diagram for stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData< TargetT >:

Public Member Functions

Succeeded set_up_before_sensitivity (shared_ptr< const TargetT > const &target_sptr) override
 set-up specifics for the derived class
 
void set_normalisation_sptr (const shared_ptr< BinNormalisation > &) override
 set_normalisation_sptr
 
void set_additive_proj_data_sptr (const shared_ptr< ExamData > &) override
 set_additive_proj_data_sptr
 
void set_input_data (const shared_ptr< ExamData > &) override
 set_input_data
 
const ListModeDataget_input_data () const override
 get input data
 
void set_max_segment_num_to_process (int)
 set maximum segment_number (in listmode data) to process
 
int get_max_segment_num_to_process () const
 get maximum segment_number (from listmode data) to process
 
- Public Member Functions inherited from stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >
void compute_sub_gradient_without_penalty (TargetT &gradient, const TargetT &current_estimate, const int subset_num) override
 Compute the subset gradient of the (unregularised) objective function.
 
virtual void compute_sub_gradient_without_penalty_plus_sensitivity (TargetT &gradient, const TargetT &current_estimate, const int subset_num)
 This should compute the subset gradient of the (unregularised) objective function plus the subset sensitivity.
 
Succeeded set_up (shared_ptr< TargetT > const &target_sptr) override
 set-up sensitivity etc if possible
 
const TargetT & get_sensitivity () const
 Get a const reference to the total sensitivity.
 
const TargetT & get_subset_sensitivity (const int subset_num) const
 Get a const reference to the sensitivity for a subset.
 
virtual void add_subset_sensitivity (TargetT &sensitivity, const int subset_num) const =0
 Add subset sensitivity to existing data.
 
bool get_use_subset_sensitivities () const
 find out if subset_sensitivities are used
 
bool get_recompute_sensitivity () const
 find current value of recompute_sensitivity
 
std::string get_sensitivity_filename () const
 get filename to read (or write) the total sensitivity
 
std::string get_subsensitivity_filenames () const
 get filename pattern to read (or write) the subset sensitivities
 
std::string get_subsensitivity_filename (const int subset_num) const
 Return the filename for a particular subset.
 
void set_recompute_sensitivity (const bool)
 
void set_subset_sensitivity_sptr (const shared_ptr< TargetT > &, const int subset_num)
 
void set_use_subset_sensitivities (const bool)
 See get_use_subset_sensitivities()
 
void set_sensitivity_filename (const std::string &)
 set filename to read (or write) the total sensitivity
 
void set_subsensitivity_filenames (const std::string &)
 set filename pattern to read (or write) the subset sensitivities
 
void fill_nonidentifiable_target_parameters (TargetT &target, const float value) const override
 
- Public Member Functions inherited from stir::GeneralisedObjectiveFunction< TargetT >
virtual TargetT * construct_target_ptr () const =0
 Creates a suitable target as determined by the parameters.
 
virtual void compute_sub_gradient (TargetT &gradient, const TargetT &current_estimate, const int subset_num)
 Compute the subset-gradient of the objective function at current_estimate.
 
virtual void compute_gradient (TargetT &gradient, const TargetT &current_estimate)
 Compute the gradient of the objective function at the current_estimate.
 
virtual void compute_gradient_without_penalty (TargetT &gradient, const TargetT &current_estimate)
 Compute the gradient of the unregularised objective function at the current_estimate.
 
virtual double compute_objective_function_without_penalty (const TargetT &current_estimate, const int subset_num)
 Compute the value of the unregularised sub-objective function at the current_estimate.
 
virtual double compute_objective_function_without_penalty (const TargetT &current_estimate)
 Compute the value of the unregularised objective function at the current_estimate.
 
double compute_penalty (const TargetT &current_estimate, const int subset_num)
 Compute the value of the sub-penalty at the current_estimate.
 
double compute_penalty (const TargetT &current_estimate)
 Compute the value of the penalty at the current_estimate.
 
double compute_objective_function (const TargetT &current_estimate, const int subset_num)
 Compute the value of the sub-objective function at the current_estimate.
 
double compute_objective_function (const TargetT &current_estimate)
 Compute the value of the objective function at the current_estimate.
 
double compute_value (const TargetT &current_estimate)
 Alias for compute_objective_function(const TargetT&)
 
Succeeded add_multiplication_with_approximate_sub_Hessian_without_penalty (TargetT &output, const TargetT &input, const int subset_num) const
 Functions that multiply the approximate (sub)Hessian with a \'vector\'.
 
Succeeded add_multiplication_with_approximate_sub_Hessian (TargetT &output, const TargetT &input, const int subset_num) const
 
Succeeded add_multiplication_with_approximate_Hessian_without_penalty (TargetT &output, const TargetT &input) const
 
Succeeded add_multiplication_with_approximate_Hessian (TargetT &output, const TargetT &input) const
 
Succeeded accumulate_Hessian_times_input (TargetT &output, const TargetT &current_image_estimate, const TargetT &input) const
 Functions that multiply the True (sub)Hessian with a \'vector\'.
 
Succeeded accumulate_Hessian_times_input_without_penalty (TargetT &output, const TargetT &current_image_estimate, const TargetT &input) const
 
Succeeded accumulate_sub_Hessian_times_input (TargetT &output, const TargetT &current_image_estimate, const TargetT &input, const int subset_num) const
 
Succeeded accumulate_sub_Hessian_times_input_without_penalty (TargetT &output, const TargetT &current_image_estimate, const TargetT &input, const int subset_num) const
 
std::string get_objective_function_values_report (const TargetT &current_estimate)
 Construct a string with info on the value of objective function with and without penalty.
 
int get_num_subsets () const
 Return the number of subsets in-use.
 
virtual std::unique_ptr< ExamInfoget_exam_info_uptr_for_target () const
 
virtual int set_num_subsets (const int num_subsets)=0
 Attempts to change the number of subsets.
 
bool subsets_are_approximately_balanced () const
 Checks of the current subset scheme is approximately balanced.
 
bool subsets_are_approximately_balanced (std::string &warning_message) const
 Checks of the current subset scheme is approximately balanced and constructs a warning message.
 
bool prior_is_zero () const
 check if the prior is set (or the penalisation factor is 0)
 
GeneralisedPrior< TargetT > *const get_prior_ptr () const
 Read-only access to the prior.
 
shared_ptr< GeneralisedPrior< TargetT > > get_prior_sptr ()
 
void set_prior_sptr (const shared_ptr< GeneralisedPrior< TargetT > > &)
 Change the prior.
 
- Public Member Functions inherited from stir::RegisteredObjectBase
virtual std::string get_registered_name () const =0
 Returns the name of the type of the object.
 
- Public Member Functions inherited from stir::ParsingObject
 ParsingObject (const ParsingObject &)
 
ParsingObjectoperator= (const ParsingObject &)
 
bool parse (std::istream &f)
 
bool parse (const char *const filename)
 
void ask_parameters ()
 
virtual std::string parameter_info ()
 

Public Attributes

TimeFrameDefinitions frame_defs
 time frame definitions
 

caching-related methods

These functions can be used to cache listmode events into memory, allowing parallelised processing.

Currently, the cached data is written to one or more files (

See also
get_cache_filename)
Warning
This code is experimental and likely to change in future versions.
When re-using an existing cache, there is no check if time-frames etc are the same as what was used when creating the cache. This is therefore quite risky.
Cache-files are written in a binary format that likely depends on compiler, endianness etc.
Todo
It should be possible to read only part of the cache in memory.
std::string frame_defs_filename
 
std::string list_mode_filename
 Filename with input projection data.
 
shared_ptr< ProjDataadditive_proj_data_sptr
 
std::string additive_projection_data_filename
 filename for additive data (only used when parsing)
 
bool reduce_memory_usage
 If true, the additive sinogram will not be read in memory.
 
shared_ptr< BinNormalisationnormalisation_sptr
 
shared_ptr< ListModeDatalist_mode_data_sptr
 Listmode pointer.
 
unsigned int current_frame_num
 
long int num_events_to_use
 This is part of some functionality I transfer from LmToProjData.
 
bool do_time_frame
 Reconstruct based on time frames.
 
ParseAndCreateFrom< TargetT, ListModeDatatarget_parameter_parser
 
unsigned long int cache_size
 This is the number of records to be cached. If this parameter is more than zero, then the flag cache_lm_file will be set to true. The listmode file up to this size will be loaded in the RAM, alongside with any additive sinograms.
 
bool cache_lm_file
 This flag is true when cache_size is more than zero.
 
bool recompute_cache
 On the first cached run, the cache will be written in the cache_path. If recompute_cache is set to zero then every consecutive reconstruction will use that cache file. If you want to create a new, either delete the previous or set this 1.
 
bool skip_lm_input_file
 This flag is set when we don't set an input lm filename and rely only on the cache file.
 
std::string cache_path
 Path to read/write the cached listmode file.
 
bool has_add
 The data set has additive corrections.
 
shared_ptr< ProjDataInfoproj_data_info_sptr
 ProjDataInfo.
 
virtual void set_cache_path (const std::string &cache_path)
 Set the directory where data will be cached.
 
virtual std::string get_cache_path () const
 Get the directory where data will be cached.
 
virtual std::string get_cache_filename (unsigned int icache) const
 Get the filename for a cache file.
 
void set_recompute_cache (bool)
 Set if existing cache files should be used or not.
 
bool get_recompute_cache () const
 
void set_skip_lm_input_file (const bool arg)
 Skip reading of listmode file entirely, only read from cache (broken!)
 
virtual void set_cache_max_size (const unsigned long int arg)
 Set maximum size (in bytes) of cache in memory.
 
virtual unsigned long int get_cache_max_size () const
 Get maximum size (in bytes) of cache.
 
void set_defaults () override
 sets any default values
 
void initialise_keymap () override
 sets keys
 
bool post_processing () override
 This will be called at the end of the parsing.
 
virtual void start_new_time_frame (const unsigned int new_frame_num)
 will be called when a new time frame starts
 

Additional Inherited Members

- Static Public Member Functions inherited from stir::RegisteredObject< GeneralisedObjectiveFunction< TargetT > >
static GeneralisedObjectiveFunction< TargetT > * 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.
 
static GeneralisedObjectiveFunction< TargetT > * ask_type_and_parameters ()
 ask the user for the type, and then calls read_registered_object(0, type)
 
static void list_registered_names (std::ostream &stream)
 List all possible registered names to the stream.
 
- Protected Types inherited from stir::RegisteredObject< GeneralisedObjectiveFunction< TargetT > >
typedef GeneralisedObjectiveFunction< TargetT > *(* 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_lessRegistryType
 The type of the registry.
 
void compute_sensitivities ()
 compute subset and total sensitivity
 
virtual void actual_compute_subset_gradient_without_penalty (TargetT &gradient, const TargetT &current_estimate, const int subset_num, const bool add_sensitivity)=0
 computes the subset gradient of the objective function without the penalty (optional: add subset sensitivity)
 
virtual bool actual_subsets_are_approximately_balanced (std::string &warning_message) const =0
 Implementation of function that checks subset balancing.
 
virtual double actual_compute_objective_function_without_penalty (const TargetT &current_estimate, const int subset_num)=0
 Implementation of function that computes the objective function for the current subset.
 
virtual Succeeded actual_add_multiplication_with_approximate_sub_Hessian_without_penalty (TargetT &output, const TargetT &input, const int subset_num) const
 Implementation of the function that multiplies the approximate sub-Hessian with a vector.
 
virtual Succeeded actual_accumulate_sub_Hessian_times_input_without_penalty (TargetT &output, const TargetT &current_image_estimate, const TargetT &input, const int subset_num) const
 Implementation of the function computes the sub-Hessian and multiplies by a vector.
 
virtual void set_key_values ()
 This will be called before parsing or parameter_info is called.
 
- Static Protected Member Functions inherited from stir::RegisteredObject< GeneralisedObjectiveFunction< TargetT > >
static RegistryTyperegistry ()
 Static function returning the registry.
 
int num_subsets
 
bool already_set_up
 
shared_ptr< GeneralisedPrior< TargetT > > prior_sptr
 
KeyParser parser
 

Detailed Description

template<typename TargetT>
class stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData< TargetT >

An objective function class appropriate for PET list mode data.

The statistics for list mode data is slightly different from having a set of counts, see the paper H. H. Barrett, L. Parra, and T. White, List-mode likelihood, J. Optical Soc. Amer. A, vol. 14, no. 11, 1997.

However, it is intuitive that the list mode likelihood can be derived from that 'binned' likelihood by taking smaller and smaller time bins. One can then see that the gradient of the list mode likelihood can be computed similar to the 'binned' case, but now with a sum over events. The sensitivity still needs a sum (or integral) over all possible detections.

At present, STIR does not contain any classes to forward/back project list mode data explictly. So, the sum over events cannot be computed for arbitrary list mode data. This is currently done in derived classes.

For list mode reconstructions, computing the sensitivity is sometimes conceptually very difficult (how to compute the integral) or impractical. So, even if we will be able to compute the sum over events in a generic way, the add_subset_sensitivity() function will have to be implemented by a derived class, specific for the measurement.

Member Function Documentation

◆ set_up_before_sensitivity()

◆ set_normalisation_sptr()

template<typename TargetT>
void stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData< TargetT >::set_normalisation_sptr ( const shared_ptr< BinNormalisation > & )
overridevirtual

set_normalisation_sptr

Author
Nikos Efthimiou

In the case the reconstruction process is called from another piece of code, the user should be able to set any additive sinogram

Implements stir::GeneralisedObjectiveFunction< TargetT >.

References set_normalisation_sptr().

Referenced by set_normalisation_sptr().

◆ set_additive_proj_data_sptr()

template<typename TargetT>
void stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData< TargetT >::set_additive_proj_data_sptr ( const shared_ptr< ExamData > & )
overridevirtual

set_additive_proj_data_sptr

Author
Nikos Efthimiou

In the case the reconstruction process is called from another piece of code, the user should be able to set any additive sinogram

Implements stir::GeneralisedObjectiveFunction< TargetT >.

References stir::error(), has_add, reduce_memory_usage, and set_additive_proj_data_sptr().

Referenced by set_additive_proj_data_sptr().

◆ set_input_data()

template<typename TargetT>
void stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData< TargetT >::set_input_data ( const shared_ptr< ExamData > & )
overridevirtual

set_input_data

Author
Nikos Efthimiou

It can be used to set the data to be reconstructed within some other code, as opposed to via parsing.

Implements stir::GeneralisedObjectiveFunction< TargetT >.

References list_mode_data_sptr, and set_input_data().

Referenced by set_input_data().

◆ get_input_data()

template<typename TargetT>
const ListModeData & stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData< TargetT >::get_input_data ( ) const
overridevirtual

◆ set_max_segment_num_to_process()

template<typename TargetT>
void stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData< TargetT >::set_max_segment_num_to_process ( int arg)

set maximum segment_number (in listmode data) to process

minimum will be -max_segment_num_to_process

Use -1 to process all.

References set_max_segment_num_to_process().

Referenced by set_max_segment_num_to_process().

◆ get_max_segment_num_to_process()

template<typename TargetT>
int stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData< TargetT >::get_max_segment_num_to_process ( ) const

◆ set_cache_path()

template<typename TargetT>
void stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData< TargetT >::set_cache_path ( const std::string & cache_path)
virtual

Set the directory where data will be cached.

Parameters
cache_pathdirectory-name (defaults to current directory). The directory has to exist.

References cache_path, and set_cache_path().

Referenced by set_cache_path().

◆ get_cache_filename()

template<typename TargetT>
std::string stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData< TargetT >::get_cache_filename ( unsigned int icache) const
virtual

Get the filename for a cache file.

Filenames are currently "my_CACHE%d.bin" with get_cache_path() prepended.

References stir::FilePath::get_as_string(), get_cache_filename(), get_cache_path(), and stir::FilePath::prepend_directory_name().

Referenced by get_cache_filename().

◆ set_recompute_cache()

template<typename TargetT>
void stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData< TargetT >::set_recompute_cache ( bool v)

Set if existing cache files should be used or not.

defaults to true

References recompute_cache, and set_recompute_cache().

Referenced by set_recompute_cache().

◆ set_skip_lm_input_file()

template<typename TargetT>
void stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData< TargetT >::set_skip_lm_input_file ( const bool arg)

Skip reading of listmode file entirely, only read from cache (broken!)

Warning
This currently aborts, as functionality is broken. (We would need to be able to read proj_data_info and exam_info).
Todo
replace with reading from a custom-listmode file (although this would have to support the additive term).
Todo
in the future the following statements should be removed.

References cache_path, stir::error(), set_skip_lm_input_file(), stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >::set_subsensitivity_filenames(), stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >::set_use_subset_sensitivities(), and skip_lm_input_file.

Referenced by set_skip_lm_input_file().

◆ set_cache_max_size()

template<typename TargetT>
void stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData< TargetT >::set_cache_max_size ( const unsigned long int arg)
virtual

Set maximum size (in bytes) of cache in memory.

When all events do not fit in the cache, several cache files will be used.

As multi-threading occurs over cached events, it is likely that better performance will be obtained with a large cache size.

References cache_size, and set_cache_max_size().

Referenced by set_cache_max_size().

◆ set_defaults()

◆ initialise_keymap()

◆ post_processing()

◆ start_new_time_frame()

template<typename TargetT>
void stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData< TargetT >::start_new_time_frame ( const unsigned int new_frame_num)
protectedvirtual

will be called when a new time frame starts

The frame numbers start from 1.

References start_new_time_frame().

Referenced by start_new_time_frame().

Member Data Documentation

◆ frame_defs

time frame definitions

Todo
This is currently used to be able to compute the gradient for one time frame. However, it probably does not belong here. For instance when fitting kinetic model parameters from list mode data, time frames are in principle irrelevant. So, we will probably shift this to the derived class.

Referenced by set_defaults(), set_up_before_sensitivity(), and stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin< TargetT >::set_up_before_sensitivity().

◆ recompute_cache

template<typename TargetT>
bool stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData< TargetT >::recompute_cache
protected

On the first cached run, the cache will be written in the cache_path. If recompute_cache is set to zero then every consecutive reconstruction will use that cache file. If you want to create a new, either delete the previous or set this 1.

Todo
multiple cache files need to be supported!

Referenced by set_defaults(), and set_recompute_cache().

◆ cache_path

template<typename TargetT>
std::string stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData< TargetT >::cache_path
protected

Path to read/write the cached listmode file.

Todo
add the ability to set a filename.

Referenced by get_cache_path(), set_cache_path(), set_defaults(), and set_skip_lm_input_file().

◆ proj_data_info_sptr


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