FIX: Attempt (and fail) to re-centre image in gantry for back-compat

src/buildblock/ProjDataInfo.cxx

@@ -44,6 +44,7 @@
 #include "stir/Coordinate3D.h"
 #include "stir/IndexRange2D.h"
 #include "stir/IndexRange3D.h"
+#include "stir/DiscretisedDensity.h"
 #include "stir/Bin.h"
 // include for ask and ask_num
 #include "stir/utilities.h"
@@ -673,6 +674,25 @@ get_vector_centre_of_first_ring_to_centre_of_gantry() const
   return CartesianCoordinate3D<float>(middle_of_first_ring_to_middle_of_last / 2.F, 0, 0);
 }
 
+CartesianCoordinate3D<float>
+ProjDataInfo::
+get_gantry_origin_in_physical_coordinates(const DiscretisedDensity<3,float>& image) const
+{
+  // TODO - STIR currently enforces that an image is centred in the z-plane.
+  // In future, this would make more sense to be implemented as:
+  // 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);
+
+  // shift by vector to centre of image in z (assume this corresponds to centre of gantry)
+  CartesianCoordinate3D<int> min_index;
+  CartesianCoordinate3D<int> max_index;
+  image.get_regular_range(min_index, max_index);
+  CartesianCoordinate3D<float> middle_of_image = image.get_physical_coordinates_for_indices(
+    (min_index + max_index) / 2.);
+  return CartesianCoordinate3D<float>(middle_of_image.z(), 0, 0);
+}
+
 CartesianCoordinate3D<float>
 ProjDataInfo::get_bed_position() const
 {

src/include/stir/ProjDataInfo.h

@@ -47,6 +47,7 @@ template <typename elemT> class Viewgram;
 template <typename elemT> class SegmentByView;
 template <typename elemT> class SegmentBySinogram;
 template <typename elemT> class RelatedViewgrams;
+template<int num_dimensions, typename elemT> class DiscretisedDensity;
 class DataSymmetriesForViewSegmentNumbers;
 class ViewSegmentNumbers;
 class Bin;
@@ -389,10 +390,12 @@ class ProjDataInfo
   //! Vector represention bed position in 3D
   CartesianCoordinate3D<float> get_bed_position() const;
 
-  // Convert coordinates from a sinogram-based
+  //! Convert coordinates from a sinogram-based
+  // STIR always centers images axially in the gantry, so unfortunately
+  // this requires knowledge of the image.
   inline CartesianCoordinate3D<float>
   get_physical_coordinates_for_gantry_coordinates
-  (const CartesianCoordinate3D<float>& coords) const;
+  (const CartesianCoordinate3D<float>& coords, const DiscretisedDensity<3,float>& image) const;
 
 protected:
   virtual bool blindly_equals(const root_type * const) const = 0;
@@ -411,6 +414,12 @@ class ProjDataInfo
   //! Vector from image frame of reference (centre of first ring) to gantry centre
   CartesianCoordinate3D<float>
   get_vector_centre_of_first_ring_to_centre_of_gantry() const;
+
+  //! Vector from image frame of reference (centre of first ring) to gantry centre
+  // STIR always centers images axially in the gantry, so unfortunately
+  // this requires knowledge of the image.
+  CartesianCoordinate3D<float>
+  get_gantry_origin_in_physical_coordinates(const DiscretisedDensity<3,float>& image) const;
 };
 
 END_NAMESPACE_STIR

src/include/stir/ProjDataInfo.inl

@@ -122,13 +122,13 @@ ProjDataInfo::get_scanner_sptr() const
 CartesianCoordinate3D<float>
 ProjDataInfo::
 get_physical_coordinates_for_gantry_coordinates
-(const CartesianCoordinate3D<float>& coords) const
+(const CartesianCoordinate3D<float>& coords, const DiscretisedDensity<3,float>& image) const
 {
   // TODO: bed postion
   return coords + get_vector_centre_of_first_ring_to_centre_of_gantry();
+  //return coords + get_gantry_origin_in_physical_coordinates(image);
 }
 
-
 int
 ProjDataInfo::get_num_sinograms() const
 {

src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx

@@ -394,13 +394,14 @@ CartesianCoordinate3D<float>
 get_point_on_lor_in_index_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 DiscretisedDensity<3, float>& density_info,
+ const DiscretisedDensity<3, float>& density,
  const ProjDataInfo& proj_data_info)
 {
-  return density_info.get_index_coordinates_for_physical_coordinates
+  return density.get_index_coordinates_for_physical_coordinates
     (proj_data_info.get_physical_coordinates_for_gantry_coordinates
      (proj_data_info.get_point_on_lor_in_gantry_coordinates
-      (s_in_mm, m_in_mm, a_in_mm, cphi, sphi, tantheta)));
+      (s_in_mm, m_in_mm, a_in_mm, cphi, sphi, tantheta),
+      density));
 }
 
 // just do 1 LOR, returns true if lor is not empty