Merge pull request #1291 from markus-jehl/issue/1290-speeding-up-scat…

…ter-simulation-for-blocks-geometry Added downsampling of BlocksOnCylindrical scanners in scatter
UCL · Aug 30, 2024 · 33a0106 · 33a0106
2 parents c51067d + f3e697d
commit 33a0106
Show file tree

Hide file tree

Showing 11 changed files with 740 additions and 97 deletions.
diff --git a/documentation/release_6.3.htm b/documentation/release_6.3.htm
@@ -0,0 +1,66 @@
+<!DOCTYPE HTML>
+<html lang="en">
+  <head>
+    <title>Summary of changes in STIR release 6.3</title>
+  </head>
+
+  <body>
+    <h1>Summary of changes in STIR release 6.3</h1>
+
+
+<h2>Overall summary</h2>
+
+
+<h2>Patch release info</h2>
+
+
+<h2> Summary for end users (also to be read by developers)</h2>
+
+
+<h3>New functionality</h3>
+<ul>
+  <li>
+    <code>ScatterSimulation</code> can now downsample the scanner transaxially (crystals per ring) for <code>BlocksOnCylindrical</code>,
+    scanners, which speeds up <code>ScatterEstimation</code> considerably. By default, downsampling the detectors per reading
+    is disabled for backwards compatibility.<br>
+    <a href=https://github.com/UCL/STIR/pull/1291>PR #1291</a>
+  </li>
+</ul>
+
+
+<h3>Changed functionality</h3>
+
+
+<h3>Bug fixes</h3>
+
+
+<h3>Build system</h3>
+
+
+<h3>Known problems</h3>
+    <p>See <a href=https://github.com/UCL/STIR/labels/bug>our issue tracker</a>.</p>
+
+
+<H2>What is new for developers (aside from what should be obvious from the above):</H2>
+
+
+<h3>Changed functionality</h3>
+
+
+<h3>Bug fixes</h3>
+
+
+<h3>Other code changes</h3>
+
+
+<h3>Test changes</h3>
+
+
+<h4>C++ tests</h4>
+
+
+<h4>recon_test_pack</h4>
+
+</body>
+
+</html>
diff --git a/src/buildblock/extend_projdata.cxx b/src/buildblock/extend_projdata.cxx
@@ -93,7 +93,7 @@ extend_segment(const SegmentBySinogram<float>& segment,
           if (extend_without_wrapping)
             {
               out[axial_pos][min_dim[2] + view_edge] = out[axial_pos][min_dim[2] + view_extension];
-              out[axial_pos][max_dim[2] - view_extension] = out[axial_pos][max_dim[2] - view_extension];
+              out[axial_pos][max_dim[2] - view_edge] = out[axial_pos][max_dim[2] - view_extension];
             }
           else if (flip_views)
             {

diff --git a/src/buildblock/interpolate_projdata.cxx b/src/buildblock/interpolate_projdata.cxx
diff --git a/src/include/stir/interpolate_projdata.h b/src/include/stir/interpolate_projdata.h
@@ -59,6 +59,8 @@ Succeeded interpolate_projdata(ProjData& proj_data_out,
                                const ProjData& proj_data_in,
                                const BasicCoordinate<3, BSpline::BSplineType>& these_types,
                                const bool remove_interleaving);
+Succeeded
+interpolate_blocks_on_cylindrical_projdata(ProjData& proj_data_out, const ProjData& proj_data_in, bool remove_interleaving);
 //@}
 
 END_NAMESPACE_STIR
diff --git a/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl b/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl
@@ -216,6 +216,7 @@ DataSymmetriesForBins_PET_CartesianGrid::find_sym_op_bin0(int segment_num, int v
       // No symmetry is implemented for generic scanner
       return new TrivialSymmetryOperation();
     }
+  return new TrivialSymmetryOperation();
 }
 
 // from symmetries
@@ -390,6 +391,8 @@ DataSymmetriesForBins_PET_CartesianGrid::find_sym_op_general_bin(int s, int segm
       // No symmetry is implemented for generic scanner
       return new TrivialSymmetryOperation();
     }
+
+  return new TrivialSymmetryOperation();
 }
 
 bool

diff --git a/src/include/stir/scatter/ScatterEstimation.h b/src/include/stir/scatter/ScatterEstimation.h
@@ -174,6 +174,12 @@ class ScatterEstimation : public ParsingObject
 
   inline void set_num_iterations(int);
 
+  inline unsigned int get_half_filter_width() const;
+  inline void set_half_filter_width(unsigned int);
+
+  inline void
+  set_downsample_scanner(bool downsample_scanner, int downsampled_number_of_rings = -1, int downsampled_detectors_per_ring = -1);
+
   void set_output_scatter_estimate_prefix(const std::string&);
   void set_export_scatter_estimates_of_each_iteration(bool);
 
