UCL · ashgillman · Sep 19, 2018 · Sep 3, 2018 · Jul 21, 2020 · Jul 29, 2020
diff --git a/src/analytic/FBP3DRP/FBP3DRPReconstruction.cxx b/src/analytic/FBP3DRP/FBP3DRPReconstruction.cxx
@@ -5,11 +5,13 @@
   \author Kris Thielemans
   \author Claire LABBE
   \author PARAPET project
+  \author Ashley Gillman
 */
 /*
     Copyright (C) 2000 PARAPET partners
     Copyright (C) 2000- 2012, Hammersmith Imanet Ltd
     Copyright (C) 2020, University College London
+    Copyright (C) 2020, CSIRO
     This file is part of STIR.
 
     SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license
@@ -19,6 +21,10 @@
 
 /*
  Modification history: (highlights in anti-chronological order)
+ AG 2020/07/21
+ - We no longer want to force images to be centre-aligned, rewrote implementation
+   of find_rmin_rmax().
+
  KT Oct 2004
  - no longer use Numerical Recipes fourier
  - option to 'stretch' the colsher filter during definition for better results
@@ -83,8 +89,9 @@
 
 #include "stir/analytic/FBP3DRP/ColsherFilter.h"
 #include "stir/display.h"
-// #include "stir/recon_buildblock/distributable.h"
-// #include "stir/FBP3DRP/process_viewgrams.h"
+//#include "stir/info.h"
+//#include "stir/recon_buildblock/distributable.h"
+//#include "stir/FBP3DRP/process_viewgrams.h"
 
 #include "stir/analytic/FBP3DRP/FBP3DRPReconstruction.h"
 #include "stir/analytic/FBP2D/FBP2DReconstruction.h"
@@ -124,50 +131,47 @@ find_rmin_rmax(int& rmin,
                const int seg_num,
                const VoxelsOnCartesianGrid<float>& image)
 {
-
+  // precomupute a few values:
+  // Radius of the FOV as per the ProjDataInfo, in mm
   const float fovrad = proj_data_info_cyl.get_s(Bin(0, 0, 0, proj_data_info_cyl.get_num_tangential_poss() / 2 - 1));
-  // Compute minimum and maximum rings of 'missing' projections
-
+  // Radius of the detector rings in mm
+  const float ringrad = proj_data_info_cyl.get_ring_radius();
+  // delta, ring difference in units of n.rings
   const float delta = proj_data_info_cyl.get_average_ring_difference(seg_num);
-
-  // find correspondence between ring coordinates and image coordinates:
-  // z = num_planes_per_virtual_ring * ring + virtual_ring_offset
-  // compute the offset by matching up the centre of the scanner
-  // in the 2 coordinate systems
-  // TODO get all this from ProjDataInfo or so
-
-  const int num_planes_per_virtual_ring
-      = (proj_data_info_cyl.get_max_ring_difference(seg_num) == proj_data_info_cyl.get_min_ring_difference(seg_num)) ? 2 : 1;
+  // if span > 1, we essentially have twice as many "virtual rings"
   const int num_virtual_rings_per_physical_ring
       = (proj_data_info_cyl.get_max_ring_difference(seg_num) == proj_data_info_cyl.get_min_ring_difference(seg_num)) ? 1 : 2;
 
-  const float virtual_ring_offset
-      = (image.get_max_z() + image.get_min_z()) / 2.F
-        - num_planes_per_virtual_ring
-              * (proj_data_info_cyl.get_max_axial_pos_num(seg_num) + num_virtual_rings_per_physical_ring * delta
-                 + proj_data_info_cyl.get_min_axial_pos_num(seg_num))
-              / 2;
-
-  // we first consider the LOR at s=0, phi=0 which goes through z=0,y=0, x=fovrad
-  // later on, we will shift it to the 'left'most edge of the FOV.
-
-  // find z position of intersection of this LOR with the detector radius
-  // (i.e. y=0, x=-ring_radius)
-  // use image coordinates first
-  float z_in_image_coordinates = -fabs(delta) * num_planes_per_virtual_ring * num_virtual_rings_per_physical_ring
-                                 * (fovrad + proj_data_info_cyl.get_ring_radius()) / (2 * proj_data_info_cyl.get_ring_radius());
-  // now shift it to the edge of the FOV
-  // (taking into account that z==get_min_z() is in the middle of the voxel)
-  z_in_image_coordinates += image.get_min_z() - .5F;
-
-  // now convert to virtual_ring_coordinates using z = num_planes_per_virtual_ring * ring + virtual_ring_offset
-  const float z_in_virtual_ring_coordinates = (z_in_image_coordinates - virtual_ring_offset) / num_planes_per_virtual_ring;
+  // This can be derived graphically by drawing a line from a top view of
+  // a scanner from ring to ring at angle theta (based on the segment)
+  // that just skims the field-of-view at the end of the ring.
+  // (Assume phi = y = 0, so no out-of-plane elements).
+  // The width of this line is delta (in units of n.rings), and height
+  // ringrad, and this triangle is similar to the smaller portion that
+  // still lies within the rings, which has width (delta - z) (where z
+  // is the overhang from the rings in units of n.rings) and height
+  // (ringrad - fovrad) / 2
+  const float lor_overhang_in_num_rings = fabs(delta) / 2 * (1 + fovrad / ringrad);
+  const float lor_overhang_in_num_virtual_rings = lor_overhang_in_num_rings * num_virtual_rings_per_physical_ring;
+
+  // rmin should just be this number of virtual rings backward from 0.
+  rmin = static_cast<int>(floor(-lor_overhang_in_num_virtual_rings));
+  // and rmax can be calculated symmetrically (we assume here the axial positions
+  // are correctly centred, this won't work otherwise.)
+  rmax = proj_data_info_cyl.get_max_axial_pos_num(seg_num) + (proj_data_info_cyl.get_min_axial_pos_num(seg_num) - rmin);
 
-  // finally find the 'ring' number
-  rmin = static_cast<int>(floor(z_in_virtual_ring_coordinates));
+  // lastly, if we were, for some reason, provided a larger-than-standard
+  // sinogram, make sure we use at least that size.
+  if (proj_data_info_cyl.get_min_axial_pos_num(seg_num) < rmin)
+    {
+      rmin = proj_data_info_cyl.get_min_axial_pos_num(seg_num);
+      rmax = proj_data_info_cyl.get_max_axial_pos_num(seg_num);
+    }
 
-  // rmax is determined by using symmetry: at both ends there are just as many missing rings
-  rmax = proj_data_info_cyl.get_max_axial_pos_num(seg_num) + (proj_data_info_cyl.get_min_axial_pos_num(seg_num) - rmin);
+  // info(boost::format("seg: %s : amin: %s, amax: %s")
+  //   % seg_num % proj_data_info_cyl.get_min_axial_pos_num(seg_num) % proj_data_info_cyl.get_max_axial_pos_num(seg_num));
+  // info(boost::format("seg: %s : rmin: %s, rmax: %s")
+  //   % seg_num % rmin % rmax);
 }
 
 const char* const FBP3DRPReconstruction::registered_name = "FBP3DRP";

diff --git a/src/buildblock/ProjDataInfo.cxx b/src/buildblock/ProjDataInfo.cxx
@@ -801,4 +801,26 @@ ProjDataInfo::operator>=(const ProjDataInfo& proj_data_info) const
   return (proj_data_info == *smaller_proj_data_info_sptr);
 }
 
+CartesianCoordinate3D<float>
+ProjDataInfo::get_point_on_lor_in_gantry_coordinates(
+    const float s_in_mm, const float m_in_mm, const float a_in_mm, const float cphi, const float sphi, const float tantheta) const
+{
+  return CartesianCoordinate3D<float>(m_in_mm - a_in_mm * tantheta,     // z
+                                      s_in_mm * sphi - a_in_mm * cphi,  // y
+                                      s_in_mm * cphi + a_in_mm * sphi); // x
+}
+
+CartesianCoordinate3D<float>
+ProjDataInfo::get_vector_centre_of_first_ring_to_centre_of_gantry() const
+{
+  float middle_of_first_ring_to_middle_of_last = (scanner_ptr->get_num_rings() - 1) * scanner_ptr->get_ring_spacing();
+  return CartesianCoordinate3D<float>(middle_of_first_ring_to_middle_of_last / 2.F, 0, 0);
+}
+
+CartesianCoordinate3D<float>
+ProjDataInfo::get_bed_position() const
+{
+  return CartesianCoordinate3D<float>(bed_position_horizontal, bed_position_vertical, 0);
+}
+
 END_NAMESPACE_STIR
diff --git a/src/buildblock/ProjDataInfoBlocksOnCylindricalNoArcCorr.cxx b/src/buildblock/ProjDataInfoBlocksOnCylindricalNoArcCorr.cxx
@@ -100,8 +100,8 @@ ProjDataInfoBlocksOnCylindricalNoArcCorr::find_scanner_coordinates_given_cartesi
 
   DetectionPosition<> det_pos1;
   DetectionPosition<> det_pos2;
-  if (get_scanner_ptr()->find_detection_position_given_cartesian_coordinate(det_pos1, c1 + this->z_shift) == Succeeded::no
-      || get_scanner_ptr()->find_detection_position_given_cartesian_coordinate(det_pos2, c2 + this->z_shift) == Succeeded::no)
+  if (get_scanner_ptr()->find_detection_position_given_cartesian_coordinate(det_pos1, c1) == Succeeded::no
+      || get_scanner_ptr()->find_detection_position_given_cartesian_coordinate(det_pos2, c2) == Succeeded::no)
     {
       return Succeeded::no;
     }

diff --git a/src/buildblock/ProjDataInfoCylindrical.cxx b/src/buildblock/ProjDataInfoCylindrical.cxx
@@ -147,6 +147,7 @@ ProjDataInfoCylindrical::initialise_ring_diff_arrays() const
     /* m_offsets are found by requiring
     get_m(..., min_axial_pos_num,...) == - get_m(..., max_axial_pos_num,...)
     */
