STIR  6.2.0
Public Member Functions | Static Public Attributes | Protected Member Functions | Protected Attributes | List of all members

A class in the GeneralisedPrior hierarchy. This implements a quadratic Gibbs prior. More...

#include "stir/recon_buildblock/QuadraticPrior.h"

Inheritance diagram for stir::QuadraticPrior< elemT >:
Inheritance graph
[legend]

Public Member Functions

 QuadraticPrior ()
 Default constructor.
 
 QuadraticPrior (const bool only_2D, float penalization_factor)
 Constructs it explicitly.
 
bool parabolic_surrogate_curvature_depends_on_argument () const override
 A function that allows skipping some computations if the curvature is independent of the current_estimate. More...
 
bool is_convex () const override
 Indicates if the prior is a smooth convex function. More...
 
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 parabolic_surrogate_curvature (DiscretisedDensity< 3, elemT > &parabolic_surrogate_curvature, const DiscretisedDensity< 3, elemT > &current_image_estimate) override
 compute the parabolic surrogate for the prior More...
 
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
 
void add_multiplication_with_approximate_Hessian (DiscretisedDensity< 3, elemT > &output, const DiscretisedDensity< 3, elemT > &input) const override
 Call accumulate_Hessian_times_input.
 
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. For the quadratic function, the hessian of the prior is 1. Therefore this will return the weights multiplied by the input.
 
Array< 3, float > get_weights () const
 get penalty weights for the neighbourhood More...
 
void set_weights (const Array< 3, float > &)
 set penalty weights for the neighbourhood More...
 
shared_ptr< const DiscretisedDensity< 3, elemT > > get_kappa_sptr () const
 get current kappa image More...
 
void set_kappa_sptr (const shared_ptr< const DiscretisedDensity< 3, elemT >> &)
 set kappa image
 
Succeeded set_up (shared_ptr< const DiscretisedDensity< 3, elemT >> const &target_sptr) override
 Has to be called before using this object.
 
template<>
const char *const registered_name
 
- Public Member Functions inherited from stir::RegisteredParsingObject< QuadraticPrior< elemT >, GeneralisedPrior< DiscretisedDensity< 3, elemT > >, PriorWithParabolicSurrogate< 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.
 
- Public Member Functions inherited from stir::GeneralisedPrior< DiscretisedDensity< 3, elemT > >
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 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::ParsingObject
 ParsingObject (const ParsingObject &)
 
ParsingObjectoperator= (const ParsingObject &)
 
void ask_parameters ()
 
bool parse (std::istream &f)
 
bool parse (const char *const filename)
 

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

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...
 
Array< 3, float > weights
 penalty weights More...
 
std::string kappa_filename
 Filename for the $\kappa$ image that will be read by post_processing()
 
- 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::RegisteredParsingObject< QuadraticPrior< elemT >, GeneralisedPrior< DiscretisedDensity< 3, elemT > >, PriorWithParabolicSurrogate< DiscretisedDensity< 3, elemT > > >
static GeneralisedPrior< DiscretisedDensity< 3, elemT > > * read_from_stream (std::istream *)
 Construct a new object (of type Derived) by parsing the istream. More...
 
- 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_lessRegistryType
 The type of the registry.
 
- Static Protected Member Functions inherited from stir::RegisteredObject< GeneralisedPrior< DiscretisedDensity< 3, elemT > > >
static RegistryTyperegistry ()
 Static function returning the registry. More...
 

Detailed Description

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

A class in the GeneralisedPrior hierarchy. This implements a quadratic Gibbs prior.

The prior is computed as follows:

\[ f = \frac{1}{4} \sum_{r,dr} w_{dr} (\lambda_r - \lambda_{r+dr})^2 * \kappa_r * \kappa_{r+dr} \]

with gradient

\[ g_r = \sum_{dr} w_{dr} (\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 $\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.

Parsing
These are the keywords that can be used in addition to the ones in GeneralPrior.
Quadratic 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}}}
; 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 Quadratic Prior Parameters:=

Member Function Documentation

◆ parabolic_surrogate_curvature_depends_on_argument()

template<typename elemT>
bool stir::QuadraticPrior< elemT >::parabolic_surrogate_curvature_depends_on_argument ( ) const
inlineoverridevirtual

A function that allows skipping some computations if the curvature is independent of the current_estimate.

Defaults to return true, but can be overloaded by the derived class.

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

◆ is_convex()

template<typename elemT >
bool stir::QuadraticPrior< 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 > >.

◆ parabolic_surrogate_curvature()

template<typename elemT>
void stir::QuadraticPrior< elemT >::parabolic_surrogate_curvature ( DiscretisedDensity< 3, elemT > &  parabolic_surrogate_curvature,
const DiscretisedDensity< 3, elemT > &  current_image_estimate 
)
overridevirtual

compute the parabolic surrogate for the prior

in the case of quadratic priors this will just be the sum of weighting coefficients

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

◆ get_weights()

template<typename elemT >
Array< 3, float > stir::QuadraticPrior< elemT >::get_weights ( ) const

get penalty weights for the neighbourhood

get penalty weights for the neigbourhood

◆ set_weights()

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

set penalty weights for the neighbourhood

set penalty weights for the neigbourhood

◆ get_kappa_sptr()

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

◆ post_processing()

template<typename elemT >
bool stir::QuadraticPrior< 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.

Member Data Documentation

◆ gradient_filename_prefix

template<typename elemT>
std::string stir::QuadraticPrior< 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.

◆ weights

template<typename elemT>
Array<3, float> stir::QuadraticPrior< 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.

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