Fix SSRB for ACFs in scatter estimation

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>
 

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

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
     {