Add "additional origin"

src/IO/InterfileHeader.cxx

@@ -441,6 +441,8 @@ InterfileImageHeader::InterfileImageHeader()
 
   add_key("first pixel offset (mm)",
 	   KeyArgument::DOUBLE, &first_pixel_offsets);
+  add_key("first pixel vendor offset (mm)",
+	   KeyArgument::DOUBLE, &first_pixel_vendor_offsets);
   add_key("number of image data types", 
     KeyArgument::INT,	(KeywordProcessor)&InterfileImageHeader::read_image_data_types,&num_image_data_types);
   add_key("index nesting level", 
@@ -466,8 +468,11 @@ read_matrix_info()
 {
   base_type::read_matrix_info();
   this->first_pixel_offsets.resize(num_dimensions);
+  this->first_pixel_vendor_offsets.resize(num_dimensions);
   std::fill(this->first_pixel_offsets.begin(), this->first_pixel_offsets.end(),
 	    base_type::double_value_not_set);
+  std::fill(this->first_pixel_vendor_offsets.begin(), this->first_pixel_vendor_offsets.end(),
+	    base_type::double_value_not_set);
 }
 
 // MJ 17/05/2000 made bool
@@ -498,6 +503,7 @@ bool InterfileImageHeader::post_processing()
       return true; 
     }
   std::vector<double>	first_pixel_offsets;
+  std::vector<double>	first_pixel_vendor_offsets;
 
   if (num_time_frames > 1 && num_image_data_types > 1)
     { 

src/IO/interfile.cxx

@@ -131,12 +131,25 @@ create_image_and_header_from(InterfileImageHeader& hdr,
 	- voxel_size * BasicCoordinate<3,float>(min_indices);
     }
 
-  return
+  VoxelsOnCartesianGrid<float>* image =
     new VoxelsOnCartesianGrid<float>
     (hdr.get_exam_info_sptr(),
      IndexRange<3>(min_indices, max_indices),
      origin,
      voxel_size);
+
+  if (hdr.first_pixel_vendor_offsets[2] != InterfileHeader::double_value_not_set)
+    {
+      // make sure that origin is such that 
+      // first_pixel_offsets =  min_indices*voxel_size + origin
+      image->set_vendor_origin
+        (make_coordinate(float(hdr.first_pixel_vendor_offsets[2]),
+                         float(hdr.first_pixel_vendor_offsets[1]),
+                         float(hdr.first_pixel_vendor_offsets[0]))
+         - voxel_size * BasicCoordinate<3,float>(min_indices));
+    }
+
+  return image;
 }
 
 VoxelsOnCartesianGrid<float> *
