From 98a2382f75320e5491051f4b41c8da7081c5576c Mon Sep 17 00:00:00 2001 From: Ashley Gillman Date: Mon, 3 Sep 2018 12:05:39 +0100 Subject: [PATCH 1/8] [REF] Use m_in_mm instead of t_in_mm in ProjMatrixByBinUsingRayTracing --- .../ProjMatrixByBinUsingRayTracing.cxx | 26 ++++++++++--------- 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx b/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx index 8560fc0f42..06deee462d 100644 --- a/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx +++ b/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx @@ -387,7 +387,7 @@ static inline int sign(const T& t) // just do 1 LOR, returns true if lor is not empty static void ray_trace_one_lor(ProjMatrixElemsForOneBin& lor, - const float s_in_mm, const float t_in_mm, + const float s_in_mm, const float m_in_mm, const float cphi, const float sphi, const float costheta, const float tantheta, const float offset_in_z, @@ -458,12 +458,14 @@ ray_trace_one_lor(ProjMatrixElemsForOneBin& lor, } //!restrict_to_cylindrical_FOV + + // start and stop point in voxel coordinates start_point.x() = (s_in_mm*cphi + max_a*sphi)/voxel_size.x(); - start_point.y() = (s_in_mm*sphi - max_a*cphi)/voxel_size.y(); - start_point.z() = (t_in_mm/costheta+offset_in_z - max_a*tantheta)/voxel_size.z(); + start_point.y() = (s_in_mm*sphi - max_a*cphi)/voxel_size.y(); + start_point.z() = (m_in_mm+offset_in_z - max_a*tantheta)/voxel_size.z(); stop_point.x() = (s_in_mm*cphi + min_a*sphi)/voxel_size.x(); - stop_point.y() = (s_in_mm*sphi - min_a*cphi)/voxel_size.y(); - stop_point.z() = (t_in_mm/costheta+offset_in_z - min_a*tantheta)/voxel_size.z(); + stop_point.y() = (s_in_mm*sphi - min_a*cphi)/voxel_size.y(); + stop_point.z() = (m_in_mm+offset_in_z - min_a*tantheta)/voxel_size.z(); #if 0 // KT 18/05/2005 this is no longer necessary @@ -595,7 +597,7 @@ calculate_proj_matrix_elems_for_one_bin( const float tantheta = proj_data_info_ptr->get_tantheta(bin); const float costheta = 1/sqrt(1+square(tantheta)); - const float t_in_mm = proj_data_info_ptr->get_t(bin); + const float m_in_mm = proj_data_info_ptr->get_m(bin); const float sampling_distance_of_adjacent_LORs_z = proj_data_info_ptr->get_sampling_in_t(bin)/costheta; @@ -656,10 +658,10 @@ calculate_proj_matrix_elems_for_one_bin( { // make sure we don't ray-trace exactly between 2 planes // z-coordinate (in voxel units) will be - // (t_in_mm+offset_in_z)/voxel_size.z(); + // (m_in_mm+offset_in_z)/voxel_size.z(); // if so, we ray trace first to the voxels at smaller z, but will add the // other plane later (in add_adjacent_z) - if (fabs(modulo((t_in_mm+offset_in_z)/voxel_size.z(),1.F)-.5)<.001) + if (fabs(modulo((m_in_mm+offset_in_z)/voxel_size.z(),1.F)-.5)<.001) offset_in_z -= .1F*voxel_size.z(); } @@ -678,7 +680,7 @@ calculate_proj_matrix_elems_for_one_bin( if (num_tangential_LORs == 1) { - ray_trace_one_lor(lor, s_in_mm, t_in_mm, + ray_trace_one_lor(lor, s_in_mm, m_in_mm, cphi, sphi, costheta, tantheta, offset_in_z, fovrad_in_mm, voxel_size, @@ -699,7 +701,7 @@ calculate_proj_matrix_elems_for_one_bin( for (int s_LOR_num=1; s_LOR_num<=num_tangential_LORs; ++s_LOR_num, current_s_in_mm+=s_inc) { ray_traced_lor.erase(); - ray_trace_one_lor(ray_traced_lor, current_s_in_mm, t_in_mm, + ray_trace_one_lor(ray_traced_lor, current_s_in_mm, m_in_mm, cphi, sphi, costheta, tantheta, offset_in_z, fovrad_in_mm, voxel_size, @@ -720,10 +722,10 @@ calculate_proj_matrix_elems_for_one_bin( origin.z()/voxel_size.z() - (max_index.z() + min_index.z())/2.F; const float left_edge_of_TOR = - (t_in_mm - sampling_distance_of_adjacent_LORs_z/2 + (m_in_mm - sampling_distance_of_adjacent_LORs_z/2 )/voxel_size.z(); const float right_edge_of_TOR = - (t_in_mm + sampling_distance_of_adjacent_LORs_z/2 + (m_in_mm + sampling_distance_of_adjacent_LORs_z/2 )/voxel_size.z(); add_adjacent_z(lor, z_of_first_voxel - left_edge_of_TOR, right_edge_of_TOR -left_edge_of_TOR); From 10d625d1b44d8ce7910bbf87675925eb601303fb Mon Sep 17 00:00:00 2001 From: Ashley Gillman Date: Mon, 3 Sep 2018 16:21:19 +0100 Subject: [PATCH 2/8] [REF] Simplify calculation and correction of phi and s_in_mm In ProjMatrixByBinUsingRayTracing --- .../ProjMatrixByBinUsingRayTracing.cxx | 27 +++++-------------- 1 file changed, 6 insertions(+), 21 deletions(-) diff --git a/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx b/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx index 06deee462d..ab294a1ace 100644 --- a/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx +++ b/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx @@ -545,25 +545,10 @@ calculate_proj_matrix_elems_for_one_bin( assert(lor.size() == 0); - float phi; + float phi = proj_data_info_ptr->get_phi(bin); float s_in_mm = proj_data_info_ptr->get_s(bin); - /* Implementation note. - KT initialised s_in_mm above instead of in the if because this meant - that gcc 3.0.1 generated identical results to the previous version of this file. - Otherwise, some pixels at the boundary appear to be treated differently - (probably due to different floating point rounding errors), at least - on Linux on x86. - A bit of a mistery that. - - TODO this is maybe solved now by having more decent handling of - start and end voxels. - */ - if (!use_actual_detector_boundaries) - { - phi = proj_data_info_ptr->get_phi(bin); - //s_in_mm = proj_data_info_ptr->get_s(bin); - } - else + + if (use_actual_detector_boundaries) { // can be static_cast later on const ProjDataInfoCylindricalNoArcCorr& proj_data_info_noarccor = @@ -580,16 +565,16 @@ calculate_proj_matrix_elems_for_one_bin( det_num2, bin.view_num(), bin.tangential_pos_num()); + + const float old_phi = phi; phi = static_cast((det_num1+det_num2)*_PI/num_detectors-_PI/2); - const float old_phi=proj_data_info_ptr->get_phi(bin); if (fabs(phi-old_phi)>2*_PI/num_detectors) warning("view %d old_phi %g new_phi %g\n",bin.view_num(), old_phi, phi); + const float old_s_in_mm = s_in_mm; s_in_mm = static_cast(ring_radius*sin((det_num1-det_num2)*_PI/num_detectors+_PI/2)); - const float old_s_in_mm=proj_data_info_ptr->get_s(bin); if (fabs(s_in_mm-old_s_in_mm)>proj_data_info_ptr->get_sampling_in_s(bin)*.0001) warning("tangential_pos_num %d old_s_in_mm %g new_s_in_mm %g\n",bin.tangential_pos_num(), old_s_in_mm, s_in_mm); - } const float cphi = cos(phi); From fa49ee19c0a05a29becc05b0f72738697de0acc3 Mon Sep 17 00:00:00 2001 From: Ashley Gillman Date: Mon, 3 Sep 2018 17:25:57 +0100 Subject: [PATCH 3/8] [REF] Remove unneeded argument to ray_trace_one_lor() --- src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx b/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx index ab294a1ace..232e2a5121 100644 --- a/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx +++ b/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx @@ -389,7 +389,7 @@ static void ray_trace_one_lor(ProjMatrixElemsForOneBin& lor, const float s_in_mm, const float m_in_mm, const float cphi, const float sphi, - const float costheta, const float tantheta, + const float tantheta, const float offset_in_z, const float fovrad_in_mm, const CartesianCoordinate3D& voxel_size, @@ -666,7 +666,7 @@ calculate_proj_matrix_elems_for_one_bin( if (num_tangential_LORs == 1) { ray_trace_one_lor(lor, s_in_mm, m_in_mm, - cphi, sphi, costheta, tantheta, + cphi, sphi, tantheta, offset_in_z, fovrad_in_mm, voxel_size, restrict_to_cylindrical_FOV, @@ -687,7 +687,7 @@ calculate_proj_matrix_elems_for_one_bin( { ray_traced_lor.erase(); ray_trace_one_lor(ray_traced_lor, current_s_in_mm, m_in_mm, - cphi, sphi, costheta, tantheta, + cphi, sphi, tantheta, offset_in_z, fovrad_in_mm, voxel_size, restrict_to_cylindrical_FOV, From 073baa1d82169c7a81c579ac29373ddbbe00cd71 Mon Sep 17 00:00:00 2001 From: Ashley Gillman Date: Thu, 18 Oct 2018 15:30:31 +1000 Subject: [PATCH 4/8] [REF] Address comments https://github.com/UCL/STIR/pull/260/files/0365f964a565fa488b278e6a5eb20fa337572b11 --- src/buildblock/ProjDataInfo.cxx | 26 +++++++++++++ src/include/stir/ProjDataInfo.h | 26 ++++++++++++- src/include/stir/ProjDataInfo.inl | 9 +++++ .../ProjMatrixByBinUsingRayTracing.cxx | 39 ++++++++++++------- 4 files changed, 84 insertions(+), 16 deletions(-) diff --git a/src/buildblock/ProjDataInfo.cxx b/src/buildblock/ProjDataInfo.cxx index 053157f8c4..38995f23d2 100644 --- a/src/buildblock/ProjDataInfo.cxx +++ b/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 +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(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 +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(middle_of_first_ring_to_middle_of_last / 2.F, 0, 0); +} + +CartesianCoordinate3D +ProjDataInfo::get_bed_position() const +{ + return CartesianCoordinate3D + (bed_position_horizontal, bed_position_vertical, 0); +} + END_NAMESPACE_STIR diff --git a/src/include/stir/ProjDataInfo.h b/src/include/stir/ProjDataInfo.h index 2bd1712155..46e94f4c4b 100644 --- a/src/include/stir/ProjDataInfo.h +++ b/src/include/stir/ProjDataInfo.h @@ -31,6 +31,7 @@ #ifndef __stir_ProjDataInfo_H__ #define __stir_ProjDataInfo_H__ +#include "stir/CartesianCoordinate3D.h" #include "stir/VectorWithOffset.h" #include "stir/Scanner.h" #include "stir/shared_ptr.h" @@ -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 + 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 get_bed_position() const; + + // Convert coordinates from a sinogram-based + inline CartesianCoordinate3D + get_relative_coordinates_for_gantry_coordinates + (const CartesianCoordinate3D& coords) const; + protected: virtual bool blindly_equals(const root_type * const) const = 0; @@ -378,7 +397,10 @@ class ProjDataInfo VectorWithOffset 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 + get_vector_centre_of_first_ring_to_centre_of_gantry() const; }; END_NAMESPACE_STIR diff --git a/src/include/stir/ProjDataInfo.inl b/src/include/stir/ProjDataInfo.inl index c01983db66..f5d88b750d 100644 --- a/src/include/stir/ProjDataInfo.inl +++ b/src/include/stir/ProjDataInfo.inl @@ -117,6 +117,15 @@ ProjDataInfo::get_scanner_sptr() const return scanner_ptr; } +CartesianCoordinate3D +ProjDataInfo:: +get_relative_coordinates_for_gantry_coordinates +(const CartesianCoordinate3D& coords) const +{ + // TODO: bed postion + return coords + get_vector_centre_of_first_ring_to_centre_of_gantry(); +} + END_NAMESPACE_STIR diff --git a/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx b/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx index 232e2a5121..e7bbe80c04 100644 --- a/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx +++ b/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx @@ -384,6 +384,19 @@ static inline int sign(const T& t) return t<0 ? -1 : 1; } +CartesianCoordinate3D +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 ProjDataInfo& proj_data_info) +{ + return density_info.get_index_coordinates_for_relative_coordinates + (proj_data_info.get_relative_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))); +} + // just do 1 LOR, returns true if lor is not empty static void ray_trace_one_lor(ProjMatrixElemsForOneBin& lor, @@ -605,11 +618,9 @@ calculate_proj_matrix_elems_for_one_bin( // find offset in z, taking into account if there are 1 or more LORs - // KT 20/06/2001 take origin.z() into account - // KT 15/05/2002 move +(max_index.z()+min_index.z())/2.F offset here instead of in formulas for Z1f,Z2f /* Here is how we find the offset of the first ray: - for only 1 ray, it is simply found by refering to the middle of the image - minus the origin.z(). + for only 1 ray, it is 0. + For multiple rays, the following reasoning is followed. First we look at oblique rays. @@ -666,11 +677,11 @@ calculate_proj_matrix_elems_for_one_bin( if (num_tangential_LORs == 1) { ray_trace_one_lor(lor, s_in_mm, m_in_mm, - cphi, sphi, tantheta, - offset_in_z, fovrad_in_mm, - voxel_size, - restrict_to_cylindrical_FOV, - num_lors_per_axial_pos); + cphi, sphi, tantheta, + offset_in_z, fovrad_in_mm, + voxel_size, + restrict_to_cylindrical_FOV, + num_lors_per_axial_pos); } else { @@ -687,11 +698,11 @@ calculate_proj_matrix_elems_for_one_bin( { ray_traced_lor.erase(); ray_trace_one_lor(ray_traced_lor, current_s_in_mm, m_in_mm, - cphi, sphi, tantheta, - offset_in_z, fovrad_in_mm, - voxel_size, - restrict_to_cylindrical_FOV, - num_lors_per_axial_pos*num_tangential_LORs); + cphi, sphi, tantheta, + offset_in_z, fovrad_in_mm, + voxel_size, + restrict_to_cylindrical_FOV, + num_lors_per_axial_pos*num_tangential_LORs); //std::cerr << "ray traced size " << ray_traced_lor.size() << std::endl; lor.merge(ray_traced_lor); } From 13708859966b09defc2745a4b334f36b8f33323b Mon Sep 17 00:00:00 2001 From: Ashley Gillman Date: Tue, 9 Oct 2018 14:01:32 +1000 Subject: [PATCH 5/8] [ENH] DiscDens Get vendor origin from ProjDataInfo, and use for LPS. --- src/buildblock/ProjDataInfo.cxx | 6 ++ src/buildblock/VoxelsOnCartesianGrid.cxx | 73 +++++++----------------- src/include/stir/DiscretisedDensity.h | 15 ++++- src/include/stir/DiscretisedDensity.inl | 58 ++++++++++++++----- src/include/stir/ProjDataInfo.h | 4 ++ src/include/stir/VoxelsOnCartesianGrid.h | 6 ++ 6 files changed, 95 insertions(+), 67 deletions(-) diff --git a/src/buildblock/ProjDataInfo.cxx b/src/buildblock/ProjDataInfo.cxx index 38995f23d2..d6f902a233 100644 --- a/src/buildblock/ProjDataInfo.cxx +++ b/src/buildblock/ProjDataInfo.cxx @@ -678,5 +678,11 @@ ProjDataInfo::get_bed_position() const (bed_position_horizontal, bed_position_vertical, 0); } +const CartesianCoordinate3D +ProjDataInfo:: +get_location_of_vendor_frame_of_reference_in_physical_coordinates() const { + return get_vector_centre_of_first_ring_to_centre_of_gantry(); +} + END_NAMESPACE_STIR diff --git a/src/buildblock/VoxelsOnCartesianGrid.cxx b/src/buildblock/VoxelsOnCartesianGrid.cxx index c36461108c..56d0463bfc 100644 --- a/src/buildblock/VoxelsOnCartesianGrid.cxx +++ b/src/buildblock/VoxelsOnCartesianGrid.cxx @@ -185,63 +185,29 @@ VoxelsOnCartesianGrid::VoxelsOnCartesianGrid(const ProjDataInfo& proj_dat const CartesianCoordinate3D& 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(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(ceil(FOVradius_in_mm / get_voxel_size().x())) + 1; - if (sizes.y()==-1) - y_size_used = 2*static_cast(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 -VoxelsOnCartesianGrid::VoxelsOnCartesianGrid(const shared_ptr < ExamInfo > & exam_info_sptr_v, - const ProjDataInfo& proj_data_info, - const float zoom, - const CartesianCoordinate3D& origin, - const CartesianCoordinate3D& sizes) +VoxelsOnCartesianGrid:: +VoxelsOnCartesianGrid(const shared_ptr& exam_info_sptr_v, + const ProjDataInfo& proj_data_info, + const float zoom, + const CartesianCoordinate3D& origin, + const CartesianCoordinate3D& 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 +void +VoxelsOnCartesianGrid:: +init_from_proj_data_info(const ProjDataInfo& proj_data_info, + const float zoom, + const CartesianCoordinate3D& origin, + const CartesianCoordinate3D& sizes) +{ int z_size = sizes.z(); // initialise to 0 to prevent compiler warnings //int z_size = 0; @@ -281,6 +247,11 @@ VoxelsOnCartesianGrid::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); } diff --git a/src/include/stir/DiscretisedDensity.h b/src/include/stir/DiscretisedDensity.h index de0eeca6f5..3a21002eb6 100644 --- a/src/include/stir/DiscretisedDensity.h +++ b/src/include/stir/DiscretisedDensity.h @@ -135,12 +135,20 @@ template const IndexRange& range, const CartesianCoordinate3D& origin); - //! Return the origin + //! Return the origin in physical coordinates inline const CartesianCoordinate3D& get_origin() const; - //! Set the origin + //! Set the origin, in physical coordinates inline void set_origin(const CartesianCoordinate3D &origin); + //! Return the vendor-defined origin in physical coordinates. Defaults to origin if unset. + // NB: could be in LPS, would this make more sense? + inline const CartesianCoordinate3D& get_vendor_origin() const; + + //! Set the vendor-defined origin, in physical coordinates + // NB: could be in LPS, would this make more sense? + inline void set_vendor_origin(const CartesianCoordinate3D &vendor_origin); + //! \name Translation between indices and physical coordinates /*! We distinguish between physical coordinates, relative coordinates (which are physical coordinates relative to the origin) and index coordinates (which run @@ -325,6 +333,9 @@ template private: CartesianCoordinate3D origin; + bool vendor_origin_set = false; + CartesianCoordinate3D vendor_origin; + static inline CartesianCoordinate3D swap_axes_based_on_orientation(const CartesianCoordinate3D& coordinates, const PatientPosition patient_position); diff --git a/src/include/stir/DiscretisedDensity.inl b/src/include/stir/DiscretisedDensity.inl index a15e494168..8071c84556 100644 --- a/src/include/stir/DiscretisedDensity.inl +++ b/src/include/stir/DiscretisedDensity.inl @@ -75,6 +75,21 @@ DiscretisedDensity:: get_origin() const { return origin; } +template +void +DiscretisedDensity:: +set_vendor_origin(const CartesianCoordinate3D &vendor_origin_v) +{ + vendor_origin = vendor_origin_v; + vendor_origin_set = true; +} + +template +const CartesianCoordinate3D& +DiscretisedDensity::get_vendor_origin() const +{ return vendor_origin_set ? vendor_origin : origin; } + + template bool @@ -252,28 +267,37 @@ swap_axes_based_on_orientation(const CartesianCoordinate3D& coords, return flip_coords; } +// The following two functions are the only ones that do actual novel work + template CartesianCoordinate3D DiscretisedDensity:: -get_LPS_coordinates_for_physical_coordinates(const CartesianCoordinate3D& coords) const +get_LPS_coordinates_for_indices +(const BasicCoordinate& indices) const { return swap_axes_based_on_orientation - (coords, this->get_exam_info().patient_position); + (get_relative_coordinates_for_indices(indices) + get_vendor_origin(), + get_exam_info().patient_position); } template -CartesianCoordinate3D +BasicCoordinate DiscretisedDensity:: -get_LPS_coordinates_for_indices(const BasicCoordinate& indices) const +get_index_coordinates_for_LPS_coordinates(const CartesianCoordinate3D& coords) const { - return this->get_LPS_coordinates_for_physical_coordinates - (this->get_physical_coordinates_for_indices(indices)); + return + this->actual_get_index_coordinates_for_relative_coordinates + (swap_axes_based_on_orientation(coords, get_exam_info().patient_position) + - this->get_vendor_origin()); } +// The rest of the LPS functions just leverage the last two + template CartesianCoordinate3D DiscretisedDensity:: -get_LPS_coordinates_for_indices(const BasicCoordinate& indices) const +get_LPS_coordinates_for_indices +(const BasicCoordinate& indices) const { return get_LPS_coordinates_for_indices(BasicCoordinate(indices)); } @@ -281,19 +305,25 @@ get_LPS_coordinates_for_indices(const BasicCoordinate& indi template CartesianCoordinate3D DiscretisedDensity:: -get_physical_coordinates_for_LPS_coordinates(const CartesianCoordinate3D& coords) const +get_LPS_coordinates_for_physical_coordinates +(const CartesianCoordinate3D& coords) const { - // operation is symmetric - return this->get_LPS_coordinates_for_physical_coordinates(coords); + // this could also be done by adding/subtracting origins and swapping + // but this is simpler to comprehend and means we only implement once above + return get_LPS_coordinates_for_indices + (get_index_coordinates_for_physical_coordinates(coords)); } template -BasicCoordinate +CartesianCoordinate3D DiscretisedDensity:: -get_index_coordinates_for_LPS_coordinates(const CartesianCoordinate3D& coords) const +get_physical_coordinates_for_LPS_coordinates +(const CartesianCoordinate3D& coords) const { - return this->get_index_coordinates_for_physical_coordinates - (this->get_physical_coordinates_for_LPS_coordinates(coords)); + // this could also be done by adding/subtracting origins and swapping + // but this is simpler to comprehend and means we only implement once above + return get_physical_coordinates_for_indices + (get_index_coordinates_for_LPS_coordinates(coords)); } template diff --git a/src/include/stir/ProjDataInfo.h b/src/include/stir/ProjDataInfo.h index 46e94f4c4b..26d36d05ac 100644 --- a/src/include/stir/ProjDataInfo.h +++ b/src/include/stir/ProjDataInfo.h @@ -384,6 +384,10 @@ class ProjDataInfo get_relative_coordinates_for_gantry_coordinates (const CartesianCoordinate3D& coords) const; + //! Get the location of the scanners frame-of-reference + const CartesianCoordinate3D + get_location_of_vendor_frame_of_reference_in_physical_coordinates() const; + protected: virtual bool blindly_equals(const root_type * const) const = 0; diff --git a/src/include/stir/VoxelsOnCartesianGrid.h b/src/include/stir/VoxelsOnCartesianGrid.h index 694deedfe6..fa27e5e377 100644 --- a/src/include/stir/VoxelsOnCartesianGrid.h +++ b/src/include/stir/VoxelsOnCartesianGrid.h @@ -190,6 +190,12 @@ void grow_z_range(const int min_z, const int max_z); //@} +private: + void + init_from_proj_data_info(const ProjDataInfo& proj_data_info, + const float zoom, + const CartesianCoordinate3D& offset, + const CartesianCoordinate3D& sizes); }; From f4285801484e06f3706f14c9758bcdfc8c34201b Mon Sep 17 00:00:00 2001 From: Ashley Gillman Date: Thu, 1 Nov 2018 11:48:22 +1000 Subject: [PATCH 6/8] [ENH] BasicCoordiante and CartesianCoordinateXD express to stringstream. BasicCoordinate in format [x1,x2,x3] CartesianCoordinate in format (x,y,z) i.e., reverse dim direction to BC --- src/buildblock/ProjDataInfo.cxx | 4 +++- src/include/stir/BasicCoordinate.h | 5 +++++ src/include/stir/BasicCoordinate.inl | 13 +++++++++++++ src/include/stir/CartesianCoordinate2D.h | 3 +++ src/include/stir/CartesianCoordinate2D.inl | 5 +++++ src/include/stir/CartesianCoordinate3D.h | 3 +++ src/include/stir/CartesianCoordinate3D.inl | 6 ++++++ 7 files changed, 38 insertions(+), 1 deletion(-) diff --git a/src/buildblock/ProjDataInfo.cxx b/src/buildblock/ProjDataInfo.cxx index d6f902a233..332be939fe 100644 --- a/src/buildblock/ProjDataInfo.cxx +++ b/src/buildblock/ProjDataInfo.cxx @@ -681,7 +681,9 @@ ProjDataInfo::get_bed_position() const const CartesianCoordinate3D ProjDataInfo:: get_location_of_vendor_frame_of_reference_in_physical_coordinates() const { - return get_vector_centre_of_first_ring_to_centre_of_gantry(); + std::cerr << "bed possition: " << get_bed_position() << std::endl; + return get_vector_centre_of_first_ring_to_centre_of_gantry() + + get_bed_position(); } END_NAMESPACE_STIR diff --git a/src/include/stir/BasicCoordinate.h b/src/include/stir/BasicCoordinate.h index fb8a8ead27..2b26ac9359 100644 --- a/src/include/stir/BasicCoordinate.h +++ b/src/include/stir/BasicCoordinate.h @@ -310,6 +310,11 @@ inline BasicCoordinate convert_int_to_float(const BasicCoordinate& cint); +//! BasicCoordinate can be converted to a string stream "[x1,x2,x3]" +/*! \ingroup Coordinate */ +template +std::ostream &operator<<(std::ostream &os, const BasicCoordinate& c); + END_NAMESPACE_STIR diff --git a/src/include/stir/BasicCoordinate.inl b/src/include/stir/BasicCoordinate.inl index 67a7ed871f..32e2f8431b 100644 --- a/src/include/stir/BasicCoordinate.inl +++ b/src/include/stir/BasicCoordinate.inl @@ -612,6 +612,19 @@ operator<(const BasicCoordinate& c1,const BasicCoordinat c1, c); \ } + +template +std::ostream &operator<<(std::ostream &os, const BasicCoordinate& c) { +{ + os << "["; + for (int i=1; i<=num_dimensions; i++) { + os << c[i]; + if (i }; +template +std::ostream &operator<<(std::ostream &os, const CartesianCoordinate2D &coord); + END_NAMESPACE_STIR #include "stir/CartesianCoordinate2D.inl" diff --git a/src/include/stir/CartesianCoordinate2D.inl b/src/include/stir/CartesianCoordinate2D.inl index 1010d98d74..aee6837822 100644 --- a/src/include/stir/CartesianCoordinate2D.inl +++ b/src/include/stir/CartesianCoordinate2D.inl @@ -92,5 +92,10 @@ CartesianCoordinate2D::x() const return this->operator[](2); } +template +std::ostream &operator<<(std::ostream &os, const CartesianCoordinate2D &coord) { + return os << "(x=" << coord.x() + << ",y=" << coord.y() << ")"; +} END_NAMESPACE_STIR diff --git a/src/include/stir/CartesianCoordinate3D.h b/src/include/stir/CartesianCoordinate3D.h index 2de936ad83..85260d3694 100644 --- a/src/include/stir/CartesianCoordinate3D.h +++ b/src/include/stir/CartesianCoordinate3D.h @@ -84,6 +84,9 @@ class CartesianCoordinate3D : public Coordinate3D }; +template +std::ostream &operator<<(std::ostream &os, const CartesianCoordinate3D &coord); + END_NAMESPACE_STIR #include "stir/CartesianCoordinate3D.inl" diff --git a/src/include/stir/CartesianCoordinate3D.inl b/src/include/stir/CartesianCoordinate3D.inl index 92a1061563..032797be50 100644 --- a/src/include/stir/CartesianCoordinate3D.inl +++ b/src/include/stir/CartesianCoordinate3D.inl @@ -131,5 +131,11 @@ CartesianCoordinate3D::x() const return this->operator[](3); } +template +std::ostream &operator<<(std::ostream &os, const CartesianCoordinate3D &coord) { + return os << "(x=" << coord.x() + << ",y=" << coord.y() + << ",z=" << coord.z() << ")"; +} END_NAMESPACE_STIR From 44b57cf8d981099a9515da232f59d41689a5f3d5 Mon Sep 17 00:00:00 2001 From: Ashley Gillman Date: Thu, 1 Nov 2018 16:38:35 +1000 Subject: [PATCH 7/8] [ENH] Account for bed position when finding vendor FoR. [FIX] Sign of correction to centre of gantry --- src/buildblock/ProjDataInfo.cxx | 13 ++++++++++--- src/buildblock/VoxelsOnCartesianGrid.cxx | 2 +- 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/src/buildblock/ProjDataInfo.cxx b/src/buildblock/ProjDataInfo.cxx index 332be939fe..78abd70384 100644 --- a/src/buildblock/ProjDataInfo.cxx +++ b/src/buildblock/ProjDataInfo.cxx @@ -675,14 +675,21 @@ CartesianCoordinate3D ProjDataInfo::get_bed_position() const { return CartesianCoordinate3D - (bed_position_horizontal, bed_position_vertical, 0); + (/*z=*/bed_position_horizontal, /*y=*/bed_position_vertical, /*x=*/0); } const CartesianCoordinate3D ProjDataInfo:: get_location_of_vendor_frame_of_reference_in_physical_coordinates() const { - std::cerr << "bed possition: " << get_bed_position() << std::endl; - return get_vector_centre_of_first_ring_to_centre_of_gantry() + 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(); } diff --git a/src/buildblock/VoxelsOnCartesianGrid.cxx b/src/buildblock/VoxelsOnCartesianGrid.cxx index 56d0463bfc..dea31c7653 100644 --- a/src/buildblock/VoxelsOnCartesianGrid.cxx +++ b/src/buildblock/VoxelsOnCartesianGrid.cxx @@ -250,7 +250,7 @@ init_from_proj_data_info(const ProjDataInfo& proj_data_info, this->set_origin(origin); this->set_vendor_origin (origin - + (proj_data_info + - (proj_data_info .get_location_of_vendor_frame_of_reference_in_physical_coordinates())); this->grow(range); From dd8bf3167599d9cca2a03a359b8b45c2afb5a967 Mon Sep 17 00:00:00 2001 From: Ashley Gillman Date: Thu, 1 Nov 2018 16:40:18 +1000 Subject: [PATCH 8/8] [ENH] Interfile images save vendor origin location. --- src/IO/InterfileHeader.cxx | 6 +++ src/IO/interfile.cxx | 47 ++++++++++++++++++-- src/include/stir/DiscretisedDensity.h | 2 - src/include/stir/IO/InterfileHeader.h | 1 + src/include/stir/IO/InterfileHeaderSiemens.h | 1 + src/include/stir/IO/interfile.h | 12 ++++- 6 files changed, 62 insertions(+), 7 deletions(-) diff --git a/src/IO/InterfileHeader.cxx b/src/IO/InterfileHeader.cxx index 2eff2674ad..d4acce35e2 100644 --- a/src/IO/InterfileHeader.cxx +++ b/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 first_pixel_offsets; + std::vector first_pixel_vendor_offsets; if (num_time_frames > 1 && num_image_data_types > 1) { diff --git a/src/IO/interfile.cxx b/src/IO/interfile.cxx index 22b466fb4b..e53d7b3712 100644 --- a/src/IO/interfile.cxx +++ b/src/IO/interfile.cxx @@ -131,12 +131,25 @@ create_image_and_header_from(InterfileImageHeader& hdr, - voxel_size * BasicCoordinate<3,float>(min_indices); } - return + VoxelsOnCartesianGrid* image = new VoxelsOnCartesianGrid (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 * @@ -531,7 +544,8 @@ write_basic_interfile_image_header(const string& header_file_name, const ByteOrder byte_order, const VectorWithOffset& scaling_factors, const VectorWithOffset& file_offsets, - const std::vector& data_type_descriptions) + const CartesianCoordinate3D& vendor_origin, + const std::vector& data_type_descriptions) { CartesianCoordinate3D min_indices; CartesianCoordinate3D 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 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 origin; origin.fill(static_cast(InterfileHeader::double_value_not_set)); + CartesianCoordinate3D vendor_origin; + vendor_origin.fill(static_cast(InterfileHeader::double_value_not_set)); return write_basic_interfile(filename, image, CartesianCoordinate3D(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& voxel_size, const CartesianCoordinate3D& origin, + const CartesianCoordinate3D& 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& voxel_size, const CartesianCoordinate3D& origin, + const CartesianCoordinate3D& 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& voxel_size, const CartesianCoordinate3D& origin, + const CartesianCoordinate3D& 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& voxel_size, const CartesianCoordinate3D& origin, + const CartesianCoordinate3D& 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& voxel_size, const CartesianCoordinate3D& origin, + const CartesianCoordinate3D& vendor_origin, const NumericType output_type, const float scale, const ByteOrder byte_order); diff --git a/src/include/stir/DiscretisedDensity.h b/src/include/stir/DiscretisedDensity.h index 3a21002eb6..3ceb977840 100644 --- a/src/include/stir/DiscretisedDensity.h +++ b/src/include/stir/DiscretisedDensity.h @@ -142,11 +142,9 @@ template inline void set_origin(const CartesianCoordinate3D &origin); //! Return the vendor-defined origin in physical coordinates. Defaults to origin if unset. - // NB: could be in LPS, would this make more sense? inline const CartesianCoordinate3D& get_vendor_origin() const; //! Set the vendor-defined origin, in physical coordinates - // NB: could be in LPS, would this make more sense? inline void set_vendor_origin(const CartesianCoordinate3D &vendor_origin); //! \name Translation between indices and physical coordinates diff --git a/src/include/stir/IO/InterfileHeader.h b/src/include/stir/IO/InterfileHeader.h index e6c1747000..ed33e5734e 100644 --- a/src/include/stir/IO/InterfileHeader.h +++ b/src/include/stir/IO/InterfileHeader.h @@ -187,6 +187,7 @@ class InterfileImageHeader : public InterfileHeader public: InterfileImageHeader(); std::vector first_pixel_offsets; + std::vector first_pixel_vendor_offsets; int num_image_data_types; std::vector index_nesting_level; std::vector image_data_type_description; diff --git a/src/include/stir/IO/InterfileHeaderSiemens.h b/src/include/stir/IO/InterfileHeaderSiemens.h index 08bfd27b06..26c7df0e3a 100644 --- a/src/include/stir/IO/InterfileHeaderSiemens.h +++ b/src/include/stir/IO/InterfileHeaderSiemens.h @@ -97,6 +97,7 @@ class InterfileImageHeader : public InterfileHeaderSiemens public: InterfileImageHeader(); std::vector first_pixel_offsets; + std::vector first_pixel_vendor_offsets; protected: virtual void read_matrix_info(); diff --git a/src/include/stir/IO/interfile.h b/src/include/stir/IO/interfile.h index 243ff631a0..fc6aedcdee 100644 --- a/src/include/stir/IO/interfile.h +++ b/src/include/stir/IO/interfile.h @@ -39,11 +39,17 @@ // has to include Succeeded.h (even if it doesn't use the return value). #include "stir/Succeeded.h" #include "stir/ByteOrder.h" +#include "stir/IO/InterfileHeader.h" #include #include START_NAMESPACE_STIR +#define CARTESIANCOORDINATE_NOTSET CartesianCoordinate3D \ + (InterfileHeader::double_value_not_set, \ + InterfileHeader::double_value_not_set, \ + InterfileHeader::double_value_not_set) + template class IndexRange; template class Array; template class DiscretisedDensity; @@ -149,7 +155,9 @@ write_basic_interfile_image_header(const std::string& header_file_name, const ByteOrder byte_order, const VectorWithOffset& scaling_factors, const VectorWithOffset& file_offsets, - const std::vector& data_type_descriptions = std::vector()); + const CartesianCoordinate3D& vendor_origin = CARTESIANCOORDINATE_NOTSET, + const std::vector& data_type_descriptions = std::vector() + ); //! a utility function that computes the file offsets of subsequent images @@ -176,6 +184,7 @@ write_basic_interfile(const std::string&filename, const Array<3,elemT>& image, const CartesianCoordinate3D& voxel_size, const CartesianCoordinate3D& origin, + const CartesianCoordinate3D& vendor_origin = CARTESIANCOORDINATE_NOTSET, const NumericType output_type = NumericType::FLOAT, const float scale= 0, const ByteOrder byte_order=ByteOrder::native); @@ -194,6 +203,7 @@ write_basic_interfile(const std::string&filename, const Array<3,elemT>& image, const CartesianCoordinate3D& voxel_size, const CartesianCoordinate3D& origin, + const CartesianCoordinate3D& vendor_origin = CARTESIANCOORDINATE_NOTSET, const NumericType output_type = NumericType::FLOAT, const float scale= 0, const ByteOrder byte_order=ByteOrder::native);