Prepare for origin shift 2020

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.
 
     This file is free software; you can redistribute it and/or modify
@@ -27,6 +29,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
@@ -88,6 +94,7 @@
 
 #include "stir/analytic/FBP3DRP/ColsherFilter.h" 
 #include "stir/display.h"
+//#include "stir/info.h"
 //#include "stir/recon_buildblock/distributable.h"
 //#include "stir/FBP3DRP/process_viewgrams.h"
 
@@ -131,55 +138,47 @@ static void find_rmin_rmax(int& rmin, int& rmax,
                            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
-
-  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;
+  // 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);
+  // 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;
 
+  // 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);
 
-   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;
+  // 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);
+  }
 
-  // finally find the 'ring' number
-  rmin = static_cast<int>(floor(z_in_virtual_ring_coordinates));
-
-
-  // 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);
 }
 
 

src/buildblock/ProjDataInfo.cxx

@@ -654,5 +654,31 @@ 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
 

src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx

@@ -363,9 +363,10 @@ get_all_det_pos_pairs_for_bin(vector<DetectionPositionPair<> >& dps,
   assert(current_dp_num == get_num_det_pos_pairs_for_bin(bin));
 }
 
+//get_det_pair_for_gantry_coordinate_pair
 Succeeded
 ProjDataInfoCylindricalNoArcCorr::
