Implementation of the relaxed Ordered Subsets Separable Paraboloidal Surrogate ( OSSPS) More...

#include "stir/OSSPS/OSSPSReconstruction.h"

Inheritance diagram for stir::OSSPSReconstruction< TargetT >:

[legend]

Public Member Functions
	OSSPSReconstruction ()
	Default constructor (calls set_defaults())

	OSSPSReconstruction (const std::string &parameter_filename)
	Constructor, initialises everything from parameter file, or (when parameter_filename == "") by calling ask_parameters().

void	ask_parameters () override
	prompts the user to enter parameter values manually

OSSPSReconstruction &	get_parameters ()
	accessor for the external parameters

const OSSPSReconstruction &	get_parameters () const
	accessor for the external parameters

std::string	method_info () const override
	gives method information

Succeeded	precompute_denominator_of_conditioner_without_penalty ()
	Precompute the data-dependent part of the denominator for the preconditioner.

Succeeded	set_up (shared_ptr< TargetT > const &target_image_ptr) override
	operations prior to the iterations

void	update_estimate (TargetT &current_image_estimate) override
	the principal operations for updating the image iterates at each iteration

Public Member Functions inherited from stir::RegisteredParsingObject< OSSPSReconstruction< TargetT >, Reconstruction< TargetT >, IterativeReconstruction< TargetT > >
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::IterativeReconstruction< TargetT >
int	get_subiteration_num () const
	accessor for the subiteration counter

int	get_subset_num ()
	accessor for finding the current subset number

virtual TargetT *	get_initial_data_ptr () const
	Gets a pointer to the initial data.

Succeeded	reconstruct () override
	executes the reconstruction

Succeeded	reconstruct (shared_ptr< TargetT > const &target_data_sptr) override
	executes the reconstruction with target_data_sptr as initial value

std::string	make_filename_prefix_subiteration_num (const std::string &filename_prefix) const
	A utility function that creates a filename_prefix by appending the current subiteration number.

std::string	make_filename_prefix_subiteration_num () const
	A utility function that creates the output filename_prefix for the current subiteration number.

GeneralisedObjectiveFunction< TargetT > const &	get_objective_function () const

shared_ptr< GeneralisedObjectiveFunction< TargetT > >	get_objective_function_sptr () const

const int	get_max_num_full_iterations () const
	the maximum allowed number of full iterations

const int	get_num_subsets () const
	the number of ordered subsets

const int	get_num_subiterations () const
	the number of subiterations

const int	get_start_subiteration_num () const
	value with which to initialize the subiteration counter

const int	get_start_subset_num () const
	the starting subset number

const int	get_save_interval () const
	subiteration interval at which data will be saved

const bool	get_randomise_subset_order () const
	signals whether to randomise the subset order in each iteration

const DataProcessor< TargetT > &	get_inter_iteration_filter () const
	inter-iteration filter

shared_ptr< DataProcessor< TargetT > >	get_inter_iteration_filter_sptr ()

const int	get_inter_iteration_filter_interval () const
	subiteration interval at which to apply inter-iteration filters

const int	get_report_objective_function_values_interval () const
	subiteration interval at which to report the values of the objective function

void	set_objective_function_sptr (const shared_ptr< GeneralisedObjectiveFunction< TargetT > > &)
	The objective function that will be optimised.

void	set_max_num_full_iterations (const int)
	the maximum allowed number of full iterations

void	set_num_subsets (const int)
	the number of ordered subsets

void	set_num_subiterations (const int)
	the number of subiterations

void	set_start_subiteration_num (const int)
	value with which to initialize the subiteration counter

void	set_start_subset_num (const int)
	the starting subset number

void	set_save_interval (const int)
	subiteration interval at which data will be saved

void	set_randomise_subset_order (const bool)
	signals whether to randomise the subset order in each iteration

void	set_inter_iteration_filter_ptr (const shared_ptr< DataProcessor< TargetT > > &)
	inter-iteration filter

void	set_inter_iteration_filter_interval (const int)
	subiteration interval at which to apply inter-iteration filters

void	set_report_objective_function_values_interval (const int)
	subiteration interval at which to report the values of the objective function

void	set_input_data (const shared_ptr< ExamData > &arg) override
	set_input_data

const ExamData &	get_input_data () const override
	get input data

Public Member Functions inherited from stir::Reconstruction< TargetT >
	Reconstruction ()
	default constructor (calls set_defaults())

	~Reconstruction () override
	virtual destructor