@@ -382,6 +388,9 @@ class ScatterEstimation : public ParsingObject
   float min_scale_value;
 
   bool downsample_scanner_bool;
+  int downsampled_number_of_rings;
+  int downsampled_detectors_per_ring;
+
   //!
   unsigned int half_filter_width;
 

diff --git a/src/include/stir/scatter/ScatterEstimation.inl b/src/include/stir/scatter/ScatterEstimation.inl
@@ -66,4 +66,26 @@ ScatterEstimation::set_num_iterations(int arg)
   this->num_scatter_iterations = arg;
 }
 
+unsigned int
+ScatterEstimation::get_half_filter_width() const
+{
+  return this->half_filter_width;
+}
+
+void
+ScatterEstimation::set_half_filter_width(unsigned int arg)
+{
+  this->half_filter_width = arg;
+}
+
+void
+ScatterEstimation::set_downsample_scanner(bool downsample_scanner,
+                                          int downsampled_number_of_rings,
+                                          int downsampled_detectors_per_ring)
+{
+  this->downsample_scanner_bool = downsample_scanner;
+  this->downsampled_number_of_rings = downsampled_number_of_rings;
+  this->downsampled_detectors_per_ring = downsampled_detectors_per_ring;
+}
+
 END_NAMESPACE_STIR
diff --git a/src/scatter_buildblock/ScatterEstimation.cxx b/src/scatter_buildblock/ScatterEstimation.cxx
@@ -23,6 +23,7 @@
 #include "stir/recon_buildblock/ChainedBinNormalisation.h"
 #include "stir/ProjDataInterfile.h"
 #include "stir/ProjDataInMemory.h"
+#include "stir/inverse_SSRB.h"
 #include "stir/ExamInfo.h"
 #include "stir/ProjDataInfo.h"
 #include "stir/ProjDataInfoCylindricalNoArcCorr.h"
@@ -77,6 +78,8 @@ ScatterEstimation::set_defaults()
   this->override_scanner_template = true;
   this->override_density_image = true;
   this->downsample_scanner_bool = true;
+  this->downsampled_number_of_rings = -1;
+  this->downsampled_detectors_per_ring = -1;
   this->remove_interleaving = true;
   this->atten_image_filename = "";
   this->atten_coeff_filename = "";
@@ -124,6 +127,8 @@ ScatterEstimation::initialise_keymap()
   this->parser.add_parsing_key("Scatter Simulation type", &this->scatter_simulation_sptr);
   this->parser.add_key("scatter simulation parameter filename", &this->scatter_sim_par_filename);
   this->parser.add_key("use scanner downsampling in scatter simulation", &this->downsample_scanner_bool);
+  this->parser.add_key("override number of downsampled rings", &this->downsampled_number_of_rings);
+  this->parser.add_key("override number of downsampled detectors per ring", &this->downsampled_detectors_per_ring);
 
   this->parser.add_key("override attenuation image", &this->override_density_image);
   this->parser.add_key("override scanner template", &this->override_scanner_template);
@@ -597,7 +602,7 @@ ScatterEstimation::set_up()
     }
 
   if (this->downsample_scanner_bool)
-    this->scatter_simulation_sptr->downsample_scanner();
+    this->scatter_simulation_sptr->downsample_scanner(this->downsampled_number_of_rings, this->downsampled_detectors_per_ring);
 
   // Check if Load a mask proj_data
 
@@ -834,7 +839,7 @@ ScatterEstimation::process_data()
   float local_min_scale_value = 0.5f;
   float local_max_scale_value = 0.5f;
 
