From 4953e27e683845ee3967895ec955af09127425eb Mon Sep 17 00:00:00 2001
From: Antonia Reeves <133121395+rmapree@users.noreply.github.com>
Date: Wed, 20 Nov 2024 17:52:42 +0000
Subject: [PATCH] Fix SSRB for ACFs in scatter estimation (#1531)

* Scatter Estimation - add boolean for 2d proj data and make sure acfs 3d are obtained by averaging

* Update release_6.3.htm

---------

Co-authored-by: Kris Thielemans <KrisThielemans@users.noreply.github.com>
---
 documentation/release_6.3.htm                | 11 ++++++++++-
 src/include/stir/scatter/ScatterEstimation.h |  3 ++-
 src/scatter_buildblock/ScatterEstimation.cxx |  8 +++++---
 3 files changed, 17 insertions(+), 5 deletions(-)

diff --git a/documentation/release_6.3.htm b/documentation/release_6.3.htm
index 52370aee3..bdcc5dce7 100644
--- a/documentation/release_6.3.htm
+++ b/documentation/release_6.3.htm
@@ -63,7 +63,16 @@ <h3>Changed functionality</h3>
     </ul>
 
 <h3>Bug fixes</h3>
-
+    <ul>
+        <li>
+            Fixed a bug in the scatter estimation code (introduced in release 5.1.0) if input data is 3D and "cylindrical" (there was no bug for "blocksoncylindrical" data).
+            The scatter estimation runs on data constructed via SSRB. However, the attenuation correction factors were incorrectly obtained with adding of oblique segments (as opposed to averaging).
+            This resulted in intermediate images that had the wrong attenuation correction which were approximately num_segments times larger. This was compensated by the tail-fitting, but resulted in unexpected scale factors (scale factors were around 1/num_segments times what was expected).
+            This means that if you used the "min/max scale factor" feature in the scatter estimate, you will have to adjust your threshold values. Expected scatter tail-fitting scale factors should now be restored to ~1-1.5 (depending on the amount of multiple and out-of-FOV scatter).
+            See <a href="https://github.com/UCL/STIR/issues/1532">Issue #1532</a> for more detail. Fixed by using averaging functionality of SSRB instead of adding segments for attenuation correction factors.
+            <a href=https://github.com/UCL/STIR/pull/1531>PR #1531</a>
+        </li>
+    </ul>
 
 <h3>Build system</h3>
 
diff --git a/src/include/stir/scatter/ScatterEstimation.h b/src/include/stir/scatter/ScatterEstimation.h
index 785cd3b25..e3ed9a722 100644
--- a/src/include/stir/scatter/ScatterEstimation.h
+++ b/src/include/stir/scatter/ScatterEstimation.h
@@ -119,7 +119,8 @@ class ScatterEstimation : public ParsingObject
 
   //! make projdata 2D shared pointer
   shared_ptr<ProjData> make_2D_projdata_sptr(const shared_ptr<ProjData> in_3d_sptr);
-  shared_ptr<ProjData> make_2D_projdata_sptr(const shared_ptr<ProjData> in_3d_sptr, string template_filename);
+  shared_ptr<ProjData>
+  make_2D_projdata_sptr(const shared_ptr<ProjData> in_3d_sptr, string template_filename, const bool do_normalisation = false);
 
   //!
   //! \brief set_up
diff --git a/src/scatter_buildblock/ScatterEstimation.cxx b/src/scatter_buildblock/ScatterEstimation.cxx
index 07772ae99..efe3e2d9b 100644
--- a/src/scatter_buildblock/ScatterEstimation.cxx
+++ b/src/scatter_buildblock/ScatterEstimation.cxx
@@ -187,7 +187,9 @@ ScatterEstimation::make_2D_projdata_sptr(const shared_ptr<ProjData> in_3d_sptr)
 }
 
 shared_ptr<ProjData>
-ScatterEstimation::make_2D_projdata_sptr(const shared_ptr<ProjData> in_3d_sptr, string template_filename)
+ScatterEstimation::make_2D_projdata_sptr(const shared_ptr<ProjData> in_3d_sptr,
+                                         string template_filename,
+                                         const bool do_normalisation)
 {
   shared_ptr<ProjData> out_2d_sptr;
   if (in_3d_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_scanner_geometry() == "Cylindrical")
@@ -198,7 +200,7 @@ ScatterEstimation::make_2D_projdata_sptr(const shared_ptr<ProjData> in_3d_sptr,
                                          this->input_projdata_2d_sptr->get_exam_info_sptr(),
                                          this->input_projdata_2d_sptr->get_proj_data_info_sptr()->create_shared_clone());
 
-      SSRB(*out_2d_sptr, *in_3d_sptr, false);
+      SSRB(*out_2d_sptr, *in_3d_sptr, do_normalisation);
     }
   else
     {
@@ -658,7 +660,7 @@ ScatterEstimation::set_up_iterative(shared_ptr<IterativeReconstruction<Discretis
       info("ScatterEstimation: Running SSRB on attenuation correction coefficients ...");
 
       std::string out_filename = "tmp_atten_sino_2d.hs";
-      atten_projdata_2d_sptr = make_2D_projdata_sptr(tmp_atten_projdata_sptr, out_filename);
+      atten_projdata_2d_sptr = make_2D_projdata_sptr(tmp_atten_projdata_sptr, out_filename, true);
     }
   else
     {