@@ -531,7 +544,8 @@ write_basic_interfile_image_header(const string& header_file_name,
 				   const ByteOrder byte_order,
 				   const VectorWithOffset<float>& scaling_factors,
 				   const VectorWithOffset<unsigned long>& file_offsets,
-                   const std::vector<std::string>& data_type_descriptions)
+				   const CartesianCoordinate3D<float>& vendor_origin,
+				   const std::vector<std::string>& data_type_descriptions)
 {
   CartesianCoordinate3D<int> min_indices;
   CartesianCoordinate3D<int> max_indices;
@@ -615,6 +629,18 @@ write_basic_interfile_image_header(const string& header_file_name,
       output_header << "first pixel offset (mm) [3] := " 
 		    << first_pixel_offsets.z() << '\n';
     }
+
+  if (vendor_origin.z() != InterfileHeader::double_value_not_set)
+    {
+      const CartesianCoordinate3D<float> first_pixel_vendor_offsets =
+	voxel_size * BasicCoordinate<3,float>(min_indices) + vendor_origin;
+      output_header << "first pixel vendor offset (mm) [1] := " 
+		    << first_pixel_vendor_offsets.x() << '\n';
+      output_header << "first pixel vendor offset (mm) [2] := " 
+		    << first_pixel_vendor_offsets.y() << '\n';
+      output_header << "first pixel vendor offset (mm) [3] := " 
+		    << first_pixel_vendor_offsets.z() << '\n';
+    }
 
   for (int i=1; i<=scaling_factors.get_length();i++)
     {
@@ -742,11 +768,14 @@ write_basic_interfile(const string& filename,
 {
   CartesianCoordinate3D<float> origin;
   origin.fill(static_cast<float>(InterfileHeader::double_value_not_set));
+  CartesianCoordinate3D<float> vendor_origin;
+  vendor_origin.fill(static_cast<float>(InterfileHeader::double_value_not_set));
   return
     write_basic_interfile(filename, 
 			  image, 
 			  CartesianCoordinate3D<float>(1,1,1), 
 			  origin,
+			  vendor_origin,
 			  output_type,
 			  scale,
 			  byte_order);
@@ -761,6 +790,7 @@ Succeeded write_basic_interfile(const string&  filename,
                 const Array<3,NUMBER>& image,
                 const CartesianCoordinate3D<float>& voxel_size,
                 const CartesianCoordinate3D<float>& origin,
+                const CartesianCoordinate3D<float>& vendor_origin,
                 const NumericType output_type,
                 const float scale,
                 const ByteOrder byte_order)
@@ -789,7 +819,8 @@ Succeeded write_basic_interfile(const string&  filename,
                          output_type,
                          byte_order,
                          scaling_factors,
-                         file_offsets);
+                         file_offsets,
+                         vendor_origin);
   #if 0
     delete[] header_name;
     delete[] data_name;
@@ -802,6 +833,7 @@ Succeeded write_basic_interfile(const string&  filename,
 				const Array<3,NUMBER>& image,
 				const CartesianCoordinate3D<float>& voxel_size,
 				const CartesianCoordinate3D<float>& origin,
+				const CartesianCoordinate3D<float>& vendor_origin,
 				const NumericType output_type,
 				const float scale,
 				const ByteOrder byte_order)
@@ -812,6 +844,7 @@ Succeeded write_basic_interfile(const string&  filename,
                                image,
                                voxel_size,
                                origin,
+                               vendor_origin,
                                output_type,
                                scale,
                                byte_order);
@@ -830,6 +863,7 @@ write_basic_interfile(const string&  filename,
 			  image, // use automatic reference to base class
 			  image.get_grid_spacing(), 
 			  image.get_origin(),
+			  image.get_vendor_origin(),
 			  output_type,
 			  scale,
 			  byte_order);
@@ -890,6 +924,7 @@ write_basic_interfile(const string& filename,
                          byte_order,
                          scaling_factors,
                          file_offsets,
+                         image.get_vendor_origin(),
                          data_type_descriptions);
   #if 0
     delete[] header_name;
@@ -934,7 +969,8 @@ write_basic_interfile(const string& filename,
                          output_type,
                          byte_order,
                          scaling_factors,
-                         file_offsets);
+                         file_offsets,
+                         image.get_density(1).get_vendor_origin());
   #if 0
     delete[] header_name;
     delete[] data_name;
@@ -1456,6 +1492,7 @@ write_basic_interfile<>(const string&  filename,
 			const Array<3,signed short>&,
 			const CartesianCoordinate3D<float>& voxel_size,
 			const CartesianCoordinate3D<float>& origin,
+			const CartesianCoordinate3D<float>& vendor_origin,
 			const NumericType output_type,
 			const float scale,
 			const ByteOrder byte_order);
@@ -1465,6 +1502,7 @@ write_basic_interfile<>(const string&  filename,
 			const Array<3,unsigned short>&,
 			const CartesianCoordinate3D<float>& voxel_size,
 			const CartesianCoordinate3D<float>& origin,
+			const CartesianCoordinate3D<float>& vendor_origin,
 			const NumericType output_type,
 			const float scale,
 			const ByteOrder byte_order);
@@ -1475,6 +1513,7 @@ write_basic_interfile<>(const string&  filename,
 			const Array<3,float>&,
 			const CartesianCoordinate3D<float>& voxel_size,
 			const CartesianCoordinate3D<float>& origin,
+			const CartesianCoordinate3D<float>& vendor_origin,
 			const NumericType output_type,
 			const float scale,
 			const ByteOrder byte_order);

src/buildblock/ProjDataInfo.cxx

@@ -652,5 +652,46 @@ 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>
+    (/*z=*/bed_position_horizontal, /*y=*/bed_position_vertical, /*x=*/0);
+}
+
+const CartesianCoordinate3D<float>
+ProjDataInfo::
+get_location_of_vendor_frame_of_reference_in_physical_coordinates() const {
+  std::cerr
+    << "vector_centre_of_first_ring_to_centre_of_gantry: "
+    << get_vector_centre_of_first_ring_to_centre_of_gantry() << std::endl
+    << "bed position: "
+    << get_bed_position() << std::endl;
+  return
+    // vendor FoR is from the centre of the gantry, STIR is from first ring
+    get_vector_centre_of_first_ring_to_centre_of_gantry()
+    // vendor FoR is bed-based
+    + get_bed_position();
+}
+
 END_NAMESPACE_STIR
 

src/buildblock/VoxelsOnCartesianGrid.cxx

@@ -185,63 +185,29 @@ VoxelsOnCartesianGrid<elemT>::VoxelsOnCartesianGrid(const ProjDataInfo& proj_dat
                                                     const CartesianCoordinate3D<int>& sizes)
 
 {
-  this->set_origin(origin);
-
-  int z_size = sizes.z();
-  // initialise to 0 to prevent compiler warnings
-  //int z_size = 0;
-  float z_sampling = 0;
-  float s_sampling = 0;
-  find_sampling_and_z_size(z_sampling, s_sampling, z_size, &proj_data_info);
-
-  this->set_grid_spacing(
-      CartesianCoordinate3D<float>(z_sampling, s_sampling/zoom, s_sampling/zoom)
-      );
-  int x_size_used = sizes.x();
-  int y_size_used = sizes.y();
-
-  if (sizes.x()==-1 || sizes.y()==-1)
-    {
-      // default it to cover full FOV by taking image_size>=2*FOVradius_in_pixs+1
-      const float FOVradius_in_mm = 
-        max(proj_data_info.get_s(Bin(0,0,0,proj_data_info.get_max_tangential_pos_num())),
-            -proj_data_info.get_s(Bin(0,0,0,proj_data_info.get_min_tangential_pos_num())));
-      if (sizes.x()==-1)
-        x_size_used = 2*static_cast<int>(ceil(FOVradius_in_mm / get_voxel_size().x())) + 1;
-      if (sizes.y()==-1)
-        y_size_used = 2*static_cast<int>(ceil(FOVradius_in_mm / get_voxel_size().y())) + 1;        
-    }
-  if (x_size_used<0)
-    error("VoxelsOnCartesianGrid: attempt to construct image with negative x_size %d\n", 
-          x_size_used);
-  if (x_size_used==0)
-    warning("VoxelsOnCartesianGrid: constructed image with x_size 0\n");
-  if (y_size_used<0)
-    error("VoxelsOnCartesianGrid: attempt to construct image with negative y_size %d\n", 
-          y_size_used);
-  if (y_size_used==0)
-    warning("VoxelsOnCartesianGrid: constructed image with y_size 0\n");
-
-  IndexRange3D range (0, z_size-1, 
-                      -(y_size_used/2), -(y_size_used/2) + y_size_used-1,
-                      -(x_size_used/2), -(x_size_used/2) + x_size_used-1);
-
-
-  this->grow(range);
+  init_from_proj_data_info(proj_data_info, zoom, origin, sizes);
 }
 
 template<class elemT>                                                 
-VoxelsOnCartesianGrid<elemT>::VoxelsOnCartesianGrid(const shared_ptr < ExamInfo > & exam_info_sptr_v,
-                                                    const ProjDataInfo& proj_data_info,
-                                                    const float zoom, 
-                                                    const CartesianCoordinate3D<float>& origin,
-                                                    const CartesianCoordinate3D<int>& sizes)
+VoxelsOnCartesianGrid<elemT>::
+VoxelsOnCartesianGrid(const shared_ptr<ExamInfo>& exam_info_sptr_v,
+                      const ProjDataInfo& proj_data_info,
+                      const float zoom,
+                      const CartesianCoordinate3D<float>& origin,
+                      const CartesianCoordinate3D<int>& sizes)
 {
   this->exam_info_sptr = exam_info_sptr_v;
-  // sadly, this code is a complete copy of the above
-  // probably avoidable in C++11
-  this->set_origin(origin);
+  init_from_proj_data_info(proj_data_info, zoom, origin, sizes);
+}
 
+template<class elemT>
+void
+VoxelsOnCartesianGrid<elemT>::
+init_from_proj_data_info(const ProjDataInfo& proj_data_info,
+                         const float zoom,
+                         const CartesianCoordinate3D<float>& origin,
+                         const CartesianCoordinate3D<int>& sizes)
+{
   int z_size = sizes.z();
   // initialise to 0 to prevent compiler warnings
   //int z_size = 0;
@@ -281,6 +247,11 @@ VoxelsOnCartesianGrid<elemT>::VoxelsOnCartesianGrid(const shared_ptr < ExamInfo
                       -(y_size_used/2), -(y_size_used/2) + y_size_used-1,
                       -(x_size_used/2), -(x_size_used/2) + x_size_used-1);
 
+  this->set_origin(origin);
+  this->set_vendor_origin
+    (origin
+     - (proj_data_info
+        .get_location_of_vendor_frame_of_reference_in_physical_coordinates()));
 
   this->grow(range);
 }

src/include/stir/BasicCoordinate.h

@@ -310,6 +310,11 @@ inline
 BasicCoordinate<num_dimensions,float> 
 convert_int_to_float(const BasicCoordinate<num_dimensions,int>& cint);
 
+//! BasicCoordinate can be converted to a string stream "[x1,x2,x3]"
+/*! \ingroup Coordinate */
+template <int num_dimensions, class coordT>
+std::ostream &operator<<(std::ostream &os, const BasicCoordinate<num_dimensions, coordT>& c);
+
 
 END_NAMESPACE_STIR
 

src/include/stir/BasicCoordinate.inl

@@ -612,6 +612,19 @@ operator<(const BasicCoordinate<num_dimensions, coordT>& c1,const BasicCoordinat
 					   c1, c); \
 }
 
