correct test on kappa characteristics in all priors

The Hessian() functions all incorrectly threw an error if kappa
was alright.

I've now moved this test into the check() functions, cleaning
code a little bit.

src/include/stir/recon_buildblock/LogcoshPrior.h

@@ -175,6 +175,9 @@ class LogcoshPrior : public RegisteredParsingObject<LogcoshPrior<elemT>,
   void initialise_keymap() override;
   bool post_processing() override;
 
+  //! Check that the prior is ready to be used
+  void check(DiscretisedDensity<3, elemT> const& current_image_estimate) const override;
+
 private:
   //! Spatially variant penalty penalty image ptr
   shared_ptr<DiscretisedDensity<3, elemT>> kappa_ptr;

src/recon_buildblock/LogcoshPrior.cxx

@@ -142,6 +142,21 @@ LogcoshPrior<elemT>::LogcoshPrior(const bool only_2D_v, float penalisation_facto
   this->scalar = scalar_v;
 }
 
+template <typename elemT>
+void
+LogcoshPrior<elemT>::check(DiscretisedDensity<3, elemT> const& current_image_estimate) const
+{
+  // Do base-class check
+  base_type::check(current_image_estimate);
+  if (!is_null_ptr(this->kappa_ptr))
+    {
+      std::string explanation;
+      if (!this->kappa_ptr->has_same_characteristics(current_image_estimate, explanation))
+        error(std::string(registered_name)
+              + ": kappa image does not have the same index range as the reconstructed image:" + explanation);
+    }
+}
+
 template <typename elemT>
 bool
 LogcoshPrior<elemT>::is_convex() const