-find_scanner_coordinates_given_cartesian_coordinates(int& det1, int& det2, int& ring1, int& ring2,
+get_det_pair_for_gantry_coordinate_pair(int& det1, int& det2, int& ring1, int& ring2,
 					             const CartesianCoordinate3D<float>& c1,
 						     const CartesianCoordinate3D<float>& c2) const
 {
@@ -406,17 +407,25 @@ find_scanner_coordinates_given_cartesian_coordinates(int& det1, int& det2, int&
   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)
+
+  // 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()/(2.*_PI/num_detectors)), num_detectors);
-  det2 = modulo(round(cyl_coords.p2().psi()/(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()/(2.*_PI/num_detectors)), num_detectors);
+  det2 = modulo(round(cyl_in_ring_coords.p2().psi()/(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
 
@@ -432,7 +441,7 @@ find_scanner_coordinates_given_cartesian_coordinates(int& det1, int& det2, int&
 
 void 
 ProjDataInfoCylindricalNoArcCorr::
-find_cartesian_coordinates_of_detection(
+get_bin_detector_locations_in_gantry_coordinates(
 					CartesianCoordinate3D<float>& coord_1,
 					CartesianCoordinate3D<float>& coord_2,
 					const Bin& bin) const
@@ -446,14 +455,14 @@ find_cartesian_coordinates_of_detection(
                        det_num_b, ring_b, bin);
 
   // find corresponding cartesian coordinates
-  find_cartesian_coordinates_given_scanner_coordinates(coord_1,coord_2,
+  get_det_pair_locations_in_gantry_coordinates(coord_1,coord_2,
     ring_a,ring_b,det_num_a,det_num_b);
 }
 
 
 void
 ProjDataInfoCylindricalNoArcCorr::
-find_cartesian_coordinates_given_scanner_coordinates (CartesianCoordinate3D<float>& coord_1,
+get_det_pair_locations_in_gantry_coordinates(CartesianCoordinate3D<float>& coord_1,
 				 CartesianCoordinate3D<float>& coord_2,
 				 const int Ring_A,const int Ring_B, 
 				 const int det1, const int det2) const
@@ -486,22 +495,26 @@ find_cartesian_coordinates_given_scanner_coordinates (CartesianCoordinate3D<floa
   assert(0<=det2);
   assert(det2<num_detectors_per_ring);
 
-  LORInCylinderCoordinates<float> cyl_coords(get_scanner_ptr()->get_effective_ring_radius());
-  cyl_coords.p1().psi() = static_cast<float>((2.*_PI/num_detectors_per_ring)*(det1));
-  cyl_coords.p2().psi() = static_cast<float>((2.*_PI/num_detectors_per_ring)*(det2));
-  cyl_coords.p1().z() = Ring_A*get_scanner_ptr()->get_ring_spacing();
-  cyl_coords.p2().z() = Ring_B*get_scanner_ptr()->get_ring_spacing();
-  LORAs2Points<float> lor(cyl_coords);  
-  coord_1 = lor.p1();
-  coord_2 = lor.p2();
+  // 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.
+  LORInCylinderCoordinates<float> cyl_in_ring_coords(get_scanner_ptr()->get_effective_ring_radius());
+  cyl_in_ring_coords.p1().psi() = static_cast<float>((2.*_PI/num_detectors_per_ring)*(det1));
+  cyl_in_ring_coords.p2().psi() = static_cast<float>((2.*_PI/num_detectors_per_ring)*(det2));
+  cyl_in_ring_coords.p1().z() = Ring_A*get_scanner_ptr()->get_ring_spacing();
+  cyl_in_ring_coords.p2().z() = Ring_B*get_scanner_ptr()->get_ring_spacing();
+  LORAs2Points<float> lor(cyl_in_ring_coords);  
+  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
 }
 
 
 void 
 ProjDataInfoCylindricalNoArcCorr::
-find_bin_given_cartesian_coordinates_of_detection(Bin& bin,
+get_bin_for_gantry_coordinate_pair(Bin& bin,
 						  const CartesianCoordinate3D<float>& coord_1,
 						  const CartesianCoordinate3D<float>& coord_2) const
 {
@@ -511,7 +524,7 @@ find_bin_given_cartesian_coordinates_of_detection(Bin& bin,
   int ring_b;
 
   // given two CartesianCoordinates find the intersection     
-  if (find_scanner_coordinates_given_cartesian_coordinates(det_num_a,det_num_b,
+  if (get_det_pair_for_gantry_coordinate_pair(det_num_a,det_num_b,
 							   ring_a, ring_b,
 							   coord_1,
 							   coord_2) ==
@@ -522,7 +535,7 @@ find_bin_given_cartesian_coordinates_of_detection(Bin& bin,
   }
 
   // check rings are in valid range
-  // this should have been done by find_scanner_coordinates_given_cartesian_coordinates
+  // this should have been done by get_det_pair_for_gantry_coordinate_pair
   assert(!(ring_a<0 ||
 	   ring_a>=get_scanner_ptr()->get_num_rings() ||
 	   ring_b<0 ||
@@ -543,8 +556,8 @@ get_bin(const LOR<float>& lor) const
   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;
@@ -554,11 +567,11 @@ get_bin(const LOR<float>& lor) const
   const int num_rings = 
     get_scanner_ptr()->get_num_rings();
 
-  const int det1 = modulo(round(cyl_coords.p1().psi()/(2.*_PI/num_detectors_per_ring)),num_detectors_per_ring);
-  const int det2 = modulo(round(cyl_coords.p2().psi()/(2.*_PI/num_detectors_per_ring)),num_detectors_per_ring);
+  const int det1 = modulo(round(cyl_in_gantry_coords.p1().psi()/(2.*_PI/num_detectors_per_ring)),num_detectors_per_ring);
+  const int det2 = modulo(round(cyl_in_gantry_coords.p2().psi()/(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);

src/buildblock/VoxelsOnCartesianGrid.cxx

@@ -372,35 +372,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
 
 /****************************

src/experimental/motion/RigidObject3DTransformation.cxx

@@ -274,7 +274,7 @@ transform_bin(Bin& bin,const ProjDataInfo& out_proj_data_info,
   CartesianCoordinate3D<float> coord_1;
   CartesianCoordinate3D<float> coord_2;
   dynamic_cast<const ProjDataInfoCylindricalNoArcCorr&>(in_proj_data_info).
-    find_cartesian_coordinates_of_detection(coord_1,coord_2,bin);
+    get_bin_detector_locations_in_gantry_coordinates(coord_1,coord_2,bin);
 
   // now do the movement
 
@@ -285,27 +285,30 @@ transform_bin(Bin& bin,const ProjDataInfo& out_proj_data_info,
     coord2transformed = transform_point(coord_2);
 
   dynamic_cast<const ProjDataInfoCylindricalNoArcCorr&>(out_proj_data_info).
-    find_bin_given_cartesian_coordinates_of_detection(bin,
+    get_bin_for_gantry_coordinate_pair(bin,
                                                       coord_1_transformed,
 					              coord2transformed);
 #else
   LORInAxialAndNoArcCorrSinogramCoordinates<float> lor;
   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?
   // 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.
   // This is pretty horrible though, as the transform_point function has no clue 
   // where the origin is
   // Note that the present shift will make this version compatible with the 
-  // version above, as find_bin_given_cartesian_coordinates_of_detection
+  // version above, as get_bin_for_gantry_coordinate_pair
   // 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();
   lor_as_points.p1().z() += z_shift;
   lor_as_points.p2().z() += z_shift;
+#endif
   LORAs2Points<float> 
     transformed_lor_as_points(transform_point(lor_as_points.p1()),
 			      transform_point(lor_as_points.p2()));