A class in the GeneralisedPrior hierarchy. This implements a Relative Difference prior. More...

#include "stir/recon_buildblock/RelativeDifferencePrior.h"

Inheritance diagram for stir::RelativeDifferencePrior< elemT >:

[legend]

Public Member Functions
	RelativeDifferencePrior ()
	Default constructor.

	RelativeDifferencePrior (const bool only_2D, float penalization_factor, float gamma, float epsilon)
	Constructs it explicitly.

double	compute_value (const DiscretisedDensity< 3, elemT > &current_image_estimate) override
	compute the value of the function

void	compute_gradient (DiscretisedDensity< 3, elemT > &prior_gradient, const DiscretisedDensity< 3, elemT > &current_image_estimate) override
	compute gradient

void	compute_Hessian (DiscretisedDensity< 3, elemT > &prior_Hessian_for_single_densel, const BasicCoordinate< 3, int > &coords, const DiscretisedDensity< 3, elemT > &current_image_estimate) const override
	This computes a single row of the Hessian.

void	add_multiplication_with_approximate_Hessian (DiscretisedDensity< 3, elemT > &output, const DiscretisedDensity< 3, elemT > &input) const override
	This should compute the multiplication of the Hessian with a vector and add it to output.

void	accumulate_Hessian_times_input (DiscretisedDensity< 3, elemT > &output, const DiscretisedDensity< 3, elemT > &current_estimate, const DiscretisedDensity< 3, elemT > &input) const override
	Compute the multiplication of the hessian of the prior multiplied by the input.

float	get_gamma () const
	get the gamma value used in RDP

void	set_gamma (float e)
	set the gamma value used in the RDP

float	get_epsilon () const
	get the epsilon value used in RDP

void	set_epsilon (float e)
	set the epsilon value used in the RDP

Array< 3, float >	get_weights () const
	get penalty weights for the neigbourhood

void	set_weights (const Array< 3, float > &)
	set penalty weights for the neigbourhood

shared_ptr< DiscretisedDensity< 3, elemT > >	get_kappa_sptr () const
	get current kappa image

void	set_kappa_sptr (const shared_ptr< DiscretisedDensity< 3, elemT > > &)
	set kappa image

virtual Succeeded	set_up (shared_ptr< const DiscretisedDensity< 3, elemT > > const &target_sptr) override
	Has to be called before using this object.

bool	is_convex () const override
	Indicates if the prior is a smooth convex function.

Public Member Functions inherited from stir::RegisteredParsingObject< RelativeDifferencePrior< elemT >, GeneralisedPrior< DiscretisedDensity< 3, elemT > >, GeneralisedPrior< DiscretisedDensity< 3, elemT > > >
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::GeneralisedPrior< DiscretisedDensity< 3, elemT > >
virtual double	compute_gradient_times_input (const DiscretisedDensity< 3, elemT > &input, const DiscretisedDensity< 3, elemT > &current_estimate)
	compute the dot product of the gradient of the log of the prior function at the current_estimate with input

virtual void	compute_Hessian_diagonal (DiscretisedDensity< 3, elemT > &Hessian_diagonal, const DiscretisedDensity< 3, elemT > &current_estimate) const
	This computes the diagonal of the Hessian of the log of the prior function at the current_estimate and stores it in Hessian_diagonal.

float	get_penalisation_factor () const

void	set_penalisation_factor (float new_penalisation_factor)

virtual double	compute_gradient_times_input (const DiscretisedDensity< 3, elemT > &input, const DiscretisedDensity< 3, elemT > &current_estimate)
	compute the dot product of the gradient of the log of the prior function at the current_estimate with input

virtual void	compute_Hessian (DiscretisedDensity< 3, elemT > &prior_Hessian_for_single_densel, const BasicCoordinate< 3, int > &coords, const DiscretisedDensity< 3, elemT > &current_image_estimate) const
	This computes a single row of the Hessian.

virtual void	compute_Hessian_diagonal (DiscretisedDensity< 3, elemT > &Hessian_diagonal, const DiscretisedDensity< 3, elemT > &current_estimate) const
	This computes the diagonal of the Hessian of the log of the prior function at the current_estimate and stores it in Hessian_diagonal.

virtual void	add_multiplication_with_approximate_Hessian (DiscretisedDensity< 3, elemT > &output, const DiscretisedDensity< 3, elemT > &input) const
	This should compute the multiplication of the Hessian with a vector and add it to output.

