Refactor build_crystal_maps and fix some issues with initial x position

src/buildblock/GeometryBlocksOnCylindrical.cxx

@@ -35,6 +35,7 @@ limitations under the License.
 #include "stir/Array.h"
 #include "stir/make_array.h"
 #include "stir/numerics/MatrixFunction.h"
+#include "stir/warning.h"
 #include <map>
 #include <iostream>
 #include <fstream>
@@ -63,86 +64,137 @@ 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_axial_buckets = scanner.get_num_axial_buckets();
+  int num_axial_blocks_per_bucket = scanner.get_num_axial_blocks_per_bucket();
   int num_transaxial_crystals_per_block = scanner.get_num_transaxial_crystals_per_block();
   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_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();
-  float axial_crystal_spacing = scanner.get_axial_crystal_spacing();
-  float transaxial_crystal_spacing = scanner.get_transaxial_crystal_spacing();
 
   det_pos_to_coord_type cartesian_coord_map_given_detection_position_keys;
-  /*Building starts from a bucket perpendicular to y axis, from its first crystal.
-          see start_x*/
-
-  // calculate start_point to build the map.
-
-  //    estimate the angle covered by half bucket, csi
-  float csi = _PI / num_transaxial_buckets;
-  float trans_blocks_gap = transaxial_block_spacing - num_transaxial_crystals_per_block * transaxial_crystal_spacing;
-  float ax_blocks_gap = axial_block_spacing - num_axial_crystals_per_block * axial_crystal_spacing;
-  float csi_minus_csiGaps = csi - (csi / transaxial_block_spacing * 2) * (transaxial_crystal_spacing / 2 + trans_blocks_gap);
-  //    distance between the center of the scannner and the first crystal in the bucket, r=Reffective/cos(csi)
-  float r = scanner.get_effective_ring_radius() / cos(csi_minus_csiGaps);
-
-  float start_z = -(axial_block_spacing * (num_axial_blocks_per_bucket)*num_axial_buckets - axial_crystal_spacing
-                    - 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
-
-  stir::CartesianCoordinate3D<float> start_point(start_z, start_y, start_x);
 
+  float csi_minus_csiGaps = get_csi_minus_csi_gaps(scanner);
+
+  // calculate first_crystal_offset to build the map
+  stir::CartesianCoordinate3D<float> first_crystal_offset(get_initial_axial_z_offset(scanner),
+                                                          -scanner.get_effective_ring_radius(),
+                                                          get_initial_axial_x_offset_for_each_bucket(scanner));
+  int radial_coord = 0;
+
+  // Lopp over axial geometry
   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_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)
-            for (int trans_crys_num = 0; trans_crys_num < num_transaxial_crystals_per_block; ++trans_crys_num)
-              {
-                // calculate detection position for a given detector
-                // note: in STIR convention, crystal(0,0,0) corresponds to card_coord(z=0,y=-r,x=0)
-                int tangential_coord;
-                tangential_coord = trans_bucket_num * num_transaxial_blocks_per_bucket * num_transaxial_crystals_per_block
-                                   + trans_block_num * num_transaxial_crystals_per_block + trans_crys_num;
-
-                if (tangential_coord < 0)
-                  tangential_coord += num_detectors_per_ring;
-
-                int axial_coord = ax_bucket_num * num_axial_blocks_per_bucket * num_axial_crystals_per_block
-                                  + ax_block_num * num_axial_crystals_per_block + ax_crys_num;
-                int radial_coord = 0;
-                stir::DetectionPosition<> det_pos(tangential_coord, axial_coord, radial_coord);
-
-                // calculate cartesian coordinate for a given detector
-                stir::CartesianCoordinate3D<float> transformation_matrix(
-                    ax_block_num * axial_block_spacing + ax_crys_num * axial_crystal_spacing,
-                    0.,
-                    trans_block_num * transaxial_block_spacing + trans_crys_num * transaxial_crystal_spacing);
-                float alpha = scanner.get_intrinsic_azimuthal_tilt() + trans_bucket_num * (2 * _PI) / num_transaxial_buckets
-                              + csi_minus_csiGaps;
-
-                stir::Array<2, float> rotation_matrix = get_rotation_matrix(alpha);
-                // to match index range of CartesianCoordinate3D, which is 1 to 3
-                rotation_matrix.set_min_index(1);
-                rotation_matrix[1].set_min_index(1);
-                rotation_matrix[2].set_min_index(1);
-                rotation_matrix[3].set_min_index(1);
-
-                stir::CartesianCoordinate3D<float> transformed_coord = start_point + transformation_matrix;
-                stir::CartesianCoordinate3D<float> cart_coord = stir::matrix_multiply(rotation_matrix, transformed_coord);
-
-                cartesian_coord_map_given_detection_position_keys[det_pos] = cart_coord;
-              }
+        {
+          int axial_coord = get_axial_coord(scanner, ax_bucket_num, ax_block_num, ax_crys_num);
+          float axial_translation = get_axial_translation(scanner, ax_bucket_num, ax_block_num, ax_crys_num);
+
+          // Loop over transaxial geometry
+          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)
+              for (int trans_crys_num = 0; trans_crys_num < num_transaxial_crystals_per_block; ++trans_crys_num)
+                {
+                  // calculate detection position for a given detector
+                  // note: in STIR convention, crystal(0,0,0) corresponds to card_coord(z=0,y=-r,x=0)
+                  int transaxial_coord = get_transaxial_coord(scanner, trans_bucket_num, trans_block_num, trans_crys_num);
+                  stir::DetectionPosition<> det_pos(transaxial_coord, axial_coord, radial_coord);
+
+                  // The translation matrix from the first crystal in the block
+                  stir::CartesianCoordinate3D<float> translation_matrix(
+                      axial_translation,
+                      0.,
+                      get_crystal_in_bucket_transaxial_translation(scanner, trans_block_num, trans_crys_num));
+
+                  stir::CartesianCoordinate3D<float> transformed_coord = first_crystal_offset + translation_matrix;
+
+                  // Calculate the rotation of the crystal
+                  float alpha = scanner.get_intrinsic_azimuthal_tilt() + trans_bucket_num * (2 * _PI) / num_transaxial_buckets
+                                + csi_minus_csiGaps;
+                  cartesian_coord_map_given_detection_position_keys[det_pos]
+                      = calculate_crystal_rotation(transformed_coord, alpha);
+                }
+        }
   set_detector_map(cartesian_coord_map_given_detection_position_keys);
 }
 