+
     for (int segment_num = get_min_segment_num(); segment_num <= get_max_segment_num(); ++segment_num)
       {
         m_offset[segment_num]
@@ -283,10 +284,7 @@ ProjDataInfoCylindrical::initialise_ring_diff_arrays() const
                                                    + (get_scanner_ptr()->get_num_rings() - 1);
               const int ring1_plus_ring2 = round(ring1_plus_ring2_float);
               // check that it was integer
-              if (get_scanner_sptr()->get_scanner_geometry() == "Cylindrical")
-                {
-                  assert(fabs(ring1_plus_ring2 - ring1_plus_ring2_float) < 1E-4);
-                }
+              assert(fabs(ring1_plus_ring2 - ring1_plus_ring2_float) < 1E-4);
               segment_axial_pos_to_ring1_plus_ring2[s_num][ax_pos_num] = ring1_plus_ring2;
             }
         }

diff --git a/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx b/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx
@@ -357,6 +357,7 @@ ProjDataInfoCylindricalNoArcCorr::get_all_det_pos_pairs_for_bin(vector<Detection
   assert(current_dp_num == get_num_det_pos_pairs_for_bin(bin, ignore_non_spatial_dimensions));
 }
 
+// find_scanner_coordinates_given_cartesian_coordinates
 Succeeded
 ProjDataInfoCylindricalNoArcCorr::find_scanner_coordinates_given_cartesian_coordinates(
     int& det1, int& det2, int& ring1, int& ring2, const CartesianCoordinate3D<float>& c1, const CartesianCoordinate3D<float>& c2)