@@ -251,9 +266,6 @@ LogcoshPrior<elemT>::compute_value(const DiscretisedDensity<3, elemT>& current_i
 
   const bool do_kappa = !is_null_ptr(kappa_ptr);
 
-  if (do_kappa && !kappa_ptr->has_same_characteristics(current_image_estimate))
-    error("LogcoshPrior: kappa image has not the same index range as the reconstructed image\n");
-
   double result = 0.;
   const int min_z = current_image_estimate.get_min_index();
   const int max_z = current_image_estimate.get_max_index();
@@ -323,8 +335,6 @@ LogcoshPrior<elemT>::compute_gradient(DiscretisedDensity<3, elemT>& prior_gradie
     }
 
   const bool do_kappa = !is_null_ptr(kappa_ptr);
-  if (do_kappa && !kappa_ptr->has_same_characteristics(current_image_estimate))
-    error("LogcoshPrior: kappa image has not the same index range as the reconstructed image\n");
 
   const int min_z = current_image_estimate.get_min_index();
   const int max_z = current_image_estimate.get_max_index();
@@ -409,9 +419,6 @@ LogcoshPrior<elemT>::compute_Hessian(DiscretisedDensity<3, elemT>& prior_Hessian
 
   const bool do_kappa = !is_null_ptr(kappa_ptr);
 
-  if (do_kappa && kappa_ptr->has_same_characteristics(current_image_estimate))
-    error("LogcoshPrior: 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];
@@ -479,9 +486,6 @@ LogcoshPrior<elemT>::parabolic_surrogate_curvature(DiscretisedDensity<3, elemT>&
 
   const bool do_kappa = !is_null_ptr(kappa_ptr);
 
-  if (do_kappa && !kappa_ptr->has_same_characteristics(current_image_estimate))
-    error("LogcoshPrior: kappa image has not the same index range as the reconstructed image\n");
-
   const int min_z = current_image_estimate.get_min_index();
   const int max_z = current_image_estimate.get_max_index();
   for (int z = min_z; z <= max_z; z++)
@@ -550,9 +554,6 @@ LogcoshPrior<elemT>::accumulate_Hessian_times_input(DiscretisedDensity<3, elemT>
 
   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++)

src/recon_buildblock/PLSPrior.cxx

@@ -121,6 +121,13 @@ PLSPrior<elemT>::check(DiscretisedDensity<3, elemT> const& current_image_estimat
   // Check anatomical and current image have same characteristics
   if (!this->anatomical_sptr->has_same_characteristics(current_image_estimate))
     error("The anatomical image must have the same charateristics as the PET image");
+  if (!is_null_ptr(this->kappa_ptr))
+    {
+      std::string explanation;
+      if (!this->kappa_ptr->has_same_characteristics(current_image_estimate, explanation))
+        error(std::string(registered_name)
+              + ": kappa image does not have the same index range as the reconstructed image:" + explanation);
+    }
 }
 
 template <typename elemT>
@@ -471,9 +478,6 @@ PLSPrior<elemT>::compute_value(const DiscretisedDensity<3, elemT>& current_image
 
   const bool do_kappa = !is_null_ptr(kappa_ptr);
 
-  if (do_kappa && !kappa_ptr->has_same_characteristics(current_image_estimate))
-    error("PLSPrior: kappa image has not the same index range as the reconstructed image\n");
-
   double result = 0.;
   const int min_z = current_image_estimate.get_min_index();
   const int max_z = current_image_estimate.get_max_index();
@@ -545,8 +549,6 @@ PLSPrior<elemT>::compute_gradient(DiscretisedDensity<3, elemT>& prior_gradient,
       *inner_product_sptr, *penalty_sptr, *pet_im_grad_z_sptr, *pet_im_grad_y_sptr, *pet_im_grad_x_sptr, current_image_estimate);
 
   const bool do_kappa = !is_null_ptr(kappa_ptr);
-  if (do_kappa && !kappa_ptr->has_same_characteristics(current_image_estimate))
-    error("PLSPrior: kappa image has not the same index range as the reconstructed image\n");
   shared_ptr<DiscretisedDensity<3, elemT>> gradient_sptr(this->anatomical_sptr->get_empty_copy());
 
   const int min_z = current_image_estimate.get_min_index();

src/recon_buildblock/QuadraticPrior.cxx

@@ -118,6 +118,13 @@ QuadraticPrior<elemT>::check(DiscretisedDensity<3, elemT> const& current_image_e
 {
   // Do base-class check
   base_type::check(current_image_estimate);
+  if (!is_null_ptr(this->kappa_ptr))
+    {
+      std::string explanation;
+      if (!this->kappa_ptr->has_same_characteristics(current_image_estimate, explanation))
+        error(std::string(registered_name)
+              + ": kappa image does not have the same index range as the reconstructed image:" + explanation);
+    }
 }
 
 template <typename elemT>
@@ -241,9 +248,6 @@ QuadraticPrior<elemT>::compute_value(const DiscretisedDensity<3, elemT>& current
 
   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");
-
   double result = 0.;
   const int min_z = current_image_estimate.get_min_index();
   const int max_z = current_image_estimate.get_max_index();
@@ -316,8 +320,6 @@ QuadraticPrior<elemT>::compute_gradient(DiscretisedDensity<3, elemT>& prior_grad
     }
 
   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 min_z = current_image_estimate.get_min_index();
   const int max_z = current_image_estimate.get_max_index();
@@ -407,9 +409,6 @@ QuadraticPrior<elemT>::compute_Hessian(DiscretisedDensity<3, elemT>& prior_Hessi
 
   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];
@@ -479,9 +478,6 @@ QuadraticPrior<elemT>::parabolic_surrogate_curvature(DiscretisedDensity<3, elemT
 
   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 min_z = current_image_estimate.get_min_index();
   const int max_z = current_image_estimate.get_max_index();
   for (int z = min_z; z <= max_z; z++)
@@ -559,9 +555,6 @@ QuadraticPrior<elemT>::add_multiplication_with_approximate_Hessian(DiscretisedDe
 
   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++)
@@ -629,9 +622,6 @@ QuadraticPrior<elemT>::accumulate_Hessian_times_input(DiscretisedDensity<3, elem
 
   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++)

src/recon_buildblock/RelativeDifferencePrior.cxx

@@ -122,6 +122,13 @@ RelativeDifferencePrior<elemT>::check(DiscretisedDensity<3, elemT> const& curren
 {
   // Do base-class check
   base_type::check(current_image_estimate);
+  if (!is_null_ptr(this->kappa_ptr))
+    {
+      std::string explanation;
+      if (!this->kappa_ptr->has_same_characteristics(current_image_estimate, explanation))
+        error(std::string(registered_name)
+              + ": kappa image does not have the same index range as the reconstructed image:" + explanation);
+    }
 }
 
 template <typename elemT>
@@ -312,9 +319,6 @@ RelativeDifferencePrior<elemT>::compute_value(const DiscretisedDensity<3, elemT>
 
   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");
-
   double result = 0.;
   const int min_z = current_image_estimate.get_min_index();
   const int max_z = current_image_estimate.get_max_index();
@@ -389,8 +393,6 @@ RelativeDifferencePrior<elemT>::compute_gradient(DiscretisedDensity<3, elemT>& p
     }
 
   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 min_z = current_image_estimate.get_min_index();
   const int max_z = current_image_estimate.get_max_index();
@@ -475,9 +477,6 @@ RelativeDifferencePrior<elemT>::compute_Hessian(DiscretisedDensity<3, elemT>& pr
 
   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];
@@ -554,9 +553,6 @@ RelativeDifferencePrior<elemT>::accumulate_Hessian_times_input(DiscretisedDensit
 
   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++)