Added downsampling of BlocksOnCylindrical scanners in scatter

src/buildblock/interpolate_projdata.cxx

@@ -249,16 +249,20 @@ 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())
+  { // for BlocksOnCylindrical, views and tangential positions are scaled by a fixed value that should be a factor 
+    // of the number of detectors per bucket, such that the intersections between buckets are handled correctly
+    auto scale_factor = (double)proj_data_out_info.get_num_tangential_poss() / (double)proj_data_in_info.get_num_tangential_poss();
+
+    if (proj_data_in_info.get_scanner_sptr()->get_num_detectors_per_ring() * (int)scale_factor != 
+        proj_data_out_info.get_scanner_sptr()->get_num_detectors_per_ring())
     {
-      error("Interpolation of BlocksOnCylindrical scanners assumes that number of tangential positions "
-            "is the same in the downsampled scanner.");
+      warning("The scanner downsampling is not done by an integer factor. This can produce issues in the interpolation due to the block geometry.");
     }
-    if (proj_data_in_info.get_num_views() != proj_data_out_info.get_num_views())
+    if (proj_data_in_info.get_scanner_sptr()->get_num_transaxial_crystals_per_bucket() * (int)scale_factor != 
+        proj_data_out_info.get_scanner_sptr()->get_num_transaxial_crystals_per_bucket())
     {
-      error("Interpolation of BlocksOnCylindrical scanners assumes that number of views "
-            "is the same in the downsampled scanner.");
+      warning("The scanner downsampling per bucket is not done by an integer factor. "
+              "This can produce issues in the interpolation due to the block geometry");
     }
 
     // only extending in axial direction - an extension of 2 was found to be sufficient
@@ -275,7 +279,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](const BasicCoordinate<3, int>& index_out) 
+    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
@@ -285,8 +289,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

@@ -172,6 +172,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);
 
@@ -385,6 +387,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

@@ -73,4 +73,12 @@ 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

@@ -78,6 +78,7 @@ 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 = "";
@@ -146,6 +147,8 @@ initialise_keymap()
                          &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);
@@ -605,7 +608,7 @@ 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

@@ -909,39 +909,40 @@ 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")
     {
+      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();
-        approx_num_non_arccorrected_bins = this->proj_data_info_sptr->get_num_tangential_poss();
-        // 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_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 new_ring_spacing=scanner_length_block/new_scanner_sptr->get_num_rings();
-//        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_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());
+      }
+
+      // 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_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 new_ring_spacing=scanner_length_block/new_scanner_sptr->get_num_rings();
+      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_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_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{
         if (new_num_dets <= 0)