From 3f0b210af1cd01e1908bf17eea4eba3a8acb2634 Mon Sep 17 00:00:00 2001 From: Markus Jehl Date: Mon, 25 Jul 2022 11:19:05 +0200 Subject: [PATCH] Revert "Merge pull request #901 from robbietuk/RDP/HessianVectorProduct" This reverts commit 4d087f3e7cf3c6d2d95b87f64e9e68bf9f1b186e, reversing changes made to f98246f96f939e6a0c725d31205a995ecb395264. --- documentation/release_5.1.htm | 10 +- .../stir/recon_buildblock/FilterRootPrior.h | 2 - .../stir/recon_buildblock/GeneralisedPrior.h | 21 +- .../stir/recon_buildblock/LogcoshPrior.h | 40 +- src/include/stir/recon_buildblock/PLSPrior.h | 2 - .../stir/recon_buildblock/QuadraticPrior.h | 27 +- .../RelativeDifferencePrior.h | 28 +- src/recon_buildblock/FilterRootPrior.cxx | 6 - .../GeneralisedObjectiveFunction.cxx | 9 +- src/recon_buildblock/GeneralisedPrior.cxx | 19 +- src/recon_buildblock/LogcoshPrior.cxx | 85 +---- src/recon_buildblock/PLSPrior.cxx | 7 - src/recon_buildblock/QuadraticPrior.cxx | 157 ++------ .../RelativeDifferencePrior.cxx | 205 +--------- src/recon_test/test_priors.cxx | 354 ++---------------- 15 files changed, 117 insertions(+), 855 deletions(-) diff --git a/documentation/release_5.1.htm b/documentation/release_5.1.htm index e76c54f41b..30783fb131 100644 --- a/documentation/release_5.1.htm +++ b/documentation/release_5.1.htm @@ -25,8 +25,7 @@

Overall summary

improvements to the documentation.

-

This release contains mainly code written by Nikos Eftimiou (UPenn and MGH), with support by Kris Thielemans (UCL). - Additional changes to penalty functions were made by Robert Twyman (UCL). +

This release contains mainly code written by Nikos Eftimiou (UPenn and MGH), with support by Kris Thielemans (UCL)

Patch release info

@@ -44,8 +43,6 @@

