From 62a55213a1fba672152bc38c25f552bf28919496 Mon Sep 17 00:00:00 2001 From: robbietuk Date: Sun, 6 Jun 2021 17:30:01 -0700 Subject: [PATCH 01/27] RDP/accumulate_Hessian_times_input * RDP accumulate_Hessian_times_input * Added RDP Hessian documentation * Add div by 0 safety * Add epsilon to Hessian function --- .../RelativeDifferencePrior.h | 26 +++++++ .../RelativeDifferencePrior.cxx | 75 +++++++++++++++++++ 2 files changed, 101 insertions(+) diff --git a/src/include/stir/recon_buildblock/RelativeDifferencePrior.h b/src/include/stir/recon_buildblock/RelativeDifferencePrior.h index b4eaf4c74e..c8e4361ba0 100644 --- a/src/include/stir/recon_buildblock/RelativeDifferencePrior.h +++ b/src/include/stir/recon_buildblock/RelativeDifferencePrior.h @@ -130,6 +130,11 @@ class RelativeDifferencePrior: public 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 Succeeded 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 @@ -193,6 +198,27 @@ class RelativeDifferencePrior: public virtual bool post_processing(); private: shared_ptr > kappa_ptr; + + //! The Hessian of the Relative Difference Prior + /*! + This function returns the hessian (second derivative) of the RDP. + It is assumed that the diagonal elements of the Hessian are 0, or the weighing is, and thus only compute the partial + derivative w.r.t \c x and \c y. + * @param x is the target voxel. + * @param y is the neighbourhood voxel. + * @param gamma is the edge preservation value controlling the transition between the quadratic and linear behaviour + * @param eps is a small non-negative value included to prevent division by zero, see epsilon. + * @return the second derivative of the Relative Difference Prior + + */ + static inline float Hessian(const float x, const float y, const float gamma, const float eps) + { + if (x > 0.0 && y > 0.0 && eps > 0.0) + return 2 * (2 * x + eps)*(2 * y + eps) / + pow(x + y + gamma * abs(x - y) + eps, 3); + else + return 0.0; + } }; diff --git a/src/recon_buildblock/RelativeDifferencePrior.cxx b/src/recon_buildblock/RelativeDifferencePrior.cxx index e7fce65ea5..5eb64f81a1 100644 --- a/src/recon_buildblock/RelativeDifferencePrior.cxx +++ b/src/recon_buildblock/RelativeDifferencePrior.cxx @@ -430,6 +430,81 @@ add_multiplication_with_approximate_Hessian(DiscretisedDensity<3,elemT>& output, return Succeeded::no; } +template +Succeeded +RelativeDifferencePrior:: +accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, + const DiscretisedDensity<3,elemT>& current_estimate, + const DiscretisedDensity<3,elemT>& input) const +{ + // TODO this function overlaps enormously with parabolic_surrogate_curvature + // the only difference is that parabolic_surrogate_curvature uses input==1 + + assert( output.has_same_characteristics(input)); + if (this->penalisation_factor==0) + { + return Succeeded::yes; + } + + DiscretisedDensityOnCartesianGrid<3,elemT>& output_cast = + dynamic_cast &>(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)) + 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 current = weights[dz][dy][dx] * + Hessian(current_estimate[z][y][x], current_estimate[z+dz][y+dy][x+dx], + this->get_gamma(), this->get_epsilon()) * + input[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 Succeeded::yes; +} + # ifdef _MSC_VER // prevent warning message on reinstantiation, // note that we get a linking error if we don't have the explicit instantiation below From eb077909f43e2fb3134a7b1604d1fbf6ee7ca74f Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Sun, 6 Jun 2021 18:10:20 -0700 Subject: [PATCH 02/27] Fix issue with logic --- src/include/stir/recon_buildblock/RelativeDifferencePrior.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/include/stir/recon_buildblock/RelativeDifferencePrior.h b/src/include/stir/recon_buildblock/RelativeDifferencePrior.h index c8e4361ba0..b027be725b 100644 --- a/src/include/stir/recon_buildblock/RelativeDifferencePrior.h +++ b/src/include/stir/recon_buildblock/RelativeDifferencePrior.h @@ -213,7 +213,7 @@ class RelativeDifferencePrior: public */ static inline float Hessian(const float x, const float y, const float gamma, const float eps) { - if (x > 0.0 && y > 0.0 && eps > 0.0) + if (x > 0.0 || y > 0.0 || eps > 0.0) return 2 * (2 * x + eps)*(2 * y + eps) / pow(x + y + gamma * abs(x - y) + eps, 3); else From b7007b4be6c62d324db2ea29869650cd35c74729 Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Tue, 8 Jun 2021 14:26:07 -0700 Subject: [PATCH 03/27] Make prior test method `test_gradient` --- src/recon_test/test_priors.cxx | 48 ++++++++++++++++++++++------------ 1 file changed, 31 insertions(+), 17 deletions(-) diff --git a/src/recon_test/test_priors.cxx b/src/recon_test/test_priors.cxx index 961afdf0a8..5a9169e840 100644 --- a/src/recon_test/test_priors.cxx +++ b/src/recon_test/test_priors.cxx @@ -79,6 +79,10 @@ class GeneralisedPriorTests : public RunTests void run_tests_for_objective_function(const std::string& test_name, GeneralisedPrior& objective_function, shared_ptr target_sptr); + + void test_gradient(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr); }; GeneralisedPriorTests:: @@ -96,6 +100,16 @@ run_tests_for_objective_function(const std::string& test_name, if (!check(objective_function.set_up(target_sptr)==Succeeded::yes, "set-up of objective function")) return; + test_gradient(test_name, objective_function, target_sptr); +} + + +void +GeneralisedPriorTests:: +test_gradient(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr) +{ // setup images target_type& target(*target_sptr); shared_ptr gradient_sptr(target.get_empty_copy()); @@ -109,32 +123,32 @@ run_tests_for_objective_function(const std::string& test_name, const double value_at_target = objective_function.compute_value(target); target_type::full_iterator target_iter=target.begin_all(); target_type::full_iterator gradient_iter=gradient_sptr->begin_all(); - target_type::full_iterator gradient_2_iter=gradient_2_sptr->begin_all(); + target_type::full_iterator gradient_2_iter=gradient_2_sptr->begin_all(); // setup perturbation response const float eps = 1e-3F; bool testOK = true; info("Computing gradient of objective function by numerical differences (this will take a while)",3); while(target_iter!=target.end_all())// && testOK) + { + const float org_image_value = *target_iter; + *target_iter += eps; // perturb current voxel + const double value_at_inc = objective_function.compute_value(target); + *target_iter = org_image_value; // restore + const float ngradient_at_iter = static_cast((value_at_inc - value_at_target)/eps); + *gradient_2_iter = ngradient_at_iter; + testOK = testOK && this->check_if_equal(ngradient_at_iter, *gradient_iter, "gradient"); + //for (int i=0; i<5 && target_iter!=target.end_all(); ++i) { - const float org_image_value = *target_iter; - *target_iter += eps; // perturb current voxel - const double value_at_inc = objective_function.compute_value(target); - *target_iter = org_image_value; // restore - const float ngradient_at_iter = static_cast((value_at_inc - value_at_target)/eps); - *gradient_2_iter = ngradient_at_iter; - testOK = testOK && this->check_if_equal(ngradient_at_iter, *gradient_iter, "gradient"); - //for (int i=0; i<5 && target_iter!=target.end_all(); ++i) - { - ++gradient_2_iter; ++target_iter; ++ gradient_iter; - } + ++gradient_2_iter; ++target_iter; ++ gradient_iter; } + } if (!testOK) - { - info("Writing diagnostic files gradient" + test_name + ".hv, numerical_gradient" + test_name + ".hv"); - write_to_file("gradient" + test_name + ".hv", *gradient_sptr); - write_to_file("numerical_gradient" + test_name + ".hv", *gradient_2_sptr); - } + { + info("Writing diagnostic files gradient" + test_name + ".hv, numerical_gradient" + test_name + ".hv"); + write_to_file("gradient" + test_name + ".hv", *gradient_sptr); + write_to_file("numerical_gradient" + test_name + ".hv", *gradient_2_sptr); + } } From 96b6ec095e7660b71ac62f8489220dce76c3a1e6 Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Tue, 8 Jun 2021 15:16:23 -0700 Subject: [PATCH 04/27] Adds test_Hessian method --- src/recon_test/test_priors.cxx | 49 ++++++++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) diff --git a/src/recon_test/test_priors.cxx b/src/recon_test/test_priors.cxx index 5a9169e840..e1b852d3e8 100644 --- a/src/recon_test/test_priors.cxx +++ b/src/recon_test/test_priors.cxx @@ -83,6 +83,10 @@ class GeneralisedPriorTests : public RunTests void test_gradient(const std::string& test_name, GeneralisedPrior& objective_function, shared_ptr target_sptr); + + void test_Hessian(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr); }; GeneralisedPriorTests:: @@ -101,6 +105,7 @@ run_tests_for_objective_function(const std::string& test_name, return; test_gradient(test_name, objective_function, target_sptr); + test_Hessian(test_name, objective_function, target_sptr); } @@ -149,7 +154,51 @@ test_gradient(const std::string& test_name, write_to_file("gradient" + test_name + ".hv", *gradient_sptr); write_to_file("numerical_gradient" + test_name + ".hv", *gradient_2_sptr); } +} + +void +GeneralisedPriorTests:: +test_Hessian(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr) +{ + + // setup images + target_type& target(*target_sptr); + shared_ptr output(target.get_empty_copy()); + shared_ptr input(target.get_empty_copy()); + +// std::cout << ". max = " << input->find_max() << ". min = " << input->find_min() << "\n"; + + { + /// Construct an input with negative values by subtracting 0.5 from the input target + target_type::full_iterator input_iter=input->begin_all(); + target_type::full_iterator target_iter=target.begin_all(); + while(input_iter!=input->end_all())// && testOK) + { + *input_iter = *target_iter - 0.5; + ++input_iter; ++target_iter; + } + } + objective_function.accumulate_Hessian_times_input(*output, target, *input); +// std::cout << "accumulate_Hessian_times_input output -> max = " << output->find_max() << ". min = " << output->find_min() << "\n"; + + { + /// Construct an input with negative values by subtracting 0.5 from the input target + target_type::full_iterator input_iter=input->begin_all(); + target_type::full_iterator output_iter=output->begin_all(); + float my_sum = 0.0; + while(input_iter!=input->end_all())// && testOK) + { + my_sum += *output_iter * *input_iter; + ++input_iter; ++output_iter; + } +// std::cout << "my_sum = " << my_sum << "\n"; + if (this->check_if_less(0, my_sum)) + info("Computation of (x^T H x < 0) " + test_name); + } +// std::cout << ". max = " << input->find_max() << ". min = " << input->find_min() << "\n"; } void From 2f418bcf02fb46c5ed6f2fa410ba3b476b06b78e Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Wed, 9 Jun 2021 14:33:23 -0700 Subject: [PATCH 05/27] Test an array of different cases for the Hessian condition --- src/recon_test/test_priors.cxx | 122 ++++++++++++++++++++++++++------- 1 file changed, 97 insertions(+), 25 deletions(-) diff --git a/src/recon_test/test_priors.cxx b/src/recon_test/test_priors.cxx index e1b852d3e8..d567e30d4d 100644 --- a/src/recon_test/test_priors.cxx +++ b/src/recon_test/test_priors.cxx @@ -53,6 +53,8 @@ START_NAMESPACE_STIR This test compares the result of GeneralisedPrior::compute_gradient() with a numerical gradient computed by using the GeneralisedPrior::compute_value() function. + Additionally, the Hessian is tested, via GeneralisedPrior::accumulate_Hessian_times_input(), + by evaluating the x^T Hx > 0 constraint. */ class GeneralisedPriorTests : public RunTests @@ -84,9 +86,26 @@ class GeneralisedPriorTests : public RunTests GeneralisedPrior& objective_function, shared_ptr target_sptr); - void test_Hessian(const std::string& test_name, - GeneralisedPrior& objective_function, - shared_ptr target_sptr); + //! Test various configurations of the prior's Hessian via accumulate_Hessian_times_input() + /*! + Tests the concave function condition + \f[ x^T \cdot H_{\lambda}x >= 0 \f] + for all non-negative \c x and non-zero \c \lambda (Relative Difference Penalty conditions). + This function constructs an array of configurations to test this condition and calls \c test_Hessian_configuration(). + */ + void test_Hessian(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr); + +private: + //! Hessian test for a particular configuration of the Hessian concave condition + void + test_Hessian_configuration(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr, + float beta, + float input_multiplication, float input_addition, + float current_image_multiplication, float current_image_addition); }; GeneralisedPriorTests:: @@ -104,7 +123,9 @@ run_tests_for_objective_function(const std::string& test_name, if (!check(objective_function.set_up(target_sptr)==Succeeded::yes, "set-up of objective function")) return; + std::cerr << "----- test " << test_name << " --> Gradient\n"; test_gradient(test_name, objective_function, target_sptr); + std::cerr << "----- test " << test_name << " --> Hessian-vector product (accumulate_Hessian_times_input)\n"; test_Hessian(test_name, objective_function, target_sptr); } @@ -162,43 +183,94 @@ test_Hessian(const std::string& test_name, GeneralisedPrior& objective_function, shared_ptr target_sptr) { + /// Construct configurations + float beta_array[] = {0.01, 1, 100}; // Penalty strength should only affect scale + // Modifications to the input image + float input_multiplication_array[] = {-100, -1, 0.01, 1, 100}; // Test negative, small and large values + float input_addition_array[] = {-10, -1, -0.5, 0.0, 1, 10}; + // Modifications to the current image (Hessian computation) + float current_image_multiplication_array[] = {0.01, 1, 100}; + float current_image_addition_array[] = {0.0, 0.5, 1, 10}; // RDP has constraint that current_image is non-negative + + float initial_beta = objective_function.get_penalisation_factor(); + for (float beta : beta_array) { + for (float input_multiplication : input_multiplication_array) { + for (float input_addition : input_addition_array) { + for (float current_image_multiplication : current_image_multiplication_array) { + for (float current_image_addition : current_image_addition_array) { + test_Hessian_configuration(test_name, objective_function, target_sptr, + beta, + input_multiplication, input_addition, + current_image_multiplication, current_image_addition); + } + } + } + } + } + /// Reset beta to original value + objective_function.set_penalisation_factor(initial_beta); +} - // setup images +void +GeneralisedPriorTests:: +test_Hessian_configuration(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr, + const float beta, + const float input_multiplication, const float input_addition, + const float current_image_multiplication, const float current_image_addition) +{ + /// setup images target_type& target(*target_sptr); shared_ptr output(target.get_empty_copy()); + shared_ptr current_image(target.get_empty_copy()); shared_ptr input(target.get_empty_copy()); - -// std::cout << ". max = " << input->find_max() << ". min = " << input->find_min() << "\n"; - + objective_function.set_penalisation_factor(beta); { - /// Construct an input with negative values by subtracting 0.5 from the input target - target_type::full_iterator input_iter=input->begin_all(); - target_type::full_iterator target_iter=target.begin_all(); - while(input_iter!=input->end_all())// && testOK) + /// Construct an current_image & input with various values that are variations of the target + target_type::full_iterator input_iter = input->begin_all(); + target_type::full_iterator current_image_iter = current_image->begin_all(); + target_type::full_iterator target_iter = target.begin_all(); + while (input_iter != input->end_all())// && testOK) { - *input_iter = *target_iter - 0.5; - ++input_iter; ++target_iter; + *input_iter = input_multiplication * *target_iter + input_addition; + *current_image_iter = current_image_multiplication * *target_iter + current_image_addition; + ++input_iter; ++target_iter; ++current_image_iter; } } - objective_function.accumulate_Hessian_times_input(*output, target, *input); -// std::cout << "accumulate_Hessian_times_input output -> max = " << output->find_max() << ". min = " << output->find_min() << "\n"; + /// Compute H x + objective_function.accumulate_Hessian_times_input(*output, *current_image, *input); + /// Compute x \cdot (H x) + float my_sum = 0.0; { - /// Construct an input with negative values by subtracting 0.5 from the input target - target_type::full_iterator input_iter=input->begin_all(); - target_type::full_iterator output_iter=output->begin_all(); - float my_sum = 0.0; - while(input_iter!=input->end_all())// && testOK) + target_type::full_iterator input_iter = input->begin_all(); + target_type::full_iterator output_iter = output->begin_all(); + while (input_iter != input->end_all())// && testOK) { my_sum += *output_iter * *input_iter; - ++input_iter; ++output_iter; + ++input_iter; + ++output_iter; } -// std::cout << "my_sum = " << my_sum << "\n"; - if (this->check_if_less(0, my_sum)) - info("Computation of (x^T H x < 0) " + test_name); } -// std::cout << ". max = " << input->find_max() << ". min = " << input->find_min() << "\n"; + + if (this->check_if_less(0, my_sum)) { +// info("PASS: Computation of x^T H x = " + std::to_string(my_sum) + " > 0"); + } else { + // print to console the FAILED configuration + info("FAIL: Computation of x^T H x = " + std::to_string(my_sum) + " < 0" + + "\ntest_name=" + test_name + + "\nbeta=" + std::to_string(beta) + + "\ninput_multiplication=" + std::to_string(input_multiplication) + + "\ninput_addition=" + std::to_string(input_addition) + + "\ncurrent_image_multiplication=" + std::to_string(current_image_multiplication) + + "\ncurrent_image_addition=" + std::to_string(current_image_addition) + + "\n >input image max=" + std::to_string(input->find_max()) + + "\n >input image min=" + std::to_string(input->find_min()) + + "\n >target image max=" + std::to_string(target.find_max()) + + "\n >target image min=" + std::to_string(target.find_min())); + } } void From fa65422ff838bb8b8ea3d18a9f3d367d8e824bb9 Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Wed, 9 Jun 2021 16:00:43 -0700 Subject: [PATCH 06/27] improve Hessian documentation and info --- src/recon_test/test_priors.cxx | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/recon_test/test_priors.cxx b/src/recon_test/test_priors.cxx index d567e30d4d..7d231ace1d 100644 --- a/src/recon_test/test_priors.cxx +++ b/src/recon_test/test_priors.cxx @@ -53,7 +53,7 @@ START_NAMESPACE_STIR This test compares the result of GeneralisedPrior::compute_gradient() with a numerical gradient computed by using the GeneralisedPrior::compute_value() function. - Additionally, the Hessian is tested, via GeneralisedPrior::accumulate_Hessian_times_input(), + Additionally, the Hessian's convexity is tested, via GeneralisedPrior::accumulate_Hessian_times_input(), by evaluating the x^T Hx > 0 constraint. */ @@ -99,13 +99,12 @@ class GeneralisedPriorTests : public RunTests private: //! Hessian test for a particular configuration of the Hessian concave condition - void - test_Hessian_configuration(const std::string& test_name, - GeneralisedPrior& objective_function, - shared_ptr target_sptr, - float beta, - float input_multiplication, float input_addition, - float current_image_multiplication, float current_image_addition); + void test_Hessian_configuration(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr, + float beta, + float input_multiplication, float input_addition, + float current_image_multiplication, float current_image_addition); }; GeneralisedPriorTests:: @@ -255,11 +254,12 @@ test_Hessian_configuration(const std::string& test_name, } } + // test for a CONVEX function if (this->check_if_less(0, my_sum)) { -// info("PASS: Computation of x^T H x = " + std::to_string(my_sum) + " > 0"); +// info("PASS: Computation of x^T H x = " + std::to_string(my_sum) + " > 0 and is therefore convex"); } else { // print to console the FAILED configuration - info("FAIL: Computation of x^T H x = " + std::to_string(my_sum) + " < 0" + + info("FAIL: Computation of x^T H x = " + std::to_string(my_sum) + " < 0 and is therefore NOT convex" + "\ntest_name=" + test_name + "\nbeta=" + std::to_string(beta) + "\ninput_multiplication=" + std::to_string(input_multiplication) + From 3e4e5f44ffe917c64bb06f1563f30f894cd78bdb Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Fri, 11 Jun 2021 13:43:51 -0700 Subject: [PATCH 07/27] Improve dot product --- src/recon_test/test_priors.cxx | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/src/recon_test/test_priors.cxx b/src/recon_test/test_priors.cxx index 7d231ace1d..3b8bebbf90 100644 --- a/src/recon_test/test_priors.cxx +++ b/src/recon_test/test_priors.cxx @@ -243,16 +243,7 @@ test_Hessian_configuration(const std::string& test_name, /// Compute x \cdot (H x) float my_sum = 0.0; - { - target_type::full_iterator input_iter = input->begin_all(); - target_type::full_iterator output_iter = output->begin_all(); - while (input_iter != input->end_all())// && testOK) - { - my_sum += *output_iter * *input_iter; - ++input_iter; - ++output_iter; - } - } + my_sum = std::inner_product(input->begin_all(), input->end_all(), output->begin_all(), my_sum); // test for a CONVEX function if (this->check_if_less(0, my_sum)) { From 101c4fc10f4d4888119c38bcec3de237ec5365f0 Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Fri, 11 Jun 2021 13:51:33 -0700 Subject: [PATCH 08/27] Adds _is_convex variable to priors. default is false --- src/include/stir/recon_buildblock/GeneralisedPrior.h | 3 +++ src/recon_buildblock/GeneralisedPrior.cxx | 1 + src/recon_buildblock/LogcoshPrior.cxx | 1 + src/recon_buildblock/PLSPrior.cxx | 1 + src/recon_buildblock/QuadraticPrior.cxx | 1 + src/recon_buildblock/RelativeDifferencePrior.cxx | 1 + 6 files changed, 8 insertions(+) diff --git a/src/include/stir/recon_buildblock/GeneralisedPrior.h b/src/include/stir/recon_buildblock/GeneralisedPrior.h index caa4bf4dab..86135e5b8d 100644 --- a/src/include/stir/recon_buildblock/GeneralisedPrior.h +++ b/src/include/stir/recon_buildblock/GeneralisedPrior.h @@ -102,6 +102,9 @@ class GeneralisedPrior: virtual void check(DataT const& current_estimate) const; bool _already_set_up; + + //! Variable to indicate that the prior is a convex function + bool _is_convex; }; END_NAMESPACE_STIR diff --git a/src/recon_buildblock/GeneralisedPrior.cxx b/src/recon_buildblock/GeneralisedPrior.cxx index 5751c0a6b5..d62c520824 100644 --- a/src/recon_buildblock/GeneralisedPrior.cxx +++ b/src/recon_buildblock/GeneralisedPrior.cxx @@ -40,6 +40,7 @@ GeneralisedPrior::set_defaults() { _already_set_up = false; this->penalisation_factor = 0; + this->_is_convex = false; } template diff --git a/src/recon_buildblock/LogcoshPrior.cxx b/src/recon_buildblock/LogcoshPrior.cxx index f14742d7b1..d03498f6e9 100644 --- a/src/recon_buildblock/LogcoshPrior.cxx +++ b/src/recon_buildblock/LogcoshPrior.cxx @@ -108,6 +108,7 @@ void LogcoshPrior::set_defaults() { base_type::set_defaults(); + this->_is_convex = true; this->only_2D = false; this->scalar = 1.0; this->kappa_ptr.reset(); diff --git a/src/recon_buildblock/PLSPrior.cxx b/src/recon_buildblock/PLSPrior.cxx index 0b782546fa..daa09ca30c 100644 --- a/src/recon_buildblock/PLSPrior.cxx +++ b/src/recon_buildblock/PLSPrior.cxx @@ -127,6 +127,7 @@ void PLSPrior::set_defaults() { base_type::set_defaults(); + this->_is_convex = true; this->only_2D = false; this->alpha=1; this->eta=1; diff --git a/src/recon_buildblock/QuadraticPrior.cxx b/src/recon_buildblock/QuadraticPrior.cxx index 904b4b0d1b..3e247ce89f 100644 --- a/src/recon_buildblock/QuadraticPrior.cxx +++ b/src/recon_buildblock/QuadraticPrior.cxx @@ -123,6 +123,7 @@ void QuadraticPrior::set_defaults() { base_type::set_defaults(); + this->_is_convex = true; this->only_2D = false; this->kappa_ptr.reset(); this->weights.recycle(); diff --git a/src/recon_buildblock/RelativeDifferencePrior.cxx b/src/recon_buildblock/RelativeDifferencePrior.cxx index 5eb64f81a1..c3aa15ce41 100644 --- a/src/recon_buildblock/RelativeDifferencePrior.cxx +++ b/src/recon_buildblock/RelativeDifferencePrior.cxx @@ -127,6 +127,7 @@ void RelativeDifferencePrior::set_defaults() { base_type::set_defaults(); + this->_is_convex = true; this->only_2D = false; this->kappa_ptr.reset(); this->weights.recycle(); From 8187715bdea7498b09829e587fa81787525dacd4 Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Fri, 11 Jun 2021 14:21:56 -0700 Subject: [PATCH 09/27] Reimplementation add_multiplication_with_approximate_Hessian --- src/recon_buildblock/QuadraticPrior.cxx | 67 ++++++++++++++++++++++++- 1 file changed, 66 insertions(+), 1 deletion(-) diff --git a/src/recon_buildblock/QuadraticPrior.cxx b/src/recon_buildblock/QuadraticPrior.cxx index 3e247ce89f..a0899491d3 100644 --- a/src/recon_buildblock/QuadraticPrior.cxx +++ b/src/recon_buildblock/QuadraticPrior.cxx @@ -582,7 +582,72 @@ QuadraticPrior:: add_multiplication_with_approximate_Hessian(DiscretisedDensity<3,elemT>& output, const DiscretisedDensity<3,elemT>& input) const { - return accumulate_Hessian_times_input(output, input, input); + // TODO this function overlaps enormously with parabolic_surrogate_curvature + // the only difference is that parabolic_surrogate_curvature uses input==1 + + assert( output.has_same_characteristics(input)); + if (this->penalisation_factor==0) + { + return Succeeded::yes; + } + + this->check(input); + + DiscretisedDensityOnCartesianGrid<3,elemT>& output_cast = + dynamic_cast &>(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)) + 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++) + { + 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 current = + weights[dz][dy][dx] * input[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 Succeeded::yes; } template From 410943e34c333480dc988bb1d7ab3eb70147c8a7 Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Fri, 11 Jun 2021 14:22:21 -0700 Subject: [PATCH 10/27] Correct math in QP accumulate_Hessian_times_input --- src/recon_buildblock/QuadraticPrior.cxx | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/src/recon_buildblock/QuadraticPrior.cxx b/src/recon_buildblock/QuadraticPrior.cxx index a0899491d3..d2309d26a7 100644 --- a/src/recon_buildblock/QuadraticPrior.cxx +++ b/src/recon_buildblock/QuadraticPrior.cxx @@ -704,13 +704,27 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, 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); + /// At this point, we have j = [z][y][x] + // The next for loops will have k = [z+dz][y+dy][x+dx] + // The following computes + //(H_{wf} y)_j = + // \sum_{k\in N_j} w_{(j,k)} f''_{d}(x_j,x_k) y_j + + // \sum_{(i \in N_j) \ne j} w_{(j,i)} f''_{od}(x_j, x_i) y_i + 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 current = - weights[dz][dy][dx] * input[z+dz][y+dy][x+dx]; + elemT current = 0.0; + if (dz != 0 || dy != 0 || dz != 0) + { + current = weights[dz][dy][dx] * + ( (1.0) * input[z][y][x] + (-1.0) * input[z+dz][y+dy][x+dx]) ; + } else { + // The j == k case + current = weights[dz][dy][dx] * (1.0) * input[z][y][x]; + } if (do_kappa) current *= From 948397f6aa8e5d3d5ec09841b401b79ef719a510 Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Fri, 11 Jun 2021 14:24:29 -0700 Subject: [PATCH 11/27] Improve QP comments --- src/recon_buildblock/QuadraticPrior.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/recon_buildblock/QuadraticPrior.cxx b/src/recon_buildblock/QuadraticPrior.cxx index d2309d26a7..6c6d27b6c4 100644 --- a/src/recon_buildblock/QuadraticPrior.cxx +++ b/src/recon_buildblock/QuadraticPrior.cxx @@ -720,7 +720,8 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, if (dz != 0 || dy != 0 || dz != 0) { current = weights[dz][dy][dx] * - ( (1.0) * input[z][y][x] + (-1.0) * input[z+dz][y+dy][x+dx]) ; + ( (1.0) * input[z][y][x] + // diagonal Hessian terms + (-1.0) * input[z+dz][y+dy][x+dx] ); // off-diagonal Hessian terms } else { // The j == k case current = weights[dz][dy][dx] * (1.0) * input[z][y][x]; From a339ac05d1d2a63ab1d15a1eeea03c4c5c99c8d8 Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Tue, 15 Jun 2021 09:42:24 -0700 Subject: [PATCH 12/27] [ci skip] Add two methods to RDP class :(off_)diagonal_second_derivative Updates the math for the RDP --- .../RelativeDifferencePrior.h | 10 +--- .../RelativeDifferencePrior.cxx | 47 +++++++++++++++++-- 2 files changed, 45 insertions(+), 12 deletions(-) diff --git a/src/include/stir/recon_buildblock/RelativeDifferencePrior.h b/src/include/stir/recon_buildblock/RelativeDifferencePrior.h index b027be725b..adff6673ff 100644 --- a/src/include/stir/recon_buildblock/RelativeDifferencePrior.h +++ b/src/include/stir/recon_buildblock/RelativeDifferencePrior.h @@ -211,14 +211,8 @@ class RelativeDifferencePrior: public * @return the second derivative of the Relative Difference Prior */ - static inline float Hessian(const float x, const float y, const float gamma, const float eps) - { - if (x > 0.0 || y > 0.0 || eps > 0.0) - return 2 * (2 * x + eps)*(2 * y + eps) / - pow(x + y + gamma * abs(x - y) + eps, 3); - else - return 0.0; - } + float diagonal_second_derivative(const float x_j, const float x_k) const; + float off_diagonal_second_derivative(const float x_j, const float x_k) const; }; diff --git a/src/recon_buildblock/RelativeDifferencePrior.cxx b/src/recon_buildblock/RelativeDifferencePrior.cxx index c3aa15ce41..57c742ad39 100644 --- a/src/recon_buildblock/RelativeDifferencePrior.cxx +++ b/src/recon_buildblock/RelativeDifferencePrior.cxx @@ -483,15 +483,28 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, 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); + /// At this point, we have j = [z][y][x] + // The next for loops will have k = [z+dz][y+dy][x+dx] + // The following computes + //(H_{wf} y)_j = + // \sum_{k\in N_j} w_{(j,k)} f''_{d}(x_j,x_k) y_j + + // \sum_{(i \in N_j) \ne j} w_{(j,i)} f''_{od}(x_j, x_i) y_i + 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 current = weights[dz][dy][dx] * - Hessian(current_estimate[z][y][x], current_estimate[z+dz][y+dy][x+dx], - this->get_gamma(), this->get_epsilon()) * - input[z+dz][y+dy][x+dx]; + elemT current = 0.0; + if (dz != 0 || dy != 0 || dz != 0) + { + current = weights[dz][dy][dx] * + (diagonal_second_derivative(current_estimate[z][y][x], current_estimate[z+dz][y+dy][x+dx]) * input[z][y][x] + + off_diagonal_second_derivative(current_estimate[z][y][x], current_estimate[z+dz][y+dy][x+dx])* input[z+dz][y+dy][x+dx] ); + } else { + // The j == k case + current = weights[dz][dy][dx] * (1.0) * input[z][y][x]; + } if (do_kappa) current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; @@ -506,6 +519,32 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, return Succeeded::yes; } +template +float +RelativeDifferencePrior:: +diagonal_second_derivative(const float x_j, const float x_k) const +{ + // THIS MATH IS WRONG!!! + if (x_j > 0.0 || x_k > 0.0 || this->epsilon > 0.0) + return 2 * (2 * x_j + this->epsilon)*(2 * x_k + this->epsilon) / + pow(x_j + x_k + this->gamma * abs(x_j - x_k) + this->epsilon, 3); + else + return 0.0; +} + +template +float +RelativeDifferencePrior:: +off_diagonal_second_derivative(const float x_j, const float x_k) const +{ + // Not sure if math is correct + if (x_j > 0.0 || x_k > 0.0 || this->epsilon > 0.0) + return 2 * (2 * x_j + this->epsilon)*(2 * x_k + this->epsilon) / + pow(x_j + x_k + this->gamma * abs(x_j - x_k) + this->epsilon, 3); + else + return 0.0; +} + # ifdef _MSC_VER // prevent warning message on reinstantiation, // note that we get a linking error if we don't have the explicit instantiation below From b6e6b54d54ab154e61856381cbde843f672c2afc Mon Sep 17 00:00:00 2001 From: robbietuk Date: Tue, 15 Jun 2021 14:42:29 -0700 Subject: [PATCH 13/27] Prior Hessian methods use second partial derivative functions (#14) * RDP Documententation * QP use (off_)diagonal_second_derivative methods and documentation * Restructure RDP and QP Hessian logic * Improve RDP and QP documentation * Implement log cosh accumulate_Hessian_times_input in terms of derivatives * Documentation --- .../stir/recon_buildblock/LogcoshPrior.h | 14 ++++++++ .../stir/recon_buildblock/QuadraticPrior.h | 14 ++++++++ .../RelativeDifferencePrior.h | 18 +++++----- src/recon_buildblock/LogcoshPrior.cxx | 30 +++++++++++++++-- src/recon_buildblock/QuadraticPrior.cxx | 33 ++++++++++++++----- .../RelativeDifferencePrior.cxx | 19 ++++++----- 6 files changed, 99 insertions(+), 29 deletions(-) diff --git a/src/include/stir/recon_buildblock/LogcoshPrior.h b/src/include/stir/recon_buildblock/LogcoshPrior.h index a3bed89478..574436300f 100644 --- a/src/include/stir/recon_buildblock/LogcoshPrior.h +++ b/src/include/stir/recon_buildblock/LogcoshPrior.h @@ -236,6 +236,20 @@ class LogcoshPrior: public const float x = d * scalar; return square((1/ cosh(x))); } + + //! The second partial derivatives of the LogCosh Prior + /*! + Diagonal refers to the second derivative w.r.t. x_j only (i.e. diagonal of the Hessian) + Off-diagonal refers to the second derivative w.r.t. x_j and x_k (i.e. off-diagonal of the Hessian) + For LogCosh, the off diagonal is the negative of the diagonal. + * @param x_j is the target voxel. + * @param x_k is the voxel in the neighbourhood. + * @return the second partial derivative of the LogCosh Prior + */ + //@{ + float diagonal_second_derivative(const float x_j, const float x_k) const; + float off_diagonal_second_derivative(const float x_j, const float x_k) const; + //@} }; diff --git a/src/include/stir/recon_buildblock/QuadraticPrior.h b/src/include/stir/recon_buildblock/QuadraticPrior.h index 4028169276..ca09af6f22 100644 --- a/src/include/stir/recon_buildblock/QuadraticPrior.h +++ b/src/include/stir/recon_buildblock/QuadraticPrior.h @@ -179,6 +179,20 @@ class QuadraticPrior: public virtual bool post_processing(); private: shared_ptr > kappa_ptr; + + //! The second partial derivatives of the Quadratic Prior + /*! + Diagonal refers to the second derivative w.r.t. x_j only (i.e. diagonal of the Hessian) + Off-diagonal refers to the second derivative w.r.t. x_j and x_k (i.e. off-diagonal of the Hessian) + For the Quadratic Prior, the off diagonal is the negative of the diagonal. + * @param x_j is the target voxel. + * @param x_k is the voxel in the neighbourhood. + * @return the second partial derivative of the Quadratic Prior + */ + //@{ + float diagonal_second_derivative(const float x_j, const float x_k) const; + float off_diagonal_second_derivative(const float x_j, const float x_k) const; + //@} }; diff --git a/src/include/stir/recon_buildblock/RelativeDifferencePrior.h b/src/include/stir/recon_buildblock/RelativeDifferencePrior.h index adff6673ff..2e93b57784 100644 --- a/src/include/stir/recon_buildblock/RelativeDifferencePrior.h +++ b/src/include/stir/recon_buildblock/RelativeDifferencePrior.h @@ -199,20 +199,18 @@ class RelativeDifferencePrior: public private: shared_ptr > kappa_ptr; - //! The Hessian of the Relative Difference Prior + //! The second partial derivatives of the Relative Difference Prior /*! - This function returns the hessian (second derivative) of the RDP. - It is assumed that the diagonal elements of the Hessian are 0, or the weighing is, and thus only compute the partial - derivative w.r.t \c x and \c y. - * @param x is the target voxel. - * @param y is the neighbourhood voxel. - * @param gamma is the edge preservation value controlling the transition between the quadratic and linear behaviour - * @param eps is a small non-negative value included to prevent division by zero, see epsilon. - * @return the second derivative of the Relative Difference Prior - + Diagonal refers to the second derivative w.r.t. x_j only (i.e. diagonal of the Hessian) + Off-diagonal refers to the second derivative w.r.t. x_j and x_k (i.e. off-diagonal of the Hessian) + * @param x_j is the target voxel. + * @param x_k is the voxel in the neighbourhood. + * @return the second partial derivative of the Relative Difference Prior */ + //@{ float diagonal_second_derivative(const float x_j, const float x_k) const; float off_diagonal_second_derivative(const float x_j, const float x_k) const; + //@} }; diff --git a/src/recon_buildblock/LogcoshPrior.cxx b/src/recon_buildblock/LogcoshPrior.cxx index d03498f6e9..43134b84e4 100644 --- a/src/recon_buildblock/LogcoshPrior.cxx +++ b/src/recon_buildblock/LogcoshPrior.cxx @@ -562,8 +562,18 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, for (int dy=min_dy;dy<=max_dy;++dy) for (int dx=min_dx;dx<=max_dx;++dx) { - elemT voxel_diff= current_estimate[z][y][x] - current_estimate[z+dz][y+dy][x+dx]; - elemT current = weights[dz][dy][dx] * Hessian( voxel_diff, this->scalar) * input[z+dz][y+dy][x+dx]; + elemT current = weights[dz][dy][dx]; + if (dz == dy == dz == 0) { + // The j == k case + current *= diagonal_second_derivative(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x]; + } else { + current *= (diagonal_second_derivative(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x] + + off_diagonal_second_derivative(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * + input[z + dz][y + dy][x + dx]); + } if (do_kappa) current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; @@ -578,6 +588,22 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, return Succeeded::yes; } +template +float +LogcoshPrior:: +diagonal_second_derivative(const float x_j, const float x_k) const +{ + return square((1/ cosh((x_j - x_k) * scalar))); +} + +template +float +LogcoshPrior:: +off_diagonal_second_derivative(const float x_j, const float x_k) const +{ + return - diagonal_second_derivative(x_j, x_k); +} + # ifdef _MSC_VER // prevent warning message on reinstantiation, // note that we get a linking error if we don't have the explicit instantiation below diff --git a/src/recon_buildblock/QuadraticPrior.cxx b/src/recon_buildblock/QuadraticPrior.cxx index 6c6d27b6c4..657b8a7be6 100644 --- a/src/recon_buildblock/QuadraticPrior.cxx +++ b/src/recon_buildblock/QuadraticPrior.cxx @@ -654,7 +654,7 @@ template Succeeded QuadraticPrior:: accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, - const DiscretisedDensity<3,elemT>& /*current_estimate*/, + const DiscretisedDensity<3,elemT>& current_estimate, const DiscretisedDensity<3,elemT>& input) const { // TODO this function overlaps enormously with parabolic_surrogate_curvature @@ -710,21 +710,24 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, //(H_{wf} y)_j = // \sum_{k\in N_j} w_{(j,k)} f''_{d}(x_j,x_k) y_j + // \sum_{(i \in N_j) \ne j} w_{(j,i)} f''_{od}(x_j, x_i) y_i + // Note the condition in the second sum that i is not equal to j 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 current = 0.0; - if (dz != 0 || dy != 0 || dz != 0) - { - current = weights[dz][dy][dx] * - ( (1.0) * input[z][y][x] + // diagonal Hessian terms - (-1.0) * input[z+dz][y+dy][x+dx] ); // off-diagonal Hessian terms - } else { + elemT current = weights[dz][dy][dx]; + if (dz == dy == dz == 0) { // The j == k case - current = weights[dz][dy][dx] * (1.0) * input[z][y][x]; + current *= diagonal_second_derivative(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x]; + } else { + current *= (diagonal_second_derivative(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x] + + off_diagonal_second_derivative(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * + input[z + dz][y + dy][x + dx]); } if (do_kappa) @@ -740,6 +743,18 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, return Succeeded::yes; } +template +float +QuadraticPrior:: +diagonal_second_derivative(const float x_j, const float x_k) const +{ return 1.0;} + +template +float +QuadraticPrior:: +off_diagonal_second_derivative(const float x_j, const float x_k) const +{ return -1.0;} + # ifdef _MSC_VER // prevent warning message on reinstantiation, // note that we get a linking error if we don't have the explicit instantiation below diff --git a/src/recon_buildblock/RelativeDifferencePrior.cxx b/src/recon_buildblock/RelativeDifferencePrior.cxx index 57c742ad39..beb9557e7c 100644 --- a/src/recon_buildblock/RelativeDifferencePrior.cxx +++ b/src/recon_buildblock/RelativeDifferencePrior.cxx @@ -489,21 +489,24 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, //(H_{wf} y)_j = // \sum_{k\in N_j} w_{(j,k)} f''_{d}(x_j,x_k) y_j + // \sum_{(i \in N_j) \ne j} w_{(j,i)} f''_{od}(x_j, x_i) y_i + // Note the condition in the second sum that i is not equal to j 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 current = 0.0; - if (dz != 0 || dy != 0 || dz != 0) - { - current = weights[dz][dy][dx] * - (diagonal_second_derivative(current_estimate[z][y][x], current_estimate[z+dz][y+dy][x+dx]) * input[z][y][x] + - off_diagonal_second_derivative(current_estimate[z][y][x], current_estimate[z+dz][y+dy][x+dx])* input[z+dz][y+dy][x+dx] ); - } else { + elemT current = weights[dz][dy][dx]; + if (dz == dy == dz == 0) { // The j == k case - current = weights[dz][dy][dx] * (1.0) * input[z][y][x]; + current *= diagonal_second_derivative(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x]; + } else { + current *= (diagonal_second_derivative(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x] + + off_diagonal_second_derivative(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * + input[z + dz][y + dy][x + dx]); } if (do_kappa) From b3c8d6f425c051d9c15964ad2e3b2fd6d98af04d Mon Sep 17 00:00:00 2001 From: robbietuk Date: Sat, 19 Jun 2021 12:41:56 -0700 Subject: [PATCH 14/27] compute_Hessian method to all priors and test improvements in general * Add method get_is_convex(), which accesses prior is_convex parameter and only test convexity of priors if convex function * compute_Hessian method added to generalised prior and errors by default but with different messages depending on is_convex * Test rename and create test_Hessian_convexity and test_Hessian_methods * Add test for compute_Hessian * RelativeDifferencePrior initialisation call set_defaults * Modernise `compute_Hessian` for QP and LogCosh. Add `compute_Hessian` for RDP * Correct RDP second derivative functions * Major changes to test_Hessian_against_numerical, which loops over each voxel for perturbation response * Add verbosity suppression to suppress gradient info calls --- .../stir/recon_buildblock/GeneralisedPrior.h | 16 +- .../stir/recon_buildblock/LogcoshPrior.h | 22 +- .../stir/recon_buildblock/QuadraticPrior.h | 9 +- .../RelativeDifferencePrior.h | 7 +- src/recon_buildblock/GeneralisedPrior.cxx | 19 ++ src/recon_buildblock/LogcoshPrior.cxx | 46 ++-- src/recon_buildblock/QuadraticPrior.cxx | 43 ++-- .../RelativeDifferencePrior.cxx | 87 ++++++- src/recon_test/test_priors.cxx | 213 ++++++++++++++---- 9 files changed, 363 insertions(+), 99 deletions(-) diff --git a/src/include/stir/recon_buildblock/GeneralisedPrior.h b/src/include/stir/recon_buildblock/GeneralisedPrior.h index 86135e5b8d..31a3249967 100644 --- a/src/include/stir/recon_buildblock/GeneralisedPrior.h +++ b/src/include/stir/recon_buildblock/GeneralisedPrior.h @@ -60,6 +60,17 @@ class GeneralisedPrior: virtual void compute_gradient(DataT& prior_gradient, const DataT ¤t_estimate) =0; + //! This should computes a single row of the Hessian + /*! Default implementation just call error(). This function needs to be overridden by the + derived class. + This method computes the Hessian of a particular voxel, indicated by \c coords, of the current image estimate. + The default method should overwrite the values in \c prior_Hessian_for_single_densel. + */ + virtual Succeeded + 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. @@ -89,6 +100,9 @@ class GeneralisedPrior: virtual Succeeded set_up(shared_ptr const& target_sptr); + //! Returns the status of the _is_convex variable + bool get_is_convex() const; + protected: float penalisation_factor; //! sets value for penalisation factor @@ -103,7 +117,7 @@ class GeneralisedPrior: bool _already_set_up; - //! Variable to indicate that the prior is a convex function + //! Variable to indicate that the prior is a convex function, should be set in defaults by the derived class bool _is_convex; }; diff --git a/src/include/stir/recon_buildblock/LogcoshPrior.h b/src/include/stir/recon_buildblock/LogcoshPrior.h index 574436300f..127e63288f 100644 --- a/src/include/stir/recon_buildblock/LogcoshPrior.h +++ b/src/include/stir/recon_buildblock/LogcoshPrior.h @@ -125,9 +125,10 @@ class LogcoshPrior: public const DiscretisedDensity<3,elemT> ¤t_image_estimate); //! compute Hessian - void compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, - const BasicCoordinate<3,int>& coords, - const DiscretisedDensity<3,elemT> ¤t_image_estimate); + virtual Succeeded + compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, + const BasicCoordinate<3,int>& coords, + const DiscretisedDensity<3,elemT> ¤t_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, @@ -224,19 +225,6 @@ class LogcoshPrior: public { return tanh(x)/x; } } - //! The Hessian of log(cosh()) is sech^2(x) = (1/cosh(x))^2 - /*! - 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 - */ - static inline float Hessian(const float d, const float scalar) - { - const float x = d * scalar; - return square((1/ cosh(x))); - } - //! The second partial derivatives of the LogCosh Prior /*! Diagonal refers to the second derivative w.r.t. x_j only (i.e. diagonal of the Hessian) @@ -244,7 +232,7 @@ class LogcoshPrior: public For LogCosh, the off diagonal is the negative of the diagonal. * @param x_j is the target voxel. * @param x_k is the voxel in the neighbourhood. - * @return the second partial derivative of the LogCosh Prior + * @return the second order partial derivatives of the LogCosh Prior */ //@{ float diagonal_second_derivative(const float x_j, const float x_k) const; diff --git a/src/include/stir/recon_buildblock/QuadraticPrior.h b/src/include/stir/recon_buildblock/QuadraticPrior.h index ca09af6f22..611358de7f 100644 --- a/src/include/stir/recon_buildblock/QuadraticPrior.h +++ b/src/include/stir/recon_buildblock/QuadraticPrior.h @@ -117,9 +117,10 @@ class QuadraticPrior: public const DiscretisedDensity<3,elemT> ¤t_image_estimate); //! compute Hessian - void compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, - const BasicCoordinate<3,int>& coords, - const DiscretisedDensity<3,elemT> ¤t_image_estimate); + virtual Succeeded + compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, + const BasicCoordinate<3,int>& coords, + const DiscretisedDensity<3,elemT> ¤t_image_estimate) const; //! Call accumulate_Hessian_times_input virtual Succeeded @@ -187,7 +188,7 @@ class QuadraticPrior: public For the Quadratic Prior, the off diagonal is the negative of the diagonal. * @param x_j is the target voxel. * @param x_k is the voxel in the neighbourhood. - * @return the second partial derivative of the Quadratic Prior + * @return the second order partial derivatives of the Quadratic Prior */ //@{ float diagonal_second_derivative(const float x_j, const float x_k) const; diff --git a/src/include/stir/recon_buildblock/RelativeDifferencePrior.h b/src/include/stir/recon_buildblock/RelativeDifferencePrior.h index 2e93b57784..d8a0707c33 100644 --- a/src/include/stir/recon_buildblock/RelativeDifferencePrior.h +++ b/src/include/stir/recon_buildblock/RelativeDifferencePrior.h @@ -125,6 +125,10 @@ class RelativeDifferencePrior: public void compute_gradient(DiscretisedDensity<3,elemT>& prior_gradient, const DiscretisedDensity<3,elemT> ¤t_image_estimate); + //! compute Hessian + virtual Succeeded compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, + const BasicCoordinate<3,int>& coords, + const DiscretisedDensity<3,elemT> ¤t_image_estimate) const; virtual Succeeded add_multiplication_with_approximate_Hessian(DiscretisedDensity<3,elemT>& output, @@ -203,9 +207,10 @@ class RelativeDifferencePrior: public /*! Diagonal refers to the second derivative w.r.t. x_j only (i.e. diagonal of the Hessian) Off-diagonal refers to the second derivative w.r.t. x_j and x_k (i.e. off-diagonal of the Hessian) + See J. Nuyts, et al., 2002, Equation 7. * @param x_j is the target voxel. * @param x_k is the voxel in the neighbourhood. - * @return the second partial derivative of the Relative Difference Prior + * @return the second order partial derivatives of the Relative Difference Prior */ //@{ float diagonal_second_derivative(const float x_j, const float x_k) const; diff --git a/src/recon_buildblock/GeneralisedPrior.cxx b/src/recon_buildblock/GeneralisedPrior.cxx index d62c520824..5ad4d81222 100644 --- a/src/recon_buildblock/GeneralisedPrior.cxx +++ b/src/recon_buildblock/GeneralisedPrior.cxx @@ -52,6 +52,20 @@ set_up(shared_ptr const&) return Succeeded::yes; } +template +Succeeded +GeneralisedPrior:: +compute_Hessian(TargetT& output, + const BasicCoordinate<3,int>& coords, + const TargetT& current_image_estimate) const +{ + if (this->get_is_convex()) + error("GeneralisedPrior:\n compute_Hessian implementation is not overloaded by your convex prior."); + else + error("GeneralisedPrior:\n compute_Hessian is not implemented because the prior is not convex."); + return Succeeded::no; +} + template Succeeded GeneralisedPrior:: @@ -82,6 +96,11 @@ void GeneralisedPrior::check(TargetT const& current_estimate) const error("The prior should already be set-up, but it's not."); } +template +bool GeneralisedPrior::get_is_convex() const { + return this->_is_convex; +} + # ifdef _MSC_VER // prevent warning message on instantiation of abstract class # pragma warning(disable:4661) diff --git a/src/recon_buildblock/LogcoshPrior.cxx b/src/recon_buildblock/LogcoshPrior.cxx index 43134b84e4..78fe8036d3 100644 --- a/src/recon_buildblock/LogcoshPrior.cxx +++ b/src/recon_buildblock/LogcoshPrior.cxx @@ -377,18 +377,21 @@ compute_gradient(DiscretisedDensity<3,elemT>& prior_gradient, } template -void +Succeeded LogcoshPrior:: compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, const BasicCoordinate<3,int>& coords, - const DiscretisedDensity<3,elemT> ¤t_image_estimate) + const DiscretisedDensity<3,elemT> ¤t_image_estimate) const { assert( prior_Hessian_for_single_densel.has_same_characteristics(current_image_estimate)); prior_Hessian_for_single_densel.fill(0); if (this->penalisation_factor==0) { - return; + return Succeeded::yes; } + + this->check(current_image_estimate); + const DiscretisedDensityOnCartesianGrid<3,elemT>& current_image_cast = dynamic_cast< const DiscretisedDensityOnCartesianGrid<3,elemT> &>(current_image_estimate); @@ -422,20 +425,37 @@ compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, for (int dy=min_dy;dy<=max_dy;++dy) for (int dx=min_dx;dx<=max_dx;++dx) { - // sech^2(x * scalar); sech(x) = 1/cosh(x) - elemT voxel_diff= current_image_estimate[z][y][x] - current_image_estimate[z+dz][y+dy][x+dx]; - elemT current = weights[dz][dy][dx] * Hessian(voxel_diff, this->scalar); - - if (do_kappa) - current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; - - diagonal += current; - prior_Hessian_for_single_densel_cast[z+dz][y+dy][x+dx] = -current*this->penalisation_factor; + elemT current = 0.0; + if (dz == 0 && dy == 0 && dx == 0) + { + // The j == k case (diagonal Hessian element), which is a sum over the neighbourhood. + for (int ddz=min_dz;ddz<=max_dz;++ddz) + for (int ddy=min_dy;ddy<=max_dy;++ddy) + for (int ddx=min_dx;ddx<=max_dx;++ddx) + { + elemT diagonal_current = weights[ddz][ddy][ddx] * + diagonal_second_derivative(current_image_estimate[z][y][x], + current_image_estimate[z+ddz][y+ddy][x+ddx]); + if (do_kappa) + diagonal_current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+ddz][y+ddy][x+ddx]; + current += diagonal_current; + } + } + else + { + // The j != k vases (off-diagonal Hessian elements) + current = weights[dz][dy][dx] * off_diagonal_second_derivative(current_image_estimate[z][y][x], + current_image_estimate[z+dz][y+dy][x+dx]); + if (do_kappa) + current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; + } + prior_Hessian_for_single_densel_cast[z+dz][y+dy][x+dx] = + current*this->penalisation_factor; } - prior_Hessian_for_single_densel[z][y][x]= diagonal * this->penalisation_factor; + return Succeeded::yes; } + template void LogcoshPrior::parabolic_surrogate_curvature(DiscretisedDensity<3,elemT>& parabolic_surrogate_curvature, diff --git a/src/recon_buildblock/QuadraticPrior.cxx b/src/recon_buildblock/QuadraticPrior.cxx index 657b8a7be6..d206f02071 100644 --- a/src/recon_buildblock/QuadraticPrior.cxx +++ b/src/recon_buildblock/QuadraticPrior.cxx @@ -431,17 +431,17 @@ compute_gradient(DiscretisedDensity<3,elemT>& prior_gradient, } template -void +Succeeded QuadraticPrior:: compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, const BasicCoordinate<3,int>& coords, - const DiscretisedDensity<3,elemT> ¤t_image_estimate) + const DiscretisedDensity<3,elemT> ¤t_image_estimate) const { assert( prior_Hessian_for_single_densel.has_same_characteristics(current_image_estimate)); prior_Hessian_for_single_densel.fill(0); if (this->penalisation_factor==0) { - return; + return Succeeded::yes; } this->check(current_image_estimate); @@ -480,19 +480,34 @@ compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, for (int dy=min_dy;dy<=max_dy;++dy) for (int dx=min_dx;dx<=max_dx;++dx) { - // dz==0,dy==0,dx==0 will have weight 0, so we can just include it in the loop - elemT current = - weights[dz][dy][dx]; - - if (do_kappa) - current *= - (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; - - diagonal += current; - prior_Hessian_for_single_densel_cast[z+dz][y+dy][x+dx] = -current*this->penalisation_factor; + elemT current = 0.0; + if (dz == 0 && dy == 0 && dx == 0) + { + // The j == k case (diagonal Hessian element), which is a sum over the neighbourhood. + for (int ddz=min_dz;ddz<=max_dz;++ddz) + for (int ddy=min_dy;ddy<=max_dy;++ddy) + for (int ddx=min_dx;ddx<=max_dx;++ddx) + { + elemT diagonal_current = weights[ddz][ddy][ddx] * + diagonal_second_derivative(current_image_estimate[z][y][x], + current_image_estimate[z+ddz][y+ddy][x+ddx]); + if (do_kappa) + diagonal_current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+ddz][y+ddy][x+ddx]; + current += diagonal_current; + } + } + else + { + // The j != k vases (off-diagonal Hessian elements) + current = weights[dz][dy][dx] * off_diagonal_second_derivative(current_image_estimate[z][y][x], + current_image_estimate[z+dz][y+dy][x+dx]); + if (do_kappa) + current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; + } + prior_Hessian_for_single_densel_cast[z+dz][y+dy][x+dx] = + current*this->penalisation_factor; } - prior_Hessian_for_single_densel[z][y][x]= diagonal * this->penalisation_factor; + return Succeeded::yes; } template diff --git a/src/recon_buildblock/RelativeDifferencePrior.cxx b/src/recon_buildblock/RelativeDifferencePrior.cxx index beb9557e7c..1fd446156b 100644 --- a/src/recon_buildblock/RelativeDifferencePrior.cxx +++ b/src/recon_buildblock/RelativeDifferencePrior.cxx @@ -179,6 +179,7 @@ template RelativeDifferencePrior::RelativeDifferencePrior(const bool only_2D_v, float penalisation_factor_v, float gamma_v, float epsilon_v) : only_2D(only_2D_v) { + set_defaults(); this->penalisation_factor = penalisation_factor_v; this->gamma = gamma_v; this->epsilon = epsilon_v; @@ -421,6 +422,86 @@ compute_gradient(DiscretisedDensity<3,elemT>& prior_gradient, } } +template +Succeeded +RelativeDifferencePrior:: +compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, + const BasicCoordinate<3,int>& coords, + const DiscretisedDensity<3,elemT> ¤t_image_estimate) const +{ + assert( prior_Hessian_for_single_densel.has_same_characteristics(current_image_estimate)); + prior_Hessian_for_single_densel.fill(0); + if (this->penalisation_factor==0) + { + return Succeeded::yes; + } + + this->check(current_image_estimate); + + const DiscretisedDensityOnCartesianGrid<3,elemT>& current_image_cast = + dynamic_cast< const DiscretisedDensityOnCartesianGrid<3,elemT> &>(current_image_estimate); + + DiscretisedDensityOnCartesianGrid<3,elemT>& prior_Hessian_for_single_densel_cast = + dynamic_cast &>(prior_Hessian_for_single_densel); + + if (weights.get_length() ==0) + { + compute_weights(weights, current_image_cast.get_grid_spacing(), this->only_2D); + } + + + 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]; + const int min_dz = max(weights.get_min_index(), prior_Hessian_for_single_densel.get_min_index()-z); + const int max_dz = min(weights.get_max_index(), prior_Hessian_for_single_densel.get_max_index()-z); + + const int min_dy = max(weights[0].get_min_index(), prior_Hessian_for_single_densel[z].get_min_index()-y); + const int max_dy = min(weights[0].get_max_index(), prior_Hessian_for_single_densel[z].get_max_index()-y); + + const int min_dx = max(weights[0][0].get_min_index(), prior_Hessian_for_single_densel[z][y].get_min_index()-x); + const int max_dx = min(weights[0][0].get_max_index(), prior_Hessian_for_single_densel[z][y].get_max_index()-x); + + elemT diagonal = 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 current = 0.0; + if (dz == 0 && dy == 0 && dx == 0) + { + // The j == k case (diagonal Hessian element), which is a sum over the neighbourhood. + for (int ddz=min_dz;ddz<=max_dz;++ddz) + for (int ddy=min_dy;ddy<=max_dy;++ddy) + for (int ddx=min_dx;ddx<=max_dx;++ddx) + { + elemT diagonal_current = weights[ddz][ddy][ddx] * + diagonal_second_derivative(current_image_estimate[z][y][x], + current_image_estimate[z+ddz][y+ddy][x+ddx]); + if (do_kappa) + diagonal_current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+ddz][y+ddy][x+ddx]; + current += diagonal_current; + } + } + else + { + // The j != k vases (off-diagonal Hessian elements) + current = weights[dz][dy][dx] * off_diagonal_second_derivative(current_image_estimate[z][y][x], + current_image_estimate[z+dz][y+dy][x+dx]); + if (do_kappa) + current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; + } + prior_Hessian_for_single_densel_cast[z+dz][y+dy][x+dx] = + current*this->penalisation_factor; + } + + return Succeeded::yes; +} + template Succeeded RelativeDifferencePrior:: @@ -527,9 +608,8 @@ float RelativeDifferencePrior:: diagonal_second_derivative(const float x_j, const float x_k) const { - // THIS MATH IS WRONG!!! if (x_j > 0.0 || x_k > 0.0 || this->epsilon > 0.0) - return 2 * (2 * x_j + this->epsilon)*(2 * x_k + this->epsilon) / + return 2 * pow(2 * x_k + this->epsilon, 2) / pow(x_j + x_k + this->gamma * abs(x_j - x_k) + this->epsilon, 3); else return 0.0; @@ -540,9 +620,8 @@ float RelativeDifferencePrior:: off_diagonal_second_derivative(const float x_j, const float x_k) const { - // Not sure if math is correct if (x_j > 0.0 || x_k > 0.0 || this->epsilon > 0.0) - return 2 * (2 * x_j + this->epsilon)*(2 * x_k + this->epsilon) / + return - 2 * (2 * x_j + this->epsilon)*(2 * x_k + this->epsilon) / pow(x_j + x_k + this->gamma * abs(x_j - x_k) + this->epsilon, 3); else return 0.0; diff --git a/src/recon_test/test_priors.cxx b/src/recon_test/test_priors.cxx index 3b8bebbf90..44ebb2366b 100644 --- a/src/recon_test/test_priors.cxx +++ b/src/recon_test/test_priors.cxx @@ -34,6 +34,7 @@ #include "stir/IO/read_from_file.h" #include "stir/IO/write_to_file.h" #include "stir/info.h" +#include "stir/Verbosity.h" #include "stir/Succeeded.h" #include "stir/num_threads.h" #include @@ -86,25 +87,34 @@ class GeneralisedPriorTests : public RunTests GeneralisedPrior& objective_function, shared_ptr target_sptr); - //! Test various configurations of the prior's Hessian via accumulate_Hessian_times_input() + //! Test various configurations of the Hessian of the prior via accumulate_Hessian_times_input() for convexity /*! - Tests the concave function condition + Tests the convexity condition: \f[ x^T \cdot H_{\lambda}x >= 0 \f] - for all non-negative \c x and non-zero \c \lambda (Relative Difference Penalty conditions). - This function constructs an array of configurations to test this condition and calls \c test_Hessian_configuration(). + for all non-negative \c x and non-zero \c \lambda (Relative Difference Prior conditions). + This function constructs an array of configurations to test this condition and calls + \c test_Hessian_convexity_configuration(). */ - void test_Hessian(const std::string& test_name, - GeneralisedPrior& objective_function, - shared_ptr target_sptr); + void test_Hessian_convexity(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr); + + //! Tests the compute_Hessian method implemented into convex priors + /*! Performs a perturbation response using compute_gradient to determine if the compute_Hessian (for a single densel) + is within tolerance. + */ + void test_Hessian_against_numerical(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr); private: //! Hessian test for a particular configuration of the Hessian concave condition - void test_Hessian_configuration(const std::string& test_name, - GeneralisedPrior& objective_function, - shared_ptr target_sptr, - float beta, - float input_multiplication, float input_addition, - float current_image_multiplication, float current_image_addition); + bool test_Hessian_convexity_configuration(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr, + const float beta, + const float input_multiplication, const float input_addition, + const float current_image_multiplication, const float current_image_addition); }; GeneralisedPriorTests:: @@ -122,10 +132,19 @@ run_tests_for_objective_function(const std::string& test_name, if (!check(objective_function.set_up(target_sptr)==Succeeded::yes, "set-up of objective function")) return; + // Gradient tests std::cerr << "----- test " << test_name << " --> Gradient\n"; test_gradient(test_name, objective_function, target_sptr); - std::cerr << "----- test " << test_name << " --> Hessian-vector product (accumulate_Hessian_times_input)\n"; - test_Hessian(test_name, objective_function, target_sptr); + + if (objective_function.get_is_convex()) + { + // The Hessian tests can only be performed on convex priors + std::cerr << "----- test " << test_name << " --> Hessian-vector product for convexity\n"; + test_Hessian_convexity(test_name, objective_function, target_sptr); + + std::cerr << "----- test " << test_name << " --> Hessian against numerical\n"; + test_Hessian_against_numerical(test_name, objective_function, target_sptr); + } } @@ -141,7 +160,10 @@ test_gradient(const std::string& test_name, shared_ptr gradient_2_sptr(target.get_empty_copy()); info("Computing gradient",3); + const int verbosity_default = Verbosity::get(); + Verbosity::set(0); objective_function.compute_gradient(*gradient_sptr, target); + Verbosity::set(verbosity_default); this->set_tolerance(std::max(fabs(double(gradient_sptr->find_min())), fabs(double(gradient_sptr->find_max()))) /1000); info("Computing objective function at target",3); @@ -170,6 +192,7 @@ test_gradient(const std::string& test_name, } if (!testOK) { + std::cerr << "Numerical gradient test failed with for " + test_name + " prior\n"; info("Writing diagnostic files gradient" + test_name + ".hv, numerical_gradient" + test_name + ".hv"); write_to_file("gradient" + test_name + ".hv", *gradient_sptr); write_to_file("numerical_gradient" + test_name + ".hv", *gradient_2_sptr); @@ -178,10 +201,12 @@ test_gradient(const std::string& test_name, void GeneralisedPriorTests:: -test_Hessian(const std::string& test_name, - GeneralisedPrior& objective_function, - shared_ptr target_sptr) +test_Hessian_convexity(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr) { + if (!objective_function.get_is_convex()) + return; /// Construct configurations float beta_array[] = {0.01, 1, 100}; // Penalty strength should only affect scale // Modifications to the input image @@ -191,35 +216,33 @@ test_Hessian(const std::string& test_name, float current_image_multiplication_array[] = {0.01, 1, 100}; float current_image_addition_array[] = {0.0, 0.5, 1, 10}; // RDP has constraint that current_image is non-negative + bool testOK = true; float initial_beta = objective_function.get_penalisation_factor(); - for (float beta : beta_array) { - for (float input_multiplication : input_multiplication_array) { - for (float input_addition : input_addition_array) { - for (float current_image_multiplication : current_image_multiplication_array) { + for (float beta : beta_array) + for (float input_multiplication : input_multiplication_array) + for (float input_addition : input_addition_array) + for (float current_image_multiplication : current_image_multiplication_array) for (float current_image_addition : current_image_addition_array) { - test_Hessian_configuration(test_name, objective_function, target_sptr, - beta, - input_multiplication, input_addition, - current_image_multiplication, current_image_addition); + if (testOK) // only compute configuration if testOK from previous tests + testOK = test_Hessian_convexity_configuration(test_name, objective_function, target_sptr, + beta, + input_multiplication, input_addition, + current_image_multiplication, current_image_addition); } - } - } - } - } /// Reset beta to original value objective_function.set_penalisation_factor(initial_beta); } -void +bool GeneralisedPriorTests:: -test_Hessian_configuration(const std::string& test_name, - GeneralisedPrior& objective_function, - shared_ptr target_sptr, - const float beta, - const float input_multiplication, const float input_addition, - const float current_image_multiplication, const float current_image_addition) +test_Hessian_convexity_configuration(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr, + const float beta, + const float input_multiplication, const float input_addition, + const float current_image_multiplication, const float current_image_addition) { - /// setup images + /// setup targets target_type& target(*target_sptr); shared_ptr output(target.get_empty_copy()); shared_ptr current_image(target.get_empty_copy()); @@ -230,7 +253,7 @@ test_Hessian_configuration(const std::string& test_name, target_type::full_iterator input_iter = input->begin_all(); target_type::full_iterator current_image_iter = current_image->begin_all(); target_type::full_iterator target_iter = target.begin_all(); - while (input_iter != input->end_all())// && testOK) + while (input_iter != input->end_all()) { *input_iter = input_multiplication * *target_iter + input_addition; *current_image_iter = current_image_multiplication * *target_iter + current_image_addition; @@ -248,6 +271,7 @@ test_Hessian_configuration(const std::string& test_name, // test for a CONVEX function if (this->check_if_less(0, my_sum)) { // info("PASS: Computation of x^T H x = " + std::to_string(my_sum) + " > 0 and is therefore convex"); + return true; } else { // print to console the FAILED configuration info("FAIL: Computation of x^T H x = " + std::to_string(my_sum) + " < 0 and is therefore NOT convex" + @@ -261,9 +285,102 @@ test_Hessian_configuration(const std::string& test_name, "\n >input image min=" + std::to_string(input->find_min()) + "\n >target image max=" + std::to_string(target.find_max()) + "\n >target image min=" + std::to_string(target.find_min())); + return false; } } + +void +GeneralisedPriorTests:: +test_Hessian_against_numerical(const std::string &test_name, + GeneralisedPrior &objective_function, + shared_ptr target_sptr) +{ + if (!objective_function.get_is_convex()) + return; + + /// Setup + const float eps = 1e-3F; + bool testOK = true; + const int verbosity_default = Verbosity::get(); + + // setup images + target_type& input(*target_sptr->get_empty_copy()); + input += *target_sptr; // make input have same values as target_sptr + shared_ptr gradient_sptr(target_sptr->get_empty_copy()); + shared_ptr pert_grad_and_numerical_Hessian_sptr(target_sptr->get_empty_copy()); + shared_ptr Hessian_sptr(target_sptr->get_empty_copy()); + + Verbosity::set(0); + objective_function.compute_gradient(*gradient_sptr, input); + Verbosity::set(verbosity_default); +// this->set_tolerance(std::max(fabs(double(gradient_sptr->find_min())), fabs(double(gradient_sptr->find_max()))) /10); + + // Setup coordinates (z,y,x) for perturbation test (Hessian will also be computed w.r.t this voxel, j) + BasicCoordinate<3, int> perturbation_coords; + + // Get min/max indices + const int min_z = input.get_min_index(); + const int max_z = input.get_max_index(); + const int min_y = input[min_z].get_min_index(); + const int max_y = input[min_z].get_max_index(); + const int min_x = input[min_z][min_y].get_min_index(); + const int max_x = input[min_z][min_y].get_max_index(); + + // Loop over each voxel j in the input and check perturbation response. + for (int z=min_z;z<= max_z;z++) + for (int y=min_y;y<= max_y;y++) + for (int x=min_x;x<= max_x;x++) + if (testOK) + { + perturbation_coords[1] = z; perturbation_coords[2] = y; perturbation_coords[3] = x; + + // Compute H(x)_j (row of the Hessian at the jth voxel) + objective_function.compute_Hessian(*Hessian_sptr, perturbation_coords, input); + this->set_tolerance(std::max(fabs(double(Hessian_sptr->find_min())), fabs(double(Hessian_sptr->find_max()))) /500); + + // Compute g(x + eps) + Verbosity::set(0); + // Perturb target at jth voxel, compute perturbed gradient, and reset voxel to original value + float perturbed_voxels_original_value = input[perturbation_coords[1]][perturbation_coords[2]][perturbation_coords[3]]; + input[perturbation_coords[1]][perturbation_coords[2]][perturbation_coords[3]] += eps; + objective_function.compute_gradient(*pert_grad_and_numerical_Hessian_sptr, input); + input[perturbation_coords[1]][perturbation_coords[2]][perturbation_coords[3]] = perturbed_voxels_original_value; + + // Now compute the numerical-Hessian = (g(x+eps) - g(x))/eps + *pert_grad_and_numerical_Hessian_sptr -= *gradient_sptr; + *pert_grad_and_numerical_Hessian_sptr /= eps; + + Verbosity::set(verbosity_default); + // Test if pert_grad_and_numerical_Hessian_sptr is all zeros. + // This can happen if the eps is too small. This is a quick test that allows for easier debugging. + if (pert_grad_and_numerical_Hessian_sptr->sum_positive() == 0.0) + { + this->everything_ok = false; + testOK = false; + info("test_Hessian_against_numerical: failed because all values are 0 in numerical Hessian"); + } + + // Loop over each of the voxels and compare the numerical-Hessian with Hessian + target_type::full_iterator numerical_Hessian_iter = pert_grad_and_numerical_Hessian_sptr->begin_all(); + target_type::full_iterator Hessian_iter = Hessian_sptr->begin_all(); + while(numerical_Hessian_iter != pert_grad_and_numerical_Hessian_sptr->end_all() && testOK) + { + testOK = testOK && this->check_if_equal(*Hessian_iter, *numerical_Hessian_iter, "Hessian"); + ++numerical_Hessian_iter; ++ Hessian_iter; + } + + if (!testOK) + { + std::cerr << "Numerical-Hessian test failed with for " + test_name + " prior\n"; + info("Writing diagnostic files `Hessian_" + test_name + ".hv` and `numerical_Hessian_" + test_name + ".hv`"); + write_to_file("Hessian_" + test_name + ".hv", *Hessian_sptr); + write_to_file("numerical_Hessian_" + test_name + ".hv", *pert_grad_and_numerical_Hessian_sptr); + write_to_file("input_" + test_name + ".hv", input); + } + } +} + void GeneralisedPriorTests:: construct_input_data(shared_ptr& density_sptr) @@ -305,19 +422,25 @@ run_tests() shared_ptr density_sptr; construct_input_data(density_sptr); - std::cerr << "Tests for QuadraticPrior\n"; + std::cerr << "\n\nTests for QuadraticPrior\n"; { QuadraticPrior objective_function(false, 1.F); this->run_tests_for_objective_function("Quadratic_no_kappa", objective_function, density_sptr); } - std::cerr << "Tests for Relative Difference Prior\n"; + std::cerr << "\n\nTests for Relative Difference Prior with epsilon = 0\n"; { - // gamma is default and epsilon is off + // gamma is default and epsilon is 0.0 RelativeDifferencePrior objective_function(false, 1.F, 2.F, 0.F); - this->run_tests_for_objective_function("RDP_no_kappa", objective_function, density_sptr); + this->run_tests_for_objective_function("RDP_no_kappa_no_eps", objective_function, density_sptr); + } + std::cerr << "\n\nTests for Relative Difference Prior with epsilon = 0.001\n"; + { + // gamma is default and epsilon is "small" + RelativeDifferencePrior objective_function(false, 1.F, 2.F, 0.0001F); + this->run_tests_for_objective_function("RDP_no_kappa_with_eps", objective_function, density_sptr); } // Disabled PLS due to known issue -// std::cerr << "Tests for PLSPrior\n"; +// std::cerr << "\n\nTests for PLSPrior\n"; // { // PLSPrior objective_function(false, 1.F); // shared_ptr > anatomical_image_sptr(density_sptr->get_empty_copy()); @@ -325,7 +448,7 @@ run_tests() // objective_function.set_anatomical_image_sptr(anatomical_image_sptr); // this->run_tests_for_objective_function("PLS_no_kappa_flat_anatomical", objective_function, density_sptr); // } - std::cerr << "Tests for Logcosh Prior\n"; + std::cerr << "\n\nTests for Logcosh Prior\n"; { // scalar is off LogcoshPrior objective_function(false, 1.F, 1.F); From d767044156b7069918e9bb104e858214362fb235 Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Sat, 19 Jun 2021 21:10:49 +0100 Subject: [PATCH 15/27] [ci skip] Add documentation to `test_gradient` --- src/recon_test/test_priors.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/src/recon_test/test_priors.cxx b/src/recon_test/test_priors.cxx index 44ebb2366b..0c39fe77d6 100644 --- a/src/recon_test/test_priors.cxx +++ b/src/recon_test/test_priors.cxx @@ -83,6 +83,7 @@ class GeneralisedPriorTests : public RunTests GeneralisedPrior& objective_function, shared_ptr target_sptr); + //! Tests the prior's gradient by comparing to the numerical gradient computed using perturbation response. void test_gradient(const std::string& test_name, GeneralisedPrior& objective_function, shared_ptr target_sptr); From c2ccf06e9eff836ce850226f222d838ad3e1ab77 Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Thu, 24 Jun 2021 23:23:51 +0100 Subject: [PATCH 16/27] Increase RDP epsilon to be 0.1, increase so well above test `eps` value --- src/recon_test/test_priors.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/recon_test/test_priors.cxx b/src/recon_test/test_priors.cxx index 0c39fe77d6..6433acfe5e 100644 --- a/src/recon_test/test_priors.cxx +++ b/src/recon_test/test_priors.cxx @@ -434,10 +434,10 @@ run_tests() RelativeDifferencePrior objective_function(false, 1.F, 2.F, 0.F); this->run_tests_for_objective_function("RDP_no_kappa_no_eps", objective_function, density_sptr); } - std::cerr << "\n\nTests for Relative Difference Prior with epsilon = 0.001\n"; + std::cerr << "\n\nTests for Relative Difference Prior with epsilon = 0.1\n"; { // gamma is default and epsilon is "small" - RelativeDifferencePrior objective_function(false, 1.F, 2.F, 0.0001F); + RelativeDifferencePrior objective_function(false, 1.F, 2.F, 0.1F); this->run_tests_for_objective_function("RDP_no_kappa_with_eps", objective_function, density_sptr); } // Disabled PLS due to known issue From f60d97cae55366200572985a00d33a5dd0eb695b Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Thu, 24 Jun 2021 23:24:37 +0100 Subject: [PATCH 17/27] Add variables and set methods for toggling tests Reenable PLS, but disable all numerical tests - just run setup. --- src/recon_test/test_priors.cxx | 86 ++++++++++++++++++++++++++++------ 1 file changed, 71 insertions(+), 15 deletions(-) diff --git a/src/recon_test/test_priors.cxx b/src/recon_test/test_priors.cxx index 6433acfe5e..d56f012b5c 100644 --- a/src/recon_test/test_priors.cxx +++ b/src/recon_test/test_priors.cxx @@ -29,7 +29,7 @@ #include "stir/recon_buildblock/QuadraticPrior.h" #include "stir/recon_buildblock/RelativeDifferencePrior.h" #include "stir/recon_buildblock/LogcoshPrior.h" -//#include "stir/recon_buildblock/PLSPrior.h" +#include "stir/recon_buildblock/PLSPrior.h" #include "stir/RunTests.h" #include "stir/IO/read_from_file.h" #include "stir/IO/write_to_file.h" @@ -73,6 +73,14 @@ class GeneralisedPriorTests : public RunTests void construct_input_data(shared_ptr& density_sptr); void run_tests(); + + //! Set methods that control which tests are run. + //@{ + void set_do_test_gradient(bool d); + void set_do_test_Hessian_convexity(bool d); + void set_do_test_Hessian_against_numerical(bool d); + //@} + protected: char const * density_filename; shared_ptr > objective_function_sptr; @@ -116,6 +124,13 @@ class GeneralisedPriorTests : public RunTests const float beta, const float input_multiplication, const float input_addition, const float current_image_multiplication, const float current_image_addition); + + //! Variables to control which tests are run, see the set methods + //@{ + bool do_test_gradient = false; + bool do_test_Hessian_convexity = false; + bool do_test_Hessian_against_numerical = false; + //@} }; GeneralisedPriorTests:: @@ -123,6 +138,27 @@ GeneralisedPriorTests(char const * const density_filename) : density_filename(density_filename) {} +void +GeneralisedPriorTests:: +set_do_test_gradient(const bool d) +{ + do_test_gradient = d; +} + +void +GeneralisedPriorTests:: +set_do_test_Hessian_convexity(const bool d) +{ + do_test_Hessian_convexity = d; +} + +void +GeneralisedPriorTests:: +set_do_test_Hessian_against_numerical(const bool d) +{ + do_test_Hessian_against_numerical = d; +} + void GeneralisedPriorTests:: run_tests_for_objective_function(const std::string& test_name, @@ -133,16 +169,20 @@ run_tests_for_objective_function(const std::string& test_name, if (!check(objective_function.set_up(target_sptr)==Succeeded::yes, "set-up of objective function")) return; - // Gradient tests - std::cerr << "----- test " << test_name << " --> Gradient\n"; - test_gradient(test_name, objective_function, target_sptr); + if (do_test_gradient) + { + std::cerr << "----- test " << test_name << " --> Gradient\n"; + test_gradient(test_name, objective_function, target_sptr); + } - if (objective_function.get_is_convex()) + if (do_test_Hessian_convexity) { - // The Hessian tests can only be performed on convex priors std::cerr << "----- test " << test_name << " --> Hessian-vector product for convexity\n"; test_Hessian_convexity(test_name, objective_function, target_sptr); + } + if (do_test_Hessian_against_numerical) + { std::cerr << "----- test " << test_name << " --> Hessian against numerical\n"; test_Hessian_against_numerical(test_name, objective_function, target_sptr); } @@ -426,33 +466,49 @@ run_tests() std::cerr << "\n\nTests for QuadraticPrior\n"; { QuadraticPrior objective_function(false, 1.F); + set_do_test_gradient(true); + set_do_test_Hessian_convexity(true); + set_do_test_Hessian_against_numerical(true); this->run_tests_for_objective_function("Quadratic_no_kappa", objective_function, density_sptr); } std::cerr << "\n\nTests for Relative Difference Prior with epsilon = 0\n"; { // gamma is default and epsilon is 0.0 RelativeDifferencePrior objective_function(false, 1.F, 2.F, 0.F); + set_do_test_gradient(true); + set_do_test_Hessian_convexity(true); + set_do_test_Hessian_against_numerical(false); // RDP, with epsilon = 0.0, will fail the numerical Hessian test this->run_tests_for_objective_function("RDP_no_kappa_no_eps", objective_function, density_sptr); } std::cerr << "\n\nTests for Relative Difference Prior with epsilon = 0.1\n"; { // gamma is default and epsilon is "small" RelativeDifferencePrior objective_function(false, 1.F, 2.F, 0.1F); + set_do_test_gradient(true); + set_do_test_Hessian_convexity(true); + set_do_test_Hessian_against_numerical(true); // With a large enough epsilon the RDP numerical test will pass this->run_tests_for_objective_function("RDP_no_kappa_with_eps", objective_function, density_sptr); } - // Disabled PLS due to known issue -// std::cerr << "\n\nTests for PLSPrior\n"; -// { -// PLSPrior objective_function(false, 1.F); -// shared_ptr > anatomical_image_sptr(density_sptr->get_empty_copy()); -// anatomical_image_sptr->fill(1.F); -// objective_function.set_anatomical_image_sptr(anatomical_image_sptr); -// this->run_tests_for_objective_function("PLS_no_kappa_flat_anatomical", objective_function, density_sptr); -// } + + std::cerr << "\n\nTests for PLSPrior\n"; + { + PLSPrior objective_function(false, 1.F); + shared_ptr > anatomical_image_sptr(density_sptr->get_empty_copy()); + anatomical_image_sptr->fill(1.F); + objective_function.set_anatomical_image_sptr(anatomical_image_sptr); + // Disabled PLS due to known issue + set_do_test_gradient(false); + set_do_test_Hessian_convexity(false); + set_do_test_Hessian_against_numerical(false); + this->run_tests_for_objective_function("PLS_no_kappa_flat_anatomical", objective_function, density_sptr); + } std::cerr << "\n\nTests for Logcosh Prior\n"; { // scalar is off LogcoshPrior objective_function(false, 1.F, 1.F); + set_do_test_gradient(true); + set_do_test_Hessian_convexity(true); + set_do_test_Hessian_against_numerical(true); this->run_tests_for_objective_function("Logcosh_no_kappa", objective_function, density_sptr); } } From 6669186592d986110fbdd7f599363f35aa1344d9 Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Thu, 24 Jun 2021 23:25:32 +0100 Subject: [PATCH 18/27] Cleanup and improve debugging checks --- src/recon_test/test_priors.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/recon_test/test_priors.cxx b/src/recon_test/test_priors.cxx index d56f012b5c..b3f7d0a1e7 100644 --- a/src/recon_test/test_priors.cxx +++ b/src/recon_test/test_priors.cxx @@ -395,7 +395,7 @@ test_Hessian_against_numerical(const std::string &test_name, Verbosity::set(verbosity_default); // Test if pert_grad_and_numerical_Hessian_sptr is all zeros. // This can happen if the eps is too small. This is a quick test that allows for easier debugging. - if (pert_grad_and_numerical_Hessian_sptr->sum_positive() == 0.0) + if (pert_grad_and_numerical_Hessian_sptr->sum_positive() == 0.0 && Hessian_sptr->sum_positive() > 0.0) { this->everything_ok = false; testOK = false; @@ -413,6 +413,7 @@ test_Hessian_against_numerical(const std::string &test_name, if (!testOK) { + // Output volumes for debug std::cerr << "Numerical-Hessian test failed with for " + test_name + " prior\n"; info("Writing diagnostic files `Hessian_" + test_name + ".hv` and `numerical_Hessian_" + test_name + ".hv`"); write_to_file("Hessian_" + test_name + ".hv", *Hessian_sptr); From 511c69b3a9445b50e36aaa790c82572ee246ea22 Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Mon, 19 Jul 2021 11:12:52 -0700 Subject: [PATCH 19/27] Remove _is_convex variable and better implement is_convex() method Make is_convex() a pure virtual method in GeneralisedPrior to be implemented in each prior. FilterRoot is not convex as it does not have 0th or 2nd order behaviour --- src/include/stir/recon_buildblock/FilterRootPrior.h | 2 ++ src/include/stir/recon_buildblock/GeneralisedPrior.h | 8 +++----- src/include/stir/recon_buildblock/LogcoshPrior.h | 2 ++ src/include/stir/recon_buildblock/PLSPrior.h | 2 ++ src/include/stir/recon_buildblock/QuadraticPrior.h | 2 ++ .../stir/recon_buildblock/RelativeDifferencePrior.h | 2 ++ src/recon_buildblock/FilterRootPrior.cxx | 6 ++++++ src/recon_buildblock/GeneralisedPrior.cxx | 8 +------- src/recon_buildblock/LogcoshPrior.cxx | 9 ++++++++- src/recon_buildblock/PLSPrior.cxx | 8 +++++++- src/recon_buildblock/QuadraticPrior.cxx | 8 +++++++- src/recon_buildblock/RelativeDifferencePrior.cxx | 10 +++++++++- src/recon_test/test_priors.cxx | 4 ++-- 13 files changed, 53 insertions(+), 18 deletions(-) diff --git a/src/include/stir/recon_buildblock/FilterRootPrior.h b/src/include/stir/recon_buildblock/FilterRootPrior.h index 979bc0c496..e2850f9ab3 100644 --- a/src/include/stir/recon_buildblock/FilterRootPrior.h +++ b/src/include/stir/recon_buildblock/FilterRootPrior.h @@ -87,6 +87,8 @@ class FilterRootPrior: public FilterRootPrior(shared_ptr >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. diff --git a/src/include/stir/recon_buildblock/GeneralisedPrior.h b/src/include/stir/recon_buildblock/GeneralisedPrior.h index 31a3249967..e4ac84b73d 100644 --- a/src/include/stir/recon_buildblock/GeneralisedPrior.h +++ b/src/include/stir/recon_buildblock/GeneralisedPrior.h @@ -100,8 +100,9 @@ class GeneralisedPrior: virtual Succeeded set_up(shared_ptr const& target_sptr); - //! Returns the status of the _is_convex variable - bool get_is_convex() const; + //! 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; @@ -116,9 +117,6 @@ class GeneralisedPrior: virtual void check(DataT const& current_estimate) const; bool _already_set_up; - - //! Variable to indicate that the prior is a convex function, should be set in defaults by the derived class - bool _is_convex; }; END_NAMESPACE_STIR diff --git a/src/include/stir/recon_buildblock/LogcoshPrior.h b/src/include/stir/recon_buildblock/LogcoshPrior.h index 127e63288f..0581ad1a66 100644 --- a/src/include/stir/recon_buildblock/LogcoshPrior.h +++ b/src/include/stir/recon_buildblock/LogcoshPrior.h @@ -112,6 +112,8 @@ class LogcoshPrior: 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> ¤t_image_estimate); diff --git a/src/include/stir/recon_buildblock/PLSPrior.h b/src/include/stir/recon_buildblock/PLSPrior.h index 91a8d83c0c..0c1458abf8 100644 --- a/src/include/stir/recon_buildblock/PLSPrior.h +++ b/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& target_sptr); + bool is_convex() const; + //! compute the value of the function double compute_value(const DiscretisedDensity<3,elemT> ¤t_image_estimate); diff --git a/src/include/stir/recon_buildblock/QuadraticPrior.h b/src/include/stir/recon_buildblock/QuadraticPrior.h index 611358de7f..649d33bc01 100644 --- a/src/include/stir/recon_buildblock/QuadraticPrior.h +++ b/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> ¤t_image_estimate); diff --git a/src/include/stir/recon_buildblock/RelativeDifferencePrior.h b/src/include/stir/recon_buildblock/RelativeDifferencePrior.h index d8a0707c33..12cfbff9a5 100644 --- a/src/include/stir/recon_buildblock/RelativeDifferencePrior.h +++ b/src/include/stir/recon_buildblock/RelativeDifferencePrior.h @@ -169,6 +169,8 @@ class RelativeDifferencePrior: public //! Has to be called before using this object virtual Succeeded set_up(shared_ptr > const& target_sptr); + bool is_convex() const; + protected: //! Create variable gamma for Relative Difference Penalty float gamma; diff --git a/src/recon_buildblock/FilterRootPrior.cxx b/src/recon_buildblock/FilterRootPrior.cxx index 90ebaa8ece..638bdcd067 100644 --- a/src/recon_buildblock/FilterRootPrior.cxx +++ b/src/recon_buildblock/FilterRootPrior.cxx @@ -41,6 +41,12 @@ FilterRootPrior(shared_ptr >const& filter_sptr, float penal this->penalisation_factor = penalisation_factor_v; } +template +bool FilterRootPrior:: +is_convex() const +{ + return false; +} template < class T> static inline int diff --git a/src/recon_buildblock/GeneralisedPrior.cxx b/src/recon_buildblock/GeneralisedPrior.cxx index 5ad4d81222..58c8b527f1 100644 --- a/src/recon_buildblock/GeneralisedPrior.cxx +++ b/src/recon_buildblock/GeneralisedPrior.cxx @@ -40,7 +40,6 @@ GeneralisedPrior::set_defaults() { _already_set_up = false; this->penalisation_factor = 0; - this->_is_convex = false; } template @@ -59,7 +58,7 @@ compute_Hessian(TargetT& output, const BasicCoordinate<3,int>& coords, const TargetT& current_image_estimate) const { - if (this->get_is_convex()) + 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 because the prior is not convex."); @@ -96,11 +95,6 @@ void GeneralisedPrior::check(TargetT const& current_estimate) const error("The prior should already be set-up, but it's not."); } -template -bool GeneralisedPrior::get_is_convex() const { - return this->_is_convex; -} - # ifdef _MSC_VER // prevent warning message on instantiation of abstract class # pragma warning(disable:4661) diff --git a/src/recon_buildblock/LogcoshPrior.cxx b/src/recon_buildblock/LogcoshPrior.cxx index 78fe8036d3..1fee3ab7a2 100644 --- a/src/recon_buildblock/LogcoshPrior.cxx +++ b/src/recon_buildblock/LogcoshPrior.cxx @@ -108,7 +108,6 @@ void LogcoshPrior::set_defaults() { base_type::set_defaults(); - this->_is_convex = true; this->only_2D = false; this->scalar = 1.0; this->kappa_ptr.reset(); @@ -145,6 +144,14 @@ LogcoshPrior::LogcoshPrior(const bool only_2D_v, float penalisation_facto this->scalar = scalar_v; } +template +bool +LogcoshPrior:: +is_convex() const +{ + return true; +} + //! get penalty weights for the neighbourhood template Array<3,float> diff --git a/src/recon_buildblock/PLSPrior.cxx b/src/recon_buildblock/PLSPrior.cxx index daa09ca30c..891d554fd5 100644 --- a/src/recon_buildblock/PLSPrior.cxx +++ b/src/recon_buildblock/PLSPrior.cxx @@ -127,7 +127,6 @@ void PLSPrior::set_defaults() { base_type::set_defaults(); - this->_is_convex = true; this->only_2D = false; this->alpha=1; this->eta=1; @@ -145,6 +144,13 @@ PLSPrior::PLSPrior() set_defaults(); } +template +bool +PLSPrior:: +is_convex() const +{ + return true; +} template PLSPrior::PLSPrior(const bool only_2D_v, float penalisation_factor_v) diff --git a/src/recon_buildblock/QuadraticPrior.cxx b/src/recon_buildblock/QuadraticPrior.cxx index d206f02071..8563c2e1f0 100644 --- a/src/recon_buildblock/QuadraticPrior.cxx +++ b/src/recon_buildblock/QuadraticPrior.cxx @@ -123,7 +123,6 @@ void QuadraticPrior::set_defaults() { base_type::set_defaults(); - this->_is_convex = true; this->only_2D = false; this->kappa_ptr.reset(); this->weights.recycle(); @@ -148,6 +147,13 @@ QuadraticPrior::QuadraticPrior(const bool only_2D_v, float penalisation_f this->penalisation_factor = penalisation_factor_v; } +template +bool +QuadraticPrior:: +is_convex() const +{ + return true; +} //! get penalty weights for the neigbourhood template diff --git a/src/recon_buildblock/RelativeDifferencePrior.cxx b/src/recon_buildblock/RelativeDifferencePrior.cxx index 1fd446156b..c9d50211b5 100644 --- a/src/recon_buildblock/RelativeDifferencePrior.cxx +++ b/src/recon_buildblock/RelativeDifferencePrior.cxx @@ -127,7 +127,7 @@ void RelativeDifferencePrior::set_defaults() { base_type::set_defaults(); - this->_is_convex = true; +// this->_is_convex = true; this->only_2D = false; this->kappa_ptr.reset(); this->weights.recycle(); @@ -146,6 +146,14 @@ RelativeDifferencePrior::RelativeDifferencePrior() set_defaults(); } +template +bool +RelativeDifferencePrior:: +is_convex() const +{ + return true; +} + // Return the value of gamma - a RDP parameter template float diff --git a/src/recon_test/test_priors.cxx b/src/recon_test/test_priors.cxx index b3f7d0a1e7..b7ce6f821d 100644 --- a/src/recon_test/test_priors.cxx +++ b/src/recon_test/test_priors.cxx @@ -246,7 +246,7 @@ test_Hessian_convexity(const std::string& test_name, GeneralisedPrior& objective_function, shared_ptr target_sptr) { - if (!objective_function.get_is_convex()) + if (!objective_function.is_convex()) return; /// Construct configurations float beta_array[] = {0.01, 1, 100}; // Penalty strength should only affect scale @@ -337,7 +337,7 @@ test_Hessian_against_numerical(const std::string &test_name, GeneralisedPrior &objective_function, shared_ptr target_sptr) { - if (!objective_function.get_is_convex()) + if (!objective_function.is_convex()) return; /// Setup From 381b6c390d56aaa713ff12da5fcbd3bac530392d Mon Sep 17 00:00:00 2001 From: robbietuk Date: Mon, 19 Jul 2021 11:18:06 -0700 Subject: [PATCH 20/27] Update src/include/stir/recon_buildblock/GeneralisedPrior.h Co-authored-by: Kris Thielemans --- src/include/stir/recon_buildblock/GeneralisedPrior.h | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/include/stir/recon_buildblock/GeneralisedPrior.h b/src/include/stir/recon_buildblock/GeneralisedPrior.h index 31a3249967..190b217427 100644 --- a/src/include/stir/recon_buildblock/GeneralisedPrior.h +++ b/src/include/stir/recon_buildblock/GeneralisedPrior.h @@ -60,11 +60,13 @@ class GeneralisedPrior: virtual void compute_gradient(DataT& prior_gradient, const DataT ¤t_estimate) =0; - //! This should computes a single row of the Hessian + //! This computes a single row of the Hessian /*! Default implementation just call error(). This function needs to be overridden by the derived class. - This method computes the Hessian of a particular voxel, indicated by \c coords, of the current image estimate. - The default method should overwrite the values in \c prior_Hessian_for_single_densel. + + 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 Succeeded compute_Hessian(DataT& prior_Hessian_for_single_densel, From 610c78457904dad109d7f1b0379002b707293bea Mon Sep 17 00:00:00 2001 From: robbietuk Date: Mon, 19 Jul 2021 11:18:20 -0700 Subject: [PATCH 21/27] Update src/include/stir/recon_buildblock/GeneralisedPrior.h Co-authored-by: Kris Thielemans --- src/include/stir/recon_buildblock/GeneralisedPrior.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/include/stir/recon_buildblock/GeneralisedPrior.h b/src/include/stir/recon_buildblock/GeneralisedPrior.h index 190b217427..46e348f9cd 100644 --- a/src/include/stir/recon_buildblock/GeneralisedPrior.h +++ b/src/include/stir/recon_buildblock/GeneralisedPrior.h @@ -102,7 +102,7 @@ class GeneralisedPrior: virtual Succeeded set_up(shared_ptr const& target_sptr); - //! Returns the status of the _is_convex variable + //! Returns if the prior is convex bool get_is_convex() const; protected: From 12f3eb691ddaaa3a0cd499287b19a55870828423 Mon Sep 17 00:00:00 2001 From: robbietuk Date: Mon, 19 Jul 2021 11:23:55 -0700 Subject: [PATCH 22/27] Update src/include/stir/recon_buildblock/LogcoshPrior.h Co-authored-by: Kris Thielemans --- src/include/stir/recon_buildblock/LogcoshPrior.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/include/stir/recon_buildblock/LogcoshPrior.h b/src/include/stir/recon_buildblock/LogcoshPrior.h index 127e63288f..f7fe81eca6 100644 --- a/src/include/stir/recon_buildblock/LogcoshPrior.h +++ b/src/include/stir/recon_buildblock/LogcoshPrior.h @@ -130,7 +130,7 @@ class LogcoshPrior: public const BasicCoordinate<3,int>& coords, const DiscretisedDensity<3,elemT> ¤t_image_estimate) const; - //! Compute the multiplication of the hessian of the prior multiplied by the input. + //! Compute the multiplication of the hessian of the prior (at \c current_estimate) and the \c input. virtual Succeeded accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, const DiscretisedDensity<3,elemT>& current_estimate, const DiscretisedDensity<3,elemT>& input) const; @@ -243,4 +243,4 @@ class LogcoshPrior: public END_NAMESPACE_STIR -#endif \ No newline at end of file +#endif From c68f4d952015c6ca95520168bae6f04bb3324454 Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Mon, 19 Jul 2021 12:59:19 -0700 Subject: [PATCH 23/27] [ci skip] Remove unnecessary documentation --- src/include/stir/recon_buildblock/LogcoshPrior.h | 1 - src/include/stir/recon_buildblock/QuadraticPrior.h | 1 - src/include/stir/recon_buildblock/RelativeDifferencePrior.h | 1 - 3 files changed, 3 deletions(-) diff --git a/src/include/stir/recon_buildblock/LogcoshPrior.h b/src/include/stir/recon_buildblock/LogcoshPrior.h index 9dcd666b63..5f4eb7ae7a 100644 --- a/src/include/stir/recon_buildblock/LogcoshPrior.h +++ b/src/include/stir/recon_buildblock/LogcoshPrior.h @@ -126,7 +126,6 @@ class LogcoshPrior: public void parabolic_surrogate_curvature(DiscretisedDensity<3,elemT>& parabolic_surrogate_curvature, const DiscretisedDensity<3,elemT> ¤t_image_estimate); - //! compute Hessian virtual Succeeded compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, const BasicCoordinate<3,int>& coords, diff --git a/src/include/stir/recon_buildblock/QuadraticPrior.h b/src/include/stir/recon_buildblock/QuadraticPrior.h index 649d33bc01..4d2f9560f8 100644 --- a/src/include/stir/recon_buildblock/QuadraticPrior.h +++ b/src/include/stir/recon_buildblock/QuadraticPrior.h @@ -118,7 +118,6 @@ class QuadraticPrior: public void parabolic_surrogate_curvature(DiscretisedDensity<3,elemT>& parabolic_surrogate_curvature, const DiscretisedDensity<3,elemT> ¤t_image_estimate); - //! compute Hessian virtual Succeeded compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, const BasicCoordinate<3,int>& coords, diff --git a/src/include/stir/recon_buildblock/RelativeDifferencePrior.h b/src/include/stir/recon_buildblock/RelativeDifferencePrior.h index 12cfbff9a5..537cd94ce5 100644 --- a/src/include/stir/recon_buildblock/RelativeDifferencePrior.h +++ b/src/include/stir/recon_buildblock/RelativeDifferencePrior.h @@ -125,7 +125,6 @@ class RelativeDifferencePrior: public void compute_gradient(DiscretisedDensity<3,elemT>& prior_gradient, const DiscretisedDensity<3,elemT> ¤t_image_estimate); - //! compute Hessian virtual Succeeded compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, const BasicCoordinate<3,int>& coords, const DiscretisedDensity<3,elemT> ¤t_image_estimate) const; From 59aed0b6c23d73011c56328affc5c643b01353bb Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Mon, 19 Jul 2021 13:19:02 -0700 Subject: [PATCH 24/27] [ci skip] Rework and rename of second derivative methods for priors --- .../stir/recon_buildblock/LogcoshPrior.h | 9 +++--- .../stir/recon_buildblock/QuadraticPrior.h | 9 +++--- .../RelativeDifferencePrior.h | 10 ++++--- src/recon_buildblock/LogcoshPrior.cxx | 30 +++++++++---------- src/recon_buildblock/QuadraticPrior.cxx | 28 ++++++++--------- .../RelativeDifferencePrior.cxx | 28 ++++++++--------- 6 files changed, 57 insertions(+), 57 deletions(-) diff --git a/src/include/stir/recon_buildblock/LogcoshPrior.h b/src/include/stir/recon_buildblock/LogcoshPrior.h index 5f4eb7ae7a..63c09153e7 100644 --- a/src/include/stir/recon_buildblock/LogcoshPrior.h +++ b/src/include/stir/recon_buildblock/LogcoshPrior.h @@ -228,16 +228,15 @@ class LogcoshPrior: public //! The second partial derivatives of the LogCosh Prior /*! - Diagonal refers to the second derivative w.r.t. x_j only (i.e. diagonal of the Hessian) - Off-diagonal refers to the second derivative w.r.t. x_j and x_k (i.e. off-diagonal of the Hessian) - For LogCosh, the off diagonal is the negative of the diagonal. + 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 */ //@{ - float diagonal_second_derivative(const float x_j, const float x_k) const; - float off_diagonal_second_derivative(const float x_j, const float x_k) const; + elemT derivative_20(const elemT x_j, const elemT x_k) const; + elemT derivative_11(const elemT x_j, const elemT x_k) const; //@} }; diff --git a/src/include/stir/recon_buildblock/QuadraticPrior.h b/src/include/stir/recon_buildblock/QuadraticPrior.h index 4d2f9560f8..09b2bef9d6 100644 --- a/src/include/stir/recon_buildblock/QuadraticPrior.h +++ b/src/include/stir/recon_buildblock/QuadraticPrior.h @@ -184,16 +184,15 @@ class QuadraticPrior: public //! The second partial derivatives of the Quadratic Prior /*! - Diagonal refers to the second derivative w.r.t. x_j only (i.e. diagonal of the Hessian) - Off-diagonal refers to the second derivative w.r.t. x_j and x_k (i.e. off-diagonal of the Hessian) - For the Quadratic Prior, the off diagonal is the negative of the diagonal. + 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 */ //@{ - float diagonal_second_derivative(const float x_j, const float x_k) const; - float off_diagonal_second_derivative(const float x_j, const float x_k) const; + elemT derivative_20(const elemT x_j, const elemT x_k) const; + elemT derivative_11(const elemT x_j, const elemT x_k) const; //@} }; diff --git a/src/include/stir/recon_buildblock/RelativeDifferencePrior.h b/src/include/stir/recon_buildblock/RelativeDifferencePrior.h index 537cd94ce5..db6185c81e 100644 --- a/src/include/stir/recon_buildblock/RelativeDifferencePrior.h +++ b/src/include/stir/recon_buildblock/RelativeDifferencePrior.h @@ -206,16 +206,18 @@ class RelativeDifferencePrior: public //! The second partial derivatives of the Relative Difference Prior /*! - Diagonal refers to the second derivative w.r.t. x_j only (i.e. diagonal of the Hessian) - Off-diagonal refers to the second derivative w.r.t. x_j and x_k (i.e. off-diagonal of the Hessian) + 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 */ //@{ - float diagonal_second_derivative(const float x_j, const float x_k) const; - float off_diagonal_second_derivative(const float x_j, const float x_k) const; + elemT derivative_20(const elemT x_j, const elemT x_k) const; + elemT derivative_11(const elemT x_j, const elemT x_k) const; //@} }; diff --git a/src/recon_buildblock/LogcoshPrior.cxx b/src/recon_buildblock/LogcoshPrior.cxx index 1fee3ab7a2..f588af7c94 100644 --- a/src/recon_buildblock/LogcoshPrior.cxx +++ b/src/recon_buildblock/LogcoshPrior.cxx @@ -441,8 +441,8 @@ compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, for (int ddx=min_dx;ddx<=max_dx;++ddx) { elemT diagonal_current = weights[ddz][ddy][ddx] * - diagonal_second_derivative(current_image_estimate[z][y][x], - current_image_estimate[z+ddz][y+ddy][x+ddx]); + derivative_20(current_image_estimate[z][y][x], + current_image_estimate[z + ddz][y + ddy][x + ddx]); if (do_kappa) diagonal_current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+ddz][y+ddy][x+ddx]; current += diagonal_current; @@ -451,8 +451,8 @@ compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, else { // The j != k vases (off-diagonal Hessian elements) - current = weights[dz][dy][dx] * off_diagonal_second_derivative(current_image_estimate[z][y][x], - current_image_estimate[z+dz][y+dy][x+dx]); + current = weights[dz][dy][dx] * derivative_11(current_image_estimate[z][y][x], + current_image_estimate[z + dz][y + dy][x + dx]); if (do_kappa) current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; } @@ -592,13 +592,13 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, elemT current = weights[dz][dy][dx]; if (dz == dy == dz == 0) { // The j == k case - current *= diagonal_second_derivative(current_estimate[z][y][x], - current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x]; + current *= derivative_20(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x]; } else { - current *= (diagonal_second_derivative(current_estimate[z][y][x], - current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x] + - off_diagonal_second_derivative(current_estimate[z][y][x], - current_estimate[z + dz][y + dy][x + dx]) * + current *= (derivative_20(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x] + + derivative_11(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z + dz][y + dy][x + dx]); } @@ -616,19 +616,19 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, } template -float +elemT LogcoshPrior:: -diagonal_second_derivative(const float x_j, const float x_k) const +derivative_20(const elemT x_j, const elemT x_k) const { return square((1/ cosh((x_j - x_k) * scalar))); } template -float +elemT LogcoshPrior:: -off_diagonal_second_derivative(const float x_j, const float x_k) const +derivative_11(const elemT x_j, const elemT x_k) const { - return - diagonal_second_derivative(x_j, x_k); + return -derivative_20(x_j, x_k); } # ifdef _MSC_VER diff --git a/src/recon_buildblock/QuadraticPrior.cxx b/src/recon_buildblock/QuadraticPrior.cxx index 8563c2e1f0..a732146c6d 100644 --- a/src/recon_buildblock/QuadraticPrior.cxx +++ b/src/recon_buildblock/QuadraticPrior.cxx @@ -495,8 +495,8 @@ compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, for (int ddx=min_dx;ddx<=max_dx;++ddx) { elemT diagonal_current = weights[ddz][ddy][ddx] * - diagonal_second_derivative(current_image_estimate[z][y][x], - current_image_estimate[z+ddz][y+ddy][x+ddx]); + derivative_20(current_image_estimate[z][y][x], + current_image_estimate[z + ddz][y + ddy][x + ddx]); if (do_kappa) diagonal_current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+ddz][y+ddy][x+ddx]; current += diagonal_current; @@ -505,8 +505,8 @@ compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, else { // The j != k vases (off-diagonal Hessian elements) - current = weights[dz][dy][dx] * off_diagonal_second_derivative(current_image_estimate[z][y][x], - current_image_estimate[z+dz][y+dy][x+dx]); + current = weights[dz][dy][dx] * derivative_11(current_image_estimate[z][y][x], + current_image_estimate[z + dz][y + dy][x + dx]); if (do_kappa) current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; } @@ -741,13 +741,13 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, elemT current = weights[dz][dy][dx]; if (dz == dy == dz == 0) { // The j == k case - current *= diagonal_second_derivative(current_estimate[z][y][x], - current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x]; + current *= derivative_20(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x]; } else { - current *= (diagonal_second_derivative(current_estimate[z][y][x], - current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x] + - off_diagonal_second_derivative(current_estimate[z][y][x], - current_estimate[z + dz][y + dy][x + dx]) * + current *= (derivative_20(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x] + + derivative_11(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z + dz][y + dy][x + dx]); } @@ -765,15 +765,15 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, } template -float +elemT QuadraticPrior:: -diagonal_second_derivative(const float x_j, const float x_k) const +derivative_20(const elemT x_j, const elemT x_k) const { return 1.0;} template -float +elemT QuadraticPrior:: -off_diagonal_second_derivative(const float x_j, const float x_k) const +derivative_11(const elemT x_j, const elemT x_k) const { return -1.0;} # ifdef _MSC_VER diff --git a/src/recon_buildblock/RelativeDifferencePrior.cxx b/src/recon_buildblock/RelativeDifferencePrior.cxx index c9d50211b5..3fcd8d6763 100644 --- a/src/recon_buildblock/RelativeDifferencePrior.cxx +++ b/src/recon_buildblock/RelativeDifferencePrior.cxx @@ -489,8 +489,8 @@ compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, for (int ddx=min_dx;ddx<=max_dx;++ddx) { elemT diagonal_current = weights[ddz][ddy][ddx] * - diagonal_second_derivative(current_image_estimate[z][y][x], - current_image_estimate[z+ddz][y+ddy][x+ddx]); + derivative_20(current_image_estimate[z][y][x], + current_image_estimate[z + ddz][y + ddy][x + ddx]); if (do_kappa) diagonal_current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+ddz][y+ddy][x+ddx]; current += diagonal_current; @@ -499,8 +499,8 @@ compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, else { // The j != k vases (off-diagonal Hessian elements) - current = weights[dz][dy][dx] * off_diagonal_second_derivative(current_image_estimate[z][y][x], - current_image_estimate[z+dz][y+dy][x+dx]); + current = weights[dz][dy][dx] * derivative_11(current_image_estimate[z][y][x], + current_image_estimate[z + dz][y + dy][x + dx]); if (do_kappa) current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; } @@ -588,13 +588,13 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, elemT current = weights[dz][dy][dx]; if (dz == dy == dz == 0) { // The j == k case - current *= diagonal_second_derivative(current_estimate[z][y][x], - current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x]; + current *= derivative_20(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x]; } else { - current *= (diagonal_second_derivative(current_estimate[z][y][x], - current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x] + - off_diagonal_second_derivative(current_estimate[z][y][x], - current_estimate[z + dz][y + dy][x + dx]) * + current *= (derivative_20(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z][y][x] + + derivative_11(current_estimate[z][y][x], + current_estimate[z + dz][y + dy][x + dx]) * input[z + dz][y + dy][x + dx]); } @@ -612,9 +612,9 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, } template -float +elemT RelativeDifferencePrior:: -diagonal_second_derivative(const float x_j, const float x_k) const +derivative_20(const elemT x_j, const elemT x_k) const { if (x_j > 0.0 || x_k > 0.0 || this->epsilon > 0.0) return 2 * pow(2 * x_k + this->epsilon, 2) / @@ -624,9 +624,9 @@ diagonal_second_derivative(const float x_j, const float x_k) const } template -float +elemT RelativeDifferencePrior:: -off_diagonal_second_derivative(const float x_j, const float x_k) const +derivative_11(const elemT x_j, const elemT x_k) const { if (x_j > 0.0 || x_k > 0.0 || this->epsilon > 0.0) return - 2 * (2 * x_j + this->epsilon)*(2 * x_k + this->epsilon) / From 15aebe74d72d14869efbde7378d46cdd7f3b4c44 Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Tue, 20 Jul 2021 10:53:42 -0700 Subject: [PATCH 25/27] All prior Hessian methods return void values This removes two checks on Succeeded::no in GeneralisedObjectiveFunction --- .../stir/recon_buildblock/GeneralisedPrior.h | 6 +++--- src/include/stir/recon_buildblock/LogcoshPrior.h | 4 ++-- .../stir/recon_buildblock/QuadraticPrior.h | 6 +++--- .../recon_buildblock/RelativeDifferencePrior.h | 6 +++--- .../GeneralisedObjectiveFunction.cxx | 9 ++------- src/recon_buildblock/GeneralisedPrior.cxx | 9 +++------ src/recon_buildblock/LogcoshPrior.cxx | 11 ++++------- src/recon_buildblock/QuadraticPrior.cxx | 16 ++++++---------- src/recon_buildblock/RelativeDifferencePrior.cxx | 14 +++++--------- 9 files changed, 31 insertions(+), 50 deletions(-) diff --git a/src/include/stir/recon_buildblock/GeneralisedPrior.h b/src/include/stir/recon_buildblock/GeneralisedPrior.h index ad35c8c166..a68bb04539 100644 --- a/src/include/stir/recon_buildblock/GeneralisedPrior.h +++ b/src/include/stir/recon_buildblock/GeneralisedPrior.h @@ -68,7 +68,7 @@ class GeneralisedPrior: 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 Succeeded + virtual void compute_Hessian(DataT& prior_Hessian_for_single_densel, const BasicCoordinate<3,int>& coords, const DataT& current_image_estimate) const; @@ -80,7 +80,7 @@ class GeneralisedPrior: 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; @@ -89,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; diff --git a/src/include/stir/recon_buildblock/LogcoshPrior.h b/src/include/stir/recon_buildblock/LogcoshPrior.h index 63c09153e7..fb9ed5aa4b 100644 --- a/src/include/stir/recon_buildblock/LogcoshPrior.h +++ b/src/include/stir/recon_buildblock/LogcoshPrior.h @@ -126,13 +126,13 @@ class LogcoshPrior: public void parabolic_surrogate_curvature(DiscretisedDensity<3,elemT>& parabolic_surrogate_curvature, const DiscretisedDensity<3,elemT> ¤t_image_estimate); - virtual Succeeded + virtual void compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, const BasicCoordinate<3,int>& coords, const DiscretisedDensity<3,elemT> ¤t_image_estimate) const; //! Compute the multiplication of the hessian of the prior (at \c current_estimate) and the \c 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; diff --git a/src/include/stir/recon_buildblock/QuadraticPrior.h b/src/include/stir/recon_buildblock/QuadraticPrior.h index 09b2bef9d6..1a328fa802 100644 --- a/src/include/stir/recon_buildblock/QuadraticPrior.h +++ b/src/include/stir/recon_buildblock/QuadraticPrior.h @@ -118,20 +118,20 @@ class QuadraticPrior: public void parabolic_surrogate_curvature(DiscretisedDensity<3,elemT>& parabolic_surrogate_curvature, const DiscretisedDensity<3,elemT> ¤t_image_estimate); - virtual Succeeded + virtual void compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, const BasicCoordinate<3,int>& coords, const DiscretisedDensity<3,elemT> ¤t_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; diff --git a/src/include/stir/recon_buildblock/RelativeDifferencePrior.h b/src/include/stir/recon_buildblock/RelativeDifferencePrior.h index db6185c81e..6df24dba0f 100644 --- a/src/include/stir/recon_buildblock/RelativeDifferencePrior.h +++ b/src/include/stir/recon_buildblock/RelativeDifferencePrior.h @@ -125,16 +125,16 @@ class RelativeDifferencePrior: public void compute_gradient(DiscretisedDensity<3,elemT>& prior_gradient, const DiscretisedDensity<3,elemT> ¤t_image_estimate); - virtual Succeeded compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, + virtual void compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, const BasicCoordinate<3,int>& coords, const DiscretisedDensity<3,elemT> ¤t_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 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; diff --git a/src/recon_buildblock/GeneralisedObjectiveFunction.cxx b/src/recon_buildblock/GeneralisedObjectiveFunction.cxx index 41b842c2ed..f3490fa5f7 100644 --- a/src/recon_buildblock/GeneralisedObjectiveFunction.cxx +++ b/src/recon_buildblock/GeneralisedObjectiveFunction.cxx @@ -263,10 +263,7 @@ add_multiplication_with_approximate_sub_Hessian(TargetT& output, { // TODO used boost:scoped_ptr shared_ptr 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 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(); diff --git a/src/recon_buildblock/GeneralisedPrior.cxx b/src/recon_buildblock/GeneralisedPrior.cxx index 58c8b527f1..d4ce84b42f 100644 --- a/src/recon_buildblock/GeneralisedPrior.cxx +++ b/src/recon_buildblock/GeneralisedPrior.cxx @@ -52,7 +52,7 @@ set_up(shared_ptr const&) } template -Succeeded +void GeneralisedPrior:: compute_Hessian(TargetT& output, const BasicCoordinate<3,int>& coords, @@ -62,22 +62,20 @@ compute_Hessian(TargetT& output, error("GeneralisedPrior:\n compute_Hessian implementation is not overloaded by your convex prior."); else error("GeneralisedPrior:\n compute_Hessian is not implemented because the prior is not convex."); - return Succeeded::no; } template -Succeeded +void GeneralisedPrior:: 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 -Succeeded +void GeneralisedPrior:: accumulate_Hessian_times_input(TargetT& output, const TargetT& current_estimate, @@ -85,7 +83,6 @@ accumulate_Hessian_times_input(TargetT& output, { error("GeneralisedPrior:\n" "accumulate_Hessian_times_input implementation is not overloaded by your prior."); - return Succeeded::no; } template diff --git a/src/recon_buildblock/LogcoshPrior.cxx b/src/recon_buildblock/LogcoshPrior.cxx index f588af7c94..5324498430 100644 --- a/src/recon_buildblock/LogcoshPrior.cxx +++ b/src/recon_buildblock/LogcoshPrior.cxx @@ -384,7 +384,7 @@ compute_gradient(DiscretisedDensity<3,elemT>& prior_gradient, } template -Succeeded +void LogcoshPrior:: compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, const BasicCoordinate<3,int>& coords, @@ -394,7 +394,7 @@ compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, prior_Hessian_for_single_densel.fill(0); if (this->penalisation_factor==0) { - return Succeeded::yes; + return; } this->check(current_image_estimate); @@ -458,8 +458,6 @@ compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, } prior_Hessian_for_single_densel_cast[z+dz][y+dy][x+dx] = + current*this->penalisation_factor; } - - return Succeeded::yes; } @@ -533,7 +531,7 @@ LogcoshPrior::parabolic_surrogate_curvature(DiscretisedDensity<3,elemT>& } template -Succeeded +void LogcoshPrior:: accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, const DiscretisedDensity<3,elemT>& current_estimate, @@ -545,7 +543,7 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, assert( output.has_same_characteristics(input)); if (this->penalisation_factor==0) { - return Succeeded::yes; + return; } DiscretisedDensityOnCartesianGrid<3,elemT>& output_cast = @@ -612,7 +610,6 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, } } } - return Succeeded::yes; } template diff --git a/src/recon_buildblock/QuadraticPrior.cxx b/src/recon_buildblock/QuadraticPrior.cxx index a732146c6d..7e1260ff41 100644 --- a/src/recon_buildblock/QuadraticPrior.cxx +++ b/src/recon_buildblock/QuadraticPrior.cxx @@ -437,7 +437,7 @@ compute_gradient(DiscretisedDensity<3,elemT>& prior_gradient, } template -Succeeded +void QuadraticPrior:: compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, const BasicCoordinate<3,int>& coords, @@ -447,7 +447,7 @@ compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, prior_Hessian_for_single_densel.fill(0); if (this->penalisation_factor==0) { - return Succeeded::yes; + return; } this->check(current_image_estimate); @@ -512,8 +512,6 @@ compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, } prior_Hessian_for_single_densel_cast[z+dz][y+dy][x+dx] = + current*this->penalisation_factor; } - - return Succeeded::yes; } template @@ -598,7 +596,7 @@ QuadraticPrior::parabolic_surrogate_curvature(DiscretisedDensity<3,elemT> } template -Succeeded +void QuadraticPrior:: add_multiplication_with_approximate_Hessian(DiscretisedDensity<3,elemT>& output, const DiscretisedDensity<3,elemT>& input) const @@ -609,7 +607,7 @@ add_multiplication_with_approximate_Hessian(DiscretisedDensity<3,elemT>& output, assert( output.has_same_characteristics(input)); if (this->penalisation_factor==0) { - return Succeeded::yes; + return; } this->check(input); @@ -668,11 +666,10 @@ add_multiplication_with_approximate_Hessian(DiscretisedDensity<3,elemT>& output, } } } - return Succeeded::yes; } template -Succeeded +void QuadraticPrior:: accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, const DiscretisedDensity<3,elemT>& current_estimate, @@ -684,7 +681,7 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, assert( output.has_same_characteristics(input)); if (this->penalisation_factor==0) { - return Succeeded::yes; + return; } this->check(input); @@ -761,7 +758,6 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, } } } - return Succeeded::yes; } template diff --git a/src/recon_buildblock/RelativeDifferencePrior.cxx b/src/recon_buildblock/RelativeDifferencePrior.cxx index 3fcd8d6763..69c443cc42 100644 --- a/src/recon_buildblock/RelativeDifferencePrior.cxx +++ b/src/recon_buildblock/RelativeDifferencePrior.cxx @@ -431,7 +431,7 @@ compute_gradient(DiscretisedDensity<3,elemT>& prior_gradient, } template -Succeeded +void RelativeDifferencePrior:: compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, const BasicCoordinate<3,int>& coords, @@ -441,7 +441,7 @@ compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, prior_Hessian_for_single_densel.fill(0); if (this->penalisation_factor==0) { - return Succeeded::yes; + return; } this->check(current_image_estimate); @@ -506,22 +506,19 @@ compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel, } prior_Hessian_for_single_densel_cast[z+dz][y+dy][x+dx] = + current*this->penalisation_factor; } - - return Succeeded::yes; } template -Succeeded +void RelativeDifferencePrior:: add_multiplication_with_approximate_Hessian(DiscretisedDensity<3,elemT>& output, const DiscretisedDensity<3,elemT>& input) const { error("add_multiplication_with_approximate_Hessian() is not implemented in Relative Difference Prior."); - return Succeeded::no; } template -Succeeded +void RelativeDifferencePrior:: accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, const DiscretisedDensity<3,elemT>& current_estimate, @@ -533,7 +530,7 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, assert( output.has_same_characteristics(input)); if (this->penalisation_factor==0) { - return Succeeded::yes; + return; } DiscretisedDensityOnCartesianGrid<3,elemT>& output_cast = @@ -608,7 +605,6 @@ accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output, } } } - return Succeeded::yes; } template From 5fc8fea8d0de5eef2d70a5b9944b78eede9e9b94 Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Thu, 30 Sep 2021 16:09:25 -0700 Subject: [PATCH 26/27] Initial commit adding a new parent class --- .../recon_buildblock/GeneralisedConvexPrior.h | 75 +++++++++++++++++++ .../stir/recon_buildblock/QuadraticPrior.h | 4 +- src/recon_buildblock/QuadraticPrior.cxx | 1 + 3 files changed, 79 insertions(+), 1 deletion(-) create mode 100644 src/include/stir/recon_buildblock/GeneralisedConvexPrior.h diff --git a/src/include/stir/recon_buildblock/GeneralisedConvexPrior.h b/src/include/stir/recon_buildblock/GeneralisedConvexPrior.h new file mode 100644 index 0000000000..43509fe872 --- /dev/null +++ b/src/include/stir/recon_buildblock/GeneralisedConvexPrior.h @@ -0,0 +1,75 @@ +// +// +/* + Copyright (C) 2000- 2007, Hammersmith Imanet Ltd + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ +/*! + \file + \ingroup priors + \brief Declaration of class stir::GeneralisedPrior + + \author Kris Thielemans + \author Sanida Mustafovic + +*/ + +#ifndef __stir_recon_buildblock_GeneralisedConvexPrior_H__ +#define __stir_recon_buildblock_GeneralisedConvexPrior_H__ + + +#include "stir/RegisteredObject.h" +#include "stir/ParsingObject.h" + +START_NAMESPACE_STIR + +class Succeeded; + +/*! + \ingroup priors + \brief + A base class for 'generalised' priors, i.e. priors for which at least + a 'gradient' is defined. + + This class exists to accomodate FilterRootPrior. Otherwise we could + just live with Prior as a base class. +*/ +template + +class GeneralisedConvexPrior + +{ +public: + +// GeneralisedConvexPrior(); + +// inline GeneralisedConvexPrior() +// { +// int x = 1; +// }; + +public: + //! compute the value of the function + /*! For derived classes where this doesn't make sense, it's recommended to return 0. + */ + virtual void + compute_Hessian(DataT& prior_Hessian_for_single_densel, + const BasicCoordinate<3,int>& coords, + const DataT& current_image_estimate) const = 0; + + + + + int my_value = 1; + +}; + +END_NAMESPACE_STIR + +//#include "stir/recon_buildblock/GeneralisedPrior.inl" + +#endif diff --git a/src/include/stir/recon_buildblock/QuadraticPrior.h b/src/include/stir/recon_buildblock/QuadraticPrior.h index 1a328fa802..7a61f2a86b 100644 --- a/src/include/stir/recon_buildblock/QuadraticPrior.h +++ b/src/include/stir/recon_buildblock/QuadraticPrior.h @@ -25,6 +25,7 @@ #include "stir/RegisteredParsingObject.h" #include "stir/recon_buildblock/PriorWithParabolicSurrogate.h" +#include "stir/recon_buildblock/GeneralisedConvexPrior.h" #include "stir/Array.h" #include "stir/DiscretisedDensity.h" #include "stir/shared_ptr.h" @@ -80,7 +81,8 @@ class QuadraticPrior: public RegisteredParsingObject< QuadraticPrior, GeneralisedPrior >, PriorWithParabolicSurrogate > - > + >, + public GeneralisedConvexPrior > { private: typedef diff --git a/src/recon_buildblock/QuadraticPrior.cxx b/src/recon_buildblock/QuadraticPrior.cxx index 7e1260ff41..58fd7813e5 100644 --- a/src/recon_buildblock/QuadraticPrior.cxx +++ b/src/recon_buildblock/QuadraticPrior.cxx @@ -19,6 +19,7 @@ */ #include "stir/recon_buildblock/QuadraticPrior.h" +#include "stir/recon_buildblock/GeneralisedConvexPrior.h" #include "stir/Succeeded.h" #include "stir/DiscretisedDensityOnCartesianGrid.h" #include "stir/IndexRange3D.h" From 30b3ee6dc31692441564c25c8962b0d33556944a Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Fri, 1 Oct 2021 17:58:36 -0700 Subject: [PATCH 27/27] Buildable version of QP using GeneralisedConvexPrior as an intermediate --- .../recon_buildblock/GeneralisedConvexPrior.h | 42 +++---- .../PriorWithParabolicSurrogate.h | 2 +- .../stir/recon_buildblock/QuadraticPrior.h | 12 +- src/iterative/OSSPS/OSSPSReconstruction.cxx | 7 +- src/recon_buildblock/CMakeLists.txt | 1 + .../GeneralisedConvexPrior.cxx | 114 ++++++++++++++++++ 6 files changed, 145 insertions(+), 33 deletions(-) create mode 100644 src/recon_buildblock/GeneralisedConvexPrior.cxx diff --git a/src/include/stir/recon_buildblock/GeneralisedConvexPrior.h b/src/include/stir/recon_buildblock/GeneralisedConvexPrior.h index 43509fe872..83e9a4d7ca 100644 --- a/src/include/stir/recon_buildblock/GeneralisedConvexPrior.h +++ b/src/include/stir/recon_buildblock/GeneralisedConvexPrior.h @@ -11,10 +11,9 @@ /*! \file \ingroup priors - \brief Declaration of class stir::GeneralisedPrior + \brief Declaration of class stir::GeneralisedConvexPrior - \author Kris Thielemans - \author Sanida Mustafovic + \author Robert Twyman */ @@ -22,8 +21,7 @@ #define __stir_recon_buildblock_GeneralisedConvexPrior_H__ -#include "stir/RegisteredObject.h" -#include "stir/ParsingObject.h" +#include "stir/recon_buildblock/GeneralisedPrior.h" START_NAMESPACE_STIR @@ -32,37 +30,37 @@ class Succeeded; /*! \ingroup priors \brief - A base class for 'generalised' priors, i.e. priors for which at least - a 'gradient' is defined. - - This class exists to accomodate FilterRootPrior. Otherwise we could - just live with Prior as a base class. + Make a brief */ template -class GeneralisedConvexPrior +class GeneralisedConvexPrior: + virtual public GeneralisedPrior { -public: - -// GeneralisedConvexPrior(); - -// inline GeneralisedConvexPrior() -// { -// int x = 1; -// }; +private: + typedef GeneralisedPrior base_type; public: - //! compute the value of the function - /*! For derived classes where this doesn't make sense, it's recommended to return 0. - */ + //! This computes a single row of the Hessian + /*! 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 = 0; +// float penalisation_factor; + + +// virtual void +// actual_compute_Hessian(DataT& prior_Hessian_for_single_densel, +// const BasicCoordinate<3,int>& coords, +// const DataT& current_image_estimate) const; int my_value = 1; diff --git a/src/include/stir/recon_buildblock/PriorWithParabolicSurrogate.h b/src/include/stir/recon_buildblock/PriorWithParabolicSurrogate.h index 57dea783d6..df6e0c8d9c 100644 --- a/src/include/stir/recon_buildblock/PriorWithParabolicSurrogate.h +++ b/src/include/stir/recon_buildblock/PriorWithParabolicSurrogate.h @@ -38,7 +38,7 @@ START_NAMESPACE_STIR */ template class PriorWithParabolicSurrogate: - public GeneralisedPrior + virtual public GeneralisedPrior { public: diff --git a/src/include/stir/recon_buildblock/QuadraticPrior.h b/src/include/stir/recon_buildblock/QuadraticPrior.h index 7a61f2a86b..281edcbb47 100644 --- a/src/include/stir/recon_buildblock/QuadraticPrior.h +++ b/src/include/stir/recon_buildblock/QuadraticPrior.h @@ -78,16 +78,16 @@ START_NAMESPACE_STIR */ template class QuadraticPrior: public - RegisteredParsingObject< QuadraticPrior, - GeneralisedPrior >, - PriorWithParabolicSurrogate > - >, - public GeneralisedConvexPrior > + RegisteredParsingObject< + QuadraticPrior, + GeneralisedPrior >, + GeneralisedConvexPrior > >, + PriorWithParabolicSurrogate > { private: typedef RegisteredParsingObject< QuadraticPrior, - GeneralisedPrior >, + GeneralisedConvexPrior >, PriorWithParabolicSurrogate > > base_type; diff --git a/src/iterative/OSSPS/OSSPSReconstruction.cxx b/src/iterative/OSSPS/OSSPSReconstruction.cxx index e73478ff91..eed4624e50 100644 --- a/src/iterative/OSSPS/OSSPSReconstruction.cxx +++ b/src/iterative/OSSPS/OSSPSReconstruction.cxx @@ -327,10 +327,10 @@ update_estimate(TargetT ¤t_image_estimate) // For the quadratic prior, this is independent of the image (only on kappa's) // And of course, it's also independent when there is no prior // TODO by default, this should be off probably (to save time). + auto* parabolic_surrogate_prior = dynamic_cast*>(this->get_prior_ptr()); const bool recompute_penalty_term_in_denominator = !this->objective_function_sptr->prior_is_zero() && - static_cast const&>(*this->get_prior_ptr()). - parabolic_surrogate_curvature_depends_on_argument(); + parabolic_surrogate_prior->parabolic_surrogate_curvature_depends_on_argument(); #ifndef PARALLEL //CPUTimer subset_timer; //subset_timer.start(); @@ -366,8 +366,7 @@ update_estimate(TargetT ¤t_image_estimate) // avoid work (or crash) when penalty is 0 if (!this->objective_function_sptr->prior_is_zero()) { - static_cast&>(*get_prior_ptr()). - parabolic_surrogate_curvature(*work_image_ptr, current_image_estimate); + parabolic_surrogate_prior->parabolic_surrogate_curvature(*work_image_ptr, current_image_estimate); //*work_image_ptr *= 2; //*work_image_ptr += *precomputed_denominator_ptr ; std::transform(work_image_ptr->begin_all(), work_image_ptr->end_all(), diff --git a/src/recon_buildblock/CMakeLists.txt b/src/recon_buildblock/CMakeLists.txt index 6ed94fb7b1..d1dbabf7ce 100644 --- a/src/recon_buildblock/CMakeLists.txt +++ b/src/recon_buildblock/CMakeLists.txt @@ -58,6 +58,7 @@ set(${dir_LIB_SOURCES} BinNormalisationFromAttenuationImage.cxx BinNormalisationSPECT.cxx GeneralisedPrior.cxx + GeneralisedConvexPrior.cxx ProjDataRebinning.cxx FourierRebinning.cxx PLSPrior.cxx diff --git a/src/recon_buildblock/GeneralisedConvexPrior.cxx b/src/recon_buildblock/GeneralisedConvexPrior.cxx new file mode 100644 index 0000000000..9ef485bfc2 --- /dev/null +++ b/src/recon_buildblock/GeneralisedConvexPrior.cxx @@ -0,0 +1,114 @@ +// +// +/* + Copyright (C) 2002- 2009, Hammersmith Imanet Ltd + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ +/*! + \file + \ingroup priors + \brief implementation of the stir::GeneralisedConvexPrior + + \author Robert Twyman +*/ + +#include "stir/recon_buildblock/GeneralisedConvexPrior.h" +//#include "stir/DiscretisedDensity.h" +//#include "stir/Succeeded.h" +#include "stir/modelling/ParametricDiscretisedDensity.h" +//#include "stir/modelling/KineticParameters.h" + +START_NAMESPACE_STIR + +//template +//void +//GeneralisedConvexPrior:: +//actual_compute_Hessian(elemT& prior_Hessian_for_single_densel, +// const BasicCoordinate<3,int>& coords, +// const elemT& current_image_estimate) const +//{ +// assert( prior_Hessian_for_single_densel.has_same_characteristics(current_image_estimate)); +// prior_Hessian_for_single_densel.fill(0); +// if (this->penalisation_factor==0) +// { +// return; +// } +// +// this->check(current_image_estimate); +// +// const DiscretisedDensityOnCartesianGrid<3,elemT>& current_image_cast = +// dynamic_cast< const DiscretisedDensityOnCartesianGrid<3,elemT> &>(current_image_estimate); +// +// DiscretisedDensityOnCartesianGrid<3,elemT>& prior_Hessian_for_single_densel_cast = +// dynamic_cast &>(prior_Hessian_for_single_densel); +// +// if (weights.get_length() ==0) +// { +// compute_weights(weights, current_image_cast.get_grid_spacing(), this->only_2D); +// } +// +// +// 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]; +// const int min_dz = max(weights.get_min_index(), prior_Hessian_for_single_densel.get_min_index()-z); +// const int max_dz = min(weights.get_max_index(), prior_Hessian_for_single_densel.get_max_index()-z); +// +// const int min_dy = max(weights[0].get_min_index(), prior_Hessian_for_single_densel[z].get_min_index()-y); +// const int max_dy = min(weights[0].get_max_index(), prior_Hessian_for_single_densel[z].get_max_index()-y); +// +// const int min_dx = max(weights[0][0].get_min_index(), prior_Hessian_for_single_densel[z][y].get_min_index()-x); +// const int max_dx = min(weights[0][0].get_max_index(), prior_Hessian_for_single_densel[z][y].get_max_index()-x); +// +// elemT diagonal = 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 current = 0.0; +// if (dz == 0 && dy == 0 && dx == 0) +// { +// // The j == k case (diagonal Hessian element), which is a sum over the neighbourhood. +// for (int ddz=min_dz;ddz<=max_dz;++ddz) +// for (int ddy=min_dy;ddy<=max_dy;++ddy) +// for (int ddx=min_dx;ddx<=max_dx;++ddx) +// { +// elemT diagonal_current = weights[ddz][ddy][ddx] * +// derivative_20(current_image_estimate[z][y][x], +// current_image_estimate[z + ddz][y + ddy][x + ddx]); +// if (do_kappa) +// diagonal_current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+ddz][y+ddy][x+ddx]; +// current += diagonal_current; +// } +// } +// else +// { +// // The j != k vases (off-diagonal Hessian elements) +// current = weights[dz][dy][dx] * derivative_11(current_image_estimate[z][y][x], +// current_image_estimate[z + dz][y + dy][x + dx]); +// if (do_kappa) +// current *= (*kappa_ptr)[z][y][x] * (*kappa_ptr)[z+dz][y+dy][x+dx]; +// } +// prior_Hessian_for_single_densel_cast[z+dz][y+dy][x+dx] = + current*this->penalisation_factor; +// } +//} + + +# ifdef _MSC_VER + // prevent warning message on instantiation of abstract class +# pragma warning(disable:4661) +# endif + + template class GeneralisedConvexPrior >; + template class GeneralisedConvexPrior; + +END_NAMESPACE_STIR