void	set_output_filename_prefix (const std::string &)
	file name for output reconstructed images

void	set_output_file_format_ptr (const shared_ptr< OutputFileFormat< TargetT > > &)
	defines the format of the output files

void	set_post_processor_sptr (const shared_ptr< DataProcessor< TargetT > > &)
	post-filter

void	set_disable_output (bool _val)
	set_disable_output

shared_ptr< TargetT >	get_target_image ()
	get_reconstructed_image

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

Public Member Functions inherited from stir::TimedObject
void	reset_timers ()
	reset all timers kept by this object

void	stop_timers () const
	stop all timers kept by this object

void	start_timers (bool do_reset=false) const
	start all timers kept by this object

double	get_CPU_timer_value () const
	get current value of the CPU timer (since first use or last reset)

double	get_wall_clock_timer_value () const
	get current value of the wall-clock timer (since first use or last reset)

Static Public Attributes
static const char *const	registered_name
	Name which will be used when parsing a ProjectorByBinPair object.

Protected Member Functions
void	set_defaults () override
	Set defaults before parsing.

void	initialise_keymap () override
	Initialise all keywords.

bool	post_processing () override
	used to check acceptable parameter ranges, etc...

virtual void	end_of_iteration_processing (TargetT &current_estimate)
	operations for the end of the iteration

virtual void	check (TargetT const &target_data) const
	do consistency checks

void	initialise (const std::string &parameter_filename)
	This function initialises all parameters, either via parsing, or by calling ask_parameters() (when parameter_filename is the empty string).

void	set_defaults () override
	Set defaults before parsing.

void	initialise_keymap () override
	Initialise all keywords.

bool	post_processing () override
	used to check acceptable parameters after parsing

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

Protected Attributes
int	enforce_initial_positivity
	determines whether non-positive values in the initial image will be set to small positive ones

double	upper_bound
	restrict values to maximum

int	write_update_image
	boolean value to determine if the update images have to be written to disk

std::string	precomputed_denominator_filename
	optional name of the file containing the "precomputed denominator" - see Erdogan & Fessler for more info

float	relaxation_parameter
	relaxation parameter used (should be around 1) (see class documentation)

float	relaxation_gamma
	parameter determining how fast relaxation goes down (see class documentation)

shared_ptr< GeneralisedObjectiveFunction< TargetT > >	objective_function_sptr

int	subiteration_num
	the subiteration counter

bool	terminate_iterations
	used to abort the loop over iterations

int	max_num_full_iterations
	the maximum allowed number of full iterations

int	num_subsets
	the number of ordered subsets

int	num_subiterations
	the number of subiterations

int	start_subiteration_num
	value with which to initialize the subiteration counter

std::string	initial_data_filename
	name of the file containing the data for intializing the reconstruction

int	start_subset_num
	the starting subset number

int	save_interval
	subiteration interval at which data will be saved

bool	randomise_subset_order
	signals whether to randomise the subset order in each iteration

shared_ptr< DataProcessor< TargetT > >	inter_iteration_filter_ptr
	inter-iteration filter

int	inter_iteration_filter_interval
	subiteration interval at which to apply inter-iteration filters

int	report_objective_function_values_interval
	subiteration interval at which to report the values of the objective function

std::string	output_filename_prefix
	file name for output reconstructed images

shared_ptr< OutputFileFormat< TargetT > >	output_file_format_ptr
	defines the format of the output files

shared_ptr< DataProcessor< TargetT > >	post_filter_sptr
	post-filter

bool	_already_set_up

shared_ptr< TargetT >	target_data_sptr
	target_data_sptr

bool	_disable_output
	_disable_output

int	_verbosity
	Verbosity level.

KeyParser	parser

Additional Inherited Members
Static Public Member Functions inherited from stir::RegisteredParsingObject< OSSPSReconstruction< TargetT >, Reconstruction< TargetT >, IterativeReconstruction< TargetT > >
static Reconstruction< TargetT > *	read_from_stream (std::istream *)
	Construct a new object (of type Derived) by parsing the istream.

static Reconstruction< TargetT > *	read_from_stream (std::istream *)
	Construct a new object (of type Derived) by parsing the istream.

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

Detailed Description

template<class TargetT>
class stir::OSSPSReconstruction< TargetT >

Implementation of the relaxed Ordered Subsets Separable Paraboloidal Surrogate ( OSSPS)