+CartesianCoordinate3D<float>
+GeometryBlocksOnCylindrical::calculate_crystal_rotation(const CartesianCoordinate3D<float>& crystal_position,
+                                                        const float alpha) const
+{
+  stir::Array<2, float> rotation_matrix = get_rotation_matrix(alpha);
+  // to match index range of CartesianCoordinate3D, which is 1 to 3
+  rotation_matrix.set_min_index(1);
+  rotation_matrix[1].set_min_index(1);
+  rotation_matrix[2].set_min_index(1);
+  rotation_matrix[3].set_min_index(1);
+  return stir::matrix_multiply(rotation_matrix, crystal_position);
+}
+
+int
+GeometryBlocksOnCylindrical::get_transaxial_coord(const Scanner& scanner,
+                                                  int transaxial_bucket_num,
+                                                  int transaxial_block_num,
+                                                  int transaxial_crystal_num)
+{
+  return transaxial_bucket_num * scanner.get_num_transaxial_blocks_per_bucket() * scanner.get_num_transaxial_crystals_per_block()
+         + transaxial_block_num * scanner.get_num_transaxial_crystals_per_block() + transaxial_crystal_num;
+}
+
+int
+GeometryBlocksOnCylindrical::get_axial_coord(const Scanner& scanner,
+                                             int axial_bucket_num,
+                                             int axial_block_num,
+                                             int axial_crystal_num)
+{
+  return axial_bucket_num * scanner.get_num_axial_blocks_per_bucket() * scanner.get_num_axial_crystals_per_block()
+         + axial_block_num * scanner.get_num_axial_crystals_per_block() + axial_crystal_num;
+}
+
+float
+GeometryBlocksOnCylindrical::get_crystal_in_bucket_transaxial_translation(const Scanner& scanner,
+                                                                          int transaxial_block_num,
+                                                                          int transaxial_crystal_num)
+{
+  // Currently, only supports 1 transaxial bucket per angle
+  return transaxial_block_num * scanner.get_transaxial_block_spacing()
+         + transaxial_crystal_num * scanner.get_transaxial_crystal_spacing();
+}
+
+float
+GeometryBlocksOnCylindrical::get_axial_translation(const Scanner& scanner,
+                                                   int axial_bucket_num,
+                                                   int axial_block_num,
+                                                   int axial_crystal_num)
+{
+  return // axial_bucket_num * scanner.get_axial_bucket_spacing() +
+      axial_block_num * scanner.get_axial_block_spacing() + axial_crystal_num * scanner.get_axial_crystal_spacing();
+}
+
+float
+GeometryBlocksOnCylindrical::get_initial_axial_z_offset(const Scanner& scanner)
+{
+  // Crystals in a block are centered, blocks in a bucket are centered, and buckets are centered in the z axis.
+  // This centers the scanner in z
+  float crystals_in_block_offset = (scanner.get_num_axial_crystals_per_block() - 1) * scanner.get_axial_crystal_spacing();
+  float blocks_in_bucket_offset = (scanner.get_num_axial_blocks_per_bucket() - 1) * scanner.get_axial_block_spacing();
+  //  float bucket_offset = (scanner.get_num_axial_buckets() - 1) * scanner.get_axial_bucket_spacing();
+  float bucket_offset = 0;
+
+  // Negative because the scanner is centered at z=0 and increases axial coordinates increase
+  // 1/2 because it is half the distance from the center to the edge of the scanner
+  return -(1.0 / 2) * (crystals_in_block_offset + blocks_in_bucket_offset + bucket_offset);
+}
+
+float
+GeometryBlocksOnCylindrical::get_initial_axial_x_offset_for_each_bucket(const Scanner& scanner)
+{
+  // This is the old method... This is probably wrong
+  //  float csi_minus_csiGaps = get_csi_minus_csi_gaps(scanner);
+  //  float r = scanner.get_effective_ring_radius() / cos(csi_minus_csiGaps);
+  //  return -1 * r * sin(csi_minus_csiGaps);
+
+  auto first_crystal_coord = get_crystal_in_bucket_transaxial_translation(scanner, 0, 0);
+  auto last_crystal_coord = get_crystal_in_bucket_transaxial_translation(
+      scanner, scanner.get_num_transaxial_blocks_per_bucket() - 1, scanner.get_num_transaxial_crystals_per_block() - 1);
+  return -(1.0 / 2) * (first_crystal_coord + last_crystal_coord);
+}
 END_NAMESPACE_STIR

