Priors/Hessian method corrections

documentation/release_5.1.htm

@@ -25,7 +25,8 @@ <h2>Overall summary</h2>
 improvements to the documentation.
 </p>
 
-<p>This release contains mainly code written by Nikos Eftimiou (UPenn and MGH), with support by Kris Thielemans (UCL)
+<p>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).
 </p>
 
 <h2>Patch release info</h2>
@@ -43,6 +44,8 @@ <h3>New functionality</h3>
   </li>    
   <li>Support for PENNPET Explorer listmode data (if proprietary libraries are found) by Nikos Efthimiou, see <a href="https://github.com/UCL/STIR/pull/1028/">PR #1028</a>. 
   </li>    
+  <li> Improvements to penalty functions including Hessian methods, see <a href="https://github.com/UCL/STIR/pull/901/">PR #901</a>.
+  </li>
 </ul>
 
 
@@ -66,6 +69,11 @@ <h3>Documentation changes</h3>
 <h3>recon_test_pack changes</h3>
 
 <h3>Other changes to tests</h3>
+  <li> Significant changes made to <code>test_priors</code> that tests the Hessian's convexity, given
+    <code>x(Hx) > 0</code>, 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).
+  </li>
 
 </body>
 

src/include/stir/recon_buildblock/FilterRootPrior.h

@@ -87,6 +87,8 @@ class FilterRootPrior:  public
   FilterRootPrior(shared_ptr<DataProcessor<DataT> >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.

src/include/stir/recon_buildblock/GeneralisedPrior.h

@@ -60,14 +60,27 @@ class GeneralisedPrior:
   virtual void compute_gradient(DataT& prior_gradient, 
 		   const DataT &current_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.
       This method assumes that the hessian of the prior is 1 and hence the function quadratic.
       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 DataT> 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

src/include/stir/recon_buildblock/LogcoshPrior.h

@@ -108,6 +108,12 @@ 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> &current_image_estimate);
@@ -120,13 +126,13 @@ class LogcoshPrior:  public
     void parabolic_surrogate_curvature(DiscretisedDensity<3,elemT>& parabolic_surrogate_curvature,
                                        const DiscretisedDensity<3,elemT> &current_image_estimate);
 
-    //! compute Hessian
-    void compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel,
-                         const BasicCoordinate<3,int>& coords,
-                         const DiscretisedDensity<3,elemT> &current_image_estimate);
+    virtual void
+    compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel,
+                    const BasicCoordinate<3,int>& coords,
+                    const DiscretisedDensity<3,elemT> &current_image_estimate) const;
 
-    //! 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;
 
@@ -220,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
+#endif

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 DiscretisedDensity<3,elemT> > const& target_sptr);
 
+  bool is_convex() const;
+
   //! compute the value of the function
   double
     compute_value(const DiscretisedDensity<3,elemT> &current_image_estimate);

src/include/stir/recon_buildblock/QuadraticPrior.h

@@ -103,6 +103,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> &current_image_estimate);
@@ -116,20 +118,20 @@ class QuadraticPrior:  public
   void parabolic_surrogate_curvature(DiscretisedDensity<3,elemT>& parabolic_surrogate_curvature, 
                         const DiscretisedDensity<3,elemT> &current_image_estimate);
 
-  //! compute Hessian 
-  void compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, 
-                const BasicCoordinate<3,int>& coords,
-                const DiscretisedDensity<3,elemT> &current_image_estimate);
+  virtual void
+  compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel,
+                  const BasicCoordinate<3,int>& coords,
+                  const DiscretisedDensity<3,elemT> &current_image_estimate) const;
 
   //! 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 +181,19 @@ class QuadraticPrior:  public
   virtual bool post_processing();
  private:
   shared_ptr<const DiscretisedDensity<3,elemT> > 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;
+  //@}
 };
 
 

src/include/stir/recon_buildblock/RelativeDifferencePrior.h

@@ -121,11 +121,19 @@ class RelativeDifferencePrior:  public
   void compute_gradient(DiscretisedDensity<3,elemT>& prior_gradient, 
                         const DiscretisedDensity<3,elemT> &current_image_estimate);
 
+  virtual void compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel,
+                                    const BasicCoordinate<3,int>& coords,
+                                    const DiscretisedDensity<3,elemT> &current_image_estimate) const;
 
-  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
@@ -156,6 +164,8 @@ class RelativeDifferencePrior:  public
   //! Has to be called before using this object
   virtual Succeeded set_up(shared_ptr<DiscretisedDensity<3,elemT> > const& target_sptr);
 
+  bool is_convex() const;
+
 protected:
   //! Create variable gamma for Relative Difference Penalty
   float gamma;
@@ -189,6 +199,22 @@ class RelativeDifferencePrior:  public
   virtual bool post_processing();
  private:
   shared_ptr<DiscretisedDensity<3,elemT> > 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;
+  //@}
 };
 
 

src/recon_buildblock/FilterRootPrior.cxx

@@ -41,6 +41,12 @@ FilterRootPrior(shared_ptr<DataProcessor<DataT> >const& filter_sptr, float penal
   this->penalisation_factor = penalisation_factor_v;
 }
 
+template <typename DataT>
+bool FilterRootPrior<DataT>::
+is_convex() const
+{
+  return false;
+}
 
 template < class T>
 static inline int

src/recon_buildblock/GeneralisedObjectiveFunction.cxx

@@ -263,10 +263,7 @@ add_multiplication_with_approximate_sub_Hessian(TargetT& output,
     {
       // TODO used boost:scoped_ptr
       shared_ptr<TargetT>  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<TargetT>  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();

src/recon_buildblock/GeneralisedPrior.cxx

@@ -52,26 +52,37 @@ set_up(shared_ptr<const TargetT> const&)
 }
 
 template <typename TargetT>
-Succeeded 
+void
+GeneralisedPrior<TargetT>::
+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 <typename TargetT>
+void
 GeneralisedPrior<TargetT>::
 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 <typename TargetT>
-Succeeded
+void
 GeneralisedPrior<TargetT>::
 accumulate_Hessian_times_input(TargetT& output,
         const TargetT& current_estimate,
         const TargetT& input) const
 {
   error("GeneralisedPrior:\n"
         "accumulate_Hessian_times_input implementation is not overloaded by your prior.");
-  return Succeeded::no;
 }
 
 template <typename TargetT>