virtual void	accumulate_Hessian_times_input (DiscretisedDensity< 3, elemT > &output, const DiscretisedDensity< 3, elemT > &current_estimate, const DiscretisedDensity< 3, elemT > &input) const
	This should compute the multiplication of the Hessian with a vector and add it to output.

float	get_penalisation_factor () const

void	set_penalisation_factor (float new_penalisation_factor)

virtual Succeeded	set_up (shared_ptr< const DiscretisedDensity< 3, elemT > > const &target_sptr)
	Has to be called before using this object.

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

Protected Member Functions
void	check (DiscretisedDensity< 3, elemT > const &current_image_estimate) const override
	Check that the prior is ready to be used.

void	set_defaults () override
	sets value for penalisation factor

void	initialise_keymap () override
	sets key for penalisation factor

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

double	value (const elemT x_j, const elemT x_k) const
	The value and partial derivatives of the Relative Difference Prior.

elemT	derivative_10 (const elemT x, const elemT y) const

elemT	derivative_20 (const elemT x_j, const elemT x_k) const

elemT	derivative_11 (const elemT x_j, const elemT x_k) const

Protected Member Functions inherited from stir::GeneralisedPrior< DiscretisedDensity< 3, elemT > >
virtual void	check (DiscretisedDensity< 3, elemT > const &current_estimate) const
	Check that the prior is ready to be used.

virtual void	set_key_values ()
	This will be called before parsing or parameter_info is called.

Protected Attributes
float	gamma
	Create variable gamma for Relative Difference Penalty.

float	epsilon
	Create variable epsilon for Relative Difference Penalty.

bool	only_2D
	can be set during parsing to restrict the weights to the 2D case

std::string	gradient_filename_prefix
	filename prefix for outputing the gradient whenever compute_gradient() is called.

Array< 3, float >	weights
	penalty weights

std::string	kappa_filename
	Filename for the $\kappa$ image that will be read by post_processing()

shared_ptr< DiscretisedDensity< 3, elemT > >	kappa_ptr

Protected Attributes inherited from stir::GeneralisedPrior< DiscretisedDensity< 3, elemT > >
float	penalisation_factor

bool	_already_set_up

float	penalisation_factor

bool	_already_set_up

KeyParser	parser

Additional Inherited Members
Static Public Member Functions inherited from stir::RegisteredParsingObject< RelativeDifferencePrior< elemT >, GeneralisedPrior< DiscretisedDensity< 3, elemT > >, GeneralisedPrior< DiscretisedDensity< 3, elemT > > >
static GeneralisedPrior< DiscretisedDensity< 3, elemT > > *	read_from_stream (std::istream *)
	Construct a new object (of type Derived) by parsing the istream.

static GeneralisedPrior< DiscretisedDensity< 3, elemT > > *	read_from_stream (std::istream *)
	Construct a new object (of type Derived) by parsing the istream.

