Refactoring for origin shift

src/buildblock/ProjDataInfo.cxx

@@ -652,5 +652,22 @@ operator>=(const ProjDataInfo& proj_data_info) const
   return (proj_data_info == *smaller_proj_data_info_sptr);
 }
 
+CartesianCoordinate3D<float>
+ProjDataInfo::
+get_bed_position() const
+{
+  return CartesianCoordinate3D<float>
+    (bed_position_horizontal, bed_position_vertical, 0);
+}
+
+CartesianCoordinate3D<float>
+ProjDataInfo::
+get_centre_of_gantry_vector_in_relative_coordinates() const
+{
+  double 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);
+}
+
 END_NAMESPACE_STIR
 

src/buildblock/ProjDataInfoCylindrical.cxx

@@ -127,14 +127,16 @@ initialise_ring_diff_arrays() const
       VectorWithOffset<float>(get_min_segment_num(),get_max_segment_num());
 
     /* m_offsets are found by requiring
-    get_m(..., min_axial_pos_num,...) == - get_m(..., max_axial_pos_num,...)
+    get_m(..., min_axial_pos_num,...) == - get_m(..., max_axial_pos_num,...),
+    and then shift from the centre to first ring and account for bed position
     */
     for (int segment_num=get_min_segment_num(); segment_num<=get_max_segment_num(); ++segment_num)
     {
       m_offset[segment_num] =
         ((get_max_axial_pos_num(segment_num) + get_min_axial_pos_num(segment_num))
 	 *get_axial_sampling(segment_num)
         )/2;
+      // ORIGINTODO: ^
     }
   }
   // initialise ax_pos_num_offset 
@@ -269,6 +271,7 @@ initialise_ring_diff_arrays() const
         const float ring1_plus_ring2_float =
           2*ax_pos_num/get_num_axial_poss_per_ring_inc(s_num)
           -2*m_offset[s_num]/ring_spacing + (get_scanner_ptr()->get_num_rings()-1);
+        // ORIGINTODO: ^
         const int ring1_plus_ring2 =
           round(ring1_plus_ring2_float);
         // check that it was integer

src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx

@@ -417,6 +417,7 @@ find_scanner_coordinates_given_cartesian_coordinates(int& det1, int& det2, int&
   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);
+  // ORIGINTODO: ^
 
 #endif
 
