An objective function class appropriate for PET emission data. More...

#include "stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.h"

Inheritance diagram for stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >:

[legend]

Public Member Functions
Functions to get parameters
Warning Be careful with changing shared pointers. If you modify the objects in one place, all objects that use the shared pointer will be affected.
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

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.

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 &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.

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 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.

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::ParsingObject
	ParsingObject (const ParsingObject &)

ParsingObject &	operator= (const ParsingObject &)

bool	parse (std::istream &f)

bool	parse (const char *const filename)

void	ask_parameters ()

Static Public Attributes
static const char *const	registered_name
	Name which will be used when parsing a GeneralisedObjectiveFunction object.

Variables for STIR_MPI
Only used when STIR_MPI is enabled. Todo move to protected area
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

	PoissonLogLikelihoodWithLinearModelForMeanAndProjData ()
	Default constructor calls set_defaults()

	~PoissonLogLikelihoodWithLinearModelForMeanAndProjData () override
	Destructor.

TargetT *	construct_target_ptr () const override
	Returns a pointer to a newly allocated target object (with 0 data).

Functions to set parameters
This can be used as alternative to the parsing mechanism. Warning After using any of these, you have to call set_up(). Be careful with changing shared pointers. If you modify the objects in one place, all objects that use the shared pointer will be affected.
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

int	max_timing_pos_num_to_process
	the maximum absolute time-of-flight bin number to use in the reconstruction

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

int	set_num_subsets (const int num_subsets) override
	Attempts to change the number of subsets.

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

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

void	set_input_data (const shared_ptr< ExamData > &) override
	set_input_data

const ProjData &	get_input_data () const override
	get input data

void	actual_compute_subset_gradient_without_penalty (TargetT &gradient, const TargetT &current_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)

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.

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 &current_estimate, const int subset_num) override
	Implementation of function that computes the objective function for the current subset.

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 &current_image_estimate, const TargetT &input, const int subset_num) const override

void	set_defaults () override
	sets any default values

void	initialise_keymap () override
	sets keys for parsing

bool	post_processing () override
	checks values after parsing

bool	actual_subsets_are_approximately_balanced (std::string &warning_message) const override
	Checks of the current subset scheme is approximately balanced.

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.

static GeneralisedObjectiveFunction< TargetT > *	read_from_stream (std::istream *)
	Construct a new object (of type Derived) by parsing the istream.

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_less >	RegistryType
	The type of the registry.

void	compute_sensitivities ()
	compute subset and total sensitivity

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 RegistryType &	registry ()
	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::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >

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.

See also: PoissonLogLikelihoodWithLinearModelForMean for notation.

Often, the probability matrix can be written as the product of a diagonal matrix and another matrix

$P = D G$

The measurement model can then be written as

$P \lambda + r = D ( F \lambda + a )$

and backprojection is obviously

$P' y = F' D y$

This can be generalised by using a different matrix to perform the backprojection.

The expression for the gradient becomes

$B D \left[ y / \left( D (F \lambda + a) \right) \right] - B D 1 = B \left[ y / \left(F \lambda + a \right) \right] - B D 1$

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 ).

Parameters for parsing

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 :=

Constructor & Destructor Documentation

◆ ~PoissonLogLikelihoodWithLinearModelForMeanAndProjData()

template<typename TargetT>

stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >::~PoissonLogLikelihoodWithLinearModelForMeanAndProjData ( )

override

Destructor.

Calls end_distributable_computation()

References stir::end_distributable_computation(), and ~PoissonLogLikelihoodWithLinearModelForMeanAndProjData().

Referenced by ~PoissonLogLikelihoodWithLinearModelForMeanAndProjData().

Member Function Documentation

◆ construct_target_ptr()

template<typename TargetT>

TargetT * stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >::construct_target_ptr ( ) const

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 >.

References construct_target_ptr(), and get_input_data().

Referenced by construct_target_ptr().

◆ set_num_subsets()

template<typename TargetT>

int stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >::set_num_subsets ( const int num_subsets )

overridevirtual

Attempts to change the number of subsets.

Returns: The number of subsets that will be used later, which is not guaranteed to be what you asked for.

Implements stir::GeneralisedObjectiveFunction< TargetT >.

References set_num_subsets().

Referenced by set_num_subsets().

◆ set_additive_proj_data_sptr()

template<typename TargetT>

void stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< 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 set_additive_proj_data_sptr().

Referenced by set_additive_proj_data_sptr().

◆ set_normalisation_sptr()

template<typename TargetT>

void stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< 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_input_data()

template<typename TargetT>

void stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< 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 proj_data_sptr, and set_input_data().

Referenced by set_input_data().

◆ get_input_data()

template<typename TargetT>

const ProjData & stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >::get_input_data ( ) const

overridevirtual

get input data

Will throw an exception if it wasn't set first

Implements stir::GeneralisedObjectiveFunction< TargetT >.

References get_input_data(), and proj_data_sptr.

Referenced by construct_target_ptr(), and get_input_data().

