From 437b1e6917bb6f4f28790244f8bf4f6be170b232 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Mon, 22 Jul 2024 22:28:25 +0100 Subject: [PATCH] [BUGFIX] remove thresholding in Hessian for Poisson...ProjData accumulate_sub_Hessian_times_input_without_penalty used divide_and_truncate, but this is incorrect as its threshold strategies are not appropriate for (measured * forward_input)/forward_current^2. Most obvious example is if forward_input contains negatives. We now use 0/x=0 (first thresholding measured data to be non-negative). https://github.com/UCL/STIR/issues/1461 --- ...LikelihoodWithLinearModelForMeanAndProjData.cxx | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx index 98fba9935..bd28fc0c1 100644 --- a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx +++ b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx @@ -27,6 +27,7 @@ #include "stir/recon_buildblock/TrivialBinNormalisation.h" #include "stir/Succeeded.h" #include "stir/RelatedViewgrams.h" +#include "stir/ArrayFunction.h" #include "stir/stream.h" #include "stir/info.h" #include "stir/warning.h" @@ -1140,11 +1141,22 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjData::actual_accumulat // final_viewgram starts as measured data RelatedViewgrams final_viewgram = this->get_proj_data().get_related_viewgrams(vg_idx_to_process[i], symmetries_sptr); { +#if 0 // Mult input_viewgram final_viewgram *= input_viewgrams_vec[i]; + // Divide final_viewgram by ybar_sq_viewgram int tmp1 = 0, tmp2 = 0; // ignore counters returned by divide_and_truncate - // Divide final_viewgeam by ybar_sq_viewgram divide_and_truncate(final_viewgram, ybar_sq_viewgram, 0, tmp1, tmp2); +#else + // use division where 0/x = 0, but truncate measured data to be >= 0 (as in the value and gradient) + constexpr auto div = [](float a, float b) { return a <= 0.F ? 0.F : a / b; }; + auto f_iter = final_viewgram.begin(); + auto yb_sq_iter = ybar_sq_viewgram.begin(); + for (; f_iter != final_viewgram.end(); ++f_iter, ++yb_sq_iter) + in_place_apply_binary_func_element_wise(*f_iter, *yb_sq_iter, div); + // Mult input_viewgram + final_viewgram *= input_viewgrams_vec[i]; +#endif } // back-project final_viewgram