A class in the GeneralisedPrior hierarchy. This implements the anatomical penalty function, Parallel Level Sets (PLS), also sometimes called Directional TV, proposed by Matthias J. Ehrhardt et. al in "PET Reconstruction With an Anatomical MRI Prior Using Parallel Level Sets", IEEE Trans. med. Imag., vol. 35, no. 9, Sept. 2016. https://doi.org/10.1109/TMI.2016.2549601 Note that PLS becomes smoothed TV when an uniform anatomical image is provided. More...

#include "stir/recon_buildblock/PLSPrior.h"

Inheritance diagram for stir::PLSPrior< elemT >:

[legend]

Public Member Functions
	PLSPrior ()
	Default constructor.

	PLSPrior (const bool only_2D, float penalization_factor)
	Constructs it explicitly.

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.

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

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

shared_ptr< const DiscretisedDensity< 3, elemT > >	get_anatomical_grad_sptr (int direction) const

shared_ptr< const DiscretisedDensity< 3, elemT > >	get_norm_sptr () const

double	get_eta () const
	get eta and alpha parameters

double	get_alpha () const

void	set_eta (const double)
	set eta parameter

void	set_alpha (const double)
	set alpha parameter

void	set_anatomical_image_sptr (const shared_ptr< const DiscretisedDensity< 3, elemT > > &)
	set anatomical pointer

shared_ptr< const DiscretisedDensity< 3, elemT > >	get_anatomical_image_sptr () const
	get anatomical pointer

void	set_anatomical_filename (const std::string &filename)
	Set anatomical filename.

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

void	set_kappa_filename (const std::string &filename)
	Set kappa filename.

void	set_only_2D (const bool arg)
	Set only 2D.

bool	get_only_2D ()
	Get only 2D.

Public Member Functions inherited from stir::RegisteredParsingObject< PLSPrior< 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 (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 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	set_defaults () override
	sets value for penalisation factor

void	initialise_keymap () override
	sets key for penalisation factor

bool	post_processing () override
	the parsing will only override any exixting kappa-image or anatomical-image if the relevant keyword is present

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

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
bool	only_2D
	can be set during parsing to restrict the gradient calculation 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()

std::string	anatomical_filename

double	eta

double	alpha

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< PLSPrior< 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::PLSPrior< elemT >

A class in the GeneralisedPrior hierarchy. This implements the anatomical penalty function, Parallel Level Sets (PLS), also sometimes called Directional TV, proposed by Matthias J. Ehrhardt et. al in "PET Reconstruction With an Anatomical MRI Prior Using Parallel Level Sets", IEEE Trans. med. Imag., vol. 35, no. 9, Sept. 2016. https://doi.org/10.1109/TMI.2016.2549601 Note that PLS becomes smoothed TV when an uniform anatomical image is provided.

The prior has 2 parameters alpha and eta. It is computed for an image as

$\phi(f) = \sqrt{\alpha^2 + |\nabla f|^2 - {\langle\nabla f,\xi\rangle}^2}$

where is the PET image, $\nabla$ is the finite difference operator (not taking voxel-sizes into account) and $\xi$ is the normalised gradient of the anatomical image calculated as follows:

$\xi = \frac{\nabla v}{\sqrt{|\nabla v|^2 + \eta^2}}$

with the anatomical image, $\alpha$ controls the edge-preservation property of PLS, and depends on the scale of the emission image, and $\eta$ avoids division by zero, and depends on the scale of the anatomical image.

A $\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.

Parsing

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

PLS Prior Parameters:=
; next defaults to 0, set to 1 for 2D images, 0 for 3D
only 2D:= 0


anatomical_filename:= file.hv; Image that provides anatomical information (i.e. CT or MR image). The
                 dimension should be the same as that of the emission image.

eta :=    ; A parameter for preventing the division by zero problem. The value depends
                 on the scale of the anatomical image.

alpha :=   ; A parameter that controls the edge-preservation property of PLS. The value
                 depends on the scale of the emission image.

; use next parameter to specify an image with penalisation factors (a la Fessler)
; see class documentation for more info
; kappa filename:=
; use next parameter to get gradient images at every subiteration
; see class documentation
gradient filename prefix:=
END PLS Prior Parameters:=

Member Function Documentation

◆ set_up()

template<typename elemT>

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

overridevirtual

Has to be called before using this object.

Todo: set the anatomical image to zero if not defined

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

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

Referenced by set_up().

◆ is_convex()

template<typename elemT>

bool stir::PLSPrior< 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().

◆ compute_value()

template<typename elemT>

double stir::PLSPrior< 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(), stir::VectorWithOffset< T >::get_max_index(), stir::VectorWithOffset< T >::get_min_index(), and only_2D.

◆ compute_gradient()

template<typename elemT>

void stir::PLSPrior< 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::VectorWithOffset< T >::get_max_index(), stir::VectorWithOffset< T >::get_min_index(), gradient_filename_prefix, stir::info(), only_2D, and stir::write_to_file().

◆ get_kappa_sptr()

template<typename elemT>

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

◆ set_defaults()

template<typename elemT>

void stir::PLSPrior< 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 only_2D, stir::ParsingObject::set_defaults(), and set_defaults().

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

◆ initialise_keymap()

template<typename elemT>

void stir::PLSPrior< 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 gradient_filename_prefix, stir::ParsingObject::initialise_keymap(), initialise_keymap(), kappa_filename, and only_2D.

Referenced by initialise_keymap().

◆ post_processing()

template<typename elemT>

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

overrideprotectedvirtual

the parsing will only override any exixting kappa-image or anatomical-image if the relevant keyword is present

Reimplemented from stir::ParsingObject.

References kappa_filename, stir::ParsingObject::post_processing(), post_processing(), set_anatomical_filename(), and set_kappa_filename().

Referenced by post_processing().

◆ check()

template<typename elemT>

void stir::PLSPrior< 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(), and compute_value().

Member Data Documentation

◆ gradient_filename_prefix

template<typename elemT>

std::string stir::PLSPrior< 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().

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

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

Public Member Functions

Static Public Attributes

Protected Member Functions

Protected Attributes

Additional Inherited Members

Detailed Description

Member Function Documentation

◆ set_up()

◆ is_convex()

◆ compute_value()

◆ compute_gradient()

◆ get_kappa_sptr()

◆ set_defaults()

◆ initialise_keymap()

◆ post_processing()

◆ check()

Member Data Documentation

◆ gradient_filename_prefix