STIR
6.2.0
|
An objective function class appropriate for PET emission data. More...
#include "stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.h"
Public Member Functions | |
PoissonLogLikelihoodWithLinearModelForMeanAndProjData () | |
Default constructor calls set_defaults() | |
~PoissonLogLikelihoodWithLinearModelForMeanAndProjData () override | |
Destructor. More... | |
TargetT * | construct_target_ptr () const override |
Returns a pointer to a newly allocated target object (with 0 data). More... | |
void | actual_compute_subset_gradient_without_penalty (TargetT &gradient, const TargetT ¤t_estimate, const int subset_num, const bool add_sensitivity) override |
computes the subset gradient of the objective function without the penalty (optional: add subset sensitivity) More... | |
std::unique_ptr< ExamInfo > | get_exam_info_uptr_for_target () const override |
void | add_subset_sensitivity (TargetT &sensitivity, const int subset_num) const override |
Add subset sensitivity to existing data. | |
Functions to get parameters | |
| |
const ProjData & | get_proj_data () const |
const shared_ptr< ProjData > & | get_proj_data_sptr () const |
const int | get_max_segment_num_to_process () const |
const int | get_max_timing_pos_num_to_process () const |
const bool | get_zero_seg0_end_planes () const |
const ProjData & | get_additive_proj_data () const |
const shared_ptr< ProjData > & | get_additive_proj_data_sptr () const |
const ProjectorByBinPair & | get_projector_pair () const |
const shared_ptr< ProjectorByBinPair > & | get_projector_pair_sptr () const |
const int | get_time_frame_num () const |
const TimeFrameDefinitions & | get_time_frame_definitions () const |
const BinNormalisation & | get_normalisation () const |
const shared_ptr< BinNormalisation > & | get_normalisation_sptr () const |
Functions to set parameters | |
This can be used as alternative to the parsing mechanism.
| |
int | set_num_subsets (const int num_subsets) override |
Attempts to change the number of subsets. More... | |
void | set_proj_data_sptr (const shared_ptr< ProjData > &) |
void | set_max_segment_num_to_process (const int) |
void | set_max_timing_pos_num_to_process (const int) |
void | set_zero_seg0_end_planes (const bool) |
void | set_additive_proj_data_sptr (const shared_ptr< ExamData > &) override |
set_additive_proj_data_sptr More... | |
void | set_projector_pair_sptr (const shared_ptr< ProjectorByBinPair > &) |
void | set_frame_num (const int) |
void | set_frame_definitions (const TimeFrameDefinitions &) |
void | set_normalisation_sptr (const shared_ptr< BinNormalisation > &) override |
set_normalisation_sptr More... | |
void | set_input_data (const shared_ptr< ExamData > &) override |
set_input_data More... | |
const ProjData & | get_input_data () const override |
get input data More... | |
Public Member Functions inherited from stir::RegisteredParsingObject< PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >, GeneralisedObjectiveFunction< TargetT >, PoissonLogLikelihoodWithLinearModelForMean< TargetT > > | |
std::string | get_registered_name () const override |
Returns Derived::registered_name. | |
std::string | parameter_info () override |
Returns a string with all parameters and their values, in a form suitable for parsing again. | |
Public Member Functions inherited from stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT > | |
void | compute_sub_gradient_without_penalty (TargetT &gradient, const TargetT ¤t_estimate, const int subset_num) override |
Compute the subset gradient of the (unregularised) objective function. More... | |
virtual void | compute_sub_gradient_without_penalty_plus_sensitivity (TargetT &gradient, const TargetT ¤t_estimate, const int subset_num) |
This should compute the subset gradient of the (unregularised) objective function plus the subset sensitivity. More... | |
Succeeded | set_up (shared_ptr< TargetT > const &target_sptr) override |
set-up sensitivity etc if possible More... | |
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. | |
bool | get_use_subset_sensitivities () const |
find out if subset_sensitivities are used More... | |
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 More... | |
std::string | get_subsensitivity_filenames () const |
get filename pattern to read (or write) the subset sensitivities More... | |
void | fill_nonidentifiable_target_parameters (TargetT &target, const float value) const override |
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 More... | |
void | set_subsensitivity_filenames (const std::string &) |
set filename pattern to read (or write) the subset sensitivities More... | |
Public Member Functions inherited from stir::GeneralisedObjectiveFunction< TargetT > | |
virtual void | compute_sub_gradient (TargetT &gradient, const TargetT ¤t_estimate, const int subset_num) |
Compute the subset-gradient of the objective function at current_estimate. More... | |
virtual void | compute_gradient (TargetT &gradient, const TargetT ¤t_estimate) |
Compute the gradient of the objective function at the current_estimate. More... | |
virtual void | compute_gradient_without_penalty (TargetT &gradient, const TargetT ¤t_estimate) |
Compute the gradient of the unregularised objective function at the current_estimate. More... | |
virtual double | compute_objective_function_without_penalty (const TargetT ¤t_estimate, const int subset_num) |
Compute the value of the unregularised sub-objective function at the current_estimate. More... | |
virtual double | compute_objective_function_without_penalty (const TargetT ¤t_estimate) |
Compute the value of the unregularised objective function at the current_estimate. More... | |
double | compute_penalty (const TargetT ¤t_estimate, const int subset_num) |
Compute the value of the sub-penalty at the current_estimate. More... | |
double | compute_penalty (const TargetT ¤t_estimate) |
Compute the value of the penalty at the current_estimate. More... | |
double | compute_objective_function (const TargetT ¤t_estimate, const int subset_num) |
Compute the value of the sub-objective function at the current_estimate. More... | |
double | compute_objective_function (const TargetT ¤t_estimate) |
Compute the value of the objective function at the current_estimate. More... | |
double | compute_value (const TargetT ¤t_estimate) |
Alias for compute_objective_function(const TargetT&) | |
std::string | get_objective_function_values_report (const TargetT ¤t_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. | |
bool | subsets_are_approximately_balanced () const |
Checks of the current subset scheme is approximately balanced. More... | |
bool | subsets_are_approximately_balanced (std::string &warning_message) const |
Checks of the current subset scheme is approximately balanced and constructs a warning message. More... | |
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. More... | |
shared_ptr< GeneralisedPrior< TargetT > > | get_prior_sptr () |
void | set_prior_sptr (const shared_ptr< GeneralisedPrior< TargetT >> &) |
Change the prior. More... | |
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\'. More... | |
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 ¤t_image_estimate, const TargetT &input) const |
Functions that multiply the True (sub)Hessian with a \'vector\'. More... | |
Succeeded | accumulate_Hessian_times_input_without_penalty (TargetT &output, const TargetT ¤t_image_estimate, const TargetT &input) const |
Succeeded | accumulate_sub_Hessian_times_input (TargetT &output, const TargetT ¤t_image_estimate, const TargetT &input, const int subset_num) const |
Succeeded | accumulate_sub_Hessian_times_input_without_penalty (TargetT &output, const TargetT ¤t_image_estimate, const TargetT &input, const int subset_num) const |
Public Member Functions inherited from stir::ParsingObject | |
ParsingObject (const ParsingObject &) | |
ParsingObject & | operator= (const ParsingObject &) |
void | ask_parameters () |
bool | parse (std::istream &f) |
bool | parse (const char *const filename) |
Public Attributes | |
Variables for STIR_MPI | |
Only used when STIR_MPI is enabled.
| |
DistributedCachingInformation * | caching_info_ptr |
points to the information object needed to support distributed caching | |
bool | distributed_cache_enabled |
enable/disable key for distributed caching | |
bool | distributed_tests_enabled |
bool | message_timings_enabled |
double | message_timings_threshold |
bool | rpc_timings_enabled |
Static Public Attributes | |
static const char *const | registered_name = "PoissonLogLikelihoodWithLinearModelForMeanAndProjData" |
Name which will be used when parsing a GeneralisedObjectiveFunction object. | |
Protected Member Functions | |
Succeeded | set_up_before_sensitivity (shared_ptr< const TargetT > const &target_sptr) override |
set-up specifics for the derived class | |
double | actual_compute_objective_function_without_penalty (const TargetT ¤t_estimate, const int subset_num) override |
Implementation of function that computes the objective function for the current subset. More... | |
Succeeded | actual_add_multiplication_with_approximate_sub_Hessian_without_penalty (TargetT &output, const TargetT &input, const int subset_num) const override |
Succeeded | actual_accumulate_sub_Hessian_times_input_without_penalty (TargetT &output, const TargetT ¤t_image_estimate, const TargetT &input, const int subset_num) const override |
void | set_defaults () override |
sets any default values More... | |
void | initialise_keymap () override |
sets keys for parsing More... | |
bool | post_processing () override |
checks values after parsing More... | |
bool | actual_subsets_are_approximately_balanced (std::string &warning_message) const override |
Checks of the current subset scheme is approximately balanced. More... | |
Protected Member Functions inherited from stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT > | |
void | compute_sensitivities () |
compute subset and total sensitivity 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::string | input_filename |
Filename with input projection data. | |
shared_ptr< ProjData > | proj_data_sptr |
points to the object for the total input projection data | |
int | max_segment_num_to_process |
the maximum absolute ring difference number to use in the reconstruction More... | |
int | max_timing_pos_num_to_process |
the maximum absolute time-of-flight bin number to use in the reconstruction More... | |
ParseAndCreateFrom< TargetT, ProjData > | target_parameter_parser |
shared_ptr< ProjectorByBinPair > | projector_pair_ptr |
Stores the projectors that are used for the computations. | |
bool | zero_seg0_end_planes |
signals whether to zero the data in the end planes of the projection data | |
bool | use_tofsens |
Triggers calculation of sensitivity using time-of-flight. | |
std::string | additive_projection_data_filename |
name of file in which additive projection data are stored | |
shared_ptr< ProjData > | additive_proj_data_sptr |
shared_ptr< BinNormalisation > | normalisation_sptr |
int | frame_num |
std::string | frame_definition_filename |
TimeFrameDefinitions | frame_defs |
int | loglikelihood_computation_interval |
subiteration interval at which the loglikelihood function is evaluated | |
bool | compute_total_loglikelihood |
indicates whether to evaluate the loglikelihood function for all bins or the current subset | |
std::string | loglikelihood_data_filename |
name of file in which loglikelihood measurements are stored | |
Protected Attributes inherited from stir::GeneralisedObjectiveFunction< TargetT > | |
int | num_subsets |
bool | already_set_up |
shared_ptr< GeneralisedPrior< TargetT > > | prior_sptr |
Protected Attributes inherited from stir::ParsingObject | |
KeyParser | parser |
Additional Inherited Members | |
Static Public Member Functions inherited from stir::RegisteredParsingObject< PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >, GeneralisedObjectiveFunction< TargetT >, PoissonLogLikelihoodWithLinearModelForMean< TargetT > > | |
static GeneralisedObjectiveFunction< TargetT > * | read_from_stream (std::istream *) |
Construct a new object (of type Derived) by parsing the istream. More... | |
Static Public Member Functions inherited from stir::RegisteredObject< GeneralisedObjectiveFunction< TargetT > > | |
static GeneralisedObjectiveFunction< TargetT > * | read_registered_object (std::istream *in, const std::string ®istered_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 GeneralisedObjectiveFunction< TargetT > * | 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 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_less > | RegistryType |
The type of the registry. | |
Static Protected Member Functions inherited from stir::RegisteredObject< GeneralisedObjectiveFunction< TargetT > > | |
static RegistryType & | registry () |
Static function returning the registry. More... | |
An objective function class appropriate for PET emission data.
Measured data is given by a ProjData object, and the linear operations necessary for computing the gradient of the objective function are performed via a ProjectorByBinPair object together with a BinNormalisation object.
Often, the probability matrix can be written as the product of a diagonal matrix and another matrix
The measurement model can then be written as
and backprojection is obviously
This can be generalised by using a different matrix to perform the backprojection.
The expression for the gradient becomes
where dropped from the first term. This was probably noticed for the first time in the context of MLEM by Hebert and Leahy, but is in fact a feature of the gradient (and indeed log-likelihood). This was first brought to Kris Thielemans' attention by Michael Zibulevsky and Arkady Nemirovsky.
Note that if and are not exact transposed operations, the gradient does no longer correspond to the gradient of a well-defined function.
In this class, the operation of multiplying with is performed by the forward projector, and with by the back projector. Multiplying with corresponds to the BinNormalisation::undo() function. Finally, the background term is stored in additive_proj_data_sptr
(note: this is not ).
PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:= ; emission projection data input file := maximum absolute segment number to process := zero end planes of segment 0 := ; see ProjectorByBinPair hierarchy for possible values Projector pair type := ; reserved value: 0 means none ; see PoissonLogLikelihoodWithLinearModelForMeanAndProjData ; class documentation additive sinogram := ; normalisation (and attenuation correction) ; time info can be used for dead-time correction ; see TimeFrameDefinitions time frame definition filename := time frame number := ; see BinNormalisation hierarchy for possible values Bin Normalisation type := End PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters :=
|
override |
|
overridevirtual |
Returns a pointer to a newly allocated target object (with 0 data).
Dimensions etc are set from the proj_data_sptr and other information set by parsing, such as zoom
, output_image_size_z
etc.
Implements stir::GeneralisedObjectiveFunction< TargetT >.
|
overridevirtual |
Attempts to change the number of subsets.
Implements stir::GeneralisedObjectiveFunction< TargetT >.
Referenced by stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::test_approximate_Hessian_concavity().
|
overridevirtual |
set_additive_proj_data_sptr
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 >.
Referenced by stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::test_approximate_Hessian_concavity().
|
overridevirtual |
set_normalisation_sptr
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 >.
Referenced by stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::test_approximate_Hessian_concavity().
|
overridevirtual |
set_input_data
It can be used to set the data to be reconstructed within some other code, as opposed to via parsing.
Implements stir::GeneralisedObjectiveFunction< TargetT >.
|
overridevirtual |
get input data
Will throw an exception if it wasn't set first
Implements stir::GeneralisedObjectiveFunction< TargetT >.
|
overridevirtual |
computes the subset gradient of the objective function without the penalty (optional: add subset sensitivity)
If add_sensitivity
is true
, this computes
(see the class general documentation). The sum will however be restricted to a subset.
However, if add_sensitivity
is false
, this function will instead compute only the gradient
Implements stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >.
References stir::error(), and stir::setup_distributable_computation().
|
overrideprotectedvirtual |
Implementation of function that computes the objective function for the current subset.
The reason we have this function is that overloading a function in a derived class, hides all functions of the same name.
Implements stir::GeneralisedObjectiveFunction< TargetT >.
References stir::error(), and stir::setup_distributable_computation().
|
overrideprotectedvirtual |
The Hessian (without penalty) is approximately given by:
where
and is the probability matrix.
This function uses the approximation of the hessian
Hence
In the above, we've used the plug-in approximation by replacing forward projection of the true image by the measured data. However, the later are noisy and this can create problems.
The LogLikelihood is a concave function and therefore the Hessian is non-positive. Thus, this (approximate-)Hessian vector product methods output volumes with negative values, if the input is non-negative.
One could smooth the data before performing the quotient. This should be done after normalisation to avoid problems with the high-frequency components in the normalisation factors:
where the probability matrix is factorised in a detection efficiency part (i.e. the normalisation factors ) times a geometric part:
It has also been suggested to use (at least if the data are still Poisson).
Reimplemented from stir::GeneralisedObjectiveFunction< TargetT >.
References stir::error().
|
overrideprotectedvirtual |
The Hessian (without penalty) is approximately given by:
where
and is the probability matrix.
Hence
This function is computationally expensive and can be approximated, see add_multiplication_with_approximate_sub_Hessian_without_penalty()
The loglikelihood is a concave function, see add_multiplication_with_approximate_sub_Hessian_without_penalty() for more details regarding Hessian methods.
Reimplemented from stir::GeneralisedObjectiveFunction< TargetT >.
References stir::error().
|
overrideprotectedvirtual |
sets any default values
Has to be called by set_defaults in the leaf-class
Reimplemented from stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >.
|
overrideprotectedvirtual |
sets keys for parsing
Has to be called by initialise_keymap in the leaf-class
Reimplemented from stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >.
|
overrideprotectedvirtual |
checks values after parsing
Has to be called by post_processing in the leaf-class
Reimplemented from stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >.
References stir::read_from_file(), and stir::warning().
|
overrideprotectedvirtual |
Checks of the current subset scheme is approximately balanced.
For this class, this means that the sub-sensitivities are approximately the same. The test simply looks at the number of views etc. It ignores unbalancing caused by normalisation_sptr (e.g. for instance when using asymmetric attenuation).
Implements stir::GeneralisedObjectiveFunction< TargetT >.
References stir::error(), stir::Array< num_dimensions, elemT >::fill(), stir::DataSymmetriesForViewSegmentNumbers::is_basic(), and stir::DataSymmetriesForViewSegmentNumbers::num_related_view_segment_numbers().
|
protected |
the maximum absolute ring difference number to use in the reconstruction
convention: if -1, use get_max_segment_num()
|
protected |
the maximum absolute time-of-flight bin number to use in the reconstruction
convention: if -1, use get_max_tof_pos_num()