@@ -399,14 +400,23 @@ ProjDataInfoCylindricalNoArcCorr::find_scanner_coordinates_given_cartesian_coord
   ring1 = round(coord_det1.z()/ring_spacing);
   ring2 = round(coord_det2.z()/ring_spacing);
 #else
-  LORInCylinderCoordinates<float> cyl_coords;
-  if (find_LOR_intersections_with_cylinder(cyl_coords, LORAs2Points<float>(c1, c2), ring_radius) == Succeeded::no)
+
+  // here we define an internal-only coord system called ring coords
+  // They just define z=0 as first ring, for ease of calculating ring no.
+  CartesianCoordinate3D<float> offset_gantry_coords_to_ring_coords = get_vector_centre_of_first_ring_to_centre_of_gantry();
+
+  LORInCylinderCoordinates<float> cyl_in_ring_coords;
+  if (find_LOR_intersections_with_cylinder(
+          cyl_in_ring_coords,
+          LORAs2Points<float>(c1 + offset_gantry_coords_to_ring_coords, c2 + offset_gantry_coords_to_ring_coords),
+          ring_radius)
+      == Succeeded::no)
     return Succeeded::no;
 
-  det1 = modulo(round((cyl_coords.p1().psi() - this->get_psi_offset()) / (2. * _PI / num_detectors)), num_detectors);
-  det2 = modulo(round((cyl_coords.p2().psi() - this->get_psi_offset()) / (2. * _PI / num_detectors)), num_detectors);
-  ring1 = round(cyl_coords.p1().z() / ring_spacing);
-  ring2 = round(cyl_coords.p2().z() / ring_spacing);
+  det1 = modulo(round((cyl_in_ring_coords.p1().psi() - this->get_psi_offset()) / (2. * _PI / num_detectors)), num_detectors);
+  det2 = modulo(round((cyl_in_ring_coords.p2().psi() - this->get_psi_offset()) / (2. * _PI / num_detectors)), num_detectors);
+  ring1 = round((cyl_in_ring_coords.p1()).z() / ring_spacing);
+  ring2 = round((cyl_in_ring_coords.p2()).z() / ring_spacing);
 
 #endif
 
