Zeljko custom branch

src/buildblock/recon_array_functions.cxx

@@ -374,12 +374,15 @@ void accumulate_loglikelihood(Viewgram<float>& projection_data,
           //  std::cerr << "Zero at " << r << ", " << b <<'\n';
 	  const float new_estimate =
 	    max(estimated_projections[r][b], 
-		projection_data[r][b]/max_quotient);
-	  if (projection_data[r][b]<=small_value)
-	    sub_result += - double(new_estimate);
-	  else
-	    sub_result += projection_data[r][b]*log(double(new_estimate)) - double(new_estimate);
-	}
+			projection_data[r][b]/max_quotient);
+
+		if (projection_data[r][b] > 0.0)
+		{
+			sub_result += - (projection_data[r][b] * log(projection_data[r][b] / double(new_estimate)) + double(new_estimate) - projection_data[r][b]);
+		} else {
+			sub_result += - (double(new_estimate));
+		}
+  }
     result += sub_result;
   }
 

src/include/stir/recon_buildblock/LogcoshPrior.h

@@ -129,6 +129,11 @@ class LogcoshPrior:  public
                          const BasicCoordinate<3,int>& coords,
                          const DiscretisedDensity<3,elemT> &current_image_estimate);
 
+    //! compute the parabolic surrogate curvature for the prior of the current image estimate and multiply by input image
+    void parabolic_surrogate_curvature_times_input(DiscretisedDensity<3,elemT>& output,
+                                                   const DiscretisedDensity<3,elemT> &current_image_estimate,
+                                                   const DiscretisedDensity<3,elemT> &input_image);
+
     //! 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,

src/include/stir/recon_buildblock/PriorWithParabolicSurrogate.h

@@ -47,7 +47,19 @@ class PriorWithParabolicSurrogate:
     parabolic_surrogate_curvature(TargetT& parabolic_surrogate_curvature, 
 				  const TargetT &current_estimate) = 0;
 
-  //! A function that allows skipping some computations if the curvature is independent of the \c current_estimate
+  //! This will calculate the parabolic surrogate curvature of the current image estimate multiplied by the input image
+  /*!
+    Function is comparable to that of accumulate_Hessian_times_input() but instead of the Hessian, we use the parabolic
+    surrogate curvature.
+    For each voxel (f_{j}) of the image, this method computes
+    \beta \sum_{i\in J} d_{i,j} f_{i}
+    where J is the neighbourhood about voxel j and d_{i,j} is the surrogate function between j and i.
+   */
+  virtual void parabolic_surrogate_curvature_times_input(TargetT& output,
+                                                         const TargetT& current_image_estimate,
+                                                         const TargetT& input_image) = 0;
+
+    //! A function that allows skipping some computations if the curvature is independent of the \c current_estimate
   /*! Defaults to return \c true, but can be overloaded by the derived class.
    */
   virtual bool

src/include/stir/recon_buildblock/QuadraticPrior.h

@@ -121,7 +121,13 @@ class QuadraticPrior:  public
                 const BasicCoordinate<3,int>& coords,
                 const DiscretisedDensity<3,elemT> &current_image_estimate);
 
-  //! Call accumulate_Hessian_times_input
+  //! compute the parabolic surrogate curvature for the prior of the current image estimate and multiply by input image
+  void parabolic_surrogate_curvature_times_input(DiscretisedDensity<3,elemT>& output,
+                                                 const DiscretisedDensity<3,elemT> &current_image_estimate,
+                                                 const DiscretisedDensity<3,elemT> &input_image);
+
+
+    //! Call accumulate_Hessian_times_input
   virtual Succeeded 
     add_multiplication_with_approximate_Hessian(DiscretisedDensity<3,elemT>& output,
                                                 const DiscretisedDensity<3,elemT>& input) const;

src/recon_buildblock/LogcoshPrior.cxx

@@ -515,6 +515,7 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output,
   // the only difference is that parabolic_surrogate_curvature uses input==1
 
   assert( output.has_same_characteristics(input));
+  output.fill(0.f);
   if (this->penalisation_factor==0)
   {
     return Succeeded::yes;
@@ -577,6 +578,79 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output,
   return Succeeded::yes;
 }
 
+template <typename elemT>
+void
+LogcoshPrior<elemT>::
+parabolic_surrogate_curvature_times_input(DiscretisedDensity<3,elemT>& output,
+                                          const DiscretisedDensity<3,elemT> &current_image_estimate,
+                                          const DiscretisedDensity<3,elemT> &input_image)
+{
+  assert( output.has_same_characteristics(input));
+  assert( output.has_same_characteristics(current_image_estimate));
+  output.fill(0.f);
+  if (this->penalisation_factor==0)
+  {
+    return;
+  }
+
+  DiscretisedDensityOnCartesianGrid<3,elemT>& output_cast =
+          dynamic_cast<DiscretisedDensityOnCartesianGrid<3,elemT> &>(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_image))
+    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);
+
+        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 voxel_diff= current_image_estimate[z][y][x] - current_image_estimate[z+dz][y+dy][x+dx];
+              elemT current = weights[dz][dy][dx] *
+                      surrogate(voxel_diff, this->scalar) * input_image[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;
+}
+
 #  ifdef _MSC_VER
     // prevent warning message on reinstantiation,
 // note that we get a linking error if we don't have the explicit instantiation below

src/recon_buildblock/QuadraticPrior.cxx

@@ -659,6 +659,78 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output,
   return Succeeded::yes;
 }
 
+template <typename elemT>
+void
+QuadraticPrior<elemT>::
+parabolic_surrogate_curvature_times_input(DiscretisedDensity<3,elemT>& output,
+                                          const DiscretisedDensity<3,elemT> &current_image_estimate,
+                                          const DiscretisedDensity<3,elemT> &input_image)
+{
+  assert( output.has_same_characteristics(input));
+  assert( output.has_same_characteristics(current_image_estimate));
+  output.fill(0.f);
+  if (this->penalisation_factor==0)
+  {
+    return;
+  }
+
+  DiscretisedDensityOnCartesianGrid<3,elemT>& output_cast =
+          dynamic_cast<DiscretisedDensityOnCartesianGrid<3,elemT> &>(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_image))
+    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);
+
+        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)
+            {
+              // The parabolic surrogate curvature of the QP is 1
+              elemT current = weights[dz][dy][dx] * input_image[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;
+      }
+    }
+  }
+}
+
+
 #  ifdef _MSC_VER
 // prevent warning message on reinstantiation, 
 // note that we get a linking error if we don't have the explicit instantiation below

src/swig/stir.i

@@ -1657,6 +1657,13 @@ stir::RegisteredParsingObject< stir::LogcoshPrior<elemT>,
         stir::PriorWithParabolicSurrogate<TargetT  > >;
 %template (LogcoshPrior3DFloat) stir::LogcoshPrior<elemT>;
 
+// Allows for access to parabolic surrogate type prior methods ( THIS IS A TEMPORARY FIX )
+%inline %{template <class T> stir::PriorWithParabolicSurrogate<T> * ToPriorWithParabolicSurrogate(stir::GeneralisedPrior<T> *b) {
+  return dynamic_cast<stir::PriorWithParabolicSurrogate<T>*>(b);
+}
+%}
+%template(ToPriorWithParabolicSurrogate) ToPriorWithParabolicSurrogate<TargetT >;
+
 %template (Reconstruction3DFloat) stir::Reconstruction<TargetT >;
 //%template () stir::Reconstruction<TargetT >;
 %template (IterativeReconstruction3DFloat) stir::IterativeReconstruction<TargetT >;