diff --git a/CHANGES.md b/CHANGES.md index e543f7e4a..680654ddf 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -10,6 +10,11 @@ - implemented basic-functionality listmode data class in C++ and Python - added objective function type for lismode reconstruction - added new demo script for the reconstruction from listmode data + - provided gradient-computing methods with return via optional argument out + in addition to the standard return, ensuring that no temorary copies of the + gradient data are created + - provided prior and objective function objects with methods for computing + the product of the Hessian and a vector ## v3.6.0 diff --git a/src/xSTIR/cSTIR/cstir.cpp b/src/xSTIR/cSTIR/cstir.cpp index bd431e98b..29680cb34 100644 --- a/src/xSTIR/cSTIR/cstir.cpp +++ b/src/xSTIR/cSTIR/cstir.cpp @@ -1,6 +1,6 @@ /* SyneRBI Synergistic Image Reconstruction Framework (SIRF) -Copyright 2017 - 2023 Rutherford Appleton Laboratory STFC +Copyright 2017 - 2024 Rutherford Appleton Laboratory STFC Copyright 2018 - 2024 University College London. This is software developed for the Collaborative Computational @@ -417,12 +417,25 @@ void* cSTIR_objectFromFile(const char* name, const char* filename) CATCH; } +typedef xSTIR_PoissonLLhLinModMeanListDataProjMatBin3DF LMObjFun; + +extern "C" +void* cSTIR_objFunListModeSetInterval(void* ptr_f, size_t ptr_data) +{ + try { + auto& objFun = objectFromHandle(ptr_f); + float* data = (float*)ptr_data; + objFun.set_time_interval((double)data[0], (double)data[1]); + return (void*)new DataHandle; + } + CATCH; +} + extern "C" void* cSTIR_setListmodeToSinogramsInterval(void* ptr_lm2s, size_t ptr_data) { try { - ListmodeToSinograms& lm2s = - objectFromHandle(ptr_lm2s); + auto& lm2s = objectFromHandle(ptr_lm2s); float *data = (float *)ptr_data; lm2s.set_time_interval((double)data[0], (double)data[1]); return (void*)new DataHandle; @@ -1249,6 +1262,50 @@ cSTIR_computeObjectiveFunctionGradientNotDivided(void* ptr_f, void* ptr_i, int s CATCH; } +extern "C" +void* +cSTIR_objectiveFunctionAccumulateHessianTimesInput + (void* ptr_fun, void* ptr_est, void* ptr_inp, int subset, void* ptr_out) +{ + try { + auto& fun = objectFromHandle(ptr_fun); + auto& est = objectFromHandle(ptr_est); + auto& inp = objectFromHandle(ptr_inp); + auto& out = objectFromHandle(ptr_out); + auto& curr_est = est.data(); + auto& input = inp.data(); + auto& output = out.data(); + if (subset >= 0) + fun.accumulate_sub_Hessian_times_input(output, curr_est, input, subset); + else { + for (int s = 0; s < fun.get_num_subsets(); s++) { + fun.accumulate_sub_Hessian_times_input(output, curr_est, input, s); + } + } + return (void*) new DataHandle; + } + CATCH; +} + +extern "C" +void* +cSTIR_objectiveFunctionComputeHessianTimesInput + (void* ptr_fun, void* ptr_est, void* ptr_inp, int subset, void* ptr_out) +{ + try { + auto& fun = objectFromHandle(ptr_fun); + auto& est = objectFromHandle(ptr_est); + auto& inp = objectFromHandle(ptr_inp); + auto& out = objectFromHandle(ptr_out); + auto& curr_est = est.data(); + auto& input = inp.data(); + auto& output = out.data(); + fun.multiply_with_Hessian(output, curr_est, input, subset); + return (void*) new DataHandle; + } + CATCH; +} + extern "C" void* cSTIR_setupPrior(void* ptr_p, void* ptr_i) @@ -1288,7 +1345,7 @@ void* cSTIR_priorGradient(void* ptr_p, void* ptr_i) { try { - Prior3DF& prior = objectFromHandle(ptr_p); + Prior3DF& prior = objectFromHandle >(ptr_p); STIRImageData& id = objectFromHandle(ptr_i); Image3DF& image = id.data(); shared_ptr sptr(new STIRImageData(image)); @@ -1299,6 +1356,42 @@ cSTIR_priorGradient(void* ptr_p, void* ptr_i) CATCH; } +extern "C" +void* +cSTIR_priorAccumulateHessianTimesInput(void* ptr_prior, void* ptr_out, void* ptr_cur, void* ptr_inp) +{ + try { + auto& prior = objectFromHandle >(ptr_prior); + auto& out = objectFromHandle(ptr_out); + auto& cur = objectFromHandle(ptr_cur); + auto& inp = objectFromHandle(ptr_inp); + auto& output = out.data(); + auto& current = cur.data(); + auto& input = inp.data(); + prior.accumulate_Hessian_times_input(output, current, input); + return (void*) new DataHandle; + } + CATCH; +} + +extern "C" +void* +cSTIR_priorComputeHessianTimesInput(void* ptr_prior, void* ptr_out, void* ptr_cur, void* ptr_inp) +{ + try { + auto& prior = objectFromHandle(ptr_prior); + auto& out = objectFromHandle(ptr_out); + auto& cur = objectFromHandle(ptr_cur); + auto& inp = objectFromHandle(ptr_inp); + auto& output = out.data(); + auto& current = cur.data(); + auto& input = inp.data(); + prior.multiply_with_Hessian(output, current, input); + return (void*) new DataHandle; + } + CATCH; +} + extern "C" void* cSTIR_computePriorGradient(void* ptr_p, void* ptr_i, void* ptr_g) diff --git a/src/xSTIR/cSTIR/cstir_p.cpp b/src/xSTIR/cSTIR/cstir_p.cpp index f1c8d056c..af5859855 100644 --- a/src/xSTIR/cSTIR/cstir_p.cpp +++ b/src/xSTIR/cSTIR/cstir_p.cpp @@ -1012,6 +1012,11 @@ sirf::cSTIR_setPoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProj else if (sirf::iequals(name, "cache_max_size")) { obj_fun.set_cache_max_size(dataFromHandle(hv)); } + else if (sirf::iequals(name, "subsensitivity_filenames")) + { + std::string s(charDataFromDataHandle(hv)); + obj_fun.set_subsensitivity_filenames(s.c_str()); + } else return parameterNotFound(name, __FILE__, __LINE__); return new DataHandle; diff --git a/src/xSTIR/cSTIR/include/sirf/STIR/cstir.h b/src/xSTIR/cSTIR/include/sirf/STIR/cstir.h index e4530c93f..7433da844 100644 --- a/src/xSTIR/cSTIR/include/sirf/STIR/cstir.h +++ b/src/xSTIR/cSTIR/include/sirf/STIR/cstir.h @@ -65,6 +65,7 @@ extern "C" { void* cSTIR_convertListmodeToSinograms(void* ptr); void* cSTIR_computeRandoms(void* ptr); void* cSTIR_lm_num_prompts_exceeds_threshold(void* ptr, const float threshold); + void* cSTIR_objFunListModeSetInterval(void* ptr_f, size_t ptr_data); // Data processor methods void* cSTIR_setupImageDataProcessor(const void* ptr_p, void* ptr_i); @@ -153,11 +154,19 @@ extern "C" { (void* ptr_f, void* ptr_i, int subset); void* cSTIR_computeObjectiveFunctionGradientNotDivided (void* ptr_f, void* ptr_i, int subset, void* ptr_g); + void* cSTIR_objectiveFunctionAccumulateHessianTimesInput + (void* ptr_fun, void* ptr_est, void* ptr_inp, int subset, void* ptr_out); + void* cSTIR_objectiveFunctionComputeHessianTimesInput + (void* ptr_fun, void* ptr_est, void* ptr_inp, int subset, void* ptr_out); // Prior methods void* cSTIR_setupPrior(void* ptr_p, void* ptr_i); void* cSTIR_priorValue(void* ptr_p, void* ptr_i); void* cSTIR_priorGradient(void* ptr_p, void* ptr_i); + void* cSTIR_priorAccumulateHessianTimesInput + (void* ptr_prior, void* ptr_out, void* ptr_curr, void* ptr_inp); + void* cSTIR_priorComputeHessianTimesInput + (void* ptr_prior, void* ptr_out, void* ptr_cur, void* ptr_inp); void* cSTIR_computePriorGradient(void* ptr_p, void* ptr_i, void* ptr_g); void* cSTIR_PLSPriorAnatomicalGradient(void* ptr_p, int dir); diff --git a/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h b/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h index 374ccd9bf..f3fd7f715 100644 --- a/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h +++ b/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h @@ -1031,6 +1031,12 @@ The actual algorithm is described in class xSTIR_GeneralisedPrior3DF : public stir::GeneralisedPrior < Image3DF > { public: + void multiply_with_Hessian(Image3DF& output, const Image3DF& curr_image_est, + const Image3DF& input) const + { + output.fill(0.0); + accumulate_Hessian_times_input(output, curr_image_est, input); + } // bool post_process() { // return post_processing(); // } @@ -1067,6 +1073,19 @@ The actual algorithm is described in class xSTIR_GeneralisedObjectiveFunction3DF : public stir::GeneralisedObjectiveFunction < Image3DF > { public: + void multiply_with_Hessian(Image3DF& output, const Image3DF& curr_image_est, + const Image3DF& input, const int subset) const + { + output.fill(0.0); + if (subset >= 0) + accumulate_sub_Hessian_times_input(output, curr_image_est, input, subset); + else { + for (int s = 0; s < get_num_subsets(); s++) { + accumulate_sub_Hessian_times_input(output, curr_image_est, input, s); + } + } + } + // bool post_process() { // return post_processing(); // } @@ -1118,6 +1137,15 @@ The actual algorithm is described in set_cache_path(filepath); } + void set_time_interval(double start, double stop) + { + std::pair interval(start, stop); + std::vector < std::pair > intervals; + intervals.push_back(interval); + frame_defs = stir::TimeFrameDefinitions(intervals); + do_time_frame = true; + } + private: //std::shared_ptr sptr_ad_; std::shared_ptr sptr_am_; diff --git a/src/xSTIR/pSTIR/STIR.py b/src/xSTIR/pSTIR/STIR.py index 7388f2aff..4d0ab9b7f 100644 --- a/src/xSTIR/pSTIR/STIR.py +++ b/src/xSTIR/pSTIR/STIR.py @@ -2333,6 +2333,23 @@ def set_up(self, image): """Sets up.""" try_calling(pystir.cSTIR_setupPrior(self.handle, image.handle)) + def accumulate_Hessian_times_input(self, current_estimate, input_, out=None): + """Computes the multiplication of the Hessian with a vector and adds it to output. + """ + if out is None or out.handle is None: + out = input_.get_uniform_copy(0.0) + try_calling(pystir.cSTIR_priorAccumulateHessianTimesInput + (self.handle, out.handle, current_estimate.handle, input_.handle)) + return out + + def multiply_with_Hessian(self, current_estimate, input_, out=None): + """Computes the multiplication of the Hessian at current_estimate with a vector. + """ + if out is None or out.handle is None: + out = input_.get_uniform_copy(0.0) + try_calling(pystir.cSTIR_priorComputeHessianTimesInput + (self.handle, current_estimate.handle, input_.handle, out.handle)) + return out class QuadraticPrior(Prior): r"""Class for the prior that is a quadratic function of the image values. @@ -2732,6 +2749,24 @@ def get_subset_gradient(self, image, subset, out=None): """ return self.gradient(image, subset, out) + def accumulate_Hessian_times_input(self, current_estimate, input_, subset=-1, out=None): + """Computes the multiplication of the Hessian at current_estimate with a vector and adds it to output. + """ + if out is None or out.handle is None: + out = input_.clone() + try_calling(pystir.cSTIR_objectiveFunctionAccumulateHessianTimesInput + (self.handle, current_estimate.handle, input_.handle, subset, out.handle)) + return out + + def multiply_with_Hessian(self, current_estimate, input_, subset=-1, out=None): + """Computes the multiplication of the Hessian at current_estimate with a vector. + """ + if out is None or out.handle is None: + out = input_.get_uniform_copy(0.0) + try_calling(pystir.cSTIR_objectiveFunctionComputeHessianTimesInput + (self.handle, current_estimate.handle, input_.handle, subset, out.handle)) + return out + @abc.abstractmethod def get_subset_sensitivity(self, subset): #print('in base class ObjectiveFunction') @@ -2854,7 +2889,8 @@ def set_acquisition_data(self, ad): self.handle, self.name, 'acquisition_data', ad.handle) -class PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin(ObjectiveFunction): +class PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin(PoissonLogLikelihoodWithLinearModelForMean): +#(ObjectiveFunction): """Class for a STIR type of Poisson loglikelihood object for listmode data. Specifically, PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin. @@ -2878,6 +2914,17 @@ def set_cache_path(self, path): def get_cache_path(self): return parms.char_par(self.handle, self.name, 'cache_path') + def set_time_interval(self, start, stop): + """Sets the time interval. + + Only data scanned during this time interval will be converted. + """ + interval = numpy.ndarray((2,), dtype=numpy.float32) + interval[0] = start + interval[1] = stop + try_calling(pystir.cSTIR_objFunListModeSetInterval( + self.handle, interval.ctypes.data)) + def set_acquisition_data(self, ad): assert_validity(ad, ListmodeData) parms.set_parameter( @@ -2918,6 +2965,9 @@ def set_cache_max_size(self, diff): def get_cache_max_size(self): return parms.int_par(self.handle, self.name, 'cache_max_size') + def set_subsensitivity_filenames(self, names): + return parms.set_char_par(self.handle, self.name, 'subsensitivity_filenames', names) + def get_subsensitivity_filenames(self): return parms.char_par(self.handle, self.name, 'subsensitivity_filenames') diff --git a/src/xSTIR/pSTIR/tests/test_ObjectiveFunction.py b/src/xSTIR/pSTIR/tests/test_ObjectiveFunction.py index 705df011e..016b42057 100644 --- a/src/xSTIR/pSTIR/tests/test_ObjectiveFunction.py +++ b/src/xSTIR/pSTIR/tests/test_ObjectiveFunction.py @@ -38,7 +38,8 @@ def setUp(self): templ = pet.AcquisitionData(os.path.join(data_path,'template_sinogram.hs')) am.set_up(templ,image) acquired_data=am.forward(image) - + am.set_background_term(acquired_data*0 + numpy.mean(acquired_data.as_array())) + am.set_up(templ,image) obj_fun = pet.make_Poisson_loglikelihood(acquired_data) obj_fun.set_acquisition_model(am) obj_fun.set_up(image) @@ -56,3 +57,21 @@ def test_Poisson_loglikelihood_call(self): numpy.testing.assert_almost_equal(a,b) + def test_Hessian(self, subset=-1, eps=1e-3): + """Checks that grad(x + dx) - grad(x) is close to H(x)*dx + """ + x = self.image + dx = x.clone() + dx *= eps/dx.norm() + dx += eps/2 + y = x + dx + gx = self.obj_fun.gradient(x, subset) + gy = self.obj_fun.gradient(y, subset) + dg = gy - gx + Hdx = self.obj_fun.multiply_with_Hessian(x, dx, subset) + q = (dg - Hdx).norm()/dg.norm() + print('norm of grad(x + dx) - grad(x): %f' % dg.norm()) + print('norm of H(x)*dx: %f' % Hdx.norm()) + print('relative difference: %f' % q) + assert q <= .002 + diff --git a/src/xSTIR/pSTIR/tests/tests_qp_lc_rdp.py b/src/xSTIR/pSTIR/tests/tests_qp_lc_rdp.py index 7de07ad71..e50f0c393 100644 --- a/src/xSTIR/pSTIR/tests/tests_qp_lc_rdp.py +++ b/src/xSTIR/pSTIR/tests/tests_qp_lc_rdp.py @@ -4,7 +4,7 @@ v{version} Usage: - tests_qp_rdp [--help | options] + tests_qp_lc_rdp [--help | options] Options: -r, --record record the measurements rather than check them @@ -14,27 +14,52 @@ {licence} """ -from sirf.STIR import * +import sirf.STIR from sirf.Utilities import runner, RE_PYEXT, __license__, examples_data_path, pTest -__version__ = "1.0.0" -__author__ = "Imraj Singh" +__version__ = "2.0.0" +__author__ = "Imraj Singh, Evgueni Ovtchinnikov, Kris Thielemans" +def test_Hessian(test, prior, x, eps=1e-3): + """Checks that grad(x + dx) - grad(x) is close to H(x)*dx + """ + if x.norm() > 0: + dx = x.clone() + dx *= eps/dx.norm() + dx += eps/2 + else: + dx = x + eps + y = x + dx + gx = prior.gradient(x.clone()) + gy = prior.gradient(y.clone()) + dg = gy - gx + Hdx = prior.multiply_with_Hessian(x, dx) + gynorm = gy.norm() + if gynorm > 0: + q = (dg - Hdx).norm()/gynorm + else: + q = (dg - Hdx).norm() + #print('norm of grad(x + dx) - grad(x): %f' % dg.norm()) + #print('norm of H(x)*dx: %f' % Hdx.norm()) + #print('relative difference: %g' % q) + test.check_if_less(q, .01*eps) + + def test_main(rec=False, verb=False, throw=True): datafile = RE_PYEXT.sub(".txt", __file__) test = pTest(datafile, rec, throw=throw) test.verbose = verb - msg_red = MessageRedirector(warn=None) + _ = sirf.STIR.MessageRedirector(warn=None) - im_thorax = ImageData(examples_data_path('PET') + '/thorax_single_slice/emission.hv') + im_thorax = sirf.STIR.ImageData(examples_data_path('PET') + '/thorax_single_slice/emission.hv') im_1 = im_thorax.get_uniform_copy(1) im_2 = im_thorax.get_uniform_copy(2) for im in [im_thorax, im_1, im_2]: for penalisation_factor in [0,1,4]: for kappa in [True, False]: - for prior in [QuadraticPrior(), LogcoshPrior(), RelativeDifferencePrior()]: + for prior in [sirf.STIR.QuadraticPrior(), sirf.STIR.LogcoshPrior(), sirf.STIR.RelativeDifferencePrior()]: if kappa: prior.set_kappa(im_thorax) # Check if the kappa is stored/returned correctly @@ -56,6 +81,10 @@ def test_main(rec=False, verb=False, throw=True): prior.set_kappa(im_thorax*2) prior.set_up(im) test.check_if_equal_within_tolerance(prior.get_gradient(im).norm(), grad_norm) + + if isinstance(prior, sirf.STIR.RelativeDifferencePrior): + prior.set_epsilon(im.max()*.01) + test_Hessian(test, prior, im, 0.03) return test.failed, test.ntest