Reverting a change that broke the blocks on cylindrical geometry. (#1462)

* Reverting a change that broke the blocks on cylindrical geometry.

* Fix for the blocks on cylindrical test by making the bin search more reliable.

* Fixing the other test by simplifying and fixing the calculation of csi using the same logic than in plot_scanner_LORs.py

---------

Co-authored-by: Markus Jehl <markus.jehl@positrigo.com>

Author: markus-jehl

documentation/release_6.2.htm

@@ -73,6 +73,12 @@ <h3>Changed functionality</h3>
 
 <h3>Bug fixes</h3>
 <ul>
+  <li>
+    There was a bug in the computation of the detector coordinates for <code>BlocksOnCylindrical</code>
+    scanners that caused the buckets to not be symmetric. <br>
+    <a href=https://github.com/UCL/STIR/pull/1462>PR #1462</a>
+  </li>
+
 </ul>
 
 <h3>Build system</h3>

src/buildblock/GeometryBlocksOnCylindrical.cxx

@@ -92,12 +92,9 @@ GeometryBlocksOnCylindrical::build_crystal_maps(const Scanner& scanner)
                     - ax_blocks_gap * (num_axial_blocks_per_bucket * num_axial_buckets - 1))
                   / 2;
   float start_y = -1 * scanner.get_effective_ring_radius();
-  float start_x
-      = -1 * r
-        * sin(csi_minus_csiGaps); //(
-                                  //								((num_transaxial_blocks_per_bucket-1)/2.)*transaxial_block_spacing
-  //							+ ((num_transaxial_crystals_per_block-1)/2.)*transaxial_crystal_spacing
-  //											); //the first crystal in the bucket
+  float start_x = -1 // the first crystal in the bucket
+                  * (((num_transaxial_blocks_per_bucket - 1) / 2.) * transaxial_block_spacing
+                     + ((num_transaxial_crystals_per_block - 1) / 2.) * transaxial_crystal_spacing);
 
   stir::CartesianCoordinate3D<float> start_point(start_z, start_y, start_x);
 

src/recon_test/test_blocks_on_cylindrical_projectors.cxx

@@ -364,27 +364,29 @@ BlocksTests::run_plane_symmetry_test(ForwardProjectorByBin& forw_projector1, For
   forw_projector2.forward_project(*projdata2, *image2_sptr);
 
   int view1_num = 0, view2_num = 0;
+  float min_diff = std::numeric_limits<float>::max();
   LORInAxialAndNoArcCorrSinogramCoordinates<float> lorB1;
   for (int i = 0; i < projdata->get_max_view_num(); i++)
     {
       Bin bin(0, i, 0, 0);
       proj_data_info_blocks_sptr->get_LOR(lorB1, bin);
-      if (std::abs(lorB1.phi() - phi1) / phi1 <= 1E-2)
+      if (std::abs(lorB1.phi() - phi1) < min_diff)
         {
+          min_diff = std::abs(lorB1.phi() - phi1);
           view1_num = i;
-          break;
         }
     }
 
   LORInAxialAndNoArcCorrSinogramCoordinates<float> lorB2;
+  min_diff = std::numeric_limits<float>::max();
   for (int i = 0; i < projdata2->get_max_view_num(); i++)
     {
       Bin bin(0, i, 0, 0);
       proj_data_info_blocks_sptr->get_LOR(lorB2, bin);
-      if (std::abs(lorB2.phi() - phi2) / phi2 <= 1E-2)
+      if (std::abs(lorB2.phi() - phi2) < min_diff)
         {
+          min_diff = std::abs(lorB2.phi() - phi2);
           view2_num = i;
-          break;
         }
     }
 
@@ -393,23 +395,25 @@ BlocksTests::run_plane_symmetry_test(ForwardProjectorByBin& forw_projector1, For
 
   //    find the tang position with the max value
   int tang1_num = 0, tang2_num = 0;
+  min_diff = std::numeric_limits<float>::max();
   for (int tang = projdata->get_min_tangential_pos_num(); tang < projdata->get_max_tangential_pos_num(); tang++)
     {
 
-      if ((max1 - projdata->get_sinogram(0, 0).at(view1_num).at(tang)) / max1 < 1E-3)
+      if ((max1 - projdata->get_sinogram(0, 0).at(view1_num).at(tang)) < min_diff)
         {
+          min_diff = (max1 - projdata->get_sinogram(0, 0).at(view1_num).at(tang));
           tang1_num = tang;
-          break;
         }
     }
 
+  min_diff = std::numeric_limits<float>::max();
   for (int tang = projdata2->get_min_tangential_pos_num(); tang < projdata2->get_max_tangential_pos_num(); tang++)
     {
 
-      if ((max2 - projdata2->get_sinogram(0, 0).at(view2_num).at(tang)) / max2 < 1E-3)
+      if ((max2 - projdata2->get_sinogram(0, 0).at(view2_num).at(tang)) < min_diff)
         {
+          min_diff = (max1 - projdata->get_sinogram(0, 0).at(view1_num).at(tang));
           tang2_num = tang;
-          break;
         }
     }
 

src/test/test_DetectorCoordinateMap.cxx

@@ -67,33 +67,13 @@ float
 DetectionPosMapTests::calculate_angle_within_half_bucket(
     const shared_ptr<Scanner> scanner_ptr, const shared_ptr<ProjDataInfoBlocksOnCylindricalNoArcCorr> proj_data_info_ptr)
 {
-  Bin bin;
-  LORInAxialAndNoArcCorrSinogramCoordinates<float> lorB;
-  float csi;
-  float C_spacing = scanner_ptr->get_transaxial_crystal_spacing();
-  float csi_crystal = std::atan((C_spacing) / scanner_ptr->get_effective_ring_radius());
-  //    float bucket_spacing=scanner_ptr->get_transaxial_block_spacing()*C_spacing;
-  //    float blocks_gap=scanner_ptr->get_transaxial_block_spacing()
-  //            -scanner_ptr->get_num_transaxial_crystals_per_block()*C_spacing;
-  //    float csi_gap=std::atan((blocks_gap)/scanner_ptr->get_effective_ring_radius());
-
-  //    get angle within half bucket
-  for (int view = 0; view <= scanner_ptr->get_max_num_views(); view++)
-    {
-      int bucket_num
-          = view / (scanner_ptr->get_num_transaxial_crystals_per_block() * scanner_ptr->get_num_transaxial_blocks_per_bucket());
-      if (bucket_num > 0)
-        break;
-
-      bin.segment_num() = 0;
-      bin.axial_pos_num() = 0;
-      bin.view_num() = view;
-      bin.tangential_pos_num() = 0;
-
-      proj_data_info_ptr->get_LOR(lorB, bin);
-      csi = lorB.phi();
-    }
-  return (csi + csi_crystal) / 2;
+  auto csi = _PI / (scanner_ptr->get_num_transaxial_blocks() / scanner_ptr->get_num_transaxial_blocks_per_bucket());
+  auto tBl_gap = scanner_ptr->get_transaxial_block_spacing()
+                 - scanner_ptr->get_num_transaxial_crystals_per_block() * scanner_ptr->get_transaxial_crystal_spacing();
+  auto csi_minus_csiGaps
+      = csi
+        - (csi / scanner_ptr->get_transaxial_block_spacing() * 2) * (scanner_ptr->get_transaxial_crystal_spacing() / 2 + tBl_gap);
+  return csi_minus_csiGaps;
 }
 
 /*!