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

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

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

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

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

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

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

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

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

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

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

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

Public Member Functions inherited from stir::ParsingObject
	ParsingObject (const ParsingObject &)

ParsingObject &	operator= (const ParsingObject &)

void	ask_parameters ()

virtual std::string	parameter_info ()

bool	parse (std::istream &f)

bool	parse (const char *const filename)

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
	Set defaults before parsing.

void	initialise_keymap () override
	Initialise all keywords.

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

Protected Member Functions inherited from stir::GeneralisedPrior< DiscretisedDensity< 3, elemT > >
void	set_defaults () override
	sets value for penalisation factor More...

void	initialise_keymap () override
	sets key for penalisation factor More...

virtual void	check (DiscretisedDensity< 3, elemT > const &current_estimate) const
	Check that the prior is ready to be used.

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

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

Protected Attributes inherited from stir::ParsingObject
KeyParser	parser

Additional Inherited Members
Static Public Member Functions inherited from stir::RegisteredObject< GeneralisedPrior< DiscretisedDensity< 3, elemT > > >
static GeneralisedPrior< DiscretisedDensity< 3, elemT > > *	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. More...

static GeneralisedPrior< DiscretisedDensity< 3, elemT > > *	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< GeneralisedPrior< DiscretisedDensity< 3, elemT > > >
typedef GeneralisedPrior< DiscretisedDensity< 3, elemT > > (	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< GeneralisedPrior< DiscretisedDensity< 3, elemT > > >
static RegistryType &	registry ()
	Static function returning the registry. More...

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 $r$ and $dr$ 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 ((x,y) = (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

◆ 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

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

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

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

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

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

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.

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

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

Public Member Functions

Protected Member Functions

Protected Attributes

Additional Inherited Members

Detailed Description

template<typename elemT, typename potentialT> class stir::GibbsPenalty< elemT, potentialT >

Member Function Documentation

◆ get_weights()

◆ set_weights()

◆ get_kappa_sptr()

◆ set_kappa_sptr()

◆ get_parsing_name()

◆ post_processing()

Member Data Documentation

◆ gradient_filename_prefix

template<typename elemT, typename potentialT>
class stir::GibbsPenalty< elemT, potentialT >