Added downsampling of BlocksOnCylindrical scanners in scatter

src/buildblock/interpolate_projdata.cxx

@@ -239,17 +239,9 @@ interpolate_projdata(ProjData& proj_data_out,
           };
         }
       else
-        { // for BlocksOnCylindrical, views and tangential positions are not subsampled and can be mapped 1:1
-          if (proj_data_in_info.get_num_tangential_poss() != proj_data_out_info.get_num_tangential_poss())
-            {
-              error("Interpolation of BlocksOnCylindrical scanners assumes that number of tangential positions "
-                    "is the same in the downsampled scanner.");
-            }
-          if (proj_data_in_info.get_num_views() != proj_data_out_info.get_num_views())
-            {
-              error("Interpolation of BlocksOnCylindrical scanners assumes that number of views "
-                    "is the same in the downsampled scanner.");
-            }
+        { // for BlocksOnCylindrical, views and tangential positions are scaled by a fixed value
+          auto scale_factor
+              = (double)proj_data_out_info.get_num_tangential_poss() / (double)proj_data_in_info.get_num_tangential_poss();
 
           // only extending in axial direction - an extension of 2 was found to be sufficient
           proj_data_interpolator.set_coef(extend_segment(segment, 0, 2, 0));
@@ -267,7 +259,7 @@ interpolate_projdata(ProjData& proj_data_out,
             }
 
           // define a function to translate indices in the output proj data to indices in input proj data
-          index_converter = [&proj_data_out_info, m_offset, m_sampling](
+          index_converter = [&proj_data_out_info, m_offset, m_sampling, scale_factor](
                                 const BasicCoordinate<3, int>& index_out) -> BasicCoordinate<3, double> {
             // translate index on output to coordinate
             auto bin
@@ -277,8 +269,8 @@ interpolate_projdata(ProjData& proj_data_out,
             // translate to indices in input proj data
             BasicCoordinate<3, double> index_in;
             index_in[1] = (out_m - m_offset) / m_sampling;
-            index_in[2] = index_out[2];
-            index_in[3] = index_out[3];
+            index_in[2] = index_out[2] / scale_factor;
+            index_in[3] = index_out[3] / scale_factor;
 
             return index_in;
           };

src/include/stir/scatter/ScatterEstimation.h

@@ -171,6 +171,8 @@ class ScatterEstimation : public ParsingObject
 
   inline void set_num_iterations(int);
 
+  inline void set_downsample_scanner(bool downsample_scanner, int downsampled_detectors_per_ring = 0);
+
   void set_output_scatter_estimate_prefix(const std::string&);
   void set_export_scatter_estimates_of_each_iteration(bool);
 
@@ -379,6 +381,8 @@ class ScatterEstimation : public ParsingObject
   float min_scale_value;
 
   bool downsample_scanner_bool;
+  int downsampled_detectors_per_ring;
+
   //!
   unsigned int half_filter_width;
 

src/include/stir/scatter/ScatterEstimation.inl

@@ -64,4 +64,11 @@ ScatterEstimation::set_num_iterations(int arg)
   this->num_scatter_iterations = arg;
 }
 
+void
+ScatterEstimation::set_downsample_scanner(bool downsample_scanner, int downsampled_detectors_per_ring)
+{
+  this->downsample_scanner_bool = downsample_scanner;
+  this->downsampled_detectors_per_ring = downsampled_detectors_per_ring;
+}
+
 END_NAMESPACE_STIR

src/scatter_buildblock/ScatterEstimation.cxx

@@ -77,6 +77,7 @@ ScatterEstimation::set_defaults()
   this->override_scanner_template = true;
   this->override_density_image = true;
   this->downsample_scanner_bool = true;
+  this->downsampled_detectors_per_ring = 0;
   this->remove_interleaving = true;
   this->atten_image_filename = "";
   this->atten_coeff_filename = "";
@@ -124,6 +125,7 @@ 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 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);
@@ -585,7 +587,7 @@ ScatterEstimation::set_up()
     }
 
   if (this->downsample_scanner_bool)
-    this->scatter_simulation_sptr->downsample_scanner();
+    this->scatter_simulation_sptr->downsample_scanner(-1, this->downsampled_detectors_per_ring);
 
   // Check if Load a mask proj_data
 

src/scatter_buildblock/ScatterSimulation.cxx

@@ -840,36 +840,39 @@ 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();
+        }
+
+      // 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 the following includes gaps
       float scanner_length_block = new_scanner_sptr->get_num_axial_buckets() * new_scanner_sptr->get_num_axial_blocks_per_bucket()
                                    * new_scanner_sptr->get_axial_block_spacing();
       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_transaxial_blocks_per_bucket(1);
 
       new_scanner_sptr->set_num_rings(new_num_rings);
-      //        float transaxial_bucket_spacing=old_scanner_ptr->get_transaxial_block_spacing()
-      //                *old_scanner_ptr->get_num_transaxial_blocks_per_bucket();
+      float transaxial_bucket_spacing
+          = old_scanner_ptr->get_transaxial_block_spacing() * old_scanner_ptr->get_num_transaxial_blocks_per_bucket();
       float new_ring_spacing = scanner_length_block / new_scanner_sptr->get_num_rings();
-      //        int num_trans_buckets=old_scanner_ptr->get_num_transaxial_buckets();
+      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
-      //        float frac,whole;
-      //        frac = std::modf(float(new_num_dets/new_scanner_sptr->get_num_transaxial_buckets()), &whole);
-      //        int newest_num_dets=whole*new_scanner_sptr->get_num_transaxial_buckets();
-      //        new_scanner_sptr->set_num_detectors_per_ring(newest_num_dets);
-      //        int new_transaxial_dets_per_bucket=newest_num_dets/num_trans_buckets;
-      //        float new_det_spacing=transaxial_bucket_spacing/new_transaxial_dets_per_bucket;
+      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_spacing / new_transaxial_dets_per_bucket;
 
       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);
 
-      //        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_det_spacing
-      //                    * new_scanner_sptr->get_num_transaxial_crystals_per_block());
+      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(transaxial_bucket_spacing);
     }
   else
     {