From 9167f1e8bc8d65ebed088c169ec17fc027e0107d Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Sun, 26 Apr 2020 10:05:50 +0000 Subject: [PATCH 1/8] adding tests for priors currently failing... --- src/recon_test/CMakeLists.txt | 1 + src/recon_test/test_priors.cxx | 212 +++++++++++++++++++++++++++++++++ 2 files changed, 213 insertions(+) create mode 100644 src/recon_test/test_priors.cxx diff --git a/src/recon_test/CMakeLists.txt b/src/recon_test/CMakeLists.txt index 9865254cc..00b6b0ae7 100644 --- a/src/recon_test/CMakeLists.txt +++ b/src/recon_test/CMakeLists.txt @@ -35,6 +35,7 @@ set(${dir_INVOLVED_TEST_EXE_SOURCES} bcktest recontest test_data_processor_projectors + test_priors ) include(stir_test_exe_targets) diff --git a/src/recon_test/test_priors.cxx b/src/recon_test/test_priors.cxx new file mode 100644 index 000000000..74a53c6ac --- /dev/null +++ b/src/recon_test/test_priors.cxx @@ -0,0 +1,212 @@ +/* + Copyright (C) 2011, Hammersmith Imanet Ltd + Copyright (C) 2020 University College London + This file is part of STIR. + + This file is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. + + This file is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + See STIR/LICENSE.txt for details +*/ +/*! + + \file + \ingroup recon_test + + \brief Test program for stir::QuadraticPrior + + \par Usage + + 
+  test_priors [ density_filename ] 
+
+ where the argument is optional. See the class documentation for more info. + + \author Kris Thielemans +*/ + +#include "stir/VoxelsOnCartesianGrid.h" +#include "stir/recon_buildblock/QuadraticPrior.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" +#include "stir/info.h" +#include "stir/Succeeded.h" +#include "stir/num_threads.h" +#include +#include +#include +#include +#include +#include + +#include "stir/IO/OutputFileFormat.h" +START_NAMESPACE_STIR + + +/*! + \ingroup test + \brief Test class for QuadraticPrior + + This test compares the result of GeneralisedPrior::compute_gradient() + with a numerical gradient computed by using the + GeneralisedPrior::compute_value() function. + +*/ +class GeneralisedPriorTests : public RunTests +{ +public: + //! Constructor that can take some input data to run the test with + /*! This makes it possible to run the test with your own data. However, beware that + it is very easy to set up a very long computation. + + \todo it would be better to parse an objective function. That would allow us to set + all parameters from the command line. + */ + GeneralisedPriorTests(char const * const density_filename = 0); + typedef DiscretisedDensity<3,float> target_type; + void construct_input_data(shared_ptr& density_sptr); + + void run_tests(); +protected: + char const * density_filename; + shared_ptr > objective_function_sptr; + + //! run the test + /*! Note that this function is not specific to a particular prior */ + void run_tests_for_objective_function(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr); +}; + +GeneralisedPriorTests:: +GeneralisedPriorTests(char const * const density_filename) + : density_filename(density_filename) +{} + +void +GeneralisedPriorTests:: +run_tests_for_objective_function(const std::string& test_name, + GeneralisedPrior& objective_function, + shared_ptr target_sptr) +{ + std::cerr << "----- test " << test_name << '\n'; + if (!check(objective_function.set_up(target_sptr)==Succeeded::yes, "set-up of objective function")) + return; + + target_type& target(*target_sptr); + + shared_ptr gradient_sptr(target.get_empty_copy()); + shared_ptr gradient_2_sptr(target.get_empty_copy()); + info("Computing gradient",3); + objective_function.compute_gradient(*gradient_sptr, target); + this->set_tolerance(std::max(fabs(double(gradient_sptr->find_min())), double(gradient_sptr->find_max()))/1000); + info("Computing objective function at target",3); + 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(); + const float eps = 1e-2F; + 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; + 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; + } + } + 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); + } + +} + +void +GeneralisedPriorTests:: +construct_input_data(shared_ptr& density_sptr) +{ + if (this->density_filename == 0) + { + // construct a small image + + shared_ptr exam_info_sptr(new ExamInfo); + exam_info_sptr->imaging_modality = ImagingModality::PT; + CartesianCoordinate3D origin (0,0,0); + CartesianCoordinate3D voxel_size(2.F,3.F,3.F); + + density_sptr.reset(new VoxelsOnCartesianGrid(exam_info_sptr, + IndexRange<3>(make_coordinate(20,19,18)), + origin, voxel_size)); + // fill with random numbers between 0 and 1 + typedef boost::mt19937 base_generator_type; + // initialize by reproducible seed + static base_generator_type generator(boost::uint32_t(42)); + static boost::uniform_01 random01(generator); + for (target_type::full_iterator iter=density_sptr->begin_all(); iter!=density_sptr->end_all(); ++iter) + *iter = static_cast(random01()); + + } + else + { + shared_ptr aptr(read_from_file(this->density_filename)); + density_sptr = aptr; + } + + return; +} + +void +GeneralisedPriorTests:: +run_tests() +{ + shared_ptr density_sptr; + construct_input_data(density_sptr); + + std::cerr << "Tests for QuadraticPrior\n"; + { + QuadraticPrior objective_function(true, 3.F); + this->run_tests_for_objective_function("Quadratic_no_kappa", objective_function, density_sptr); + } + std::cerr << "Tests for PLSPrior\n"; + { + PLSPrior objective_function(true, 3.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); + } +} + +END_NAMESPACE_STIR + + +USING_NAMESPACE_STIR + +int main(int argc, char **argv) +{ + set_default_num_threads(); + + GeneralisedPriorTests tests(argc>1? argv[1] : 0); + tests.run_tests(); + return tests.main_return_value(); +} From 8f610c0e3ae05fc59cb6a0a232149725cf126ae4 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Wed, 23 Sep 2020 08:43:53 +0100 Subject: [PATCH 2/8] Update test_priors.cxx reduced value of epsilon --- src/recon_test/test_priors.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/recon_test/test_priors.cxx b/src/recon_test/test_priors.cxx index 74a53c6ac..89ea3cdef 100644 --- a/src/recon_test/test_priors.cxx +++ b/src/recon_test/test_priors.cxx @@ -114,7 +114,7 @@ run_tests_for_objective_function(const std::string& test_name, 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(); - const float eps = 1e-2F; + const float eps = 5e-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) From cd128426fa6750e8c8be63f1c2e6830d81a9a313 Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Fri, 18 Sep 2020 10:11:58 +0100 Subject: [PATCH 3/8] Add RDP test, add additional comments --- src/recon_test/test_priors.cxx | 24 +++++++++++++++++------- 1 file changed, 17 insertions(+), 7 deletions(-) diff --git a/src/recon_test/test_priors.cxx b/src/recon_test/test_priors.cxx index 89ea3cdef..8078a556c 100644 --- a/src/recon_test/test_priors.cxx +++ b/src/recon_test/test_priors.cxx @@ -34,6 +34,7 @@ #include "stir/VoxelsOnCartesianGrid.h" #include "stir/recon_buildblock/QuadraticPrior.h" +#include "stir/recon_buildblock/RelativeDifferencePrior.h" #include "stir/recon_buildblock/PLSPrior.h" #include "stir/RunTests.h" #include "stir/IO/read_from_file.h" @@ -102,31 +103,34 @@ 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; + // setup images target_type& target(*target_sptr); - shared_ptr gradient_sptr(target.get_empty_copy()); shared_ptr gradient_2_sptr(target.get_empty_copy()); + info("Computing gradient",3); objective_function.compute_gradient(*gradient_sptr, target); this->set_tolerance(std::max(fabs(double(gradient_sptr->find_min())), double(gradient_sptr->find_max()))/1000); + info("Computing objective function at target",3); 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(); - const float eps = 5e-3F; + + // setup perturbation response + const float eps = 1e-2F; 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; + *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"); + 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; @@ -147,7 +151,7 @@ construct_input_data(shared_ptr& density_sptr) { if (this->density_filename == 0) { - // construct a small image + // construct a small image with random voxel values between 0 and 1 shared_ptr exam_info_sptr(new ExamInfo); exam_info_sptr->imaging_modality = ImagingModality::PT; @@ -168,10 +172,10 @@ construct_input_data(shared_ptr& density_sptr) } else { + // load image from file shared_ptr aptr(read_from_file(this->density_filename)); density_sptr = aptr; } - return; } @@ -187,6 +191,12 @@ run_tests() QuadraticPrior objective_function(true, 3.F); this->run_tests_for_objective_function("Quadratic_no_kappa", objective_function, density_sptr); } + std::cerr << "Tests for Relative Difference Prior\n"; + { + // gamma and epsilon are default + RelativeDifferencePrior objective_function(true, 3.F, 2.F, 0.0001); + this->run_tests_for_objective_function("RDP_no_kappa", objective_function, density_sptr); + } std::cerr << "Tests for PLSPrior\n"; { PLSPrior objective_function(true, 3.F); From 46a3a06100b6d61aa0369cef2cfae0ce221dcf38 Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Wed, 23 Sep 2020 13:27:55 +0100 Subject: [PATCH 4/8] Change number of voxels in image, no need for that large --- src/recon_test/test_priors.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/recon_test/test_priors.cxx b/src/recon_test/test_priors.cxx index 8078a556c..5912e6983 100644 --- a/src/recon_test/test_priors.cxx +++ b/src/recon_test/test_priors.cxx @@ -159,7 +159,7 @@ construct_input_data(shared_ptr& density_sptr) CartesianCoordinate3D voxel_size(2.F,3.F,3.F); density_sptr.reset(new VoxelsOnCartesianGrid(exam_info_sptr, - IndexRange<3>(make_coordinate(20,19,18)), + IndexRange<3>(make_coordinate(10,9,8)), origin, voxel_size)); // fill with random numbers between 0 and 1 typedef boost::mt19937 base_generator_type; From f2b4f2293d972752a3864179f9371f0ccbd9cf7c Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Wed, 23 Sep 2020 13:28:48 +0100 Subject: [PATCH 5/8] Test priors in 3D --- src/recon_test/test_priors.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/recon_test/test_priors.cxx b/src/recon_test/test_priors.cxx index 5912e6983..2541d48c2 100644 --- a/src/recon_test/test_priors.cxx +++ b/src/recon_test/test_priors.cxx @@ -188,18 +188,18 @@ run_tests() std::cerr << "Tests for QuadraticPrior\n"; { - QuadraticPrior objective_function(true, 3.F); + 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"; { // gamma and epsilon are default - RelativeDifferencePrior objective_function(true, 3.F, 2.F, 0.0001); + RelativeDifferencePrior objective_function(false, 1.F, 2.F, 0.0001); this->run_tests_for_objective_function("RDP_no_kappa", objective_function, density_sptr); } std::cerr << "Tests for PLSPrior\n"; { - PLSPrior objective_function(true, 3.F); + 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); From 1b5d1314eef9ecd5eb0a9644ef0614dcf204c042 Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Wed, 23 Sep 2020 13:29:03 +0100 Subject: [PATCH 6/8] Small fixes --- 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 2541d48c2..fc6f0f40f 100644 --- a/src/recon_test/test_priors.cxx +++ b/src/recon_test/test_priors.cxx @@ -110,7 +110,7 @@ run_tests_for_objective_function(const std::string& test_name, info("Computing gradient",3); objective_function.compute_gradient(*gradient_sptr, target); - this->set_tolerance(std::max(fabs(double(gradient_sptr->find_min())), double(gradient_sptr->find_max()))/1000); + 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); const double value_at_target = objective_function.compute_value(target); @@ -119,7 +119,7 @@ run_tests_for_objective_function(const std::string& test_name, target_type::full_iterator gradient_2_iter=gradient_2_sptr->begin_all(); // setup perturbation response - const float eps = 1e-2F; + 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) From 8c362dabfc8e5b5a5188cab386d71a3c4ba14849 Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Thu, 24 Sep 2020 11:53:53 +0100 Subject: [PATCH 7/8] make test actually runs test_priors --- src/recon_test/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/recon_test/CMakeLists.txt b/src/recon_test/CMakeLists.txt index 00b6b0ae7..c8d6e777d 100644 --- a/src/recon_test/CMakeLists.txt +++ b/src/recon_test/CMakeLists.txt @@ -27,6 +27,7 @@ set(${dir_SIMPLE_TEST_EXE_SOURCES} test_FBP2D test_FBP3DRP test_OSMAPOSL + test_priors ) @@ -35,7 +36,6 @@ set(${dir_INVOLVED_TEST_EXE_SOURCES} bcktest recontest test_data_processor_projectors - test_priors ) include(stir_test_exe_targets) From 52f4fa82937906f7f47cd59c3a983eef9a635711 Mon Sep 17 00:00:00 2001 From: Robert Twyman Date: Fri, 18 Dec 2020 17:00:18 +0000 Subject: [PATCH 8/8] Add logcosh prior and disable PLS --- src/recon_test/test_priors.cxx | 32 ++++++++++++++++++++------------ 1 file changed, 20 insertions(+), 12 deletions(-) diff --git a/src/recon_test/test_priors.cxx b/src/recon_test/test_priors.cxx index fc6f0f40f..0abce86ea 100644 --- a/src/recon_test/test_priors.cxx +++ b/src/recon_test/test_priors.cxx @@ -20,7 +20,7 @@ \file \ingroup recon_test - \brief Test program for stir::QuadraticPrior + \brief Test program for stir::QuadraticPrior, RelativeDifferencePrior, and LogcoshPrior \par Usage @@ -30,12 +30,14 @@ where the argument is optional. See the class documentation for more info. \author Kris Thielemans + \author Robert Twyman */ #include "stir/VoxelsOnCartesianGrid.h" #include "stir/recon_buildblock/QuadraticPrior.h" #include "stir/recon_buildblock/RelativeDifferencePrior.h" -#include "stir/recon_buildblock/PLSPrior.h" +#include "stir/recon_buildblock/LogcoshPrior.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" @@ -49,13 +51,12 @@ #include #include -#include "stir/IO/OutputFileFormat.h" START_NAMESPACE_STIR /*! \ingroup test - \brief Test class for QuadraticPrior + \brief Test class for QuadraticPrior, RelativeDifferencePrior, and LogcoshPrior This test compares the result of GeneralisedPrior::compute_gradient() with a numerical gradient computed by using the @@ -193,17 +194,24 @@ run_tests() } std::cerr << "Tests for Relative Difference Prior\n"; { - // gamma and epsilon are default - RelativeDifferencePrior objective_function(false, 1.F, 2.F, 0.0001); + // gamma is default and epsilon is off + RelativeDifferencePrior objective_function(false, 1.F, 2.F, 0.F); this->run_tests_for_objective_function("RDP_no_kappa", objective_function, density_sptr); } - std::cerr << "Tests for PLSPrior\n"; + // Disabled PLS due to known issue +// std::cerr << "Tests 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 << "Tests for Logcosh Prior\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); + // scalar is off + LogcoshPrior objective_function(false, 1.F, 1.F); + this->run_tests_for_objective_function("Logcosh_no_kappa", objective_function, density_sptr); } }