-  stir::BSpline::BSplineType spline_type = stir::BSpline::quadratic;
+  stir::BSpline::BSplineType spline_type = stir::BSpline::linear;
 
   // This has been set to 2D or 3D in the set_up()
   shared_ptr<ProjData> unscaled_est_projdata_sptr(
@@ -1026,16 +1031,13 @@ ScatterEstimation::process_data()
               shared_ptr<BinNormalisation> normalisation_factors_3d_sptr
                   = this->get_normalisation_object_sptr(this->multiplicative_binnorm_sptr);
 
-              upsample_and_fit_scatter_estimate(*scatter_estimate_sptr,
-                                                *this->input_projdata_sptr,
-                                                *temp_projdata,
-                                                *normalisation_factors_3d_sptr,
-                                                *this->input_projdata_sptr,
-                                                1.0f,
-                                                1.0f,
-                                                1,
-                                                spline_type,
-                                                false);
+              ProjDataInMemory interpolated_scatter(this->input_projdata_sptr->get_exam_info_sptr(),
+                                                    this->input_projdata_sptr->get_proj_data_info_sptr()->create_shared_clone());
+              inverse_SSRB(interpolated_scatter, *temp_projdata);
+              normalisation_factors_3d_sptr->set_up(this->input_projdata_sptr->get_exam_info_sptr(),
+                                                    this->input_projdata_sptr->get_proj_data_info_sptr()->create_shared_clone());
+              normalisation_factors_3d_sptr->undo(interpolated_scatter);
+              scatter_estimate_sptr->fill(interpolated_scatter);
             }
           else
             {

diff --git a/src/scatter_buildblock/ScatterSimulation.cxx b/src/scatter_buildblock/ScatterSimulation.cxx
@@ -841,19 +841,51 @@ ScatterSimulation::downsample_scanner(int new_num_rings, int new_num_dets)
   float approx_num_non_arccorrected_bins;
   if (new_scanner_sptr->get_scanner_geometry() != "Cylindrical")
     {
-      new_num_dets = this->proj_data_info_sptr->get_scanner_ptr()->get_num_detectors_per_ring();
-      approx_num_non_arccorrected_bins = this->proj_data_info_sptr->get_num_tangential_poss();
+      if (new_num_dets <= 0)
+        { // by default, do not downsample the detectors per ring for BlocksOnCylindrical
+          new_num_dets = this->proj_data_info_sptr->get_scanner_ptr()->get_num_detectors_per_ring();
+        }
 
+      // STIR does not like block spacings that are smaller than number of crystals times crystal spacing,
+      // therefore add a 1% extension on top for the downsampled scanner, to avoid running into floating point issues
+      const float block_spacing_factor = 1.01;
+
+      // extend the bins by a small amount to avoid edge-effects, at the expense of longer computation time
+      approx_num_non_arccorrected_bins = ceil(this->proj_data_info_sptr->get_num_tangential_poss() * float(new_num_dets)
+                                              / old_scanner_ptr->get_num_detectors_per_ring())
+                                         + 1;
       // preserve the length of the scanner, but place crystals equidistantly
       float scanner_length = old_scanner_ptr->get_axial_length();
       float new_ring_spacing = scanner_length / (new_num_rings - 1);
-
       new_scanner_sptr->set_num_axial_blocks_per_bucket(1);
+      new_scanner_sptr->set_num_transaxial_blocks_per_bucket(1);
+
       new_scanner_sptr->set_num_rings(new_num_rings);
+      new_scanner_sptr->set_num_detectors_per_ring(new_num_dets);
+
       new_scanner_sptr->set_axial_crystal_spacing(new_ring_spacing);
       new_scanner_sptr->set_ring_spacing(new_ring_spacing);
       new_scanner_sptr->set_num_axial_crystals_per_block(new_num_rings);
-      new_scanner_sptr->set_axial_block_spacing(new_ring_spacing * new_scanner_sptr->get_num_axial_crystals_per_block());
+      new_scanner_sptr->set_axial_block_spacing(new_ring_spacing * new_num_rings * block_spacing_factor);
+
+      float transaxial_bucket_width
+          = old_scanner_ptr->get_transaxial_block_spacing() * (old_scanner_ptr->get_num_transaxial_blocks_per_bucket() - 1)
+            + old_scanner_ptr->get_transaxial_crystal_spacing() * (old_scanner_ptr->get_num_transaxial_crystals_per_block() - 1);
+
+      int num_trans_buckets = old_scanner_ptr->get_num_transaxial_buckets();
+      // get a new number of detectors that is a multiple of the number of buckets to preserve scanner shape
+      new_scanner_sptr->set_num_detectors_per_ring(new_num_dets);
+      int new_transaxial_dets_per_bucket = new_num_dets / num_trans_buckets;
+      float new_det_spacing = transaxial_bucket_width / (new_transaxial_dets_per_bucket - 1);
+
+      new_scanner_sptr->set_num_transaxial_blocks_per_bucket(1);
+      new_scanner_sptr->set_num_transaxial_crystals_per_block(new_transaxial_dets_per_bucket);
+      new_scanner_sptr->set_transaxial_crystal_spacing(new_det_spacing);
+      new_scanner_sptr->set_transaxial_block_spacing(new_transaxial_dets_per_bucket * new_det_spacing * block_spacing_factor);
+      // avoid problems with Scanner checks by setting singles_units to 1 crystal
+      // (only used for dead-time processing in ECAY norm)
+      new_scanner_sptr->set_num_axial_crystals_per_singles_unit(1);
+      new_scanner_sptr->set_num_transaxial_crystals_per_singles_unit(1);
     }
   else
     {

diff --git a/src/test/CMakeLists.txt b/src/test/CMakeLists.txt
@@ -74,7 +74,7 @@ set(buildblock_simple_tests
 include(stir_test_exe_targets)
 
 foreach(source ${buildblock_simple_tests})
-  create_stir_test(${source} "buildblock;IO;buildblock;numerics_buildblock;display;IO;recon_buildblock;Shape_buildblock" "")
+  create_stir_test(${source} "buildblock;IO;buildblock;numerics_buildblock;display;IO;recon_buildblock;Shape_buildblock;scatter_buildblock" "")
 endforeach()
 #