Refactoring for origin shift

src/buildblock/ProjDataInfo.cxx

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

@@ -491,6 +491,7 @@ find_cartesian_coordinates_given_scanner_coordinates (CartesianCoordinate3D<floa
   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();
+  // ORIGINTODO: ^
   LORAs2Points<float> lor(cyl_coords);  
   coord_1 = lor.p1();
   coord_2 = lor.p2();

src/buildblock/zoom.cxx

@@ -283,6 +283,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

@@ -42,6 +42,7 @@
 #include "stir/CartesianCoordinate3D.h"
 #include "stir/Array.h"
 #include "stir/ExamData.h"
+#include "stir/ProjDataInfo.h"
 #include "stir/shared_ptr.h"
 #include <string>
 

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>
@@ -255,6 +256,16 @@ class ProjDataInfo
 	    const Bin&) const = 0;
   //@}
 
+  //! Get a point in gantry space along an LOR
+  /*!
+    The point is parameterised by s, a, m, cos(phi), sin(phi) and tan(theta).
+    Gantry space is defined w.r.t. the center of the gantry.
+  */
+  virtual CartesianCoordinate3D<float>
+  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;
+
   //! \name Functions that return info on the sampling in the different coordinates
   //@{
   //! Get sampling distance in the \c t coordinate
@@ -364,7 +375,15 @@ 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;
+
+  // Convert coordinates from a sinogram-based
+  inline CartesianCoordinate3D<float>
+  get_physical_coordinates_for_gantry_coordinates
+  (const CartesianCoordinate3D<float>& coords) const;
+
 protected:
   virtual bool blindly_equals(const root_type * const) const = 0;
 
@@ -378,7 +397,10 @@ class ProjDataInfo
   VectorWithOffset<int> max_axial_pos_per_seg;
   float bed_position_horizontal;
   float bed_position_vertical;
-
+
+  //! 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;
 };
 
 END_NAMESPACE_STIR

src/include/stir/ProjDataInfo.inl

@@ -117,6 +117,15 @@ ProjDataInfo::get_scanner_sptr() const
   return scanner_ptr;
 }
 
+CartesianCoordinate3D<float>
+ProjDataInfo::
+get_physical_coordinates_for_gantry_coordinates
+(const CartesianCoordinate3D<float>& coords) const
+{
+  // TODO: bed postion
+  return coords + get_vector_centre_of_first_ring_to_centre_of_gantry();
+}
+
 
 END_NAMESPACE_STIR
 

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

@@ -166,7 +166,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

@@ -111,6 +111,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(