src/include/stir/GeometryBlocksOnCylindrical.h

@@ -52,6 +52,45 @@ class GeometryBlocksOnCylindrical : public DetectorCoordinateMap
 public:
   GeometryBlocksOnCylindrical(const Scanner& scanner);
 
+  //! Calculates the transaxial coordinate of a crystal given a scanner and the crystal's indices
+  static int
+  get_transaxial_coord(const Scanner& scanner, int transaxial_bucket_num, int transaxial_block_num, int transaxial_crystal_num);
+
+  //! Calculates the axial coordinate of a crystal given a scanner and the crystal's indices
+  static int get_axial_coord(const Scanner& scanner, int axial_bucket_num, int axial_block_num, int axial_crystal_num);
+
+  //! Calculates the transaxial translation of a crystal given a scanner and the crystal's indices
+  static float
+  get_crystal_in_bucket_transaxial_translation(const Scanner& scanner, int transaxial_block_num, int transaxial_crystal_num);
+
+  //! Calculates the axial translation of a crystal given a scanner and the crystal's indices
+  static float get_axial_translation(const Scanner& scanner, int axial_bucket_num, int axial_block_num, int axial_crystal_num);
+
+  //! Calculate the initial axial z offset to center the scanner on 0,0,0
+  static float get_initial_axial_z_offset(const Scanner& scanner);
+
+  //! Calculate the initial transaxial x offset to center the scanner on 0,0,0
+  static float get_initial_axial_x_offset_for_each_bucket(const Scanner& scanner);
+
+  static float get_csi_minus_csi_gaps(const Scanner& scanner)
+  {
+    //! Calculate the CSI, angle covered by half a bucket
+    // 2 * PI / num_transaxial_buckets / 2 (simplified)
+    float csi = _PI / scanner.get_num_transaxial_buckets(); // TODO, this assumes 1 transaxial bucket per angle
+
+    // The difference between the transaxial block spacing and the sum of all transaxial crystal spacing's in the block
+    float trans_blocks_gap = scanner.get_transaxial_block_spacing()
+                             - scanner.get_num_transaxial_crystals_per_block() * scanner.get_transaxial_crystal_spacing();
+    // Calculate the angle covered by the gaps between the blocks
+    float csi_gaps
+        = 2 * csi * (scanner.get_transaxial_crystal_spacing() / 2 + trans_blocks_gap) / scanner.get_transaxial_block_spacing();
+    return csi - csi_gaps;
+  };
+
+  //!
+  CartesianCoordinate3D<float> calculate_crystal_rotation(const CartesianCoordinate3D<float>& crystal_position,
+                                                          const float alpha) const;
+
 private:
   //! Get rotation matrix for a given angle around z axis
   stir::Array<2, float> get_rotation_matrix(float alpha) const;