Merge pull request #696 from robbietuk/RDP_grad_fix

Fix an issue with the gradient of the relative difference prior (also change default values)

examples/samples/OSMAPOSL_RDP.par

@@ -0,0 +1,59 @@
+OSMAPOSLParameters :=
+; example file for using a relative difference prior (with One Step Late and subsets)
+; this is a minimal file. check other sample files and doc for more options
+
+objective function type:= PoissonLogLikelihoodWithLinearModelForMeanAndProjData
+PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:=
+
+input file := test.hs
+; if disabled, defaults to maximum segment number in the file
+maximum absolute segment number to process := 4
+; see User's Guide to see when you need this
+zero end planes of segment 0:= 1
+
+; if the next parameter is disabled, 
+; the sensitivity will be computed
+; use %d where you want the subset-number (a la printf)
+subset sensitivity filenames:= sens_%d.hv
+; specify additive projection data to handle randoms or so
+; see User's Guide for more info
+additive sinogram := 0
+
+
+; here comes the prior stuff
+
+
+prior type := Relative Difference Prior
+  Relative Difference Prior Parameters:=
+  penalisation factor := 1.0
+  ; next can be used to set weights explicitly. Needs to be a 3D array (of floats). value of only_2D is ignored
+  ; following example uses 2D 'nearest neighbour' penalty
+  ; weights:={{{0,1,0},{1,0,1},{0,1,0}}}
+  ; use next parameter to specify an image with penalisation factors (a la Fessler)
+  ; see class documentation for more info
+  ; kappa filename:=
+  ; gamma value is a RDP specific parameter and controls the edge preservation. Default value is 2
+  ; gamma value := 
+  ; epsilon value - a small non-negative float to smooth the prior at very small values. Default value is 0
+  ; epsilon value :=
+END Relative Difference Prior Parameters:=
+
+end PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:=
+
+initial estimate:= some_image
+; enable this when you read an initial estimate with negative data
+enforce initial positivity condition:=0
+
+output filename prefix := test_RDP
+number of subsets:= 12
+number of subiterations:= 24
+save estimates at subiteration intervals:= 12
+
+; This is the default verbosity, just putting it here as an example
+; lower value means less output, 0 being the minimum
+verbosity := 2
+
+; enable this for multiplicative form of OSMAPOSL (see User's Guide)
+;MAP model := multiplicative
+
+END :=

src/include/stir/recon_buildblock/RelativeDifferencePrior.h

@@ -59,8 +59,9 @@ START_NAMESPACE_STIR
 
   where \f$\lambda\f$ is the image and \f$r\f$ and \f$dr\f$ are indices and the sum
   is over the neighbourhood where the weights \f$w_{dr}\f$ are non-zero. \f$\gamma\f$ is
-  a smoothing scalar term and the \f$\epsilon\f$ is a small positive value included to prevent division by zero.
-  For more details, see: <em> J. Nuyts, D. Bequ&eacute;, P. Dupont, and L. Mortelmans,
+  a smoothing scalar term and the \f$\epsilon\f$ is a small non-negative value included to prevent division by zero.
+  Please note that the RDP is only well defined for non-negative voxel values.
+  For more details, see: <em> J. Nuyts, D. Beque, P. Dupont, and L. Mortelmans,
   “A Concave Prior Penalizing Relative Differences for Maximum-a-Posteriori Reconstruction in Emission Tomography,”
   vol. 49, no. 1, pp. 56–60, 2002. </em>
 
@@ -69,7 +70,7 @@ START_NAMESPACE_STIR
   penalty is computed. If \f$\kappa\f$ is not set, this class will effectively
   use 1 for all \f$\kappa\f$'s.
 
-  By default, a 3x3 or 3x3x3 neigbourhood is used where the weights are set to 
+  By default, a 3x3 or 3x3x3 neighbourhood is used where the weights are set to
   x-voxel_size divided by the Euclidean distance between the points.
 
 \par Parsing

src/recon_buildblock/RelativeDifferencePrior.cxx

@@ -1,4 +1,4 @@
-//
+//
 //
 /*
     Copyright (C) 2000- 2019, Hammersmith Imanet Ltd
@@ -138,8 +138,8 @@ RelativeDifferencePrior<elemT>::set_defaults()
   this->only_2D = false;
   this->kappa_ptr.reset();  
   this->weights.recycle();
-  this->gamma = 1;
-  this->epsilon = 0.0001;
+  this->gamma = 2;
+  this->epsilon = 0.0;
 }
 
 template <>
@@ -192,14 +192,14 @@ RelativeDifferencePrior<elemT>::RelativeDifferencePrior(const bool only_2D_v, fl
 }
 
 
-  //! get penalty weights for the neigbourhood
+  //! get penalty weights for the neighbourhood
 template <typename elemT>
 Array<3,float>  
 RelativeDifferencePrior<elemT>::
 get_weights() const
 { return this->weights; }
 
-  //! set penalty weights for the neigbourhood
+  //! set penalty weights for the neighbourhood
 template <typename elemT>
 void 
 RelativeDifferencePrior<elemT>::
@@ -208,7 +208,7 @@ set_weights(const Array<3,float>& w)
 
   //! get current kappa image
   /*! \warning As this function returns a shared_ptr, this is dangerous. You should not
-      modify the image by manipulating the image refered to by this pointer.
+      modify the image by manipulating the image referred to by this pointer.
       Unpredictable results will occur.
   */
 template <typename elemT>
@@ -313,10 +313,15 @@ compute_value(const DiscretisedDensity<3,elemT> &current_image_estimate)
                     for (int dx=min_dx;dx<=max_dx;++dx)
                       {
                         elemT current;
+                        if (this->epsilon ==0.0 && current_image_estimate[z][y][x] == 0.0 && current_image_estimate[z+dz][y+dy][x+dx] == 0.0){
+                          // handle the undefined nature of the function
+                          current = 0.0;
+                        } else {
                           current = weights[dz][dy][dx] * 0.5 *
                                   (pow(current_image_estimate[z][y][x]-current_image_estimate[z+dz][y+dy][x+dx],2)/
                                   (current_image_estimate[z][y][x]+current_image_estimate[z+dz][y+dy][x+dx]
                                   + this->gamma * abs(current_image_estimate[z][y][x]-current_image_estimate[z+dz][y+dy][x+dx]) + this->epsilon ));
+                        }
                         if (do_kappa)
                           current *= 
                             (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx];
@@ -388,12 +393,17 @@ compute_gradient(DiscretisedDensity<3,elemT>& prior_gradient,
                       {
 
                         elemT current;
+                        if (this->epsilon ==0.0 && current_image_estimate[z][y][x] == 0.0 && current_image_estimate[z+dz][y+dy][x+dx] == 0.0){
+                          // handle the undefined nature of the gradient
+                          current = 0.0;
+                        } else {
                             current = weights[dz][dy][dx] *
                                     (((current_image_estimate[z][y][x] - current_image_estimate[z+dz][y+dy][x+dx]) *
                                       (this->gamma * abs(current_image_estimate[z][y][x] - current_image_estimate[z+dz][y+dy][x+dx]) +
                                        current_image_estimate[z][y][x] + 3 * current_image_estimate[z+dz][y+dy][x+dx] + 2 * this->epsilon))/
                                      (square((current_image_estimate[z][y][x] + current_image_estimate[z+dz][y+dy][x+dx]) +
-                                      this->gamma * abs(current_image_estimate[z][y][x] - current_image_estimate[z+dz][y+dy][x+dx])) + this->epsilon));
+                                      this->gamma * abs(current_image_estimate[z][y][x] - current_image_estimate[z+dz][y+dy][x+dx]) + this->epsilon)));
+                        }
                         if (do_kappa)
                           current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx];