Axial bucket projection issue

src/buildblock/GeometryBlocksOnCylindrical.cxx

@@ -44,6 +44,7 @@ START_NAMESPACE_STIR
 
 GeometryBlocksOnCylindrical::GeometryBlocksOnCylindrical(const Scanner& scanner)
 {
+  check_scanner_configuration(scanner);
   build_crystal_maps(scanner);
 }
 
@@ -55,17 +56,32 @@ GeometryBlocksOnCylindrical::get_rotation_matrix(float alpha) const
                           stir::make_1d_array(0.F, -1 * std::sin(alpha), std::cos(alpha)));
 }
 
+void
+GeometryBlocksOnCylindrical::check_scanner_configuration(const Scanner& scanner)
+{
+  { // Ensure the number of axial crystals per bucket is a multiple of the number of rings
+    int num_axial_crystals_per_bucket = scanner.get_num_axial_crystals_per_block() * scanner.get_num_axial_blocks_per_bucket();
+    if (scanner.get_num_rings() % (num_axial_crystals_per_bucket) != 0)
+      {
+        error("Error in GeometryBlocksOnCylindrical: number of rings (%d) is not a multiple of the num_axial_crystals_per_bucket "
+              "(%d) = num_axial_crystals_per_block (%d) * num_axial_blocks_per_bucket (%d)",
+              scanner.get_num_rings(),
+              num_axial_crystals_per_bucket,
+              scanner.get_num_axial_crystals_per_block(),
+              scanner.get_num_axial_blocks_per_bucket());
+      }
+  }
+}
 void
 GeometryBlocksOnCylindrical::build_crystal_maps(const Scanner& scanner)
 {
   // local variables to describe scanner
   int num_axial_crystals_per_block = scanner.get_num_axial_crystals_per_block();
   int num_transaxial_crystals_per_block = scanner.get_num_transaxial_crystals_per_block();
-  int num_axial_blocks = scanner.get_num_axial_blocks();
   int num_transaxial_blocks_per_bucket = scanner.get_num_transaxial_blocks_per_bucket();
   int num_axial_blocks_per_bucket = scanner.get_num_axial_blocks_per_bucket();
-  int num_transaxial_buckets = scanner.get_num_transaxial_blocks() / num_transaxial_blocks_per_bucket;
-  int num_axial_buckets = scanner.get_num_axial_blocks() / num_axial_blocks_per_bucket;
+  int num_transaxial_buckets = scanner.get_num_transaxial_buckets();
+  int num_axial_buckets = scanner.get_num_axial_buckets();
   int num_detectors_per_ring = scanner.get_num_detectors_per_ring();
   float axial_block_spacing = scanner.get_axial_block_spacing();
   float transaxial_block_spacing = scanner.get_transaxial_block_spacing();
@@ -100,7 +116,7 @@ GeometryBlocksOnCylindrical::build_crystal_maps(const Scanner& scanner)
   stir::CartesianCoordinate3D<float> start_point(start_z, start_y, start_x);
 
   for (int ax_bucket_num = 0; ax_bucket_num < num_axial_buckets; ++ax_bucket_num)
-    for (int ax_block_num = 0; ax_block_num < num_axial_blocks; ++ax_block_num)
+    for (int ax_block_num = 0; ax_block_num < num_axial_blocks_per_bucket; ++ax_block_num)
       for (int ax_crys_num = 0; ax_crys_num < num_axial_crystals_per_block; ++ax_crys_num)
         for (int trans_bucket_num = 0; trans_bucket_num < num_transaxial_buckets; ++trans_bucket_num)
           for (int trans_block_num = 0; trans_block_num < num_transaxial_blocks_per_bucket; ++trans_block_num)

src/include/stir/GeometryBlocksOnCylindrical.h

@@ -56,6 +56,9 @@ class GeometryBlocksOnCylindrical : public DetectorCoordinateMap
   //! Get rotation matrix for a given angle around z axis
   stir::Array<2, float> get_rotation_matrix(float alpha) const;
 
+  //! Checks the scanner configuration is valid for a crystal map
+  static void check_scanner_configuration(const Scanner& scanner) ;
+
   //! Build crystal map in cartesian coordinate
   void build_crystal_maps(const Scanner& scanner);
 };

src/recon_test/test_blocks_on_cylindrical_projectors.cxx

@@ -53,6 +53,7 @@
 #include "stir/recon_buildblock/BackProjectorByBinUsingProjMatrixByBin.h"
 #include "stir/IO/write_to_file.h"
 #include "stir/VoxelsOnCartesianGrid.h"
+#include "stir/PixelsOnCartesianGrid.h"
 //#include "stir/Shape/Shape3D.h"
 #include <cmath>
 
@@ -80,6 +81,7 @@ class BlocksTests : public RunTests
   void run_map_orientation_test(ForwardProjectorByBin& forw_projector1, ForwardProjectorByBin& forw_projector2);
   void run_projection_test(ForwardProjectorByBin& forw_projector1, ForwardProjectorByBin& forw_projector2);
   void run_intersection_with_cylinder_test();
+  void run_back_projection_test_with_axial_buckets();
 };
 
 /*! The following is a function to allow a projdata_info BlocksOnCylindrical to be created from the scanner. */