@@ -489,8 +490,9 @@ find_cartesian_coordinates_given_scanner_coordinates (CartesianCoordinate3D<floa
   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.p1().z() = Ring_A*get_scanner_ptr()->get_ring_spacing(); 
   cyl_coords.p2().z() = Ring_B*get_scanner_ptr()->get_ring_spacing();
+  // ORIGINTODO: ^
   LORAs2Points<float> lor(cyl_coords);  
   coord_1 = lor.p1();
   coord_2 = lor.p2();

src/buildblock/zoom.cxx

@@ -280,6 +280,7 @@ construct_new_image_from_zoom_parameters(const VoxelsOnCartesianGrid<float> &ima
      new_image.get_physical_coordinates_for_indices(new_max_indices))/2;
   // now make sure that these are shifted as required
   new_image.set_origin(offsets_in_mm + middle - new_middle);
+  // ORIGINTODO: ^
   // check
   {
     const BasicCoordinate<3,float> final_middle =

src/experimental/recon_buildblock/DataSymmetriesForDensels_PET_CartesianGrid.cxx

@@ -105,6 +105,7 @@ find_relation_between_coordinate_systems(int& num_planes_per_scanner_ring,
        *(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;
+    // ORIGINDOTO: ^
   }
 }
 

src/include/stir/DiscretisedDensity.h

@@ -41,6 +41,7 @@
 #include "stir/CartesianCoordinate3D.h"
 #include "stir/Array.h"
 #include "stir/IO/ExamData.h"
+#include "stir/ProjDataInfo.h"
 #include "stir/shared_ptr.h"
 #include <string>
 
@@ -253,6 +254,16 @@ template<int num_dimensions, typename elemT>
   inline BasicCoordinate<num_dimensions,int>
   get_indices_closest_to_LPS_coordinates(const CartesianCoordinate3D<float>& coords) const;
 
+  inline CartesianCoordinate3D<float>
+  get_relative_coordinates_for_gantry_coordinates
+  (const CartesianCoordinate3D<float>& coords,
+   const shared_ptr<ProjDataInfo> proj_data_info_sptr) const;
+
+  inline BasicCoordinate<num_dimensions, float>
+  get_index_coordinates_for_gantry_coordinates
+  (const CartesianCoordinate3D<float>& coords,
+   const shared_ptr<ProjDataInfo> proj_data_info_sptr) const;
+
   //@}
 
   //! Allocate a new DiscretisedDensity object with same characteristics as the current one.

src/include/stir/DiscretisedDensity.inl

@@ -304,4 +304,29 @@ get_indices_closest_to_LPS_coordinates(const CartesianCoordinate3D<float>& coord
   return round(this->get_index_coordinates_for_LPS_coordinates(coords));
 }
 
+template<int num_dimensions, typename elemT>
+CartesianCoordinate3D<float>
+DiscretisedDensity<num_dimensions, elemT>::
+get_relative_coordinates_for_gantry_coordinates
+(const CartesianCoordinate3D<float>& coords,
+ const shared_ptr<ProjDataInfo> proj_data_info_ptr) const
+{
+  // gantry coordinates are like relative coordinates, but the (0, 0, 0)
+  // point is at the center of the gantry
+  return coords
+      + proj_data_info_ptr->get_centre_of_gantry_vector_in_relative_coordinates();
+}
+
+template<int num_dimensions, typename elemT>
+BasicCoordinate<num_dimensions,float>
+DiscretisedDensity<num_dimensions, elemT>::
+get_index_coordinates_for_gantry_coordinates
+(const CartesianCoordinate3D<float>& coords,
+ const shared_ptr<ProjDataInfo> proj_data_info_ptr) const
+{
+  return this->get_index_coordinates_for_relative_coordinates
+    (this->get_relative_coordinates_for_gantry_coordinates
+     (coords, proj_data_info_ptr));
+}
+
 END_NAMESPACE_STIR

src/include/stir/ProjDataInfo.h

@@ -33,6 +33,7 @@
 
 #include "stir/VectorWithOffset.h"
 #include "stir/Scanner.h"
+#include "stir/CartesianCoordinate3D.h"
 #include "stir/shared_ptr.h"
 #include "stir/unique_ptr.h"
 #include <string>
@@ -364,7 +365,14 @@ class ProjDataInfo
 
   //! Get vertical bed position
   float get_bed_position_vertical() const { return bed_position_vertical; }
-
+
+  //! Vector represention bed position in 3D
+  CartesianCoordinate3D<float> get_bed_position() const;
+
+  //! Vector from image frame of reference (centre of first ring) to gantry centre
+  CartesianCoordinate3D<float>
+  get_centre_of_gantry_vector_in_relative_coordinates() const;
+
 protected:
   virtual bool blindly_equals(const root_type * const) const = 0;
 

src/include/stir/recon_buildblock/ProjMatrixByBinUsingRayTracing.h

@@ -199,7 +199,8 @@ public :
   CartesianCoordinate3D<int> min_index;
   CartesianCoordinate3D<int> max_index;
 
-  shared_ptr<ProjDataInfo> proj_data_info_ptr;
+  shared_ptr<ProjDataInfo> proj_data_info_sptr;
+  shared_ptr< DiscretisedDensity<3,float> > density_info_sptr;
 
 
   virtual void 

src/include/stir/zoom.h

@@ -164,7 +164,7 @@ zoom_image(const VoxelsOnCartesianGrid<float> &image,
 
   The index range of the new image is according to the standard STIR
   conventions (z starts from 0, but x and y from -(new_size/2)).
-  The origin is then chosen such that the geometric centres of the
+  The origin is then chosen such that the geometric centres of the // ORIGINTODO
   images satisfy
   \verbatim
   offsets_in_mm == new_middle - old_middle

src/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.cxx

@@ -107,6 +107,7 @@ find_relation_between_coordinate_systems(int& num_planes_per_scanner_ring,
        *(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;
+    // ORIGINTODO: ^
   }
 }
 

src/recon_buildblock/ProjMatrixByBinUsingInterpolation.cxx

@@ -127,6 +127,7 @@ set_up(
   origin = image_info_ptr->get_origin();
   const float z_to_middle =
     (densel_range.get_max_index() + densel_range.get_min_index())*voxel_size.z()/2.F;
+  // ORIGINTODO: ^
   origin.z() -= z_to_middle;
 
   symmetries_sptr.reset(