A base class for Gibbs type penalties in the GeneralisedPrior hierarchy. More...

#include "stir/recon_buildblock/GibbsPenalty.h"

Inheritance diagram for stir::GibbsPenalty< elemT, potentialT >:

[legend]

Public Member Functions
	GibbsPenalty ()
	Default constructor.

	GibbsPenalty (const bool only_2D, float penalization_factor)
	Explicit Constructor with 2D/3D option and penalization factor.

double	compute_value (const DiscretisedDensity< 3, elemT > &current_image_estimate) override
	Compute the value of the prior for the current image estimate.

void	compute_gradient (DiscretisedDensity< 3, elemT > &prior_gradient, const DiscretisedDensity< 3, elemT > &current_image_estimate) override
	Compute the gradient of the prior for the current image estimate.

double	compute_gradient_times_input (const DiscretisedDensity< 3, elemT > &input, const DiscretisedDensity< 3, elemT > &current_image_estimate) override
	Compute the dot product of the prior gradient and an input image.

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
	Compute the Hessian row of the prior at (coords).

void	compute_Hessian_diagonal (DiscretisedDensity< 3, elemT > &Hessian_diagonal, const DiscretisedDensity< 3, elemT > &current_estimate) const override
	Compute the diagonal of the Hessian matrix.

void	accumulate_Hessian_times_input (DiscretisedDensity< 3, elemT > &output, const DiscretisedDensity< 3, elemT > &current_estimate, const DiscretisedDensity< 3, elemT > &input) const override
	Accumulate Hessian times input image into output.

const Array< 3, float > &	get_weights () const
	Get the current weights array.

virtual void	set_weights (const Array< 3, float > &)
	Set the weights array for the prior.

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

virtual void	set_kappa_sptr (const shared_ptr< const DiscretisedDensity< 3, elemT > > &)
	Set the kappa image (spatially-varying penalty factors).

Succeeded	set_up (shared_ptr< const DiscretisedDensity< 3, elemT > > const &target_sptr) override
	Set up the prior for a target image. Must be called before use.

virtual std::string	get_parsing_name () const
	Getter method to retrieve the parsing name.

bool	is_convex () const override
	Return whether the prior is convex or not.

Public Member Functions inherited from stir::GeneralisedPrior< DiscretisedDensity< 3, elemT > >
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.

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::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 &)

ParsingObject &	operator= (const ParsingObject &)

bool	parse (std::istream &f)

bool	parse (const char *const filename)

void	ask_parameters ()

virtual std::string	parameter_info ()

Protected Member Functions
void	compute_default_weights (const CartesianCoordinate3D< float > &grid_spacing, bool only_2D)
	Compute default weights for the prior.

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.

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
Array< 3, float >	weights
	The weights for the neighbourhood.

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.

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

potentialT	potential
	Gibbs Potential Function.

shared_ptr< const DiscretisedDensity< 3, elemT > >	kappa_ptr
	The kappa image (spatially-varying penalty factors).

Image and weight boundary indices
CartesianCoordinate3D< int >	image_dim
	Image dimensions.

CartesianCoordinate3D< int >	image_max_indices
	Maximum image indices.

CartesianCoordinate3D< int >	image_min_indices
	Minimum image indices.

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::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, typename potentialT>
class stir::GibbsPenalty< elemT, potentialT >

A base class for Gibbs type penalties in the GeneralisedPrior hierarchy.

The prior is computed as follows:

$f = \sum_{r,dr} w_{r,dr} \phi(\lambda_r , \lambda_{r+dr}) * \kappa_r * \kappa_{r+dr}$

with gradient given by:

$g_r = 2 * \sum_{dr} w_{r,dr} \phi'(\lambda_r , \lambda_{r+dr}) * \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.

The $\phi$ function is the potential function, which is provided via the template parameter PotentialFun. The potential function needs to be symmetric ( $\phi(x,y) = \phi(y,x)$ ). Currently the potential function is implemented in the header of the derived classes (check GibbsQuadraticPrior.h or GibbsRelativeDifferencePrior.h to see examples).

Potential Function Requirements: The potential function class (PotentialT) must implement the following methods:

value(elemT val_center, elemT val_neigh, int z, int y, int x) - Returns the value of the potential function for the two voxel values
derivative_10(elemT val_center, elemT val_neigh, int z, int y, int x) - First derivative with respect to the center voxel
derivative_11(elemT val_center, elemT val_neigh, int z, int y, int x) - Second mixed derivative
derivative_20(elemT val_center, elemT val_neigh, int z, int y, int x) - Second derivative with respect to the center voxel
static bool is_convex() - Returns whether the potential function is convex
void initialise_keymap(KeyParser& parser) - Sets up parsing for any potential-specific parameters

These methods should be declared with __host__ __device__ qualifiers, except for is_convex() and initialise_keymap() which are only used on the host. The coordinate parameters (z, y, x) may be used by the potential function for position-dependent behavior. Even if the potential has no parameters to parse, the initialise_keymap method must be implemented (possibly with an empty body).

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 neigbourhood is used where the weights are set to x-voxel_size divided by the Euclidean distance between the points. Custom weights can be set using the method set_weights, the general form of the weights is NxNxN with N odd.

Warning: Currently only symmetric weights are supported (w_{i,j} = w_{j,i}).

Parsing

These are the keywords that can be used in all priors derived from GibbsPrior:

; 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).
; weights:={{{0,1,0},{1,0,1},{0,1,0}}}
; 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
gradient filename prefix:=

Member Function Documentation

◆ compute_value()

template<typename elemT, typename PotentialT>

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

overridevirtual

Compute the value of the prior for the current image estimate.

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

References check(), stir::error(), image_max_indices, image_min_indices, kappa_ptr, potential, and weights.

◆ compute_gradient()

template<typename elemT, typename PotentialT>

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

overridevirtual

Compute the gradient of the prior for the current image estimate.

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

References check(), stir::error(), stir::Array< num_dimensions, elemT >::fill(), stir::DiscretisedDensity< num_dimensions, elemT >::has_same_characteristics(), image_max_indices, image_min_indices, kappa_ptr, potential, and weights.

◆ compute_gradient_times_input()

template<typename elemT, typename PotentialT>

double stir::GibbsPenalty< elemT, PotentialT >::compute_gradient_times_input	(	const DiscretisedDensity< 3, elemT > &	input,
		const DiscretisedDensity< 3, elemT > &	current_image_estimate )

overridevirtual

Compute the dot product of the prior gradient and an input image.

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

References check(), stir::error(), stir::DiscretisedDensity< num_dimensions, elemT >::has_same_characteristics(), image_max_indices, image_min_indices, kappa_ptr, potential, and weights.

◆ compute_Hessian()

template<typename elemT, typename PotentialT>

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

overridevirtual

Compute the Hessian row of the prior at (coords).

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

References check(), stir::error(), stir::Array< num_dimensions, elemT >::fill(), stir::DiscretisedDensity< num_dimensions, elemT >::has_same_characteristics(), image_max_indices, image_min_indices, kappa_ptr, potential, and weights.

◆ compute_Hessian_diagonal()

template<typename elemT, typename PotentialT>

void stir::GibbsPenalty< elemT, PotentialT >::compute_Hessian_diagonal	(	DiscretisedDensity< 3, elemT > &	Hessian_diagonal,
		const DiscretisedDensity< 3, elemT > &	current_estimate ) const

overridevirtual

Compute the diagonal of the Hessian matrix.

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

References check(), stir::error(), stir::Array< num_dimensions, elemT >::fill(), stir::DiscretisedDensity< num_dimensions, elemT >::has_same_characteristics(), image_max_indices, image_min_indices, kappa_ptr, potential, and weights.

◆ accumulate_Hessian_times_input()

template<typename elemT, typename PotentialT>

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

overridevirtual

Accumulate Hessian times input image into output.

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

References check(), stir::error(), stir::DiscretisedDensity< num_dimensions, elemT >::has_same_characteristics(), image_max_indices, image_min_indices, kappa_ptr, potential, and weights.

◆ get_weights()

template<typename elemT, typename PotentialT>

const Array< 3, float > & stir::GibbsPenalty< elemT, PotentialT >::get_weights ( ) const

