diff --git a/src/buildblock/ProjDataInfo.cxx b/src/buildblock/ProjDataInfo.cxx index 053157f8c4..de4afa0504 100644 --- a/src/buildblock/ProjDataInfo.cxx +++ b/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 +ProjDataInfo:: +get_bed_position() const +{ + return CartesianCoordinate3D + (bed_position_horizontal, bed_position_vertical, 0); +} + +CartesianCoordinate3D +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(middle_of_first_ring_to_middle_of_last / 2.F, 0, 0); +} + END_NAMESPACE_STIR diff --git a/src/buildblock/ProjDataInfoCylindrical.cxx b/src/buildblock/ProjDataInfoCylindrical.cxx index ef6bb8ee2e..cbaac57091 100644 --- a/src/buildblock/ProjDataInfoCylindrical.cxx +++ b/src/buildblock/ProjDataInfoCylindrical.cxx @@ -127,7 +127,8 @@ initialise_ring_diff_arrays() const VectorWithOffset(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) { @@ -135,6 +136,7 @@ initialise_ring_diff_arrays() const ((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 diff --git a/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx b/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx index 86d5cddee4..a79b971af5 100644 --- a/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx +++ b/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 cyl_coords(get_scanner_ptr()->get_effective_ring_radius()); cyl_coords.p1().psi() = static_cast((2.*_PI/num_detectors_per_ring)*(det1)); cyl_coords.p2().psi() = static_cast((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 lor(cyl_coords); coord_1 = lor.p1(); coord_2 = lor.p2(); diff --git a/src/buildblock/zoom.cxx b/src/buildblock/zoom.cxx index cd80afc760..549dda432a 100644 --- a/src/buildblock/zoom.cxx +++ b/src/buildblock/zoom.cxx @@ -280,6 +280,7 @@ construct_new_image_from_zoom_parameters(const VoxelsOnCartesianGrid &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 = diff --git a/src/experimental/recon_buildblock/DataSymmetriesForDensels_PET_CartesianGrid.cxx b/src/experimental/recon_buildblock/DataSymmetriesForDensels_PET_CartesianGrid.cxx index 79758c05b4..ea59a2a933 100644 --- a/src/experimental/recon_buildblock/DataSymmetriesForDensels_PET_CartesianGrid.cxx +++ b/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: ^ } } diff --git a/src/include/stir/DiscretisedDensity.h b/src/include/stir/DiscretisedDensity.h index 4c3a2f640f..df7f07edb0 100644 --- a/src/include/stir/DiscretisedDensity.h +++ b/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 @@ -253,6 +254,16 @@ template inline BasicCoordinate get_indices_closest_to_LPS_coordinates(const CartesianCoordinate3D& coords) const; + inline CartesianCoordinate3D + get_relative_coordinates_for_gantry_coordinates + (const CartesianCoordinate3D& coords, + const shared_ptr proj_data_info_sptr) const; + + inline BasicCoordinate + get_index_coordinates_for_gantry_coordinates + (const CartesianCoordinate3D& coords, + const shared_ptr proj_data_info_sptr) const; + //@} //! Allocate a new DiscretisedDensity object with same characteristics as the current one. diff --git a/src/include/stir/DiscretisedDensity.inl b/src/include/stir/DiscretisedDensity.inl index a15e494168..1e1f1bf2c5 100644 --- a/src/include/stir/DiscretisedDensity.inl +++ b/src/include/stir/DiscretisedDensity.inl @@ -304,4 +304,29 @@ get_indices_closest_to_LPS_coordinates(const CartesianCoordinate3D& coord return round(this->get_index_coordinates_for_LPS_coordinates(coords)); } +template +CartesianCoordinate3D +DiscretisedDensity:: +get_relative_coordinates_for_gantry_coordinates +(const CartesianCoordinate3D& coords, + const shared_ptr 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 +BasicCoordinate +DiscretisedDensity:: +get_index_coordinates_for_gantry_coordinates +(const CartesianCoordinate3D& coords, + const shared_ptr 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 diff --git a/src/include/stir/ProjDataInfo.h b/src/include/stir/ProjDataInfo.h index 2bd1712155..baee3bfeda 100644 --- a/src/include/stir/ProjDataInfo.h +++ b/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 @@ -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 get_bed_position() const; + + //! Vector from image frame of reference (centre of first ring) to gantry centre + CartesianCoordinate3D + get_centre_of_gantry_vector_in_relative_coordinates() const; + protected: virtual bool blindly_equals(const root_type * const) const = 0; diff --git a/src/include/stir/recon_buildblock/ProjMatrixByBinUsingRayTracing.h b/src/include/stir/recon_buildblock/ProjMatrixByBinUsingRayTracing.h index dface8d895..c23ea1f33d 100644 --- a/src/include/stir/recon_buildblock/ProjMatrixByBinUsingRayTracing.h +++ b/src/include/stir/recon_buildblock/ProjMatrixByBinUsingRayTracing.h @@ -199,7 +199,8 @@ public : CartesianCoordinate3D min_index; CartesianCoordinate3D max_index; - shared_ptr proj_data_info_ptr; + shared_ptr proj_data_info_sptr; + shared_ptr< DiscretisedDensity<3,float> > density_info_sptr; virtual void diff --git a/src/include/stir/zoom.h b/src/include/stir/zoom.h index c056e0bf91..f25f213808 100644 --- a/src/include/stir/zoom.h +++ b/src/include/stir/zoom.h @@ -164,7 +164,7 @@ zoom_image(const VoxelsOnCartesianGrid &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 diff --git a/src/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.cxx b/src/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.cxx index 2fcfbf9fe6..ed1a26a630 100644 --- a/src/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.cxx +++ b/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: ^ } } diff --git a/src/recon_buildblock/ProjMatrixByBinUsingInterpolation.cxx b/src/recon_buildblock/ProjMatrixByBinUsingInterpolation.cxx index fe4bf6e3b6..f10b63d793 100644 --- a/src/recon_buildblock/ProjMatrixByBinUsingInterpolation.cxx +++ b/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( diff --git a/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx b/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx index 8560fc0f42..1d8e313d63 100644 --- a/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx +++ b/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx @@ -267,33 +267,35 @@ static bool is_multiple(const float a, const float b) void ProjMatrixByBinUsingRayTracing:: set_up( - const shared_ptr& proj_data_info_ptr_v, - const shared_ptr >& density_info_ptr // TODO should be Info only + const shared_ptr& proj_data_info_sptr_v, + const shared_ptr >& density_info_sptr_v // TODO should be Info only ) { - ProjMatrixByBin::set_up(proj_data_info_ptr_v, density_info_ptr); + ProjMatrixByBin::set_up(proj_data_info_sptr_v, density_info_sptr_v); - proj_data_info_ptr= proj_data_info_ptr_v; - const VoxelsOnCartesianGrid * image_info_ptr = - dynamic_cast*> (density_info_ptr.get()); + proj_data_info_sptr = proj_data_info_sptr_v; + density_info_sptr = density_info_sptr_v; - if (image_info_ptr == NULL) + const VoxelsOnCartesianGrid * image_info_sptr = + dynamic_cast*> (density_info_sptr.get()); + + if (image_info_sptr == NULL) error("ProjMatrixByBinUsingRayTracing initialised with a wrong type of DiscretisedDensity\n"); - voxel_size = image_info_ptr->get_voxel_size(); - origin = image_info_ptr->get_origin(); - image_info_ptr->get_regular_range(min_index, max_index); + voxel_size = image_info_sptr->get_voxel_size(); + origin = image_info_sptr->get_origin(); + image_info_sptr->get_regular_range(min_index, max_index); symmetries_sptr.reset( - new DataSymmetriesForBins_PET_CartesianGrid(proj_data_info_ptr, - density_info_ptr, + new DataSymmetriesForBins_PET_CartesianGrid(proj_data_info_sptr, + density_info_sptr, do_symmetry_90degrees_min_phi, do_symmetry_180degrees_min_phi, do_symmetry_swap_segment, do_symmetry_swap_s, do_symmetry_shift_z)); const float sampling_distance_of_adjacent_LORs_xy = - proj_data_info_ptr->get_sampling_in_s(Bin(0,0,0,0)); + proj_data_info_sptr->get_sampling_in_s(Bin(0,0,0,0)); if(sampling_distance_of_adjacent_LORs_xy/num_tangential_LORs > voxel_size.x() + 1.E-3 || sampling_distance_of_adjacent_LORs_xy/num_tangential_LORs > voxel_size.y() + 1.E-3) @@ -308,8 +310,9 @@ set_up( if (use_actual_detector_boundaries) { + // Is this a memory leak? const ProjDataInfoCylindricalNoArcCorr * proj_data_info_cyl_ptr = - dynamic_cast(proj_data_info_ptr.get()); + dynamic_cast(proj_data_info_sptr.get()); if (proj_data_info_cyl_ptr== 0) { warning("ProjMatrixByBinUsingRayTracing: use_actual_detector_boundaries" @@ -344,6 +347,7 @@ set_up( // test if our 2D code does not have problems { // currently 2D code relies on the LOR falling in the middle of a voxel (in z-direction) + // If reactivated, update to use proper index -> physical space conversion const float z_shift = - origin.z()/voxel_size.z() +(max_index.z()+min_index.z())/2.F; if (fabs(z_shift - round(z_shift)) > .01) @@ -387,21 +391,22 @@ 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 tantheta, const float offset_in_z, const float fovrad_in_mm, const CartesianCoordinate3D& voxel_size, const bool restrict_to_cylindrical_FOV, - const int num_LORs) + const int num_LORs, + const shared_ptr< DiscretisedDensity<3,float> >& + density_info_sptr, + const shared_ptr proj_data_info_sptr) { assert(lor.size() == 0); /* Find Intersection points of LOR and image FOV (assuming infinitely long scanner)*/ /* (in voxel units) */ - CartesianCoordinate3D start_point; - CartesianCoordinate3D stop_point; { /* parametrisation of LOR is X= s*cphi + a*sphi, @@ -452,18 +457,54 @@ ray_trace_one_lor(ProjMatrixElemsForOneBin& lor, (fovrad_in_mm*sign(cphi) + s_in_mm*sphi)/cphi); min_a = max((-fovrad_in_mm*sign(sphi) - s_in_mm*cphi)/sphi, (-fovrad_in_mm*sign(cphi) + s_in_mm*sphi)/cphi); + if (min_a > max_a - 1.E-3*voxel_size.x()) return; } } //!restrict_to_cylindrical_FOV - 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(); - 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(); + + // start_point_gantry = (s_in_mm*cphi + max_a*sphi); + // start_point_bed = tx_gantry_to_bed(start_point_gantry); + // start_point_vx = tx_bed_to_idx(start_point_bed); + + + CartesianCoordinate3D start_point; + CartesianCoordinate3D stop_point; + + // // 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() = (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() = (m_in_mm+offset_in_z - min_a*tantheta)/voxel_size.z(); + // std::cerr << "LOR reference " << start_point << "->" << stop_point << std::endl; + + // start and stop point in gantry coordinates + // These are in mm, where (0, 0, 0) is the centre of the gantry + start_point.x() = s_in_mm*cphi + max_a*sphi; + start_point.y() = s_in_mm*sphi - max_a*cphi; + start_point.z() = m_in_mm+offset_in_z - max_a*tantheta; + stop_point.x() = s_in_mm*cphi + min_a*sphi; + stop_point.y() = s_in_mm*sphi - min_a*cphi; + stop_point.z() = m_in_mm+offset_in_z - min_a*tantheta; + // std::cerr << "LOR gantry " << start_point << "->" << stop_point << std::endl; + + // Shift to bed coordinates + start_point = density_info_sptr + ->get_index_coordinates_for_gantry_coordinates(start_point, proj_data_info_sptr); + stop_point = density_info_sptr + ->get_index_coordinates_for_gantry_coordinates(stop_point, proj_data_info_sptr); + // std::cerr << "LOR bedpos " << start_point << "->" << stop_point << std::endl; + + // // Shift to Image Index coordinates + // start_point = density_info_sptr + // ->get_index_coordinates_for_physical_coordinates(start_point); + // stop_point = density_info_sptr + // ->get_index_coordinates_for_physical_coordinates(stop_point); + // // std::cerr << "LOR index " << start_point << "->" << stop_point << std::endl; #if 0 // KT 18/05/2005 this is no longer necessary @@ -538,39 +579,24 @@ calculate_proj_matrix_elems_for_one_bin( } const Bin bin = lor.get_bin(); - assert(bin.segment_num() >= proj_data_info_ptr->get_min_segment_num()); - assert(bin.segment_num() <= proj_data_info_ptr->get_max_segment_num()); + assert(bin.segment_num() >= proj_data_info_sptr->get_min_segment_num()); + assert(bin.segment_num() <= proj_data_info_sptr->get_max_segment_num()); assert(lor.size() == 0); - float phi; - 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 + float phi = proj_data_info_sptr->get_phi(bin); + float s_in_mm = proj_data_info_sptr->get_s(bin); + + if (use_actual_detector_boundaries) { // can be static_cast later on const ProjDataInfoCylindricalNoArcCorr& proj_data_info_noarccor = - dynamic_cast(*proj_data_info_ptr); + dynamic_cast(*proj_data_info_sptr); // TODO check on 180 degrees for views const int num_detectors = - proj_data_info_ptr->get_scanner_ptr()->get_num_detectors_per_ring(); + proj_data_info_sptr->get_scanner_sptr()->get_num_detectors_per_ring(); const float ring_radius = - proj_data_info_ptr->get_scanner_ptr()->get_effective_ring_radius(); + proj_data_info_sptr->get_scanner_sptr()->get_effective_ring_radius(); int det_num1=0, det_num2=0; proj_data_info_noarccor. @@ -578,28 +604,27 @@ 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) + if (fabs(s_in_mm-old_s_in_mm)>proj_data_info_sptr->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); const float sphi = sin(phi); - const float tantheta = proj_data_info_ptr->get_tantheta(bin); + const float tantheta = proj_data_info_sptr->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_sptr->get_m(bin); const float sampling_distance_of_adjacent_LORs_z = - proj_data_info_ptr->get_sampling_in_t(bin)/costheta; - + proj_data_info_sptr->get_sampling_in_t(bin)/costheta; // find number of LORs we have to take, such that we don't miss voxels // we have to subtract a tiny amount from the quotient, to avoid having too many LORs @@ -646,20 +671,17 @@ calculate_proj_matrix_elems_for_one_bin( */ const float z_position_of_first_LOR_wrt_centre_of_TOR = (-sampling_distance_of_adjacent_LORs_z/(2*num_lors_per_axial_pos)* - (num_lors_per_axial_pos-1)) - - origin.z(); - float offset_in_z = - z_position_of_first_LOR_wrt_centre_of_TOR - +(max_index.z()+min_index.z())/2.F * voxel_size.z(); + (num_lors_per_axial_pos-1)); + float offset_in_z = z_position_of_first_LOR_wrt_centre_of_TOR; if (tantheta==0) { // 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(); } @@ -667,23 +689,32 @@ calculate_proj_matrix_elems_for_one_bin( // use FOV which is slightly 'inside' the image to avoid // index out of range #ifdef STIR_PMRT_LARGER_FOV - const float fovrad_in_mm = - min((min(max_index.x(), -min_index.x())+.45F)*voxel_size.x(), - (min(max_index.y(), -min_index.y())+.45F)*voxel_size.y()); + CartesianCoordinate3D max_pos = + density_info_sptr + ->get_relative_coordinates_for_indices(max_index + 0.45F); + CartesianCoordinate3D min_pos = + density_info_sptr + ->get_relative_coordinates_for_indices(min_index - 0.45F); #else - const float fovrad_in_mm = - min((min(max_index.x(), -min_index.x()))*voxel_size.x(), - (min(max_index.y(), -min_index.y()))*voxel_size.y()); + CartesianCoordinate3D max_pos = + density_info_sptr + ->get_relative_coordinates_for_indices(max_index); + CartesianCoordinate3D min_pos = + density_info_sptr + ->get_relative_coordinates_for_indices(min_index); #endif + const float fovrad_in_mm = + min(min(max_pos.x(), -min_pos.x()), + min(max_pos.y(), -min_pos.y())); if (num_tangential_LORs == 1) { - ray_trace_one_lor(lor, s_in_mm, t_in_mm, - cphi, sphi, costheta, tantheta, - offset_in_z, fovrad_in_mm, - voxel_size, - restrict_to_cylindrical_FOV, - num_lors_per_axial_pos); + 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, density_info_sptr, proj_data_info_sptr); } else { @@ -693,18 +724,19 @@ calculate_proj_matrix_elems_for_one_bin( // interleaved case has a sampling which is twice as high const float s_inc = (!use_actual_detector_boundaries ? 1 : 2) * - proj_data_info_ptr->get_sampling_in_s(bin)/num_tangential_LORs; + proj_data_info_sptr->get_sampling_in_s(bin)/num_tangential_LORs; float current_s_in_mm = s_in_mm - s_inc*(num_tangential_LORs-1)/2.F; 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, - cphi, sphi, costheta, tantheta, - offset_in_z, fovrad_in_mm, - voxel_size, - restrict_to_cylindrical_FOV, - num_lors_per_axial_pos*num_tangential_LORs); + 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, + density_info_sptr, proj_data_info_sptr); //std::cerr << "ray traced size " << ray_traced_lor.size() << std::endl; lor.merge(ray_traced_lor); } @@ -714,19 +746,38 @@ calculate_proj_matrix_elems_for_one_bin( if (lor.size()>0) { if (tantheta==0 ) - { - const float z_of_first_voxel= - lor.begin()->coord1() + - 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 - )/voxel_size.z(); - const float right_edge_of_TOR = - (t_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); + { + CartesianCoordinate3D + first_point_of_first_lor(lor.begin()->coord1(), + lor.begin()->coord2(), + lor.begin()->coord3()); + // std::cerr << "first_point_of_first_lor: " << first_point_of_first_lor << std::endl; + const float z_of_first_voxel_in_index_space + = first_point_of_first_lor.z(); + + // In gantry space + CartesianCoordinate3D + start_of_tor_projected_to_axis_without_offset + (m_in_mm - sampling_distance_of_adjacent_LORs_z/2, 0, 0); + CartesianCoordinate3D + end_of_tor_projected_to_axis_without_offset + (m_in_mm + sampling_distance_of_adjacent_LORs_z/2, 0, 0); + // Shift to Image Index coordinates + start_of_tor_projected_to_axis_without_offset = density_info_sptr + ->get_index_coordinates_for_gantry_coordinates + (start_of_tor_projected_to_axis_without_offset, proj_data_info_sptr); + end_of_tor_projected_to_axis_without_offset = density_info_sptr + ->get_index_coordinates_for_gantry_coordinates + (end_of_tor_projected_to_axis_without_offset, proj_data_info_sptr); + + const float first_lor_relative_to_tor + = z_of_first_voxel_in_index_space + - start_of_tor_projected_to_axis_without_offset.z(); + const float end_of_tor_relative_to_tor + = end_of_tor_projected_to_axis_without_offset.z() + - start_of_tor_projected_to_axis_without_offset.z(); + + add_adjacent_z(lor, first_lor_relative_to_tor, end_of_tor_relative_to_tor); } else if (num_lors_per_axial_pos>1) { @@ -843,6 +894,7 @@ add_adjacent_z(ProjMatrixElemsForOneBin& lor, const float overlap_of_voxel_with_TOR = std::min(right_edge_of_TOR, z_of_first_voxel + .5F) - std::max(0.F, z_of_first_voxel - .5F); + // std::cerr << "overlap_of_voxel_with_TOR: " << overlap_of_voxel_with_TOR << std::endl; assert (overlap_of_voxel_with_TOR>0); assert(overlap_of_voxel_with_TOR < 1.0001); if (overlap_of_voxel_with_TOR<.9999) // test if it is 1 diff --git a/src/scatter_buildblock/ScatterEstimationByBin.cxx b/src/scatter_buildblock/ScatterEstimationByBin.cxx index d654abde91..c0f214927b 100644 --- a/src/scatter_buildblock/ScatterEstimationByBin.cxx +++ b/src/scatter_buildblock/ScatterEstimationByBin.cxx @@ -302,6 +302,7 @@ process_data() #endif this->shift_detector_coordinates_to_origin = CartesianCoordinate3D(this->proj_data_info_ptr->get_m(Bin(0,0,0,0)),0, 0); + // ^ ORIGINTODO float total_scatter = 0 ; diff --git a/src/scatter_buildblock/sample_scatter_points.cxx b/src/scatter_buildblock/sample_scatter_points.cxx index efa28dd143..7f481659f8 100644 --- a/src/scatter_buildblock/sample_scatter_points.cxx +++ b/src/scatter_buildblock/sample_scatter_points.cxx @@ -64,6 +64,7 @@ sample_scatter_points() // TODO use class function once it exists const float z_to_middle = (image.get_max_index() + image.get_min_index())*voxel_size.z()/2.F; + // ORIGINTODO: ^ origin.z() -= z_to_middle; this->scatter_volume = voxel_size[1]*voxel_size[2]*voxel_size[3]; diff --git a/src/scatter_buildblock/scatter_detection_modelling.cxx b/src/scatter_buildblock/scatter_detection_modelling.cxx index dbecbe0ce5..72fdcbb850 100644 --- a/src/scatter_buildblock/scatter_detection_modelling.cxx +++ b/src/scatter_buildblock/scatter_detection_modelling.cxx @@ -73,6 +73,7 @@ find_detectors(unsigned& det_num_A, unsigned& det_num_B, const Bin& bin) const det_num_A = this->find_in_detection_points_vector(detector_coord_A + this->shift_detector_coordinates_to_origin); + // ORIGINTODO: ^ det_num_B = this->find_in_detection_points_vector(detector_coord_B + this->shift_detector_coordinates_to_origin); diff --git a/src/scatter_buildblock/single_scatter_integrals.cxx b/src/scatter_buildblock/single_scatter_integrals.cxx index 4410f308b1..5f0d7e81f0 100644 --- a/src/scatter_buildblock/single_scatter_integrals.cxx +++ b/src/scatter_buildblock/single_scatter_integrals.cxx @@ -100,6 +100,7 @@ integral_between_2_points(const DiscretisedDensity<3,float>& density, const float z_to_middle = (image.get_max_index() + image.get_min_index())*voxel_size.z()/2.F; origin.z() -= z_to_middle; + // ORIGINTODO: ^ /* TODO replace with image.get_index_coordinates_for_physical_coordinates */ ProjMatrixElemsForOneBin lor; RayTraceVoxelsOnCartesianGrid(lor,