From df5058a395b65d67c415843318743539f839a000 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Sun, 12 May 2024 10:00:56 +0100 Subject: [PATCH] minor clean-up and clarifications --- .gitignore | 1 + .../stir/recon_buildblock/ProjMatrixByBin.inl | 2 +- .../ProjMatrixElemsForOneBin.h | 8 ++-- src/recon_buildblock/distributable.cxx | 48 +++++++++---------- 4 files changed, 29 insertions(+), 30 deletions(-) diff --git a/.gitignore b/.gitignore index d3689c4f86..075f775d2a 100644 --- a/.gitignore +++ b/.gitignore @@ -57,6 +57,7 @@ install_manifest.txt ~$* *~ *.bak +\#*# # Ignore Visual Studio User-specific files *.suo diff --git a/src/include/stir/recon_buildblock/ProjMatrixByBin.inl b/src/include/stir/recon_buildblock/ProjMatrixByBin.inl index 7ff4ac1530..8f041306ec 100644 --- a/src/include/stir/recon_buildblock/ProjMatrixByBin.inl +++ b/src/include/stir/recon_buildblock/ProjMatrixByBin.inl @@ -149,7 +149,7 @@ ProjMatrixByBin::get_tof_value(const float d1, const float d2) const if ((d1_n >= 4.f && d2_n >= 4.f) || (d1_n <= -4.f && d2_n <= -4.f)) return 0.F; else - return 0.5f * (erf_interpolation(d2_n) - erf_interpolation(d1_n)); + return static_cast(0.5 * (erf_interpolation(d2_n) - erf_interpolation(d1_n))); } END_NAMESPACE_STIR diff --git a/src/include/stir/recon_buildblock/ProjMatrixElemsForOneBin.h b/src/include/stir/recon_buildblock/ProjMatrixElemsForOneBin.h index a006b02f4e..a619096ee4 100644 --- a/src/include/stir/recon_buildblock/ProjMatrixElemsForOneBin.h +++ b/src/include/stir/recon_buildblock/ProjMatrixElemsForOneBin.h @@ -178,14 +178,14 @@ class ProjMatrixElemsForOneBin //******************** projection operations ********************// - //! back project a single bin + //! back project a single bin (accumulates) void back_project(DiscretisedDensity<3, float>&, const Bin&) const; - //! forward project into a single bin + //! forward project into a single bin (accumulates) void forward_project(Bin&, const DiscretisedDensity<3, float>&) const; - //! back project related bins + //! back project related bins (accumulates) void back_project(DiscretisedDensity<3, float>&, const RelatedBins&) const; - //! forward project related bins + //! forward project related bins (accumulates) void forward_project(RelatedBins&, const DiscretisedDensity<3, float>&) const; private: diff --git a/src/recon_buildblock/distributable.cxx b/src/recon_buildblock/distributable.cxx index 6eb092c6b5..c2f776a531 100644 --- a/src/recon_buildblock/distributable.cxx +++ b/src/recon_buildblock/distributable.cxx @@ -685,10 +685,10 @@ LM_distributable_computation_template(const shared_ptr PM_sptr, PM_sptr->get_proj_matrix_elems_for_one_bin(local_row[thread_num], measured_bin); call_back(*local_output_image_sptrs[thread_num], - measured_bin, - *input_image_ptr, local_row[thread_num], - has_add ? record.my_corr : 0); + has_add ? record.my_corr : 0.F, + measured_bin, + *input_image_ptr); } } #ifdef STIR_OPENMP @@ -710,25 +710,25 @@ LM_distributable_computation_template(const shared_ptr PM_sptr, constexpr float max_quotient = 10000.F; +/* gradient without the sensitivity term + +\sum_e A_e^t (y_e/(A_e lambda+ c)) +*/ inline void LM_gradient(DiscretisedDensity<3, float>& output_image, - const Bin& measured_bin, - const DiscretisedDensity<3, float>& input_image, const ProjMatrixElemsForOneBin& row, - const float add_term) + const float add_term, + const Bin& measured_bin, + const DiscretisedDensity<3, float>& input_image) { Bin fwd_bin = measured_bin; fwd_bin.set_bin_value(0.0f); row.forward_project(fwd_bin, input_image); + const auto fwd = fwd_bin.get_bin_value() + add_term; - if (add_term) - fwd_bin.set_bin_value(fwd_bin.get_bin_value() + add_term); - - float measured_div_fwd = 0.F; - if (measured_bin.get_bin_value() <= max_quotient * fwd_bin.get_bin_value()) - measured_div_fwd = measured_bin.get_bin_value() / fwd_bin.get_bin_value(); - else - return; + if (measured_bin.get_bin_value() > max_quotient * fwd) + return; // cancel singularity + const auto measured_div_fwd = measured_bin.get_bin_value() / fwd; fwd_bin.set_bin_value(measured_div_fwd); row.back_project(output_image, fwd_bin); @@ -736,30 +736,28 @@ LM_gradient(DiscretisedDensity<3, float>& output_image, /* Hessian -\sum_e A_e^t (y_e/(A_e lambda+ c)^2 A_e x) +\sum_e A_e^t (y_e/(A_e lambda+ c)^2 A_e rhs) */ inline void LM_Hessian(DiscretisedDensity<3, float>& output_image, + const ProjMatrixElemsForOneBin& row, + const float add_term, const Bin& measured_bin, const DiscretisedDensity<3, float>& input_image, - const DiscretisedDensity<3, float>& rhs, - const ProjMatrixElemsForOneBin& row, - const float add_term) + const DiscretisedDensity<3, float>& rhs) { Bin fwd_bin = measured_bin; fwd_bin.set_bin_value(0.0f); row.forward_project(fwd_bin, input_image); + const auto fwd = fwd_bin.get_bin_value() + add_term; - if (add_term) - fwd_bin.set_bin_value(fwd_bin.get_bin_value() + add_term); - - if (measured_bin.get_bin_value() > max_quotient * fwd_bin.get_bin_value()) - return; - float measured_div_fwd2 = measured_bin.get_bin_value() / square(fwd_bin.get_bin_value()); + if (measured_bin.get_bin_value() > max_quotient * fwd) + return; // cancel singularity + const auto measured_div_fwd2 = measured_bin.get_bin_value() / square(fwd); + // forward project rhs fwd_bin.set_bin_value(0.0f); row.forward_project(fwd_bin, rhs); - if (fwd_bin.get_bin_value() == 0) return;