See S. Ahn and J. A. Fessler, Globally convergent image reconstruction for emission tomography using relaxed ordered subsets algorithms, IEEE Trans. Med. Imag., vol. 22, no. 3, pp. 613-626, May 2003.

OSSPS is a relaxed preconditioned sub-gradient descent algorithm:

$\lambda^\mathrm{new} = \lambda + \zeta D \nabla \Psi$

with $\lambda$ the parameters to be estimated, $\Psi$ the objective function (see class GeneralisedObjectiveFunction) a diagonal matrix (the preconditioner) and $\zeta$ an iteration-dependent number called the relaxation parameter (see below).

depends on $\Psi$ . The data-dependent term in the preconditioner suggested in Ahn and Fessler turns out to be equal to the product of the Hessian of $\Psi$ summed over all columns (while using the data plug-in approach for approximating the Hessian close to the solution). Therefore, this OSSPS implementation uses GeneralisedObjectiveFunction::add_multiplication_with_approximate_Hessian_without_penalty. It is thus completely generic and can be applied to any objective function implemented in STIR.

Note that the exact paraboloidal surrogate algorithm by Erdogan and Fessler is currently not implemented.

When no prior info is specified, this should converge to the Maximum Likelihood solution (but in a very different way from MLEM!).

Note that this implementation probably assumes 'balanced subsets', i.e.

$\sum_{b \in \mathrm{subset}} P_{bv} = \sum_b P_{bv} \over \mathrm{numsubsets}$

Relaxation scheme

The relaxation value for the additive update follows the suggestion from Ahn and Fessler:

$\zeta= { \alpha \over 1+ \gamma n }$

with the (full) iteration number. The parameter $\alpha$ corresponds to the class member relaxation_parameter, and $\gamma$ to relaxation_gamma. Ahn and Fessler recommend to set $\alpha \approx 1$ and $\gamma$ small (e.g. 0.1).

Warning: This class should be the last in the Reconstruction hierarchy.

Todo: split into a preconditioned subgradient descent class and something that computes the preconditioner.

Member Function Documentation

◆ ask_parameters()

template<class TargetT>

void stir::OSSPSReconstruction< TargetT >::ask_parameters ( )

overridevirtual

prompts the user to enter parameter values manually

Reimplemented from stir::IterativeReconstruction< TargetT >.

References stir::ask(), stir::ask_filename_with_extension(), stir::ask_num(), ask_parameters(), stir::ask_string(), enforce_initial_positivity, stir::error(), stir::RegisteredObject< GeneralisedObjectiveFunction< TargetT > >::list_registered_names(), stir::max_filename_length, precomputed_denominator_filename, stir::RegisteredObject< GeneralisedObjectiveFunction< TargetT > >::read_registered_object(), upper_bound, and write_update_image.

Referenced by ask_parameters().

◆ method_info()

template<class TargetT>

std::string stir::OSSPSReconstruction< TargetT >::method_info ( ) const

overridevirtual

gives method information

Implements stir::Reconstruction< TargetT >.

References stir::IterativeReconstruction< TargetT >::inter_iteration_filter_interval, method_info(), and stir::IterativeReconstruction< TargetT >::num_subsets.

Referenced by method_info().

◆ precompute_denominator_of_conditioner_without_penalty()

template<class TargetT>

Succeeded stir::OSSPSReconstruction< TargetT >::precompute_denominator_of_conditioner_without_penalty ( )

Precompute the data-dependent part of the denominator for the preconditioner.

This function precomputes the denominator for the SPS (separable paraboloidal surrogates) algorithm. It calls GeneralisedObjectiveFunction::add_multiplication_with_approximate_Hessian_without_penalty on a vector filled with ones. For emission and transmission tomography, this corresponds to Erdogan and Fessler's approximations.

This method assumes that the objective function is concave and the output of add_multiplication_with_approximate_Hessian_without_penalty is non-positive.

The computed denominator is saved to file as output_filename_prefix plus "_precomputed_denominator".

References stir::info(), stir::Reconstruction< TargetT >::output_file_format_ptr, stir::Reconstruction< TargetT >::output_filename_prefix, precompute_denominator_of_conditioner_without_penalty(), stir::Timer::reset(), stir::Timer::start(), stir::Timer::stop(), stir::Timer::value(), and stir::warning().

Referenced by precompute_denominator_of_conditioner_without_penalty(), and set_up().