◆ actual_compute_subset_gradient_without_penalty()

template<typename TargetT>

void stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >::actual_compute_subset_gradient_without_penalty	(	TargetT &	gradient,
		const TargetT &	current_estimate,
		const int	subset_num,
		const bool	add_sensitivity )

overridevirtual

computes the subset gradient of the objective function without the penalty (optional: add subset sensitivity)

If add_sensitivity is true, this computes

${\partial L \over \partial \lambda_v} + P_v = \sum_b P_{bv} {y_b \over Y_b}$

(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

${\partial L \over \partial \lambda_v} = \sum_b P_{bv} ({y_b \over Y_b} - 1)$

Implements stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >.

References actual_compute_subset_gradient_without_penalty(), caching_info_ptr, distributed_cache_enabled, stir::error(), proj_data_sptr, projector_pair_ptr, stir::setup_distributable_computation(), and zero_seg0_end_planes.

Referenced by actual_compute_subset_gradient_without_penalty().

◆ get_exam_info_uptr_for_target()

template<typename TargetT>

std::unique_ptr< ExamInfo > stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >::get_exam_info_uptr_for_target ( ) const

overridevirtual

Reimplemented from stir::GeneralisedObjectiveFunction< TargetT >.

◆ add_subset_sensitivity()

template<typename TargetT>

void stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >::add_subset_sensitivity	(	TargetT &	sensitivity,
		const int	subset_num ) const

overridevirtual

Add subset sensitivity to existing data.

Implements stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >.

References add_subset_sensitivity(), caching_info_ptr, distributed_cache_enabled, stir::error(), max_segment_num_to_process, max_timing_pos_num_to_process, proj_data_sptr, projector_pair_ptr, stir::setup_distributable_computation(), use_tofsens, and zero_seg0_end_planes.

Referenced by add_subset_sensitivity().

◆ set_up_before_sensitivity()

template<typename TargetT>

Succeeded stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >::set_up_before_sensitivity ( shared_ptr< const TargetT > const & target_sptr )

overrideprotectedvirtual

set-up specifics for the derived class

Implements stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >.

References caching_info_ptr, distributed_cache_enabled, stir::error(), stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >::get_recompute_sensitivity(), stir::info(), max_segment_num_to_process, max_timing_pos_num_to_process, distributed::num_processors, proj_data_sptr, projector_pair_ptr, set_up_before_sensitivity(), and use_tofsens.

Referenced by set_up_before_sensitivity().

◆ actual_compute_objective_function_without_penalty()

template<typename TargetT>

double stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >::actual_compute_objective_function_without_penalty	(	const TargetT &	current_estimate,
		const int	subset_num )

overrideprotectedvirtual

Implementation of function that computes the objective function for the current subset.

See also: compute_objective_function_without_penalty(const Target&,const int)

Developer\'s note

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 actual_compute_objective_function_without_penalty(), distributed_cache_enabled, stir::error(), proj_data_sptr, projector_pair_ptr, stir::setup_distributable_computation(), and zero_seg0_end_planes.

Referenced by actual_compute_objective_function_without_penalty().

◆ actual_add_multiplication_with_approximate_sub_Hessian_without_penalty()

template<typename TargetT>

Succeeded stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >::actual_add_multiplication_with_approximate_sub_Hessian_without_penalty	(	TargetT &	output,
		const TargetT &	input,
		const int	subset_num ) const

overrideprotectedvirtual

The Hessian (without penalty) is approximately given by:

$H_{jk} = - \sum_i P_{ij} h_i^{''}(y_i) P_{ik}$

where

$h_i(l) = y_i log (l) - l; h_i^{''}(y_i) = y_i / ((P \lambda)_i + a_i)^2;$

and $P_{ij}$ is the probability matrix.

This function uses the approximation of the hessian

$h_i^{''}(y_i) \approx -1/y_i$

Hence

$H_{jk} = - \sum_i P_{ij}(1/y_i) P_{ik}$

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.

Todo: Two work-arounds for the noisy estimate of the Hessian are listed below, but they are currently not implemented.

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:

$H_{jk} = \sum_i G_{ij}{1 \over n_i \mathrm{smooth}( n_i y_i)} G_{ik}$

where the probability matrix is factorised in a detection efficiency part (i.e. the normalisation factors ) times a geometric part:

$P_{ij} = {1 \over n_i } G_{ij}$

It has also been suggested to use $1 \over y_i+1$ (at least if the data are still Poisson).

Reimplemented from stir::GeneralisedObjectiveFunction< TargetT >.

References actual_add_multiplication_with_approximate_sub_Hessian_without_penalty(), stir::divide_and_truncate(), stir::error(), stir::RelatedViewgrams< elemT >::find_min(), stir::GeneralisedObjectiveFunction< TargetT >::get_num_subsets(), and stir::info().

Referenced by actual_add_multiplication_with_approximate_sub_Hessian_without_penalty().

◆ actual_accumulate_sub_Hessian_times_input_without_penalty()

template<typename TargetT>

Succeeded stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >::actual_accumulate_sub_Hessian_times_input_without_penalty	(	TargetT &	output,
		const TargetT &	current_image_estimate,
		const TargetT &	input,
		const int	subset_num ) const

overrideprotectedvirtual

The Hessian (without penalty) is approximately given by:

$H_{jk} = - \sum_i P_{ij} h_i^{''}(y_i) P_{ik}$

where

$h_i(l) = y_i log (l) - l; h_i^{''}(y_i) = - y_i / ((P \lambda)_i + a_i)^2;$

and $P_{ij}$ is the probability matrix.

Hence

$H_{jk} = - \sum_i P_{ij}(y_i / ((P \lambda)_i + a_i)^2) P_{ik}$

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 actual_accumulate_sub_Hessian_times_input_without_penalty(), stir::divide_and_truncate(), stir::error(), stir::GeneralisedObjectiveFunction< TargetT >::get_num_subsets(), and stir::info().

Referenced by actual_accumulate_sub_Hessian_times_input_without_penalty().

◆ set_defaults()

template<typename TargetT>

void stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >::set_defaults ( )

overrideprotectedvirtual

sets any default values

Has to be called by set_defaults in the leaf-class

Reimplemented from stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >.

References additive_projection_data_filename, distributed_cache_enabled, input_filename, max_segment_num_to_process, max_timing_pos_num_to_process, proj_data_sptr, projector_pair_ptr, stir::ParsingObject::set_defaults(), set_defaults(), use_tofsens, and zero_seg0_end_planes.

Referenced by PoissonLogLikelihoodWithLinearModelForMeanAndProjData(), and set_defaults().

◆ initialise_keymap()

template<typename TargetT>

void stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >::initialise_keymap ( )

overrideprotectedvirtual

sets keys for parsing

Has to be called by initialise_keymap in the leaf-class

Reimplemented from stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >.

References additive_projection_data_filename, stir::ParsingObject::initialise_keymap(), initialise_keymap(), input_filename, max_segment_num_to_process, projector_pair_ptr, use_tofsens, and zero_seg0_end_planes.

Referenced by initialise_keymap().

◆ post_processing()

template<typename TargetT>

bool stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >::post_processing ( )

overrideprotectedvirtual

checks values after parsing

Has to be called by post_processing in the leaf-class

Reimplemented from stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >.

References stir::error(), input_filename, stir::ParsingObject::post_processing(), post_processing(), proj_data_sptr, stir::ProjData::read_from_file(), and stir::warning().

Referenced by post_processing().

◆ actual_subsets_are_approximately_balanced()

template<typename TargetT>

bool stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >::actual_subsets_are_approximately_balanced ( std::string & warning_message ) const

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 actual_subsets_are_approximately_balanced(), stir::Array< num_dimensions, elemT >::fill(), stir::DataSymmetriesForViewSegmentNumbers::is_basic(), max_segment_num_to_process, stir::DataSymmetriesForViewSegmentNumbers::num_related_view_segment_numbers(), proj_data_sptr, and projector_pair_ptr.

Referenced by actual_subsets_are_approximately_balanced().

Member Data Documentation

◆ max_segment_num_to_process

template<typename TargetT>

int stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >::max_segment_num_to_process

protected

the maximum absolute ring difference number to use in the reconstruction

convention: if -1, use get_max_segment_num()

Referenced by actual_subsets_are_approximately_balanced(), add_subset_sensitivity(), initialise_keymap(), set_defaults(), and set_up_before_sensitivity().

◆ max_timing_pos_num_to_process

template<typename TargetT>

int stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >::max_timing_pos_num_to_process

protected

the maximum absolute time-of-flight bin number to use in the reconstruction

convention: if -1, use get_max_tof_pos_num()

Referenced by add_subset_sensitivity(), set_defaults(), and set_up_before_sensitivity().

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

/home/kris/devel/STIRdistrib/STIR/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.h
/home/kris/devel/STIRdistrib/STIR/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx

Public Member Functions

Static Public Attributes

Variables for STIR_MPI

Functions to set parameters

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ ~PoissonLogLikelihoodWithLinearModelForMeanAndProjData()

Member Function Documentation

◆ construct_target_ptr()

◆ set_num_subsets()

◆ set_additive_proj_data_sptr()

◆ set_normalisation_sptr()

◆ set_input_data()

◆ get_input_data()

◆ actual_compute_subset_gradient_without_penalty()

◆ get_exam_info_uptr_for_target()

◆ add_subset_sensitivity()

◆ set_up_before_sensitivity()

◆ actual_compute_objective_function_without_penalty()

◆ actual_add_multiplication_with_approximate_sub_Hessian_without_penalty()

◆ actual_accumulate_sub_Hessian_times_input_without_penalty()

◆ set_defaults()

◆ initialise_keymap()

◆ post_processing()

◆ actual_subsets_are_approximately_balanced()

Member Data Documentation

◆ max_segment_num_to_process

◆ max_timing_pos_num_to_process