@@ -427,10 +437,6 @@ ProjDataInfoCylindricalNoArcCorr::find_cartesian_coordinates_of_detection(Cartes
   DetectionPositionPair<> dpp;
   get_det_pos_pair_for_bin(dpp, bin);
 
-  /* TODO
-   best to use Scanner::get_coordinate_for_det_pos().
-   Sadly, the latter is not yet implemented for Cylindrical scanners.
-  */
   // find corresponding cartesian coordinates
   find_cartesian_coordinates_given_scanner_coordinates(coord_1,
                                                        coord_2,
@@ -497,8 +503,9 @@ ProjDataInfoCylindricalNoArcCorr::find_cartesian_coordinates_given_scanner_coord
   cyl_coords.p1().z() = r1 * get_scanner_ptr()->get_ring_spacing();
   cyl_coords.p2().z() = r2 * get_scanner_ptr()->get_ring_spacing();
   LORAs2Points<float> lor(cyl_coords);
-  coord_1 = lor.p1();
-  coord_2 = lor.p2();
+  CartesianCoordinate3D<float> offset_gantry_coords_to_ring_coords = get_vector_centre_of_first_ring_to_centre_of_gantry();
+  coord_1 = lor.p1() - offset_gantry_coords_to_ring_coords;
+  coord_2 = lor.p2() - offset_gantry_coords_to_ring_coords;
 
 #endif
   if (tpos < 0)
@@ -538,28 +545,29 @@ ProjDataInfoCylindricalNoArcCorr::get_bin(const LOR<float>& lor, const double de
   Bin bin;
 #ifndef STIR_DEVEL
   // find nearest bin by going to nearest detectors first
-  LORInCylinderCoordinates<float> cyl_coords;
-  if (lor.change_representation(cyl_coords, get_ring_radius()) == Succeeded::no)
+  LORInCylinderCoordinates<float> cyl_in_gantry_coords;
+  if (lor.change_representation(cyl_in_gantry_coords, get_ring_radius()) == Succeeded::no)
     {
       bin.set_bin_value(-1);
       return bin;
     }
   const int num_detectors_per_ring = get_scanner_ptr()->get_num_detectors_per_ring();
   const int num_rings = get_scanner_ptr()->get_num_rings();
 
-  const int det1 = modulo(round((cyl_coords.p1().psi() - this->get_psi_offset()) / (2. * _PI / num_detectors_per_ring)),
+  const int det1 = modulo(round((cyl_in_gantry_coords.p1().psi() - this->get_psi_offset()) / (2. * _PI / num_detectors_per_ring)),
                           num_detectors_per_ring);
-  const int det2 = modulo(round((cyl_coords.p2().psi() - this->get_psi_offset()) / (2. * _PI / num_detectors_per_ring)),
+  const int det2 = modulo(round((cyl_in_gantry_coords.p2().psi() - this->get_psi_offset()) / (2. * _PI / num_detectors_per_ring)),
                           num_detectors_per_ring);
   // TODO WARNING LOR coordinates are w.r.t. centre of scanner, but the rings are numbered with the first ring at 0
-  const int ring1 = round(cyl_coords.p1().z() / get_ring_spacing() + (num_rings - 1) / 2.F);
-  const int ring2 = round(cyl_coords.p2().z() / get_ring_spacing() + (num_rings - 1) / 2.F);
+  const int ring1 = round(cyl_in_gantry_coords.p1().z() / get_ring_spacing() + (num_rings - 1) / 2.F);
+  const int ring2 = round(cyl_in_gantry_coords.p2().z() / get_ring_spacing() + (num_rings - 1) / 2.F);
 
   assert(det1 >= 0 && det1 < num_detectors_per_ring);
   assert(det2 >= 0 && det2 < num_detectors_per_ring);
 
   if (ring1 >= 0 && ring1 < num_rings && ring2 >= 0 && ring2 < num_rings
-      && get_bin_for_det_pair(bin, det1, ring1, det2, ring2, (cyl_coords.is_swapped() ? -1 : 1) * get_tof_bin(delta_time))
+      && get_bin_for_det_pair(
+             bin, det1, ring1, det2, ring2, (cyl_in_gantry_coords.is_swapped() ? -1 : 1) * get_tof_bin(delta_time))
              == Succeeded::yes
       && bin.tangential_pos_num() >= get_min_tangential_pos_num() && bin.tangential_pos_num() <= get_max_tangential_pos_num())
     {

diff --git a/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx b/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx
@@ -54,9 +54,9 @@ ProjDataInfoGenericNoArcCorr::ProjDataInfoGenericNoArcCorr(const shared_ptr<Scan
     error("ProjDataInfoGenericNoArcCorr: view mashing is not supported");
 
   // find shift between "new" centre-of-scanner and "old" centre-of-first-ring coordinate system
-  this->z_shift.z() = this->get_scanner_ptr()->get_coordinate_for_det_pos(DetectionPosition<>(0, 0, 0)).z();
-  this->z_shift.y() = 0;
-  this->z_shift.x() = 0;
+  // this->z_shift.z() = this->get_scanner_ptr()->get_coordinate_for_det_pos(DetectionPosition<>(0, 0, 0)).z();
+  // this->z_shift.y() = 0;
+  // this->z_shift.x() = 0;
 }
 
 ProjDataInfoGenericNoArcCorr*
@@ -122,8 +122,8 @@ ProjDataInfoGenericNoArcCorr::get_LOR(LORInAxialAndNoArcCorrSinogramCoordinates<
   CartesianCoordinate3D<float> _p2;
   find_cartesian_coordinates_of_detection(_p1, _p2, bin);
 
-  _p1.z() += z_shift.z();
-  _p2.z() += z_shift.z();
+  // _p1.z() += z_shift.z();
+  // _p2.z() += z_shift.z();
 
   LORAs2Points<float> lor_as_2_points(_p1, _p2);
   const double R = sqrt(std::max(square(_p1.x()) + square(_p1.y()), square(_p2.x()) + square(_p2.y())));
@@ -176,8 +176,8 @@ ProjDataInfoGenericNoArcCorr::find_cartesian_coordinates_given_scanner_coordinat
 
   coord_1 = get_scanner_ptr()->get_coordinate_for_det_pos(det_pos1);
   coord_2 = get_scanner_ptr()->get_coordinate_for_det_pos(det_pos2);
-  coord_1.z() -= z_shift.z();
-  coord_2.z() -= z_shift.z();
+  // coord_1.z() -= z_shift.z();
+  // coord_2.z() -= z_shift.z();
 
   if (timing_pos_num < 0)
     {

diff --git a/src/buildblock/VoxelsOnCartesianGrid.cxx b/src/buildblock/VoxelsOnCartesianGrid.cxx
@@ -371,32 +371,6 @@ VoxelsOnCartesianGrid<elemT>::grow_z_range(const int min_z, const int max_z)
   this->grow(IndexRange<3>(min_indices, max_indices));
 }
 
-template <class elemT>
-BasicCoordinate<3, int>
-VoxelsOnCartesianGrid<elemT>::get_lengths() const
-{
-  return make_coordinate(this->get_z_size(), this->get_y_size(), this->get_x_size());
-}
-
-template <class elemT>
-BasicCoordinate<3, int>
-VoxelsOnCartesianGrid<elemT>::get_min_indices() const
-{
-  CartesianCoordinate3D<int> min_indices;
-  CartesianCoordinate3D<int> max_indices;
-  this->get_regular_range(min_indices, max_indices);
-  return min_indices;
-}
-
-template <class elemT>
-BasicCoordinate<3, int>
-VoxelsOnCartesianGrid<elemT>::get_max_indices() const
-{
-  CartesianCoordinate3D<int> min_indices;
-  CartesianCoordinate3D<int> max_indices;
-  this->get_regular_range(min_indices, max_indices);
-  return max_indices;
-}
 #if 0
 
 /****************************

diff --git a/src/experimental/motion/RigidObject3DTransformation.cxx b/src/experimental/motion/RigidObject3DTransformation.cxx
@@ -290,6 +290,7 @@ RigidObject3DTransformation::transform_bin(Bin& bin,
   in_proj_data_info.get_LOR(lor, bin);
   LORAs2Points<float> lor_as_points;
   lor.get_intersections_with_cylinder(lor_as_points, lor.radius());
+#  if 0 // AG: No longer needed? so commented out
   // TODO origin
   // currently, the origin used for  proj_data_info is in the centre of the scanner,
   // while for standard images it is in the centre of the first ring.
@@ -300,6 +301,9 @@ RigidObject3DTransformation::transform_bin(Bin& bin,
   // also uses an origin in the centre of the first ring
   const float z_shift = (in_proj_data_info.get_scanner_ptr()->get_num_rings() - 1) / 2.F
                         * in_proj_data_info.get_scanner_ptr()->get_ring_spacing();
+#  else
+  constexpr float z_shift = 0.F;
+#  endif
   lor_as_points.p1().z() += z_shift;
   lor_as_points.p2().z() += z_shift;
   LORAs2Points<float> transformed_lor_as_points(transform_point(lor_as_points.p1()), transform_point(lor_as_points.p2()));

diff --git a/src/experimental/recon_buildblock/DataSymmetriesForDensels_PET_CartesianGrid.cxx b/src/experimental/recon_buildblock/DataSymmetriesForDensels_PET_CartesianGrid.cxx
@@ -85,15 +85,11 @@ find_relation_between_coordinate_systems(int& num_planes_per_scanner_ring,
 
       const float delta = proj_data_info_cyl_ptr->get_average_ring_difference(segment_num);
 
-      // KT 20/06/2001 take origin.z() into account
-      axial_pos_to_z_offset[segment_num]
-          = (cartesian_grid_info_ptr->get_max_index() + cartesian_grid_info_ptr->get_min_index()) / 2.F
-            - cartesian_grid_info_ptr->get_origin().z() / image_plane_spacing
-            - (num_planes_per_axial_pos[segment_num]
-                   * (proj_data_info_cyl_ptr->get_max_axial_pos_num(segment_num)
-                      + proj_data_info_cyl_ptr->get_min_axial_pos_num(segment_num))
-               + num_planes_per_scanner_ring * delta)
-                  / 2;
+      // AG: Each segment is offset because of the angle, we adjust using delta
+      const float segment_offset_in_z_gantry_coords = proj_data_info_cyl_ptr->get_m(first_bin) - ring_spacing * delta / 2;
+      axial_pos_to_z_offset[segment_num] = cartesian_grid_info_ptr->get_index_coordinates_for_physical_coordinates(
+          proj_data_info_cyl_ptr->get_physical_coordinates_for_gantry_coordinates(
+              CartesianCoordinate3D<float>(segment_offset_in_z_gantry_coords, 0, 0)))[1]; // z
     }
 }
 

diff --git a/src/experimental/test/test_proj_data_info_LOR.cxx b/src/experimental/test/test_proj_data_info_LOR.cxx
@@ -159,8 +159,7 @@ main(int argc, char* argv[])
           ring1_90,
           ring2_90,
           det1_90,
-          det2_90,
-          0); // set timing_pos_num=0 as test-code is pre-TOF
+          det2_90); // set timing_pos_num=0 as test-code is pre-TOF
 
 #if 0
        cout << coord_1_0<<endl;