diff --git a/src/include/stir/recon_buildblock/FilterRootPrior.h b/src/include/stir/recon_buildblock/FilterRootPrior.h index 979bc0c496..e2850f9ab3 100644 --- a/src/include/stir/recon_buildblock/FilterRootPrior.h +++ b/src/include/stir/recon_buildblock/FilterRootPrior.h @@ -87,6 +87,8 @@ class FilterRootPrior: public FilterRootPrior(shared_ptr >const&, const float penalization_factor); + bool is_convex() const; + //! compute the value of the function /*! \warning Generally there is no function associated to this prior, so we just return 0 and write a warning the first time it's called. diff --git a/src/include/stir/recon_buildblock/GeneralisedConvexPrior.h b/src/include/stir/recon_buildblock/GeneralisedConvexPrior.h new file mode 100644 index 0000000000..83e9a4d7ca --- /dev/null +++ b/src/include/stir/recon_buildblock/GeneralisedConvexPrior.h @@ -0,0 +1,73 @@ +// +// +/* + Copyright (C) 2000- 2007, Hammersmith Imanet Ltd + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ +/*! + \file + \ingroup priors + \brief Declaration of class stir::GeneralisedConvexPrior + + \author Robert Twyman + +*/ + +#ifndef __stir_recon_buildblock_GeneralisedConvexPrior_H__ +#define __stir_recon_buildblock_GeneralisedConvexPrior_H__ + + +#include "stir/recon_buildblock/GeneralisedPrior.h" + +START_NAMESPACE_STIR + +class Succeeded; + +/*! + \ingroup priors + \brief + Make a brief +*/ +template + +class GeneralisedConvexPrior: + virtual public GeneralisedPrior + +{ +private: + typedef GeneralisedPrior base_type; + +public: + //! This computes a single row of the Hessian + /*! The method computes a row (i.e. at a densel/voxel, indicated by \c coords) of the Hessian at \c current_estimate. + Note that a row corresponds to an object of `DataT`. + The method (as implemented in derived classes) should store the result in \c prior_Hessian_for_single_densel. + */ + virtual void + compute_Hessian(DataT& prior_Hessian_for_single_densel, + const BasicCoordinate<3,int>& coords, + const DataT& current_image_estimate) const = 0; + +// float penalisation_factor; + + + + +// virtual void +// actual_compute_Hessian(DataT& prior_Hessian_for_single_densel, +// const BasicCoordinate<3,int>& coords, +// const DataT& current_image_estimate) const; + + int my_value = 1; + +}; + +END_NAMESPACE_STIR + +//#include "stir/recon_buildblock/GeneralisedPrior.inl" + +#endif diff --git a/src/include/stir/recon_buildblock/GeneralisedPrior.h b/src/include/stir/recon_buildblock/GeneralisedPrior.h index caa4bf4dab..a68bb04539 100644 --- a/src/include/stir/recon_buildblock/GeneralisedPrior.h +++ b/src/include/stir/recon_buildblock/GeneralisedPrior.h @@ -60,6 +60,19 @@ class GeneralisedPrior: virtual void compute_gradient(DataT& prior_gradient, const DataT ¤t_estimate) =0; + //! This computes a single row of the Hessian + /*! Default implementation just call error(). This function needs to be overridden by the + derived class. + + The method computes a row (i.e. at a densel/voxel, indicated by \c coords) of the Hessian at \c current_estimate. + Note that a row corresponds to an object of `DataT`. + The method (as implemented in derived classes) should store the result in \c prior_Hessian_for_single_densel. + */ + virtual void + compute_Hessian(DataT& prior_Hessian_for_single_densel, + const BasicCoordinate<3,int>& coords, + const DataT& current_image_estimate) const; + //! This should compute the multiplication of the Hessian with a vector and add it to \a output /*! Default implementation just call error(). This function needs to be overridden by the derived class. @@ -67,7 +80,7 @@ class GeneralisedPrior: Instead, accumulate_Hessian_times_input() should be used. This method remains for backwards comparability. \warning The derived class should accumulate in \a output. */ - virtual Succeeded + virtual void add_multiplication_with_approximate_Hessian(DataT& output, const DataT& input) const; @@ -76,7 +89,7 @@ class GeneralisedPrior: derived class. \warning The derived class should accumulate in \a output. */ - virtual Succeeded + virtual void accumulate_Hessian_times_input(DataT& output, const DataT& current_estimate, const DataT& input) const; @@ -89,6 +102,10 @@ class GeneralisedPrior: virtual Succeeded set_up(shared_ptr const& target_sptr); + //! Indicates if the prior is a smooth convex function + /*! If true, the prior is expected to have 0th, 1st and 2nd order behaviour implemented.*/ + virtual bool is_convex() const = 0; + protected: float penalisation_factor; //! sets value for penalisation factor diff --git a/src/include/stir/recon_buildblock/LogcoshPrior.h b/src/include/stir/recon_buildblock/LogcoshPrior.h index a3bed89478..fb9ed5aa4b 100644 --- a/src/include/stir/recon_buildblock/LogcoshPrior.h +++ b/src/include/stir/recon_buildblock/LogcoshPrior.h @@ -112,6 +112,8 @@ class LogcoshPrior: public parabolic_surrogate_curvature_depends_on_argument() const { return false; } + bool is_convex() const; + //! compute the value of the function double compute_value(const DiscretisedDensity<3,elemT> ¤t_image_estimate); @@ -124,13 +126,13 @@ class LogcoshPrior: public void parabolic_surrogate_curvature(DiscretisedDensity<3,elemT>& parabolic_surrogate_curvature, const DiscretisedDensity<3,elemT> ¤t_image_estimate); - //! compute Hessian - void compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, - const BasicCoordinate<3,int>& coords, - const DiscretisedDensity<3,elemT> ¤t_image_estimate); + virtual void + compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, + const BasicCoordinate<3,int>& coords, + const DiscretisedDensity<3,elemT> ¤t_image_estimate) const; - //! Compute the multiplication of the hessian of the prior multiplied by the input. - virtual Succeeded accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, + //! Compute the multiplication of the hessian of the prior (at \c current_estimate) and the \c input. + virtual void accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, const DiscretisedDensity<3,elemT>& current_estimate, const DiscretisedDensity<3,elemT>& input) const; @@ -224,21 +226,21 @@ class LogcoshPrior: public { return tanh(x)/x; } } - //! The Hessian of log(cosh()) is sech^2(x) = (1/cosh(x))^2 + //! The second partial derivatives of the LogCosh Prior /*! - This function returns the hessian of the logcosh function - * @param d the difference between the ith and jth voxel. - * @param scalar is the logcosh scalar value controlling the priors transition between the quadratic and linear behaviour - * @return the second derivative of the log-cosh function + derivative_20 refers to the second derivative w.r.t. x_j only (i.e. diagonal elements of the Hessian) + derivative_11 refers to the second derivative w.r.t. x_j and x_k (i.e. off-diagonal elements of the Hessian) + * @param x_j is the target voxel. + * @param x_k is the voxel in the neighbourhood. + * @return the second order partial derivatives of the LogCosh Prior */ - static inline float Hessian(const float d, const float scalar) - { - const float x = d * scalar; - return square((1/ cosh(x))); - } + //@{ + elemT derivative_20(const elemT x_j, const elemT x_k) const; + elemT derivative_11(const elemT x_j, const elemT x_k) const; + //@} }; END_NAMESPACE_STIR -#endif \ No newline at end of file +#endif diff --git a/src/include/stir/recon_buildblock/PLSPrior.h b/src/include/stir/recon_buildblock/PLSPrior.h index 91a8d83c0c..0c1458abf8 100644 --- a/src/include/stir/recon_buildblock/PLSPrior.h +++ b/src/include/stir/recon_buildblock/PLSPrior.h @@ -121,6 +121,8 @@ class PLSPrior: public /*! \todo set the anatomical image to zero if not defined */ virtual Succeeded set_up(shared_ptr > const& target_sptr); + bool is_convex() const; + //! compute the value of the function double compute_value(const DiscretisedDensity<3,elemT> ¤t_image_estimate); diff --git a/src/include/stir/recon_buildblock/PriorWithParabolicSurrogate.h b/src/include/stir/recon_buildblock/PriorWithParabolicSurrogate.h index 57dea783d6..df6e0c8d9c 100644 --- a/src/include/stir/recon_buildblock/PriorWithParabolicSurrogate.h +++ b/src/include/stir/recon_buildblock/PriorWithParabolicSurrogate.h @@ -38,7 +38,7 @@ START_NAMESPACE_STIR */ template class PriorWithParabolicSurrogate: - public GeneralisedPrior + virtual public GeneralisedPrior { public: diff --git a/src/include/stir/recon_buildblock/QuadraticPrior.h b/src/include/stir/recon_buildblock/QuadraticPrior.h index 4028169276..281edcbb47 100644 --- a/src/include/stir/recon_buildblock/QuadraticPrior.h +++ b/src/include/stir/recon_buildblock/QuadraticPrior.h @@ -25,6 +25,7 @@ #include "stir/RegisteredParsingObject.h" #include "stir/recon_buildblock/PriorWithParabolicSurrogate.h" +#include "stir/recon_buildblock/GeneralisedConvexPrior.h" #include "stir/Array.h" #include "stir/DiscretisedDensity.h" #include "stir/shared_ptr.h" @@ -77,15 +78,16 @@ START_NAMESPACE_STIR */ template class QuadraticPrior: public - RegisteredParsingObject< QuadraticPrior, - GeneralisedPrior >, - PriorWithParabolicSurrogate > - > + RegisteredParsingObject< + QuadraticPrior, + GeneralisedPrior >, + GeneralisedConvexPrior > >, + PriorWithParabolicSurrogate > { private: typedef RegisteredParsingObject< QuadraticPrior, - GeneralisedPrior >, + GeneralisedConvexPrior >, PriorWithParabolicSurrogate > > base_type; @@ -103,6 +105,8 @@ class QuadraticPrior: public parabolic_surrogate_curvature_depends_on_argument() const { return false; } + bool is_convex() const; + //! compute the value of the function double compute_value(const DiscretisedDensity<3,elemT> ¤t_image_estimate); @@ -116,20 +120,20 @@ class QuadraticPrior: public void parabolic_surrogate_curvature(DiscretisedDensity<3,elemT>& parabolic_surrogate_curvature, const DiscretisedDensity<3,elemT> ¤t_image_estimate); - //! compute Hessian - void compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, - const BasicCoordinate<3,int>& coords, - const DiscretisedDensity<3,elemT> ¤t_image_estimate); + virtual void + compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, + const BasicCoordinate<3,int>& coords, + const DiscretisedDensity<3,elemT> ¤t_image_estimate) const; //! Call accumulate_Hessian_times_input - virtual Succeeded + virtual void add_multiplication_with_approximate_Hessian(DiscretisedDensity<3,elemT>& output, const DiscretisedDensity<3,elemT>& input) const; //! 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. - virtual Succeeded accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, + virtual void accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, const DiscretisedDensity<3,elemT>& current_estimate, const DiscretisedDensity<3,elemT>& input) const; @@ -179,6 +183,19 @@ class QuadraticPrior: public virtual bool post_processing(); private: shared_ptr > kappa_ptr; + + //! The second partial derivatives of the Quadratic Prior + /*! + derivative_20 refers to the second derivative w.r.t. x_j (i.e. diagonal elements of the Hessian) + derivative_11 refers to the second derivative w.r.t. x_j and x_k (i.e. off-diagonal elements of the Hessian) + * @param x_j is the target voxel. + * @param x_k is the voxel in the neighbourhood. + * @return the second order partial derivatives of the Quadratic Prior + */ + //@{ + elemT derivative_20(const elemT x_j, const elemT x_k) const; + elemT derivative_11(const elemT x_j, const elemT x_k) const; + //@} }; diff --git a/src/include/stir/recon_buildblock/RelativeDifferencePrior.h b/src/include/stir/recon_buildblock/RelativeDifferencePrior.h index b4eaf4c74e..6df24dba0f 100644 --- a/src/include/stir/recon_buildblock/RelativeDifferencePrior.h +++ b/src/include/stir/recon_buildblock/RelativeDifferencePrior.h @@ -125,11 +125,19 @@ class RelativeDifferencePrior: public void compute_gradient(DiscretisedDensity<3,elemT>& prior_gradient, const DiscretisedDensity<3,elemT> ¤t_image_estimate); + virtual void compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, + const BasicCoordinate<3,int>& coords, + const DiscretisedDensity<3,elemT> ¤t_image_estimate) const; - virtual Succeeded + virtual void add_multiplication_with_approximate_Hessian(DiscretisedDensity<3,elemT>& output, const DiscretisedDensity<3,elemT>& input) const; + //! Compute the multiplication of the hessian of the prior multiplied by the input. + virtual void accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, + const DiscretisedDensity<3,elemT>& current_estimate, + const DiscretisedDensity<3,elemT>& input) const; + //! get the gamma value used in RDP float get_gamma() const; //! set the gamma value used in the RDP @@ -160,6 +168,8 @@ class RelativeDifferencePrior: public //! Has to be called before using this object virtual Succeeded set_up(shared_ptr > const& target_sptr); + bool is_convex() const; + protected: //! Create variable gamma for Relative Difference Penalty float gamma; @@ -193,6 +203,22 @@ class RelativeDifferencePrior: public virtual bool post_processing(); private: shared_ptr > kappa_ptr; + + //! The second partial derivatives of the Relative Difference Prior + /*! + derivative_20 refers to the second derivative w.r.t. x_j only (i.e. diagonal elements of the Hessian) + derivative_11 refers to the second derivative w.r.t. x_j and x_k (i.e. off-diagonal elements of the Hessian) + See J. Nuyts, et al., 2002, Equation 7. + In the instance x_j, x_k and epsilon equal 0.0, these functions return 0.0 to prevent returning an undefined value + due to 0/0 computation. This is a "reasonable" solution to this issue. + * @param x_j is the target voxel. + * @param x_k is the voxel in the neighbourhood. + * @return the second order partial derivatives of the Relative Difference Prior + */ + //@{ + elemT derivative_20(const elemT x_j, const elemT x_k) const; + elemT derivative_11(const elemT x_j, const elemT x_k) const; + //@} }; diff --git a/src/iterative/OSSPS/OSSPSReconstruction.cxx b/src/iterative/OSSPS/OSSPSReconstruction.cxx index e73478ff91..eed4624e50 100644 --- a/src/iterative/OSSPS/OSSPSReconstruction.cxx +++ b/src/iterative/OSSPS/OSSPSReconstruction.cxx @@ -327,10 +327,10 @@ update_estimate(TargetT ¤t_image_estimate) // For the quadratic prior, this is independent of the image (only on kappa's) // And of course, it's also independent when there is no prior // TODO by default, this should be off probably (to save time). + auto* parabolic_surrogate_prior = dynamic_cast*>(this->get_prior_ptr()); const bool recompute_penalty_term_in_denominator = !this->objective_function_sptr->prior_is_zero() && - static_cast const&>(*this->get_prior_ptr()). - parabolic_surrogate_curvature_depends_on_argument(); + parabolic_surrogate_prior->parabolic_surrogate_curvature_depends_on_argument(); #ifndef PARALLEL //CPUTimer subset_timer; //subset_timer.start(); @@ -366,8 +366,7 @@ update_estimate(TargetT ¤t_image_estimate) // avoid work (or crash) when penalty is 0 if (!this->objective_function_sptr->prior_is_zero()) { - static_cast&>(*get_prior_ptr()). - parabolic_surrogate_curvature(*work_image_ptr, current_image_estimate); + parabolic_surrogate_prior->parabolic_surrogate_curvature(*work_image_ptr, current_image_estimate); //*work_image_ptr *= 2; //*work_image_ptr += *precomputed_denominator_ptr ; std::transform(work_image_ptr->begin_all(), work_image_ptr->end_all(), diff --git a/src/recon_buildblock/CMakeLists.txt b/src/recon_buildblock/CMakeLists.txt index 6ed94fb7b1..d1dbabf7ce 100644 --- a/src/recon_buildblock/CMakeLists.txt +++ b/src/recon_buildblock/CMakeLists.txt @@ -58,6 +58,7 @@ set(${dir_LIB_SOURCES} BinNormalisationFromAttenuationImage.cxx BinNormalisationSPECT.cxx GeneralisedPrior.cxx + GeneralisedConvexPrior.cxx ProjDataRebinning.cxx FourierRebinning.cxx PLSPrior.cxx diff --git a/src/recon_buildblock/FilterRootPrior.cxx b/src/recon_buildblock/FilterRootPrior.cxx index 90ebaa8ece..638bdcd067 100644 --- a/src/recon_buildblock/FilterRootPrior.cxx +++ b/src/recon_buildblock/FilterRootPrior.cxx @@ -41,6 +41,12 @@ FilterRootPrior(shared_ptr >const& filter_sptr, float penal this->penalisation_factor = penalisation_factor_v; } +template +bool FilterRootPrior:: +is_convex() const +{ + return false; +} template < class T> static inline int diff --git a/src/recon_buildblock/GeneralisedConvexPrior.cxx b/src/recon_buildblock/GeneralisedConvexPrior.cxx new file mode 100644 index 0000000000..9ef485bfc2 --- /dev/null +++ b/src/recon_buildblock/GeneralisedConvexPrior.cxx @@ -0,0 +1,114 @@ +// +// +/* + Copyright (C) 2002- 2009, Hammersmith Imanet Ltd + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ +/*! + \file + \ingroup priors + \brief implementation of the stir::GeneralisedConvexPrior + + \author Robert Twyman +*/ + +#include "stir/recon_buildblock/GeneralisedConvexPrior.h" +//#include "stir/DiscretisedDensity.h" +//#include "stir/Succeeded.h" +#include "stir/modelling/ParametricDiscretisedDensity.h" +//#include "stir/modelling/KineticParameters.h" + +START_NAMESPACE_STIR + +//template +//void +//GeneralisedConvexPrior:: +//actual_compute_Hessian(elemT& prior_Hessian_for_single_densel, +// const BasicCoordinate<3,int>& coords, +// const elemT& current_image_estimate) const +//{ +// assert( prior_Hessian_for_single_densel.has_same_characteristics(current_image_estimate)); +// prior_Hessian_for_single_densel.fill(0); +// if (this->penalisation_factor==0) +// { +// return; +// } +// +// this->check(current_image_estimate); +// +// const DiscretisedDensityOnCartesianGrid<3,elemT>& current_image_cast = +// dynamic_cast< const DiscretisedDensityOnCartesianGrid<3,elemT> &>(current_image_estimate); +// +// DiscretisedDensityOnCartesianGrid<3,elemT>& prior_Hessian_for_single_densel_cast = +// dynamic_cast &>(prior_Hessian_for_single_densel); +// +// if (weights.get_length() ==0) +// { +// compute_weights(weights, current_image_cast.get_grid_spacing(), this->only_2D); +// } +// +// +// const bool do_kappa = !is_null_ptr(kappa_ptr); +// +// if (do_kappa && kappa_ptr->has_same_characteristics(current_image_estimate)) +// error("QuadraticPrior: kappa image has not the same index range as the reconstructed image\n"); +// +// const int z = coords[1]; +// const int y = coords[2]; +// const int x = coords[3]; +// const int min_dz = max(weights.get_min_index(), prior_Hessian_for_single_densel.get_min_index()-z); +// const int max_dz = min(weights.get_max_index(), prior_Hessian_for_single_densel.get_max_index()-z); +// +// const int min_dy = max(weights[0].get_min_index(), prior_Hessian_for_single_densel[z].get_min_index()-y); +// const int max_dy = min(weights[0].get_max_index(), prior_Hessian_for_single_densel[z].get_max_index()-y); +// +// const int min_dx = max(weights[0][0].get_min_index(), prior_Hessian_for_single_densel[z][y].get_min_index()-x); +// const int max_dx = min(weights[0][0].get_max_index(), prior_Hessian_for_single_densel[z][y].get_max_index()-x); +// +// elemT diagonal = 0; +// for (int dz=min_dz;dz<=max_dz;++dz) +// for (int dy=min_dy;dy<=max_dy;++dy) +// for (int dx=min_dx;dx<=max_dx;++dx) +// { +// elemT current = 0.0; +// if (dz == 0 && dy == 0 && dx == 0) +// { +// // The j == k case (diagonal Hessian element), which is a sum over the neighbourhood. +// for (int ddz=min_dz;ddz<=max_dz;++ddz) +// for (int ddy=min_dy;ddy<=max_dy;++ddy) +// for (int ddx=min_dx;ddx<=max_dx;++ddx) +// { +// elemT diagonal_current = weights[ddz][ddy][ddx] * +// derivative_20(current_image_estimate[z][y][x], +// current_image_estimate[z + ddz][y + ddy][x + ddx]); +// if (do_kappa) +// diagonal_current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+ddz][y+ddy][x+ddx]; +// current += diagonal_current; +// } +// } +// else +// { +// // The j != k vases (off-diagonal Hessian elements) +// current = weights[dz][dy][dx] * derivative_11(current_image_estimate[z][y][x], +// current_image_estimate[z + dz][y + dy][x + dx]); +// if (do_kappa) +// current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; +// } +// prior_Hessian_for_single_densel_cast[z+dz][y+dy][x+dx] = + current*this->penalisation_factor; +// } +//} + + +# ifdef _MSC_VER + // prevent warning message on instantiation of abstract class +# pragma warning(disable:4661) +# endif + + template class GeneralisedConvexPrior >; + template class GeneralisedConvexPrior; + +END_NAMESPACE_STIR diff --git a/src/recon_buildblock/GeneralisedObjectiveFunction.cxx b/src/recon_buildblock/GeneralisedObjectiveFunction.cxx index 41b842c2ed..f3490fa5f7 100644 --- a/src/recon_buildblock/GeneralisedObjectiveFunction.cxx +++ b/src/recon_buildblock/GeneralisedObjectiveFunction.cxx @@ -263,10 +263,7 @@ add_multiplication_with_approximate_sub_Hessian(TargetT& output, { // TODO used boost:scoped_ptr shared_ptr prior_output_sptr(output.get_empty_copy()); - if (this->prior_sptr->add_multiplication_with_approximate_Hessian(*prior_output_sptr, output) == - Succeeded::no) - return Succeeded::no; - + this->prior_sptr->add_multiplication_with_approximate_Hessian(*prior_output_sptr, output); // (*prior_output_sptr)/= num_subsets; // output -= *prior_output_sptr; @@ -383,9 +380,7 @@ accumulate_sub_Hessian_times_input(TargetT& output, { // TODO used boost:scoped_ptr shared_ptr prior_output_sptr(output.get_empty_copy()); - if (this->prior_sptr->accumulate_Hessian_times_input(*prior_output_sptr, current_image_estimate, output) == - Succeeded::no) - return Succeeded::no; + this->prior_sptr->accumulate_Hessian_times_input(*prior_output_sptr, current_image_estimate, output); typename TargetT::const_full_iterator prior_output_iter = prior_output_sptr->begin_all_const(); const typename TargetT::const_full_iterator end_prior_output_iter = prior_output_sptr->end_all_const(); diff --git a/src/recon_buildblock/GeneralisedPrior.cxx b/src/recon_buildblock/GeneralisedPrior.cxx index 5751c0a6b5..d4ce84b42f 100644 --- a/src/recon_buildblock/GeneralisedPrior.cxx +++ b/src/recon_buildblock/GeneralisedPrior.cxx @@ -52,18 +52,30 @@ set_up(shared_ptr const&) } template -Succeeded +void +GeneralisedPrior:: +compute_Hessian(TargetT& output, + const BasicCoordinate<3,int>& coords, + const TargetT& current_image_estimate) const +{ + if (this->is_convex()) + error("GeneralisedPrior:\n compute_Hessian implementation is not overloaded by your convex prior."); + else + error("GeneralisedPrior:\n compute_Hessian is not implemented because the prior is not convex."); +} + +template +void GeneralisedPrior:: add_multiplication_with_approximate_Hessian(TargetT& output, const TargetT& input) const { error("GeneralisedPrior:\n" "add_multiplication_with_approximate_Hessian implementation is not overloaded by your prior."); - return Succeeded::no; } template -Succeeded +void GeneralisedPrior:: accumulate_Hessian_times_input(TargetT& output, const TargetT& current_estimate, @@ -71,7 +83,6 @@ accumulate_Hessian_times_input(TargetT& output, { error("GeneralisedPrior:\n" "accumulate_Hessian_times_input implementation is not overloaded by your prior."); - return Succeeded::no; } template diff --git a/src/recon_buildblock/LogcoshPrior.cxx b/src/recon_buildblock/LogcoshPrior.cxx index f14742d7b1..5324498430 100644 --- a/src/recon_buildblock/LogcoshPrior.cxx +++ b/src/recon_buildblock/LogcoshPrior.cxx @@ -144,6 +144,14 @@ LogcoshPrior::LogcoshPrior(const bool only_2D_v, float penalisation_facto this->scalar = scalar_v; } +template +bool +LogcoshPrior:: +is_convex() const +{ + return true; +} + //! get penalty weights for the neighbourhood template Array<3,float> @@ -380,7 +388,7 @@ void LogcoshPrior:: compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, const BasicCoordinate<3,int>& coords, - const DiscretisedDensity<3,elemT> ¤t_image_estimate) + const DiscretisedDensity<3,elemT> ¤t_image_estimate) const { assert( prior_Hessian_for_single_densel.has_same_characteristics(current_image_estimate)); prior_Hessian_for_single_densel.fill(0); @@ -388,6 +396,9 @@ compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, { return; } + + this->check(current_image_estimate); + const DiscretisedDensityOnCartesianGrid<3,elemT>& current_image_cast = dynamic_cast< const DiscretisedDensityOnCartesianGrid<3,elemT> &>(current_image_estimate); @@ -421,20 +432,35 @@ compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, for (int dy=min_dy;dy<=max_dy;++dy) for (int dx=min_dx;dx<=max_dx;++dx) { - // sech^2(x * scalar); sech(x) = 1/cosh(x) - elemT voxel_diff= current_image_estimate[z][y][x] - current_image_estimate[z+dz][y+dy][x+dx]; - elemT current = weights[dz][dy][dx] * Hessian(voxel_diff, this->scalar); - - if (do_kappa) - current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; - - diagonal += current; - prior_Hessian_for_single_densel_cast[z+dz][y+dy][x+dx] = -current*this->penalisation_factor; + elemT current = 0.0; + if (dz == 0 && dy == 0 && dx == 0) + { + // The j == k case (diagonal Hessian element), which is a sum over the neighbourhood. + for (int ddz=min_dz;ddz<=max_dz;++ddz) + for (int ddy=min_dy;ddy<=max_dy;++ddy) + for (int ddx=min_dx;ddx<=max_dx;++ddx) + { + elemT diagonal_current = weights[ddz][ddy][ddx] * + derivative_20(current_image_estimate[z][y][x], + current_image_estimate[z + ddz][y + ddy][x + ddx]); + if (do_kappa) + diagonal_current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+ddz][y+ddy][x+ddx]; + current += diagonal_current; + } + } + else + { + // The j != k vases (off-diagonal Hessian elements) + current = weights[dz][dy][dx] * derivative_11(current_image_estimate[z][y][x], + current_image_estimate[z + dz][y + dy][x + dx]); + if (do_kappa) + current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; + } + prior_Hessian_for_single_densel_cast[z+dz][y+dy][x+dx] = + current*this->penalisation_factor; } - - prior_Hessian_for_single_densel[z][y][x]= diagonal * this->penalisation_factor; } + template void LogcoshPrior::parabolic_surrogate_curvature(DiscretisedDensity<3,elemT>& parabolic_surrogate_curvature, @@ -505,7 +531,7 @@ LogcoshPrior::parabolic_surrogate_curvature(DiscretisedDensity<3,elemT>& } template -Succeeded +void LogcoshPrior:: accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, const DiscretisedDensity<3,elemT>& current_estimate, @@ -517,7 +543,7 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, assert( output.has_same_characteristics(input)); if (this->penalisation_factor==0) { - return Succeeded::yes; + return; } DiscretisedDensityOnCartesianGrid<3,elemT>& output_cast = @@ -561,8 +587,18 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, for (int dy=min_dy;dy<=max_dy;++dy) for (int dx=min_dx;dx<=max_dx;++dx) { - elemT voxel_diff= current_estimate[z][y][x] - current_estimate[z+dz][y+dy][x+dx]; - elemT current = weights[dz][dy][dx] * Hessian( voxel_diff, this->scalar) * input[z+dz][y+dy][x+dx]; + elemT current = weights[dz][dy][dx]; + if (dz == dy == dz == 0) { + // The j == k case + current *= derivative_20(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x]; + } else { + current *= (derivative_20(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x] + + derivative_11(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * + input[z + dz][y + dy][x + dx]); + } if (do_kappa) current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; @@ -574,7 +610,22 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, } } } - return Succeeded::yes; +} + +template +elemT +LogcoshPrior:: +derivative_20(const elemT x_j, const elemT x_k) const +{ + return square((1/ cosh((x_j - x_k) * scalar))); +} + +template +elemT +LogcoshPrior:: +derivative_11(const elemT x_j, const elemT x_k) const +{ + return -derivative_20(x_j, x_k); } # ifdef _MSC_VER diff --git a/src/recon_buildblock/PLSPrior.cxx b/src/recon_buildblock/PLSPrior.cxx index 0b782546fa..891d554fd5 100644 --- a/src/recon_buildblock/PLSPrior.cxx +++ b/src/recon_buildblock/PLSPrior.cxx @@ -144,6 +144,13 @@ PLSPrior::PLSPrior() set_defaults(); } +template +bool +PLSPrior:: +is_convex() const +{ + return true; +} template PLSPrior::PLSPrior(const bool only_2D_v, float penalisation_factor_v) diff --git a/src/recon_buildblock/QuadraticPrior.cxx b/src/recon_buildblock/QuadraticPrior.cxx index 904b4b0d1b..58fd7813e5 100644 --- a/src/recon_buildblock/QuadraticPrior.cxx +++ b/src/recon_buildblock/QuadraticPrior.cxx @@ -19,6 +19,7 @@ */ #include "stir/recon_buildblock/QuadraticPrior.h" +#include "stir/recon_buildblock/GeneralisedConvexPrior.h" #include "stir/Succeeded.h" #include "stir/DiscretisedDensityOnCartesianGrid.h" #include "stir/IndexRange3D.h" @@ -147,6 +148,13 @@ QuadraticPrior::QuadraticPrior(const bool only_2D_v, float penalisation_f this->penalisation_factor = penalisation_factor_v; } +template +bool +QuadraticPrior:: +is_convex() const +{ + return true; +} //! get penalty weights for the neigbourhood template @@ -430,11 +438,11 @@ compute_gradient(DiscretisedDensity<3,elemT>& prior_gradient, } template -void +void QuadraticPrior:: compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, const BasicCoordinate<3,int>& coords, - const DiscretisedDensity<3,elemT> ¤t_image_estimate) + const DiscretisedDensity<3,elemT> ¤t_image_estimate) const { assert( prior_Hessian_for_single_densel.has_same_characteristics(current_image_estimate)); prior_Hessian_for_single_densel.fill(0); @@ -479,19 +487,32 @@ compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, for (int dy=min_dy;dy<=max_dy;++dy) for (int dx=min_dx;dx<=max_dx;++dx) { - // dz==0,dy==0,dx==0 will have weight 0, so we can just include it in the loop - elemT current = - weights[dz][dy][dx]; - - if (do_kappa) - current *= - (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; - - diagonal += current; - prior_Hessian_for_single_densel_cast[z+dz][y+dy][x+dx] = -current*this->penalisation_factor; + elemT current = 0.0; + if (dz == 0 && dy == 0 && dx == 0) + { + // The j == k case (diagonal Hessian element), which is a sum over the neighbourhood. + for (int ddz=min_dz;ddz<=max_dz;++ddz) + for (int ddy=min_dy;ddy<=max_dy;++ddy) + for (int ddx=min_dx;ddx<=max_dx;++ddx) + { + elemT diagonal_current = weights[ddz][ddy][ddx] * + derivative_20(current_image_estimate[z][y][x], + current_image_estimate[z + ddz][y + ddy][x + ddx]); + if (do_kappa) + diagonal_current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+ddz][y+ddy][x+ddx]; + current += diagonal_current; + } + } + else + { + // The j != k vases (off-diagonal Hessian elements) + current = weights[dz][dy][dx] * derivative_11(current_image_estimate[z][y][x], + current_image_estimate[z + dz][y + dy][x + dx]); + if (do_kappa) + current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; + } + prior_Hessian_for_single_densel_cast[z+dz][y+dy][x+dx] = + current*this->penalisation_factor; } - - prior_Hessian_for_single_densel[z][y][x]= diagonal * this->penalisation_factor; } template @@ -576,19 +597,83 @@ QuadraticPrior::parabolic_surrogate_curvature(DiscretisedDensity<3,elemT> } template -Succeeded +void QuadraticPrior:: add_multiplication_with_approximate_Hessian(DiscretisedDensity<3,elemT>& output, const DiscretisedDensity<3,elemT>& input) const { - return accumulate_Hessian_times_input(output, input, input); + // TODO this function overlaps enormously with parabolic_surrogate_curvature + // the only difference is that parabolic_surrogate_curvature uses input==1 + + assert( output.has_same_characteristics(input)); + if (this->penalisation_factor==0) + { + return; + } + + this->check(input); + + DiscretisedDensityOnCartesianGrid<3,elemT>& output_cast = + dynamic_cast &>(output); + + if (weights.get_length() ==0) + { + compute_weights(weights, output_cast.get_grid_spacing(), this->only_2D); + } + + const bool do_kappa = !is_null_ptr(kappa_ptr); + + if (do_kappa && !kappa_ptr->has_same_characteristics(input)) + error("QuadraticPrior: kappa image has not the same index range as the reconstructed image\n"); + + const int min_z = output.get_min_index(); + const int max_z = output.get_max_index(); + for (int z=min_z; z<=max_z; z++) + { + const int min_dz = max(weights.get_min_index(), min_z-z); + const int max_dz = min(weights.get_max_index(), max_z-z); + + const int min_y = output[z].get_min_index(); + const int max_y = output[z].get_max_index(); + + for (int y=min_y;y<= max_y;y++) + { + const int min_dy = max(weights[0].get_min_index(), min_y-y); + const int max_dy = min(weights[0].get_max_index(), max_y-y); + + const int min_x = output[z][y].get_min_index(); + const int max_x = output[z][y].get_max_index(); + + for (int x=min_x;x<= max_x;x++) + { + const int min_dx = max(weights[0][0].get_min_index(), min_x-x); + const int max_dx = min(weights[0][0].get_max_index(), max_x-x); + + elemT result = 0; + for (int dz=min_dz;dz<=max_dz;++dz) + for (int dy=min_dy;dy<=max_dy;++dy) + for (int dx=min_dx;dx<=max_dx;++dx) + { + elemT current = + weights[dz][dy][dx] * input[z+dz][y+dy][x+dx]; + + if (do_kappa) + current *= + (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; + result += current; + } + + output[z][y][x] += result * this->penalisation_factor; + } + } + } } template -Succeeded +void QuadraticPrior:: accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, - const DiscretisedDensity<3,elemT>& /*current_estimate*/, + const DiscretisedDensity<3,elemT>& current_estimate, const DiscretisedDensity<3,elemT>& input) const { // TODO this function overlaps enormously with parabolic_surrogate_curvature @@ -597,7 +682,7 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, assert( output.has_same_characteristics(input)); if (this->penalisation_factor==0) { - return Succeeded::yes; + return; } this->check(input); @@ -638,13 +723,31 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, const int min_dx = max(weights[0][0].get_min_index(), min_x-x); const int max_dx = min(weights[0][0].get_max_index(), max_x-x); + /// At this point, we have j = [z][y][x] + // The next for loops will have k = [z+dz][y+dy][x+dx] + // The following computes + //(H_{wf} y)_j = + // \sum_{k\in N_j} w_{(j,k)} f''_{d}(x_j,x_k) y_j + + // \sum_{(i \in N_j) \ne j} w_{(j,i)} f''_{od}(x_j, x_i) y_i + // Note the condition in the second sum that i is not equal to j + elemT result = 0; for (int dz=min_dz;dz<=max_dz;++dz) for (int dy=min_dy;dy<=max_dy;++dy) for (int dx=min_dx;dx<=max_dx;++dx) { - elemT current = - weights[dz][dy][dx] * input[z+dz][y+dy][x+dx]; + elemT current = weights[dz][dy][dx]; + if (dz == dy == dz == 0) { + // The j == k case + current *= derivative_20(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x]; + } else { + current *= (derivative_20(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x] + + derivative_11(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * + input[z + dz][y + dy][x + dx]); + } if (do_kappa) current *= @@ -656,9 +759,20 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, } } } - return Succeeded::yes; } +template +elemT +QuadraticPrior:: +derivative_20(const elemT x_j, const elemT x_k) const +{ return 1.0;} + +template +elemT +QuadraticPrior:: +derivative_11(const elemT x_j, const elemT x_k) const +{ return -1.0;} + # ifdef _MSC_VER // prevent warning message on reinstantiation, // note that we get a linking error if we don't have the explicit instantiation below diff --git a/src/recon_buildblock/RelativeDifferencePrior.cxx b/src/recon_buildblock/RelativeDifferencePrior.cxx index e7fce65ea5..69c443cc42 100644 --- a/src/recon_buildblock/RelativeDifferencePrior.cxx +++ b/src/recon_buildblock/RelativeDifferencePrior.cxx @@ -127,6 +127,7 @@ void RelativeDifferencePrior::set_defaults() { base_type::set_defaults(); +// this->_is_convex = true; this->only_2D = false; this->kappa_ptr.reset(); this->weights.recycle(); @@ -145,6 +146,14 @@ RelativeDifferencePrior::RelativeDifferencePrior() set_defaults(); } +template +bool +RelativeDifferencePrior:: +is_convex() const +{ + return true; +} + // Return the value of gamma - a RDP parameter template float @@ -178,6 +187,7 @@ template RelativeDifferencePrior::RelativeDifferencePrior(const bool only_2D_v, float penalisation_factor_v, float gamma_v, float epsilon_v) : only_2D(only_2D_v) { + set_defaults(); this->penalisation_factor = penalisation_factor_v; this->gamma = gamma_v; this->epsilon = epsilon_v; @@ -421,13 +431,204 @@ compute_gradient(DiscretisedDensity<3,elemT>& prior_gradient, } template -Succeeded +void +RelativeDifferencePrior:: +compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, + const BasicCoordinate<3,int>& coords, + const DiscretisedDensity<3,elemT> ¤t_image_estimate) const +{ + assert( prior_Hessian_for_single_densel.has_same_characteristics(current_image_estimate)); + prior_Hessian_for_single_densel.fill(0); + if (this->penalisation_factor==0) + { + return; + } + + this->check(current_image_estimate); + + const DiscretisedDensityOnCartesianGrid<3,elemT>& current_image_cast = + dynamic_cast< const DiscretisedDensityOnCartesianGrid<3,elemT> &>(current_image_estimate); + + DiscretisedDensityOnCartesianGrid<3,elemT>& prior_Hessian_for_single_densel_cast = + dynamic_cast &>(prior_Hessian_for_single_densel); + + if (weights.get_length() ==0) + { + compute_weights(weights, current_image_cast.get_grid_spacing(), this->only_2D); + } + + + const bool do_kappa = !is_null_ptr(kappa_ptr); + + if (do_kappa && kappa_ptr->has_same_characteristics(current_image_estimate)) + error("RelativeDifferencePrior: kappa image has not the same index range as the reconstructed image\n"); + + const int z = coords[1]; + const int y = coords[2]; + const int x = coords[3]; + const int min_dz = max(weights.get_min_index(), prior_Hessian_for_single_densel.get_min_index()-z); + const int max_dz = min(weights.get_max_index(), prior_Hessian_for_single_densel.get_max_index()-z); + + const int min_dy = max(weights[0].get_min_index(), prior_Hessian_for_single_densel[z].get_min_index()-y); + const int max_dy = min(weights[0].get_max_index(), prior_Hessian_for_single_densel[z].get_max_index()-y); + + const int min_dx = max(weights[0][0].get_min_index(), prior_Hessian_for_single_densel[z][y].get_min_index()-x); + const int max_dx = min(weights[0][0].get_max_index(), prior_Hessian_for_single_densel[z][y].get_max_index()-x); + + elemT diagonal = 0; + for (int dz=min_dz;dz<=max_dz;++dz) + for (int dy=min_dy;dy<=max_dy;++dy) + for (int dx=min_dx;dx<=max_dx;++dx) + { + elemT current = 0.0; + if (dz == 0 && dy == 0 && dx == 0) + { + // The j == k case (diagonal Hessian element), which is a sum over the neighbourhood. + for (int ddz=min_dz;ddz<=max_dz;++ddz) + for (int ddy=min_dy;ddy<=max_dy;++ddy) + for (int ddx=min_dx;ddx<=max_dx;++ddx) + { + elemT diagonal_current = weights[ddz][ddy][ddx] * + derivative_20(current_image_estimate[z][y][x], + current_image_estimate[z + ddz][y + ddy][x + ddx]); + if (do_kappa) + diagonal_current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+ddz][y+ddy][x+ddx]; + current += diagonal_current; + } + } + else + { + // The j != k vases (off-diagonal Hessian elements) + current = weights[dz][dy][dx] * derivative_11(current_image_estimate[z][y][x], + current_image_estimate[z + dz][y + dy][x + dx]); + if (do_kappa) + current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; + } + prior_Hessian_for_single_densel_cast[z+dz][y+dy][x+dx] = + current*this->penalisation_factor; + } +} + +template +void RelativeDifferencePrior:: add_multiplication_with_approximate_Hessian(DiscretisedDensity<3,elemT>& output, const DiscretisedDensity<3,elemT>& input) const { error("add_multiplication_with_approximate_Hessian() is not implemented in Relative Difference Prior."); - return Succeeded::no; +} + +template +void +RelativeDifferencePrior:: +accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, + const DiscretisedDensity<3,elemT>& current_estimate, + const DiscretisedDensity<3,elemT>& input) const +{ + // TODO this function overlaps enormously with parabolic_surrogate_curvature + // the only difference is that parabolic_surrogate_curvature uses input==1 + + assert( output.has_same_characteristics(input)); + if (this->penalisation_factor==0) + { + return; + } + + DiscretisedDensityOnCartesianGrid<3,elemT>& output_cast = + dynamic_cast &>(output); + + if (weights.get_length() ==0) + { + compute_weights(weights, output_cast.get_grid_spacing(), this->only_2D); + } + + const bool do_kappa = !is_null_ptr(kappa_ptr); + + if (do_kappa && !kappa_ptr->has_same_characteristics(input)) + error("LogcoshPrior: kappa image has not the same index range as the reconstructed image\n"); + + const int min_z = output.get_min_index(); + const int max_z = output.get_max_index(); + for (int z=min_z; z<=max_z; z++) + { + const int min_dz = max(weights.get_min_index(), min_z-z); + const int max_dz = min(weights.get_max_index(), max_z-z); + + const int min_y = output[z].get_min_index(); + const int max_y = output[z].get_max_index(); + + for (int y=min_y;y<= max_y;y++) + { + const int min_dy = max(weights[0].get_min_index(), min_y-y); + const int max_dy = min(weights[0].get_max_index(), max_y-y); + + const int min_x = output[z][y].get_min_index(); + const int max_x = output[z][y].get_max_index(); + + for (int x=min_x;x<= max_x;x++) + { + const int min_dx = max(weights[0][0].get_min_index(), min_x-x); + const int max_dx = min(weights[0][0].get_max_index(), max_x-x); + + /// At this point, we have j = [z][y][x] + // The next for loops will have k = [z+dz][y+dy][x+dx] + // The following computes + //(H_{wf} y)_j = + // \sum_{k\in N_j} w_{(j,k)} f''_{d}(x_j,x_k) y_j + + // \sum_{(i \in N_j) \ne j} w_{(j,i)} f''_{od}(x_j, x_i) y_i + // Note the condition in the second sum that i is not equal to j + + elemT result = 0; + for (int dz=min_dz;dz<=max_dz;++dz) + for (int dy=min_dy;dy<=max_dy;++dy) + for (int dx=min_dx;dx<=max_dx;++dx) + { + elemT current = weights[dz][dy][dx]; + if (dz == dy == dz == 0) { + // The j == k case + current *= derivative_20(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x]; + } else { + current *= (derivative_20(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x] + + derivative_11(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * + input[z + dz][y + dy][x + dx]); + } + + if (do_kappa) + current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; + + result += current; + } + + output[z][y][x] += result * this->penalisation_factor; + } + } + } +} + +template +elemT +RelativeDifferencePrior:: +derivative_20(const elemT x_j, const elemT x_k) const +{ + if (x_j > 0.0 || x_k > 0.0 || this->epsilon > 0.0) + return 2 * pow(2 * x_k + this->epsilon, 2) / + pow(x_j + x_k + this->gamma * abs(x_j - x_k) + this->epsilon, 3); + else + return 0.0; +} + +template +elemT +RelativeDifferencePrior:: +derivative_11(const elemT x_j, const elemT x_k) const +{ + if (x_j > 0.0 || x_k > 0.0 || this->epsilon > 0.0) + return - 2 * (2 * x_j + this->epsilon)*(2 * x_k + this->epsilon) / + pow(x_j + x_k + this->gamma * abs(x_j - x_k) + this->epsilon, 3); + else + return 0.0; } # ifdef _MSC_VER diff --git a/src/recon_test/test_priors.cxx b/src/recon_test/test_priors.cxx index 961afdf0a8..b7ce6f821d 100644 --- a/src/recon_test/test_priors.cxx +++ b/src/recon_test/test_priors.cxx @@ -29,11 +29,12 @@ #include "stir/recon_buildblock/QuadraticPrior.h" #include "stir/recon_buildblock/RelativeDifferencePrior.h" #include "stir/recon_buildblock/LogcoshPrior.h" -//#include "stir/recon_buildblock/PLSPrior.h" +#include "stir/recon_buildblock/PLSPrior.h" #include "stir/RunTests.h" #include "stir/IO/read_from_file.h" #include "stir/IO/write_to_file.h" #include "stir/info.h" +#include "stir/Verbosity.h" #include "stir/Succeeded.h" #include "stir/num_threads.h" #include @@ -53,6 +54,8 @@ START_NAMESPACE_STIR This test compares the result of GeneralisedPrior::compute_gradient() with a numerical gradient computed by using the GeneralisedPrior::compute_value() function. + Additionally, the Hessian's convexity is tested, via GeneralisedPrior::accumulate_Hessian_times_input(), + by evaluating the x^T Hx > 0 constraint. */ class GeneralisedPriorTests : public RunTests @@ -70,6 +73,14 @@ class GeneralisedPriorTests : public RunTests void construct_input_data(shared_ptr& density_sptr); void run_tests(); + + //! Set methods that control which tests are run. + //@{ + void set_do_test_gradient(bool d); + void set_do_test_Hessian_convexity(bool d); + void set_do_test_Hessian_against_numerical(bool d); + //@} + protected: char const * density_filename; shared_ptr > objective_function_sptr; @@ -79,6 +90,47 @@ class GeneralisedPriorTests : public RunTests void run_tests_for_objective_function(const std::string& test_name, GeneralisedPrior& objective_function, shared_ptr target_sptr); + + //! Tests the prior's gradient by comparing to the numerical gradient computed using perturbation response. + void test_gradient(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr); + + //! Test various configurations of the Hessian of the prior via accumulate_Hessian_times_input() for convexity + /*! + Tests the convexity condition: + \f[ x^T \cdot H_{\lambda}x >= 0 \f] + for all non-negative \c x and non-zero \c \lambda (Relative Difference Prior conditions). + This function constructs an array of configurations to test this condition and calls + \c test_Hessian_convexity_configuration(). + */ + void test_Hessian_convexity(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr); + + //! Tests the compute_Hessian method implemented into convex priors + /*! Performs a perturbation response using compute_gradient to determine if the compute_Hessian (for a single densel) + is within tolerance. + */ + void test_Hessian_against_numerical(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr); + +private: + //! Hessian test for a particular configuration of the Hessian concave condition + bool test_Hessian_convexity_configuration(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr, + const float beta, + const float input_multiplication, const float input_addition, + const float current_image_multiplication, const float current_image_addition); + + //! Variables to control which tests are run, see the set methods + //@{ + bool do_test_gradient = false; + bool do_test_Hessian_convexity = false; + bool do_test_Hessian_against_numerical = false; + //@} }; GeneralisedPriorTests:: @@ -86,6 +138,27 @@ GeneralisedPriorTests(char const * const density_filename) : density_filename(density_filename) {} +void +GeneralisedPriorTests:: +set_do_test_gradient(const bool d) +{ + do_test_gradient = d; +} + +void +GeneralisedPriorTests:: +set_do_test_Hessian_convexity(const bool d) +{ + do_test_Hessian_convexity = d; +} + +void +GeneralisedPriorTests:: +set_do_test_Hessian_against_numerical(const bool d) +{ + do_test_Hessian_against_numerical = d; +} + void GeneralisedPriorTests:: run_tests_for_objective_function(const std::string& test_name, @@ -96,46 +169,258 @@ run_tests_for_objective_function(const std::string& test_name, if (!check(objective_function.set_up(target_sptr)==Succeeded::yes, "set-up of objective function")) return; + if (do_test_gradient) + { + std::cerr << "----- test " << test_name << " --> Gradient\n"; + test_gradient(test_name, objective_function, target_sptr); + } + + if (do_test_Hessian_convexity) + { + std::cerr << "----- test " << test_name << " --> Hessian-vector product for convexity\n"; + test_Hessian_convexity(test_name, objective_function, target_sptr); + } + + if (do_test_Hessian_against_numerical) + { + std::cerr << "----- test " << test_name << " --> Hessian against numerical\n"; + test_Hessian_against_numerical(test_name, objective_function, target_sptr); + } +} + + +void +GeneralisedPriorTests:: +test_gradient(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr) +{ // setup images target_type& target(*target_sptr); shared_ptr gradient_sptr(target.get_empty_copy()); shared_ptr gradient_2_sptr(target.get_empty_copy()); info("Computing gradient",3); + const int verbosity_default = Verbosity::get(); + Verbosity::set(0); objective_function.compute_gradient(*gradient_sptr, target); + Verbosity::set(verbosity_default); this->set_tolerance(std::max(fabs(double(gradient_sptr->find_min())), fabs(double(gradient_sptr->find_max()))) /1000); info("Computing objective function at target",3); const double value_at_target = objective_function.compute_value(target); target_type::full_iterator target_iter=target.begin_all(); target_type::full_iterator gradient_iter=gradient_sptr->begin_all(); - target_type::full_iterator gradient_2_iter=gradient_2_sptr->begin_all(); + target_type::full_iterator gradient_2_iter=gradient_2_sptr->begin_all(); // setup perturbation response const float eps = 1e-3F; bool testOK = true; info("Computing gradient of objective function by numerical differences (this will take a while)",3); while(target_iter!=target.end_all())// && testOK) + { + const float org_image_value = *target_iter; + *target_iter += eps; // perturb current voxel + const double value_at_inc = objective_function.compute_value(target); + *target_iter = org_image_value; // restore + const float ngradient_at_iter = static_cast((value_at_inc - value_at_target)/eps); + *gradient_2_iter = ngradient_at_iter; + testOK = testOK && this->check_if_equal(ngradient_at_iter, *gradient_iter, "gradient"); + //for (int i=0; i<5 && target_iter!=target.end_all(); ++i) { - const float org_image_value = *target_iter; - *target_iter += eps; // perturb current voxel - const double value_at_inc = objective_function.compute_value(target); - *target_iter = org_image_value; // restore - const float ngradient_at_iter = static_cast((value_at_inc - value_at_target)/eps); - *gradient_2_iter = ngradient_at_iter; - testOK = testOK && this->check_if_equal(ngradient_at_iter, *gradient_iter, "gradient"); - //for (int i=0; i<5 && target_iter!=target.end_all(); ++i) - { - ++gradient_2_iter; ++target_iter; ++ gradient_iter; - } + ++gradient_2_iter; ++target_iter; ++ gradient_iter; } + } if (!testOK) + { + std::cerr << "Numerical gradient test failed with for " + test_name + " prior\n"; + info("Writing diagnostic files gradient" + test_name + ".hv, numerical_gradient" + test_name + ".hv"); + write_to_file("gradient" + test_name + ".hv", *gradient_sptr); + write_to_file("numerical_gradient" + test_name + ".hv", *gradient_2_sptr); + } +} + +void +GeneralisedPriorTests:: +test_Hessian_convexity(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr) +{ + if (!objective_function.is_convex()) + return; + /// Construct configurations + float beta_array[] = {0.01, 1, 100}; // Penalty strength should only affect scale + // Modifications to the input image + float input_multiplication_array[] = {-100, -1, 0.01, 1, 100}; // Test negative, small and large values + float input_addition_array[] = {-10, -1, -0.5, 0.0, 1, 10}; + // Modifications to the current image (Hessian computation) + float current_image_multiplication_array[] = {0.01, 1, 100}; + float current_image_addition_array[] = {0.0, 0.5, 1, 10}; // RDP has constraint that current_image is non-negative + + bool testOK = true; + float initial_beta = objective_function.get_penalisation_factor(); + for (float beta : beta_array) + for (float input_multiplication : input_multiplication_array) + for (float input_addition : input_addition_array) + for (float current_image_multiplication : current_image_multiplication_array) + for (float current_image_addition : current_image_addition_array) { + if (testOK) // only compute configuration if testOK from previous tests + testOK = test_Hessian_convexity_configuration(test_name, objective_function, target_sptr, + beta, + input_multiplication, input_addition, + current_image_multiplication, current_image_addition); + } + /// Reset beta to original value + objective_function.set_penalisation_factor(initial_beta); +} + +bool +GeneralisedPriorTests:: +test_Hessian_convexity_configuration(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr, + const float beta, + const float input_multiplication, const float input_addition, + const float current_image_multiplication, const float current_image_addition) +{ + /// setup targets + target_type& target(*target_sptr); + shared_ptr output(target.get_empty_copy()); + shared_ptr current_image(target.get_empty_copy()); + shared_ptr input(target.get_empty_copy()); + objective_function.set_penalisation_factor(beta); + { + /// Construct an current_image & input with various values that are variations of the target + target_type::full_iterator input_iter = input->begin_all(); + target_type::full_iterator current_image_iter = current_image->begin_all(); + target_type::full_iterator target_iter = target.begin_all(); + while (input_iter != input->end_all()) { - info("Writing diagnostic files gradient" + test_name + ".hv, numerical_gradient" + test_name + ".hv"); - write_to_file("gradient" + test_name + ".hv", *gradient_sptr); - write_to_file("numerical_gradient" + test_name + ".hv", *gradient_2_sptr); + *input_iter = input_multiplication * *target_iter + input_addition; + *current_image_iter = current_image_multiplication * *target_iter + current_image_addition; + ++input_iter; ++target_iter; ++current_image_iter; } + } + + /// Compute H x + objective_function.accumulate_Hessian_times_input(*output, *current_image, *input); + + /// Compute x \cdot (H x) + float my_sum = 0.0; + my_sum = std::inner_product(input->begin_all(), input->end_all(), output->begin_all(), my_sum); + + // test for a CONVEX function + if (this->check_if_less(0, my_sum)) { +// info("PASS: Computation of x^T H x = " + std::to_string(my_sum) + " > 0 and is therefore convex"); + return true; + } else { + // print to console the FAILED configuration + info("FAIL: Computation of x^T H x = " + std::to_string(my_sum) + " < 0 and is therefore NOT convex" + + "\ntest_name=" + test_name + + "\nbeta=" + std::to_string(beta) + + "\ninput_multiplication=" + std::to_string(input_multiplication) + + "\ninput_addition=" + std::to_string(input_addition) + + "\ncurrent_image_multiplication=" + std::to_string(current_image_multiplication) + + "\ncurrent_image_addition=" + std::to_string(current_image_addition) + + "\n >input image max=" + std::to_string(input->find_max()) + + "\n >input image min=" + std::to_string(input->find_min()) + + "\n >target image max=" + std::to_string(target.find_max()) + + "\n >target image min=" + std::to_string(target.find_min())); + return false; + } +} + + +void +GeneralisedPriorTests:: +test_Hessian_against_numerical(const std::string &test_name, + GeneralisedPrior &objective_function, + shared_ptr target_sptr) +{ + if (!objective_function.is_convex()) + return; + /// Setup + const float eps = 1e-3F; + bool testOK = true; + const int verbosity_default = Verbosity::get(); + + // setup images + target_type& input(*target_sptr->get_empty_copy()); + input += *target_sptr; // make input have same values as target_sptr + shared_ptr gradient_sptr(target_sptr->get_empty_copy()); + shared_ptr pert_grad_and_numerical_Hessian_sptr(target_sptr->get_empty_copy()); + shared_ptr Hessian_sptr(target_sptr->get_empty_copy()); + + Verbosity::set(0); + objective_function.compute_gradient(*gradient_sptr, input); + Verbosity::set(verbosity_default); +// this->set_tolerance(std::max(fabs(double(gradient_sptr->find_min())), fabs(double(gradient_sptr->find_max()))) /10); + + // Setup coordinates (z,y,x) for perturbation test (Hessian will also be computed w.r.t this voxel, j) + BasicCoordinate<3, int> perturbation_coords; + + // Get min/max indices + const int min_z = input.get_min_index(); + const int max_z = input.get_max_index(); + const int min_y = input[min_z].get_min_index(); + const int max_y = input[min_z].get_max_index(); + const int min_x = input[min_z][min_y].get_min_index(); + const int max_x = input[min_z][min_y].get_max_index(); + + // Loop over each voxel j in the input and check perturbation response. + for (int z=min_z;z<= max_z;z++) + for (int y=min_y;y<= max_y;y++) + for (int x=min_x;x<= max_x;x++) + if (testOK) + { + perturbation_coords[1] = z; perturbation_coords[2] = y; perturbation_coords[3] = x; + + // Compute H(x)_j (row of the Hessian at the jth voxel) + objective_function.compute_Hessian(*Hessian_sptr, perturbation_coords, input); + this->set_tolerance(std::max(fabs(double(Hessian_sptr->find_min())), fabs(double(Hessian_sptr->find_max()))) /500); + + // Compute g(x + eps) + Verbosity::set(0); + // Perturb target at jth voxel, compute perturbed gradient, and reset voxel to original value + float perturbed_voxels_original_value = input[perturbation_coords[1]][perturbation_coords[2]][perturbation_coords[3]]; + input[perturbation_coords[1]][perturbation_coords[2]][perturbation_coords[3]] += eps; + objective_function.compute_gradient(*pert_grad_and_numerical_Hessian_sptr, input); + input[perturbation_coords[1]][perturbation_coords[2]][perturbation_coords[3]] = perturbed_voxels_original_value; + + // Now compute the numerical-Hessian = (g(x+eps) - g(x))/eps + *pert_grad_and_numerical_Hessian_sptr -= *gradient_sptr; + *pert_grad_and_numerical_Hessian_sptr /= eps; + + Verbosity::set(verbosity_default); + // Test if pert_grad_and_numerical_Hessian_sptr is all zeros. + // This can happen if the eps is too small. This is a quick test that allows for easier debugging. + if (pert_grad_and_numerical_Hessian_sptr->sum_positive() == 0.0 && Hessian_sptr->sum_positive() > 0.0) + { + this->everything_ok = false; + testOK = false; + info("test_Hessian_against_numerical: failed because all values are 0 in numerical Hessian"); + } + + // Loop over each of the voxels and compare the numerical-Hessian with Hessian + target_type::full_iterator numerical_Hessian_iter = pert_grad_and_numerical_Hessian_sptr->begin_all(); + target_type::full_iterator Hessian_iter = Hessian_sptr->begin_all(); + while(numerical_Hessian_iter != pert_grad_and_numerical_Hessian_sptr->end_all() && testOK) + { + testOK = testOK && this->check_if_equal(*Hessian_iter, *numerical_Hessian_iter, "Hessian"); + ++numerical_Hessian_iter; ++ Hessian_iter; + } + + if (!testOK) + { + // Output volumes for debug + std::cerr << "Numerical-Hessian test failed with for " + test_name + " prior\n"; + info("Writing diagnostic files `Hessian_" + test_name + ".hv` and `numerical_Hessian_" + test_name + ".hv`"); + write_to_file("Hessian_" + test_name + ".hv", *Hessian_sptr); + write_to_file("numerical_Hessian_" + test_name + ".hv", *pert_grad_and_numerical_Hessian_sptr); + write_to_file("input_" + test_name + ".hv", input); + } + } } void @@ -179,30 +464,52 @@ run_tests() shared_ptr density_sptr; construct_input_data(density_sptr); - std::cerr << "Tests for QuadraticPrior\n"; + std::cerr << "\n\nTests for QuadraticPrior\n"; { QuadraticPrior objective_function(false, 1.F); + set_do_test_gradient(true); + set_do_test_Hessian_convexity(true); + set_do_test_Hessian_against_numerical(true); this->run_tests_for_objective_function("Quadratic_no_kappa", objective_function, density_sptr); } - std::cerr << "Tests for Relative Difference Prior\n"; + std::cerr << "\n\nTests for Relative Difference Prior with epsilon = 0\n"; { - // gamma is default and epsilon is off + // gamma is default and epsilon is 0.0 RelativeDifferencePrior objective_function(false, 1.F, 2.F, 0.F); - this->run_tests_for_objective_function("RDP_no_kappa", objective_function, density_sptr); + set_do_test_gradient(true); + set_do_test_Hessian_convexity(true); + set_do_test_Hessian_against_numerical(false); // RDP, with epsilon = 0.0, will fail the numerical Hessian test + this->run_tests_for_objective_function("RDP_no_kappa_no_eps", objective_function, density_sptr); + } + std::cerr << "\n\nTests for Relative Difference Prior with epsilon = 0.1\n"; + { + // gamma is default and epsilon is "small" + RelativeDifferencePrior objective_function(false, 1.F, 2.F, 0.1F); + set_do_test_gradient(true); + set_do_test_Hessian_convexity(true); + set_do_test_Hessian_against_numerical(true); // With a large enough epsilon the RDP numerical test will pass + this->run_tests_for_objective_function("RDP_no_kappa_with_eps", objective_function, density_sptr); + } + + std::cerr << "\n\nTests for PLSPrior\n"; + { + PLSPrior objective_function(false, 1.F); + shared_ptr > anatomical_image_sptr(density_sptr->get_empty_copy()); + anatomical_image_sptr->fill(1.F); + objective_function.set_anatomical_image_sptr(anatomical_image_sptr); + // Disabled PLS due to known issue + set_do_test_gradient(false); + set_do_test_Hessian_convexity(false); + set_do_test_Hessian_against_numerical(false); + this->run_tests_for_objective_function("PLS_no_kappa_flat_anatomical", objective_function, density_sptr); } - // Disabled PLS due to known issue -// std::cerr << "Tests for PLSPrior\n"; -// { -// PLSPrior objective_function(false, 1.F); -// shared_ptr > anatomical_image_sptr(density_sptr->get_empty_copy()); -// anatomical_image_sptr->fill(1.F); -// objective_function.set_anatomical_image_sptr(anatomical_image_sptr); -// this->run_tests_for_objective_function("PLS_no_kappa_flat_anatomical", objective_function, density_sptr); -// } - std::cerr << "Tests for Logcosh Prior\n"; + std::cerr << "\n\nTests for Logcosh Prior\n"; { // scalar is off LogcoshPrior objective_function(false, 1.F, 1.F); + set_do_test_gradient(true); + set_do_test_Hessian_convexity(true); + set_do_test_Hessian_against_numerical(true); this->run_tests_for_objective_function("Logcosh_no_kappa", objective_function, density_sptr); } }