Static Public Member Functions inherited from stir::RegisteredObject< Root >
static Root *	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 Root *	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< Root >
typedef Root (	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< Root >
static RegistryType &	registry ()
	Static function returning the registry.

Detailed Description

template<typename elemT>
class stir::RelativeDifferencePrior< elemT >

A class in the GeneralisedPrior hierarchy. This implements a Relative Difference prior.

The value of the prior is computed as follows:

$f= \sum_{r,dr} \frac{w_{dr}}{2} \frac{(\lambda_r - \lambda_{r+dr})^2}{(\lambda_r+ \lambda_{r+dr} + \gamma |\lambda_r - \lambda_{r+dr}| + \epsilon)} * \kappa_r * \kappa_{r+dr}$

where $\lambda$ is the image and and are indices and the sum is over the neighbourhood where the weights $w_{dr}$ are non-zero. $\gamma$ is a smoothing scalar term and the $\epsilon$ is a small non-negative value included to prevent division by zero. Please note that the RDP is only well defined for non-negative voxel values. For more details, see: J. Nuyts, D. Beque, P. Dupont, and L. Mortelmans, "A Concave Prior Penalizing Relative Differences for Maximum-a-Posteriori Reconstruction in Emission Tomography," vol. 49, no. 1, pp. 56-60, 2002.

If $\epsilon=0$ , we attempt to resolve 0/0 at $\lambda_r = \lambda_{r+dr}=0$ by using the limit. Note that the Hessian explodes to infinity when both voxel values approach 0, and we currently return INFINITY. Also, as the RDP is not differentiable at this point, we have chosen to return 0 for the gradient (such that a zero background is not modified).

Warning: the default value for $\epsilon$ is zero, which can be problematic for gradient-based algorithms.

The $\kappa$ image can be used to have spatially-varying penalties such as in Jeff Fessler's papers. It should have identical dimensions to the image for which the penalty is computed. If $\kappa$ is not set, this class will effectively use 1 for all $\kappa$ 's.

By default, a 3x3 or 3x3x3 neighbourhood is used where the weights are set to x-voxel_size divided by the Euclidean distance between the points.

The prior computation excludes voxel-pairs where one voxel is outside the volume. This is effectively the same as extending the volume by replicating the edges (which is different from zero boundary conditions).

Parsing

These are the keywords that can be used in addition to the ones in GeneralPrior.

  Relative Difference Prior Parameters:=
  ; next defaults to 0, set to 1 for 2D inverse Euclidean weights, 0 for 3D
  only 2D:= 0
  ; next can be used to set weights explicitly. Needs to be a 3D array (of floats).
  ' value of only_2D is ignored
  ; following example uses 2D 'nearest neighbour' penalty
  ; weights:={{{0,1,0},{1,0,1},{0,1,0}}}
  ; gamma value :=
  ; epsilon value :=
  ; see class documentation for more info
  ; use next parameter to specify an image with penalisation factors (a la Fessler)
  ; kappa filename:=
  ; use next parameter to get gradient images at every subiteration
  ; see class documentation
  gradient filename prefix:=
  END Relative Difference Prior Parameters:=

Member Function Documentation

◆ compute_value()

template<typename elemT>

double stir::RelativeDifferencePrior< elemT >::compute_value ( const DiscretisedDensity< 3, elemT > & current_image_estimate )

overridevirtual

compute the value of the function

Implements stir::GeneralisedPrior< DiscretisedDensity< 3, elemT > >.

References check(), epsilon, stir::DiscretisedDensityOnCartesianGrid< num_dimensions, elemT >::get_grid_spacing(), stir::VectorWithOffset< T >::get_length(), stir::VectorWithOffset< T >::get_max_index(), stir::VectorWithOffset< T >::get_min_index(), value(), and weights.

◆ compute_gradient()

template<typename elemT>

void stir::RelativeDifferencePrior< elemT >::compute_gradient	(	DiscretisedDensity< 3, elemT > &	prior_gradient,
		const DiscretisedDensity< 3, elemT > &	current_image_estimate )

overridevirtual

compute gradient

Implements stir::GeneralisedPrior< DiscretisedDensity< 3, elemT > >.

References check(), stir::Array< num_dimensions, elemT >::fill(), stir::Array< num_dimensions, elemT >::find_max(), stir::Array< num_dimensions, elemT >::find_min(), stir::DiscretisedDensityOnCartesianGrid< num_dimensions, elemT >::get_grid_spacing(), stir::VectorWithOffset< T >::get_length(), stir::VectorWithOffset< T >::get_max_index(), stir::VectorWithOffset< T >::get_min_index(), gradient_filename_prefix, stir::DiscretisedDensity< num_dimensions, elemT >::has_same_characteristics(), stir::info(), weights, and stir::write_to_file().

◆ compute_Hessian()

template<typename elemT>

void stir::RelativeDifferencePrior< elemT >::compute_Hessian	(	DiscretisedDensity< 3, elemT > &	prior_Hessian_for_single_densel,
		const BasicCoordinate< 3, int > &	coords,
		const DiscretisedDensity< 3, elemT > &	current_image_estimate ) const

overridevirtual

This computes a single row of the Hessian.

Default implementation just call error(). This function needs to be overridden by the derived class.

The method computes a row (i.e. at a densel/voxel, indicated by coords) of the Hessian at current_estimate. Note that a row corresponds to an object of DataT. The method (as implemented in derived classes) should store the result in prior_Hessian_for_single_densel.

Reimplemented from stir::GeneralisedPrior< DiscretisedDensity< 3, elemT > >.

References check(), stir::Array< num_dimensions, elemT >::fill(), stir::DiscretisedDensityOnCartesianGrid< num_dimensions, elemT >::get_grid_spacing(), stir::VectorWithOffset< T >::get_max_index(), stir::VectorWithOffset< T >::get_min_index(), stir::DiscretisedDensity< num_dimensions, elemT >::has_same_characteristics(), and weights.

◆ add_multiplication_with_approximate_Hessian()

template<typename elemT>

void stir::RelativeDifferencePrior< elemT >::add_multiplication_with_approximate_Hessian	(	DiscretisedDensity< 3, elemT > &	output,
		const DiscretisedDensity< 3, elemT > &	input ) const

overridevirtual

This should compute the multiplication of the Hessian with a vector and add it to output.

Default implementation just call error(). This function needs to be overridden by the derived class. This method assumes that the hessian of the prior is 1 and hence the function quadratic. Instead, accumulate_Hessian_times_input() should be used. This method remains for backwards comparability.

Warning: The derived class should accumulate in output.

Reimplemented from stir::GeneralisedPrior< DiscretisedDensity< 3, elemT > >.

References stir::error().

◆ accumulate_Hessian_times_input()

template<typename elemT>

void stir::RelativeDifferencePrior< elemT >::accumulate_Hessian_times_input	(	DiscretisedDensity< 3, elemT > &	output,
		const DiscretisedDensity< 3, elemT > &	current_estimate,
		const DiscretisedDensity< 3, elemT > &	input ) const

overridevirtual

Compute the multiplication of the hessian of the prior multiplied by the input.

Reimplemented from stir::GeneralisedPrior< DiscretisedDensity< 3, elemT > >.

References stir::DiscretisedDensityOnCartesianGrid< num_dimensions, elemT >::get_grid_spacing(), stir::VectorWithOffset< T >::get_max_index(), stir::VectorWithOffset< T >::get_min_index(), stir::DiscretisedDensity< num_dimensions, elemT >::has_same_characteristics(), and weights.

◆ get_weights()

template<typename elemT>

Array< 3, float > stir::RelativeDifferencePrior< elemT >::get_weights ( ) const

get penalty weights for the neigbourhood

get penalty weights for the neighbourhood

References weights.

◆ set_weights()

template<typename elemT>

void stir::RelativeDifferencePrior< elemT >::set_weights ( const Array< 3, float > & w )

set penalty weights for the neigbourhood

set penalty weights for the neighbourhood

References weights.

◆ get_kappa_sptr()

template<typename elemT>

shared_ptr< DiscretisedDensity< 3, elemT > > stir::RelativeDifferencePrior< elemT >::get_kappa_sptr ( ) const

get current kappa image

Warning: As this function returns a shared_ptr, this is dangerous. You should not modify the image by manipulating the image refered to by this pointer. Unpredictable results will occur.; As this function returns a shared_ptr, this is dangerous. You should not modify the image by manipulating the image referred to by this pointer. Unpredictable results will occur.

◆ set_up()

template<typename elemT>

Succeeded stir::RelativeDifferencePrior< elemT >::set_up ( shared_ptr< const DiscretisedDensity< 3, elemT > > const & target_sptr )

overridevirtual

Has to be called before using this object.

Reimplemented from stir::GeneralisedPrior< DiscretisedDensity< 3, elemT > >.

Reimplemented in stir::CudaRelativeDifferencePrior< elemT >.

References stir::GeneralisedPrior< DiscretisedDensity< 3, elemT > >::set_up(), and set_up().

Referenced by set_up().

◆ is_convex()

template<typename elemT>

bool stir::RelativeDifferencePrior< elemT >::is_convex ( ) const

overridevirtual

Indicates if the prior is a smooth convex function.

If true, the prior is expected to have 0th, 1st and 2nd order behaviour implemented.

Implements stir::GeneralisedPrior< DiscretisedDensity< 3, elemT > >.

References is_convex().

Referenced by is_convex().

◆ check()

template<typename elemT>

void stir::RelativeDifferencePrior< elemT >::check ( DiscretisedDensity< 3, elemT > const & current_image_estimate ) const

overrideprotectedvirtual

Check that the prior is ready to be used.

Reimplemented from stir::GeneralisedPrior< DiscretisedDensity< 3, elemT > >.

References stir::GeneralisedPrior< DiscretisedDensity< 3, elemT > >::check(), check(), stir::error(), and registered_name.

Referenced by check(), compute_gradient(), compute_Hessian(), and compute_value().

◆ set_defaults()

template<typename elemT>

void stir::RelativeDifferencePrior< elemT >::set_defaults ( )

overrideprotectedvirtual

sets value for penalisation factor

Has to be called by set_defaults in the leaf-class

Reimplemented from stir::GeneralisedPrior< DiscretisedDensity< 3, elemT > >.

References epsilon, gamma, only_2D, stir::ParsingObject::set_defaults(), set_defaults(), and weights.

Referenced by RelativeDifferencePrior(), RelativeDifferencePrior(), and set_defaults().

◆ initialise_keymap()

template<typename elemT>

void stir::RelativeDifferencePrior< elemT >::initialise_keymap ( )

overrideprotectedvirtual

sets key for penalisation factor

Has to be called by initialise_keymap in the leaf-class

Reimplemented from stir::GeneralisedPrior< DiscretisedDensity< 3, elemT > >.

References epsilon, gamma, gradient_filename_prefix, stir::ParsingObject::initialise_keymap(), initialise_keymap(), kappa_filename, only_2D, and weights.

Referenced by initialise_keymap().

◆ post_processing()

template<typename elemT>

bool stir::RelativeDifferencePrior< elemT >::post_processing ( )

overrideprotectedvirtual

This will be called at the end of the parsing.

Returns: false if everything OK, true if not

Reimplemented from stir::ParsingObject.

References kappa_filename, stir::ParsingObject::post_processing(), post_processing(), stir::read_from_file(), stir::warning(), and weights.

Referenced by post_processing().

◆ value()

template<typename elemT>

double stir::RelativeDifferencePrior< elemT >::value	(	const elemT	x_j,
		const elemT	x_k ) const

protected

The value and partial derivatives of the Relative Difference Prior.

value refers to its value derivative_10 refers to the derivative w.r.t. x_j only derivative_20 refers to the second derivative w.r.t. x_j only (i.e. diagonal elements of the Hessian) derivative_11 refers to the second derivative w.r.t. x_j and x_k (i.e. off-diagonal elements of the Hessian)

See J. Nuyts, et al., 2002, Equation 7 etc, but here an epsilon is added.

In the instance x_j, x_k and epsilon equal 0.0, these functions get ill-defined due to 0/0 computations. We currently return 0 for the value, 0 for derivative_10, and INFINITY for the second order derivatives. These follow by taking the limit, i.e. by assuming continuity. Note however that when epsilon=0, derivative_10(x,0) limits to $1/(1+\gamma)$ , while derivative_10(x,x) limits to .

Parameters

x_j	is the target voxel.
x_k	is the voxel in the neighbourhood.

Returns: the second order partial derivatives of the Relative Difference Prior

References epsilon, gamma, and stir::square().

Referenced by compute_value().

Member Data Documentation

◆ gradient_filename_prefix

template<typename elemT>

std::string stir::RelativeDifferencePrior< elemT >::gradient_filename_prefix

protected

filename prefix for outputing the gradient whenever compute_gradient() is called.

An internal counter is used to keep track of the number of times the gradient is computed. The filename will be constructed by concatenating gradient_filename_prefix and the counter.

Referenced by compute_gradient(), and initialise_keymap().

◆ weights

template<typename elemT>

Array<3, float> stir::RelativeDifferencePrior< elemT >::weights

mutableprotected

penalty weights

Todo: This member is mutable at present because some const functions initialise it. That initialisation should be moved to a new set_up() function.

Referenced by accumulate_Hessian_times_input(), compute_gradient(), compute_Hessian(), compute_value(), get_weights(), initialise_keymap(), post_processing(), set_defaults(), and set_weights().

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

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

Public Member Functions

Static Public Attributes

Protected Member Functions

Protected Attributes

Additional Inherited Members

Detailed Description

Member Function Documentation

◆ compute_value()

◆ compute_gradient()

◆ compute_Hessian()

◆ add_multiplication_with_approximate_Hessian()

◆ accumulate_Hessian_times_input()

◆ get_weights()

◆ set_weights()

◆ get_kappa_sptr()

◆ set_up()

◆ is_convex()

◆ check()

◆ set_defaults()

◆ initialise_keymap()

◆ post_processing()

◆ value()

Member Data Documentation

◆ gradient_filename_prefix

◆ weights