+
+template <int num_dimensions, class coordT>
+std::ostream &operator<<(std::ostream &os, const BasicCoordinate<num_dimensions, coordT>& c) {
+{
+  os << "[";
+  for (int i=1; i<=num_dimensions; i++) {
+    os << c[i];
+    if (i<num_dimensions)
+      os << ",";
+  }
+  return os << "]";
+}
+
 DEFINE_OPERATOR_LESS(2)
 DEFINE_OPERATOR_LESS(3)
 DEFINE_OPERATOR_LESS(4)

src/include/stir/CartesianCoordinate2D.h

@@ -70,6 +70,9 @@ class CartesianCoordinate2D : public Coordinate2D<coordT>
 
 };
 
+template <class coordT>
+std::ostream &operator<<(std::ostream &os, const CartesianCoordinate2D<coordT> &coord);
+
 END_NAMESPACE_STIR
 
 #include "stir/CartesianCoordinate2D.inl"

src/include/stir/CartesianCoordinate2D.inl

@@ -92,5 +92,10 @@ CartesianCoordinate2D<coordT>::x() const
   return this->operator[](2);
 }
 
+template <class coordT>
+std::ostream &operator<<(std::ostream &os, const CartesianCoordinate2D<coordT> &coord) {
+  return os << "(x=" << coord.x()
+            << ",y=" << coord.y() << ")";
+}
 
 END_NAMESPACE_STIR

src/include/stir/CartesianCoordinate3D.h

@@ -84,6 +84,9 @@ class CartesianCoordinate3D : public Coordinate3D<coordT>
 
 };
 
+template <class coordT>
+std::ostream &operator<<(std::ostream &os, const CartesianCoordinate3D<coordT> &coord);
+
 END_NAMESPACE_STIR
 
 #include "stir/CartesianCoordinate3D.inl"

src/include/stir/CartesianCoordinate3D.inl

@@ -131,5 +131,11 @@ CartesianCoordinate3D<coordT>::x() const
   return this->operator[](3);
 }
 
+template <class coordT>
+std::ostream &operator<<(std::ostream &os, const CartesianCoordinate3D<coordT> &coord) {
+  return os << "(x=" << coord.x()
+            << ",y=" << coord.y()
+            << ",z=" << coord.z() << ")";
+}
 
 END_NAMESPACE_STIR