New functionality

  • Support for PENNPET Explorer listmode data (if proprietary libraries are found) by Nikos Efthimiou, see PR #1028.
    • -
  • Improvements to penalty functions including Hessian methods, see PR #901. -
    • @@ -69,11 +66,6 @@

    Documentation changes

    recon_test_pack changes

    Other changes to tests

    -
  • Significant changes made to test_priors that tests the Hessian's convexity, given - x(Hx) > 0, and a perturbation response, using gradients, was added to determine the Hessian - (for a single densel) is within tolerance. - Tests the Quadratic, Relative Difference (in two configurations) and Log-Cosh Penalties (Robert Twyman, UCL). -
    • diff --git a/src/include/stir/recon_buildblock/FilterRootPrior.h b/src/include/stir/recon_buildblock/FilterRootPrior.h index e2850f9ab3..979bc0c496 100644 --- a/src/include/stir/recon_buildblock/FilterRootPrior.h +++ b/src/include/stir/recon_buildblock/FilterRootPrior.h @@ -87,8 +87,6 @@ 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/GeneralisedPrior.h b/src/include/stir/recon_buildblock/GeneralisedPrior.h index a68bb04539..caa4bf4dab 100644 --- a/src/include/stir/recon_buildblock/GeneralisedPrior.h +++ b/src/include/stir/recon_buildblock/GeneralisedPrior.h @@ -60,19 +60,6 @@ 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. @@ -80,7 +67,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 void + virtual Succeeded add_multiplication_with_approximate_Hessian(DataT& output, const DataT& input) const; @@ -89,7 +76,7 @@ class GeneralisedPrior: derived class. \warning The derived class should accumulate in \a output. */ - virtual void + virtual Succeeded accumulate_Hessian_times_input(DataT& output, const DataT& current_estimate, const DataT& input) const; @@ -102,10 +89,6 @@ 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 fb9ed5aa4b..1593681ee9 100644 --- a/src/include/stir/recon_buildblock/LogcoshPrior.h +++ b/src/include/stir/recon_buildblock/LogcoshPrior.h @@ -108,12 +108,6 @@ class LogcoshPrior: public //! Constructs it explicitly with scalar LogcoshPrior(const bool only_2D, float penalization_factor, const float scalar); - virtual bool - 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); @@ -126,13 +120,13 @@ class LogcoshPrior: public void parabolic_surrogate_curvature(DiscretisedDensity<3,elemT>& parabolic_surrogate_curvature, 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 Hessian + void compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, + const BasicCoordinate<3,int>& coords, + const DiscretisedDensity<3,elemT> ¤t_image_estimate); - //! 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, + //! Compute the multiplication of the hessian of the prior multiplied by the input. + virtual Succeeded accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, const DiscretisedDensity<3,elemT>& current_estimate, const DiscretisedDensity<3,elemT>& input) const; @@ -226,21 +220,21 @@ class LogcoshPrior: public { return tanh(x)/x; } } - //! The second partial derivatives of the LogCosh Prior + //! The Hessian of log(cosh()) is sech^2(x) = (1/cosh(x))^2 /*! - 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 + 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 */ - //@{ - elemT derivative_20(const elemT x_j, const elemT x_k) const; - elemT derivative_11(const elemT x_j, const elemT x_k) const; - //@} + static inline float Hessian(const float d, const float scalar) + { + const float x = d * scalar; + return square((1/ cosh(x))); + } }; END_NAMESPACE_STIR -#endif +#endif \ No newline at end of file diff --git a/src/include/stir/recon_buildblock/PLSPrior.h b/src/include/stir/recon_buildblock/PLSPrior.h index 0c1458abf8..91a8d83c0c 100644 --- a/src/include/stir/recon_buildblock/PLSPrior.h +++ b/src/include/stir/recon_buildblock/PLSPrior.h @@ -121,8 +121,6 @@ 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/QuadraticPrior.h b/src/include/stir/recon_buildblock/QuadraticPrior.h index 1a328fa802..4028169276 100644 --- a/src/include/stir/recon_buildblock/QuadraticPrior.h +++ b/src/include/stir/recon_buildblock/QuadraticPrior.h @@ -103,8 +103,6 @@ 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); @@ -118,20 +116,20 @@ class QuadraticPrior: public void parabolic_surrogate_curvature(DiscretisedDensity<3,elemT>& parabolic_surrogate_curvature, 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 Hessian + void compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, + const BasicCoordinate<3,int>& coords, + const DiscretisedDensity<3,elemT> ¤t_image_estimate); //! Call accumulate_Hessian_times_input - virtual void + virtual Succeeded 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 void accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, + virtual Succeeded accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, const DiscretisedDensity<3,elemT>& current_estimate, const DiscretisedDensity<3,elemT>& input) const; @@ -181,19 +179,6 @@ 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 83a3fb08a0..947f1fa78f 100644 --- a/src/include/stir/recon_buildblock/RelativeDifferencePrior.h +++ b/src/include/stir/recon_buildblock/RelativeDifferencePrior.h @@ -121,19 +121,11 @@ 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 void + virtual Succeeded 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 @@ -164,8 +156,6 @@ 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; @@ -199,22 +189,6 @@ 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/recon_buildblock/FilterRootPrior.cxx b/src/recon_buildblock/FilterRootPrior.cxx index 638bdcd067..90ebaa8ece 100644 --- a/src/recon_buildblock/FilterRootPrior.cxx +++ b/src/recon_buildblock/FilterRootPrior.cxx @@ -41,12 +41,6 @@ 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/GeneralisedObjectiveFunction.cxx b/src/recon_buildblock/GeneralisedObjectiveFunction.cxx index f3490fa5f7..41b842c2ed 100644 --- a/src/recon_buildblock/GeneralisedObjectiveFunction.cxx +++ b/src/recon_buildblock/GeneralisedObjectiveFunction.cxx @@ -263,7 +263,10 @@ add_multiplication_with_approximate_sub_Hessian(TargetT& output, { // TODO used boost:scoped_ptr shared_ptr prior_output_sptr(output.get_empty_copy()); - this->prior_sptr->add_multiplication_with_approximate_Hessian(*prior_output_sptr, output); + if (this->prior_sptr->add_multiplication_with_approximate_Hessian(*prior_output_sptr, output) == + Succeeded::no) + return Succeeded::no; + // (*prior_output_sptr)/= num_subsets; // output -= *prior_output_sptr; @@ -380,7 +383,9 @@ accumulate_sub_Hessian_times_input(TargetT& output, { // TODO used boost:scoped_ptr shared_ptr prior_output_sptr(output.get_empty_copy()); - this->prior_sptr->accumulate_Hessian_times_input(*prior_output_sptr, current_image_estimate, output); + if (this->prior_sptr->accumulate_Hessian_times_input(*prior_output_sptr, current_image_estimate, output) == + Succeeded::no) + return Succeeded::no; 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 b1e5b44f71..5751c0a6b5 100644 --- a/src/recon_buildblock/GeneralisedPrior.cxx +++ b/src/recon_buildblock/GeneralisedPrior.cxx @@ -52,30 +52,18 @@ set_up(shared_ptr const&) } template -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 for this (non-convex) prior."); -} - -template -void +Succeeded 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 -void +Succeeded GeneralisedPrior:: accumulate_Hessian_times_input(TargetT& output, const TargetT& current_estimate, @@ -83,6 +71,7 @@ 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 5324498430..f14742d7b1 100644 --- a/src/recon_buildblock/LogcoshPrior.cxx +++ b/src/recon_buildblock/LogcoshPrior.cxx @@ -144,14 +144,6 @@ 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> @@ -388,7 +380,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 + const DiscretisedDensity<3,elemT> ¤t_image_estimate) { assert( prior_Hessian_for_single_densel.has_same_characteristics(current_image_estimate)); prior_Hessian_for_single_densel.fill(0); @@ -396,9 +388,6 @@ 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); @@ -432,34 +421,19 @@ 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) { - 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; + // 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; } -} + prior_Hessian_for_single_densel[z][y][x]= diagonal * this->penalisation_factor; +} template void @@ -531,7 +505,7 @@ LogcoshPrior::parabolic_surrogate_curvature(DiscretisedDensity<3,elemT>& } template -void +Succeeded LogcoshPrior:: accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, const DiscretisedDensity<3,elemT>& current_estimate, @@ -543,7 +517,7 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, assert( output.has_same_characteristics(input)); if (this->penalisation_factor==0) { - return; + return Succeeded::yes; } DiscretisedDensityOnCartesianGrid<3,elemT>& output_cast = @@ -587,18 +561,8 @@ 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 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]); - } + 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]; if (do_kappa) current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; @@ -610,22 +574,7 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, } } } -} - -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); + return Succeeded::yes; } # ifdef _MSC_VER diff --git a/src/recon_buildblock/PLSPrior.cxx b/src/recon_buildblock/PLSPrior.cxx index 891d554fd5..0b782546fa 100644 --- a/src/recon_buildblock/PLSPrior.cxx +++ b/src/recon_buildblock/PLSPrior.cxx @@ -144,13 +144,6 @@ 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 7e1260ff41..904b4b0d1b 100644 --- a/src/recon_buildblock/QuadraticPrior.cxx +++ b/src/recon_buildblock/QuadraticPrior.cxx @@ -147,13 +147,6 @@ 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 @@ -437,11 +430,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 + const DiscretisedDensity<3,elemT> ¤t_image_estimate) { assert( prior_Hessian_for_single_densel.has_same_characteristics(current_image_estimate)); prior_Hessian_for_single_densel.fill(0); @@ -486,32 +479,19 @@ 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) { - 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; + // 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; } + + prior_Hessian_for_single_densel[z][y][x]= diagonal * this->penalisation_factor; } template @@ -596,83 +576,19 @@ QuadraticPrior::parabolic_surrogate_curvature(DiscretisedDensity<3,elemT> } template -void +Succeeded QuadraticPrior:: add_multiplication_with_approximate_Hessian(DiscretisedDensity<3,elemT>& output, 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; - } - - 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; - } - } - } + return accumulate_Hessian_times_input(output, input, input); } template -void +Succeeded 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 @@ -681,7 +597,7 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, assert( output.has_same_characteristics(input)); if (this->penalisation_factor==0) { - return; + return Succeeded::yes; } this->check(input); @@ -722,31 +638,13 @@ 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]; - 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]); - } + elemT current = + weights[dz][dy][dx] * input[z+dz][y+dy][x+dx]; if (do_kappa) current *= @@ -758,20 +656,9 @@ 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 679a701ba6..e7fce65ea5 100644 --- a/src/recon_buildblock/RelativeDifferencePrior.cxx +++ b/src/recon_buildblock/RelativeDifferencePrior.cxx @@ -127,7 +127,6 @@ void RelativeDifferencePrior::set_defaults() { base_type::set_defaults(); -// this->_is_convex = true; this->only_2D = false; this->kappa_ptr.reset(); this->weights.recycle(); @@ -146,14 +145,6 @@ RelativeDifferencePrior::RelativeDifferencePrior() set_defaults(); } -template -bool -RelativeDifferencePrior:: -is_convex() const -{ - return true; -} - // Return the value of gamma - a RDP parameter template float @@ -187,7 +178,6 @@ 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; @@ -431,204 +421,13 @@ compute_gradient(DiscretisedDensity<3,elemT>& prior_gradient, } template -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 cases (off-diagonal Hessian elements), no summing over neighbourhood - 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 +Succeeded 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."); -} - -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("RelativeDifferencePrior: 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; + return Succeeded::no; } # ifdef _MSC_VER diff --git a/src/recon_test/test_priors.cxx b/src/recon_test/test_priors.cxx index 2a27ca0d69..961afdf0a8 100644 --- a/src/recon_test/test_priors.cxx +++ b/src/recon_test/test_priors.cxx @@ -29,12 +29,11 @@ #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 @@ -42,6 +41,7 @@ #include #include #include +#include START_NAMESPACE_STIR @@ -53,8 +53,6 @@ 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 @@ -67,15 +65,11 @@ class GeneralisedPriorTests : public RunTests \todo it would be better to parse an objective function. That would allow us to set all parameters from the command line. */ - explicit GeneralisedPriorTests(char const * density_filename = nullptr); + GeneralisedPriorTests(char const * const density_filename = 0); typedef DiscretisedDensity<3,float> target_type; void construct_input_data(shared_ptr& density_sptr); - void run_tests() override; - - //! Set methods that control which tests are run. - void configure_prior_tests(bool gradient, bool Hessian_convexity, bool Hessian_numerical); - + void run_tests(); protected: char const * density_filename; shared_ptr > objective_function_sptr; @@ -84,48 +78,7 @@ class GeneralisedPriorTests : public RunTests /*! Note that this function is not specific to a particular prior */ void run_tests_for_objective_function(const std::string& test_name, GeneralisedPrior& objective_function, - const 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, - const 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, - const 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, - const 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, - const shared_ptr& target_sptr, - float beta, - float input_multiplication, float input_addition, - float current_image_multiplication, 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; - //@} + shared_ptr target_sptr); }; GeneralisedPriorTests:: @@ -133,283 +86,63 @@ GeneralisedPriorTests(char const * const density_filename) : density_filename(density_filename) {} -void -GeneralisedPriorTests::configure_prior_tests(const bool gradient, const bool Hessian_convexity, const bool Hessian_numerical) -{ - do_test_gradient = gradient; - do_test_Hessian_convexity = Hessian_convexity; - do_test_Hessian_against_numerical = Hessian_numerical; -} - void GeneralisedPriorTests:: run_tests_for_objective_function(const std::string& test_name, GeneralisedPrior& objective_function, - const shared_ptr& target_sptr) + shared_ptr target_sptr) { std::cerr << "----- test " << test_name << '\n'; 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, - const 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 auto 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; + 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; + } } - } 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, - const 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, - const 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()) { - *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; + 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); } - } - - /// 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, - const 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 && 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 GeneralisedPriorTests:: construct_input_data(shared_ptr& density_sptr) { - if (this->density_filename == nullptr) + if (this->density_filename == 0) { // construct a small image with random voxel values between 0 and 1 @@ -436,6 +169,7 @@ construct_input_data(shared_ptr& density_sptr) shared_ptr aptr(read_from_file(this->density_filename)); density_sptr = aptr; } + return; } void @@ -445,42 +179,30 @@ run_tests() shared_ptr density_sptr; construct_input_data(density_sptr); - std::cerr << "\n\nTests for QuadraticPrior\n"; + std::cerr << "Tests for QuadraticPrior\n"; { QuadraticPrior objective_function(false, 1.F); - this->configure_prior_tests(true, true, true); this->run_tests_for_objective_function("Quadratic_no_kappa", objective_function, density_sptr); } - std::cerr << "\n\nTests for Relative Difference Prior with epsilon = 0\n"; + std::cerr << "Tests for Relative Difference Prior\n"; { - // gamma is default and epsilon is 0.0 + // gamma is default and epsilon is off RelativeDifferencePrior objective_function(false, 1.F, 2.F, 0.F); - this->configure_prior_tests(true, true, 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); - this->configure_prior_tests(true, true, true); // With a large enough epsilon the RDP Hessian 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 - this->configure_prior_tests(false, false, false); - this->run_tests_for_objective_function("PLS_no_kappa_flat_anatomical", objective_function, density_sptr); + this->run_tests_for_objective_function("RDP_no_kappa", objective_function, density_sptr); } - std::cerr << "\n\nTests for Logcosh Prior\n"; + // 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"; { // scalar is off LogcoshPrior objective_function(false, 1.F, 1.F); - this->configure_prior_tests(true, true, true); this->run_tests_for_objective_function("Logcosh_no_kappa", objective_function, density_sptr); } } @@ -494,7 +216,7 @@ int main(int argc, char **argv) { set_default_num_threads(); - GeneralisedPriorTests tests(argc>1? argv[1] : nullptr); + GeneralisedPriorTests tests(argc>1? argv[1] : 0); tests.run_tests(); return tests.main_return_value(); }