Get the current weights array.

get penalty weights for the neigbourhood

References weights.

◆ set_weights()

template<typename elemT, typename PotentialT>

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

virtual

Set the weights array for the prior.

set penalty weights for the neigbourhood

Reimplemented in stir::CudaGibbsPenalty< elemT, PotentialT >.

References weights.

◆ get_kappa_sptr()

template<typename elemT, typename PotentialT>

shared_ptr< const DiscretisedDensity< 3, elemT > > stir::GibbsPenalty< elemT, PotentialT >::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.

References kappa_ptr.

◆ set_kappa_sptr()

template<typename elemT, typename PotentialT>

void stir::GibbsPenalty< elemT, PotentialT >::set_kappa_sptr ( const shared_ptr< const DiscretisedDensity< 3, elemT > > & k )

virtual

Set the kappa image (spatially-varying penalty factors).

set kappa image

Reimplemented in stir::CudaGibbsPenalty< elemT, PotentialT >.

References kappa_ptr.

◆ set_up()

template<typename elemT, typename potentialT>

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

overridevirtual

Set up the prior for a target image. Must be called before use.

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

References set_up().

Referenced by set_up().

◆ get_parsing_name()

template<typename elemT, typename potentialT>

std::string stir::GibbsPenalty< elemT, potentialT >::get_parsing_name ( ) const

virtual

Getter method to retrieve the parsing name.

Default implementation just uses get_registered_name(), and appends " Parameters".

References get_parsing_name(), and stir::RegisteredObjectBase::get_registered_name().

Referenced by get_parsing_name(), and initialise_keymap().

◆ is_convex()

template<typename elemT, typename potentialT>

bool stir::GibbsPenalty< elemT, potentialT >::is_convex ( ) const

overridevirtual

Return whether the prior is convex or not.

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

References is_convex().

Referenced by is_convex().

◆ check()

template<typename elemT, typename potentialT>

void stir::GibbsPenalty< elemT, potentialT >::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 check().

Referenced by accumulate_Hessian_times_input(), check(), compute_gradient(), compute_gradient_times_input(), compute_Hessian(), compute_Hessian_diagonal(), and compute_value().

◆ set_defaults()

template<typename elemT, typename potentialT>

void stir::GibbsPenalty< elemT, potentialT >::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 set_defaults().

Referenced by GibbsPenalty(), and set_defaults().

◆ initialise_keymap()

template<typename elemT, typename potentialT>

void stir::GibbsPenalty< elemT, potentialT >::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 get_parsing_name(), gradient_filename_prefix, stir::GeneralisedPrior< DiscretisedDensity< 3, elemT > >::initialise_keymap(), initialise_keymap(), kappa_filename, only_2D, potential, and weights.

Referenced by initialise_keymap().

◆ post_processing()

template<typename elemT, typename potentialT>

bool stir::GibbsPenalty< elemT, potentialT >::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, kappa_ptr, post_processing(), stir::ParsingObject::post_processing(), stir::read_from_file(), stir::warning(), and weights.

Referenced by post_processing().

Member Data Documentation

◆ gradient_filename_prefix

template<typename elemT, typename potentialT>

std::string stir::GibbsPenalty< elemT, potentialT >::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 initialise_keymap().

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

/home/kris/devel/STIRdistrib/STIR/src/include/stir/recon_buildblock/GibbsPenalty.h
/home/kris/devel/STIRdistrib/STIR/src/include/stir/recon_buildblock/GibbsPenalty.inl

Public Member Functions

Protected Member Functions

Protected Attributes

Additional Inherited Members

Detailed Description

Member Function Documentation

◆ compute_value()

◆ compute_gradient()

◆ compute_gradient_times_input()

◆ compute_Hessian()

◆ compute_Hessian_diagonal()

◆ accumulate_Hessian_times_input()

◆ get_weights()

◆ set_weights()

◆ get_kappa_sptr()

◆ set_kappa_sptr()

◆ set_up()

◆ get_parsing_name()

◆ is_convex()

◆ check()

◆ set_defaults()

◆ initialise_keymap()

◆ post_processing()

Member Data Documentation

◆ gradient_filename_prefix