diff --git a/src/include/stir/recon_buildblock/LogcoshPrior.h b/src/include/stir/recon_buildblock/LogcoshPrior.h index 0b2dfb45d..5b5715264 100644 --- a/src/include/stir/recon_buildblock/LogcoshPrior.h +++ b/src/include/stir/recon_buildblock/LogcoshPrior.h @@ -175,6 +175,9 @@ class LogcoshPrior : public RegisteredParsingObject, 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> kappa_ptr; diff --git a/src/recon_buildblock/LogcoshPrior.cxx b/src/recon_buildblock/LogcoshPrior.cxx index 0876407ea..68c3a602a 100644 --- a/src/recon_buildblock/LogcoshPrior.cxx +++ b/src/recon_buildblock/LogcoshPrior.cxx @@ -142,6 +142,21 @@ LogcoshPrior::LogcoshPrior(const bool only_2D_v, float penalisation_facto this->scalar = scalar_v; } +template +void +LogcoshPrior::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 bool LogcoshPrior::is_convex() const @@ -251,9 +266,6 @@ LogcoshPrior::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::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::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::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::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++) diff --git a/src/recon_buildblock/PLSPrior.cxx b/src/recon_buildblock/PLSPrior.cxx index 6ff213b1a..505a2c397 100644 --- a/src/recon_buildblock/PLSPrior.cxx +++ b/src/recon_buildblock/PLSPrior.cxx @@ -121,6 +121,13 @@ PLSPrior::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 @@ -471,9 +478,6 @@ PLSPrior::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::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> gradient_sptr(this->anatomical_sptr->get_empty_copy()); const int min_z = current_image_estimate.get_min_index(); diff --git a/src/recon_buildblock/QuadraticPrior.cxx b/src/recon_buildblock/QuadraticPrior.cxx index 2a070a294..69322c832 100644 --- a/src/recon_buildblock/QuadraticPrior.cxx +++ b/src/recon_buildblock/QuadraticPrior.cxx @@ -118,6 +118,13 @@ QuadraticPrior::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 @@ -241,9 +248,6 @@ QuadraticPrior::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::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::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::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::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::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++) diff --git a/src/recon_buildblock/RelativeDifferencePrior.cxx b/src/recon_buildblock/RelativeDifferencePrior.cxx index 3961aee7a..a25e44507 100644 --- a/src/recon_buildblock/RelativeDifferencePrior.cxx +++ b/src/recon_buildblock/RelativeDifferencePrior.cxx @@ -122,6 +122,13 @@ RelativeDifferencePrior::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 @@ -312,9 +319,6 @@ RelativeDifferencePrior::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::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::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::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++)