@@ -871,6 +873,90 @@ BlocksTests::run_intersection_with_cylinder_test()
     }
 }
 
+void
+BlocksTests::run_back_projection_test_with_axial_buckets()
+
+{
+  //    create projadata info
+  auto scanner_sptr = std::make_shared<Scanner>(Scanner::SAFIRDualRingPrototype);
+  {
+    scanner_sptr->set_average_depth_of_interaction(5);
+    scanner_sptr->set_num_axial_crystals_per_block(1);
+    scanner_sptr->set_axial_block_spacing(scanner_sptr->get_axial_crystal_spacing()
+                                          * scanner_sptr->get_num_axial_crystals_per_block());
+    scanner_sptr->set_transaxial_block_spacing(scanner_sptr->get_transaxial_crystal_spacing()
+                                               * scanner_sptr->get_num_transaxial_crystals_per_block());
+    scanner_sptr->set_num_axial_blocks_per_bucket(2);
+    scanner_sptr->set_num_rings(2 * 5);  // More than 1 bucket
+  }
+  auto segments = scanner_sptr->get_num_rings() - 1;
+  VectorWithOffset<int> num_axial_pos_per_segment(-segments, segments);
+  VectorWithOffset<int> min_ring_diff_v(-segments, segments);
+  VectorWithOffset<int> max_ring_diff_v(-segments, segments);
+  for (int i = -segments; i <= segments; i++)
+    {
+      min_ring_diff_v[i] = i;
+      max_ring_diff_v[i] = i;
+      num_axial_pos_per_segment[i] = scanner_sptr->get_num_rings() - abs(i);
+    }
+
+  scanner_sptr->set_scanner_geometry("BlocksOnCylindrical");
+  scanner_sptr->set_up();
+  auto projdata_info_sptr
+      = std::make_shared<ProjDataInfoBlocksOnCylindricalNoArcCorr>(scanner_sptr,
+                                                                   num_axial_pos_per_segment,
+                                                                   min_ring_diff_v,
+                                                                   max_ring_diff_v,
+                                                                   scanner_sptr->get_max_num_views(),
+                                                                   scanner_sptr->get_max_num_non_arccorrected_bins());
+  //  scanner_sptr->set_scanner_geometry("Cylindrical");
+  //  scanner_sptr->set_up();
+  //  auto projdata_info_sptr = std::make_shared<ProjDataInfoCylindricalNoArcCorr>(scanner_sptr,
+  //                                                                               num_axial_pos_per_segment,
+  //                                                                               min_ring_diff_v,
+  //                                                                               max_ring_diff_v,
+  //                                                                               scanner_sptr->get_max_num_views(),
+  //                                                                               scanner_sptr->get_max_num_non_arccorrected_bins());
+
+  auto projdata_sptr = std::make_shared<ProjDataInMemory>(std::make_shared<ExamInfo>(ImagingModality::PT), projdata_info_sptr);
+
+  auto origin = CartesianCoordinate3D<float>(0, 0, 0);
+  auto volume_dimensions = CartesianCoordinate3D<int>(-1, -1, -1);
+  auto volume = std::make_shared<VoxelsOnCartesianGrid<float>>(*projdata_info_sptr, 1, origin, volume_dimensions);
+
+  // Now run the test
+  volume->fill(0.0);
+  projdata_sptr->fill(1.0);
+
+  auto PM = std::make_shared<ProjMatrixByBinUsingRayTracing>();
+  PM->enable_cache(false);
+  auto back_projector = std::make_shared<BackProjectorByBinUsingProjMatrixByBin>(PM);
+  back_projector->set_up(projdata_info_sptr, volume);
+  back_projector->back_project(*volume, *projdata_sptr, 0, 144);
+
+  bool test_ok = true;
+  std::ostringstream oss;
+  oss << "BlocksTests::run_back_projection_test_with_axial_buckets: Central voxel values are:" << std::endl;
+
+  auto center_axial_values = std::vector<float>(volume->get_z_size());
+  for (int z = volume->get_min_z(); z <= volume->get_max_z(); z++)
+    {
+      center_axial_values[z] = volume->get_plane(z).at(0).at(0);
+      oss << "\tz = " << z << "/" << volume->get_max_z() << " is " << center_axial_values[z] << std::endl;
+      if (test_ok)
+        {
+          test_ok = check(center_axial_values[z] > 0,
+                          "BlocksTests::run_back_projection_test_with_axial_buckets: Central voxel value <= 0");
+          everything_ok = everything_ok && test_ok;
+        }
+    }
+
+  if (test_ok)
+    return;
+  // Something went wrong, log the error
+  std::cerr << oss.str();
+}
+
 void
 BlocksTests::run_tests()
 {
@@ -905,6 +991,8 @@ BlocksTests::run_tests()
   print_time("map orientation test took: ");
   run_intersection_with_cylinder_test();
   print_time("intersection with cylinder test took: ");
+  run_back_projection_test_with_axial_buckets();
+  print_time("back projection test with axial buckets took: ");
 
 #ifdef STIR_WITH_Parallelproj_PROJECTOR
   // run the same tests with parallelproj, if available