diff --git a/src/analytic/FBP3DRP/FBP3DRPReconstruction.cxx b/src/analytic/FBP3DRP/FBP3DRPReconstruction.cxx index e531117129..e738ebfad7 100644 --- a/src/analytic/FBP3DRP/FBP3DRPReconstruction.cxx +++ b/src/analytic/FBP3DRP/FBP3DRPReconstruction.cxx @@ -5,11 +5,13 @@ \author Kris Thielemans \author Claire LABBE \author PARAPET project + \author Ashley Gillman */ /* Copyright (C) 2000 PARAPET partners Copyright (C) 2000- 2012, Hammersmith Imanet Ltd Copyright (C) 2020, University College London + Copyright (C) 2020, CSIRO This file is part of STIR. SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license @@ -19,6 +21,10 @@ /* Modification history: (highlights in anti-chronological order) + AG 2020/07/21 + - We no longer want to force images to be centre-aligned, rewrote implementation + of find_rmin_rmax(). + KT Oct 2004 - no longer use Numerical Recipes fourier - option to 'stretch' the colsher filter during definition for better results @@ -83,8 +89,9 @@ #include "stir/analytic/FBP3DRP/ColsherFilter.h" #include "stir/display.h" -// #include "stir/recon_buildblock/distributable.h" -// #include "stir/FBP3DRP/process_viewgrams.h" +//#include "stir/info.h" +//#include "stir/recon_buildblock/distributable.h" +//#include "stir/FBP3DRP/process_viewgrams.h" #include "stir/analytic/FBP3DRP/FBP3DRPReconstruction.h" #include "stir/analytic/FBP2D/FBP2DReconstruction.h" @@ -124,50 +131,47 @@ find_rmin_rmax(int& rmin, const int seg_num, const VoxelsOnCartesianGrid& image) { - + // precomupute a few values: + // Radius of the FOV as per the ProjDataInfo, in mm const float fovrad = proj_data_info_cyl.get_s(Bin(0, 0, 0, proj_data_info_cyl.get_num_tangential_poss() / 2 - 1)); - // Compute minimum and maximum rings of 'missing' projections - + // Radius of the detector rings in mm + const float ringrad = proj_data_info_cyl.get_ring_radius(); + // delta, ring difference in units of n.rings const float delta = proj_data_info_cyl.get_average_ring_difference(seg_num); - - // find correspondence between ring coordinates and image coordinates: - // z = num_planes_per_virtual_ring * ring + virtual_ring_offset - // compute the offset by matching up the centre of the scanner - // in the 2 coordinate systems - // TODO get all this from ProjDataInfo or so - - const int num_planes_per_virtual_ring - = (proj_data_info_cyl.get_max_ring_difference(seg_num) == proj_data_info_cyl.get_min_ring_difference(seg_num)) ? 2 : 1; + // if span > 1, we essentially have twice as many "virtual rings" const int num_virtual_rings_per_physical_ring = (proj_data_info_cyl.get_max_ring_difference(seg_num) == proj_data_info_cyl.get_min_ring_difference(seg_num)) ? 1 : 2; - const float virtual_ring_offset - = (image.get_max_z() + image.get_min_z()) / 2.F - - num_planes_per_virtual_ring - * (proj_data_info_cyl.get_max_axial_pos_num(seg_num) + num_virtual_rings_per_physical_ring * delta - + proj_data_info_cyl.get_min_axial_pos_num(seg_num)) - / 2; - - // we first consider the LOR at s=0, phi=0 which goes through z=0,y=0, x=fovrad - // later on, we will shift it to the 'left'most edge of the FOV. - - // find z position of intersection of this LOR with the detector radius - // (i.e. y=0, x=-ring_radius) - // use image coordinates first - float z_in_image_coordinates = -fabs(delta) * num_planes_per_virtual_ring * num_virtual_rings_per_physical_ring - * (fovrad + proj_data_info_cyl.get_ring_radius()) / (2 * proj_data_info_cyl.get_ring_radius()); - // now shift it to the edge of the FOV - // (taking into account that z==get_min_z() is in the middle of the voxel) - z_in_image_coordinates += image.get_min_z() - .5F; - - // now convert to virtual_ring_coordinates using z = num_planes_per_virtual_ring * ring + virtual_ring_offset - const float z_in_virtual_ring_coordinates = (z_in_image_coordinates - virtual_ring_offset) / num_planes_per_virtual_ring; + // This can be derived graphically by drawing a line from a top view of + // a scanner from ring to ring at angle theta (based on the segment) + // that just skims the field-of-view at the end of the ring. + // (Assume phi = y = 0, so no out-of-plane elements). + // The width of this line is delta (in units of n.rings), and height + // ringrad, and this triangle is similar to the smaller portion that + // still lies within the rings, which has width (delta - z) (where z + // is the overhang from the rings in units of n.rings) and height + // (ringrad - fovrad) / 2 + const float lor_overhang_in_num_rings = fabs(delta) / 2 * (1 + fovrad / ringrad); + const float lor_overhang_in_num_virtual_rings = lor_overhang_in_num_rings * num_virtual_rings_per_physical_ring; + + // rmin should just be this number of virtual rings backward from 0. + rmin = static_cast(floor(-lor_overhang_in_num_virtual_rings)); + // and rmax can be calculated symmetrically (we assume here the axial positions + // are correctly centred, this won't work otherwise.) + rmax = proj_data_info_cyl.get_max_axial_pos_num(seg_num) + (proj_data_info_cyl.get_min_axial_pos_num(seg_num) - rmin); - // finally find the 'ring' number - rmin = static_cast(floor(z_in_virtual_ring_coordinates)); + // lastly, if we were, for some reason, provided a larger-than-standard + // sinogram, make sure we use at least that size. + if (proj_data_info_cyl.get_min_axial_pos_num(seg_num) < rmin) + { + rmin = proj_data_info_cyl.get_min_axial_pos_num(seg_num); + rmax = proj_data_info_cyl.get_max_axial_pos_num(seg_num); + } - // rmax is determined by using symmetry: at both ends there are just as many missing rings - rmax = proj_data_info_cyl.get_max_axial_pos_num(seg_num) + (proj_data_info_cyl.get_min_axial_pos_num(seg_num) - rmin); + // info(boost::format("seg: %s : amin: %s, amax: %s") + // % seg_num % proj_data_info_cyl.get_min_axial_pos_num(seg_num) % proj_data_info_cyl.get_max_axial_pos_num(seg_num)); + // info(boost::format("seg: %s : rmin: %s, rmax: %s") + // % seg_num % rmin % rmax); } const char* const FBP3DRPReconstruction::registered_name = "FBP3DRP"; diff --git a/src/buildblock/ProjDataInfo.cxx b/src/buildblock/ProjDataInfo.cxx index b72804f35b..6bf46996de 100644 --- a/src/buildblock/ProjDataInfo.cxx +++ b/src/buildblock/ProjDataInfo.cxx @@ -801,4 +801,26 @@ ProjDataInfo::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/buildblock/ProjDataInfoBlocksOnCylindricalNoArcCorr.cxx b/src/buildblock/ProjDataInfoBlocksOnCylindricalNoArcCorr.cxx index 6c2db040e4..811923b5cb 100644 --- a/src/buildblock/ProjDataInfoBlocksOnCylindricalNoArcCorr.cxx +++ b/src/buildblock/ProjDataInfoBlocksOnCylindricalNoArcCorr.cxx @@ -100,8 +100,8 @@ ProjDataInfoBlocksOnCylindricalNoArcCorr::find_scanner_coordinates_given_cartesi DetectionPosition<> det_pos1; DetectionPosition<> det_pos2; - if (get_scanner_ptr()->find_detection_position_given_cartesian_coordinate(det_pos1, c1 + this->z_shift) == Succeeded::no - || get_scanner_ptr()->find_detection_position_given_cartesian_coordinate(det_pos2, c2 + this->z_shift) == Succeeded::no) + if (get_scanner_ptr()->find_detection_position_given_cartesian_coordinate(det_pos1, c1) == Succeeded::no + || get_scanner_ptr()->find_detection_position_given_cartesian_coordinate(det_pos2, c2) == Succeeded::no) { return Succeeded::no; } diff --git a/src/buildblock/ProjDataInfoCylindrical.cxx b/src/buildblock/ProjDataInfoCylindrical.cxx index 412e24c63b..7a7f9c2290 100644 --- a/src/buildblock/ProjDataInfoCylindrical.cxx +++ b/src/buildblock/ProjDataInfoCylindrical.cxx @@ -147,6 +147,7 @@ ProjDataInfoCylindrical::initialise_ring_diff_arrays() const /* m_offsets are found by requiring get_m(..., min_axial_pos_num,...) == - get_m(..., max_axial_pos_num,...) */ + for (int segment_num = get_min_segment_num(); segment_num <= get_max_segment_num(); ++segment_num) { m_offset[segment_num] @@ -283,10 +284,7 @@ ProjDataInfoCylindrical::initialise_ring_diff_arrays() const + (get_scanner_ptr()->get_num_rings() - 1); const int ring1_plus_ring2 = round(ring1_plus_ring2_float); // check that it was integer - if (get_scanner_sptr()->get_scanner_geometry() == "Cylindrical") - { - assert(fabs(ring1_plus_ring2 - ring1_plus_ring2_float) < 1E-4); - } + assert(fabs(ring1_plus_ring2 - ring1_plus_ring2_float) < 1E-4); segment_axial_pos_to_ring1_plus_ring2[s_num][ax_pos_num] = ring1_plus_ring2; } } diff --git a/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx b/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx index 8932574c74..47721803a1 100644 --- a/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx +++ b/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx @@ -357,6 +357,7 @@ ProjDataInfoCylindricalNoArcCorr::get_all_det_pos_pairs_for_bin(vector& c1, const CartesianCoordinate3D& c2) @@ -399,14 +400,23 @@ ProjDataInfoCylindricalNoArcCorr::find_scanner_coordinates_given_cartesian_coord ring1 = round(coord_det1.z()/ring_spacing); ring2 = round(coord_det2.z()/ring_spacing); #else - LORInCylinderCoordinates cyl_coords; - if (find_LOR_intersections_with_cylinder(cyl_coords, LORAs2Points(c1, c2), ring_radius) == Succeeded::no) + + // here we define an internal-only coord system called ring coords + // They just define z=0 as first ring, for ease of calculating ring no. + CartesianCoordinate3D offset_gantry_coords_to_ring_coords = get_vector_centre_of_first_ring_to_centre_of_gantry(); + + LORInCylinderCoordinates cyl_in_ring_coords; + if (find_LOR_intersections_with_cylinder( + cyl_in_ring_coords, + LORAs2Points(c1 + offset_gantry_coords_to_ring_coords, c2 + offset_gantry_coords_to_ring_coords), + ring_radius) + == Succeeded::no) return Succeeded::no; - det1 = modulo(round((cyl_coords.p1().psi() - this->get_psi_offset()) / (2. * _PI / num_detectors)), num_detectors); - det2 = modulo(round((cyl_coords.p2().psi() - this->get_psi_offset()) / (2. * _PI / num_detectors)), num_detectors); - ring1 = round(cyl_coords.p1().z() / ring_spacing); - ring2 = round(cyl_coords.p2().z() / ring_spacing); + det1 = modulo(round((cyl_in_ring_coords.p1().psi() - this->get_psi_offset()) / (2. * _PI / num_detectors)), num_detectors); + det2 = modulo(round((cyl_in_ring_coords.p2().psi() - this->get_psi_offset()) / (2. * _PI / num_detectors)), num_detectors); + ring1 = round((cyl_in_ring_coords.p1()).z() / ring_spacing); + ring2 = round((cyl_in_ring_coords.p2()).z() / ring_spacing); #endif @@ -427,10 +437,6 @@ ProjDataInfoCylindricalNoArcCorr::find_cartesian_coordinates_of_detection(Cartes DetectionPositionPair<> dpp; get_det_pos_pair_for_bin(dpp, bin); - /* TODO - best to use Scanner::get_coordinate_for_det_pos(). - Sadly, the latter is not yet implemented for Cylindrical scanners. - */ // find corresponding cartesian coordinates find_cartesian_coordinates_given_scanner_coordinates(coord_1, coord_2, @@ -497,8 +503,9 @@ ProjDataInfoCylindricalNoArcCorr::find_cartesian_coordinates_given_scanner_coord cyl_coords.p1().z() = r1 * get_scanner_ptr()->get_ring_spacing(); cyl_coords.p2().z() = r2 * get_scanner_ptr()->get_ring_spacing(); LORAs2Points lor(cyl_coords); - coord_1 = lor.p1(); - coord_2 = lor.p2(); + CartesianCoordinate3D offset_gantry_coords_to_ring_coords = get_vector_centre_of_first_ring_to_centre_of_gantry(); + coord_1 = lor.p1() - offset_gantry_coords_to_ring_coords; + coord_2 = lor.p2() - offset_gantry_coords_to_ring_coords; #endif if (tpos < 0) @@ -538,8 +545,8 @@ ProjDataInfoCylindricalNoArcCorr::get_bin(const LOR& lor, const double de Bin bin; #ifndef STIR_DEVEL // find nearest bin by going to nearest detectors first - LORInCylinderCoordinates cyl_coords; - if (lor.change_representation(cyl_coords, get_ring_radius()) == Succeeded::no) + LORInCylinderCoordinates cyl_in_gantry_coords; + if (lor.change_representation(cyl_in_gantry_coords, get_ring_radius()) == Succeeded::no) { bin.set_bin_value(-1); return bin; @@ -547,19 +554,20 @@ ProjDataInfoCylindricalNoArcCorr::get_bin(const LOR& lor, const double de const int num_detectors_per_ring = get_scanner_ptr()->get_num_detectors_per_ring(); const int num_rings = get_scanner_ptr()->get_num_rings(); - const int det1 = modulo(round((cyl_coords.p1().psi() - this->get_psi_offset()) / (2. * _PI / num_detectors_per_ring)), + const int det1 = modulo(round((cyl_in_gantry_coords.p1().psi() - this->get_psi_offset()) / (2. * _PI / num_detectors_per_ring)), num_detectors_per_ring); - const int det2 = modulo(round((cyl_coords.p2().psi() - this->get_psi_offset()) / (2. * _PI / num_detectors_per_ring)), + const int det2 = modulo(round((cyl_in_gantry_coords.p2().psi() - this->get_psi_offset()) / (2. * _PI / num_detectors_per_ring)), num_detectors_per_ring); // TODO WARNING LOR coordinates are w.r.t. centre of scanner, but the rings are numbered with the first ring at 0 - const int ring1 = round(cyl_coords.p1().z() / get_ring_spacing() + (num_rings - 1) / 2.F); - const int ring2 = round(cyl_coords.p2().z() / get_ring_spacing() + (num_rings - 1) / 2.F); + const int ring1 = round(cyl_in_gantry_coords.p1().z() / get_ring_spacing() + (num_rings - 1) / 2.F); + const int ring2 = round(cyl_in_gantry_coords.p2().z() / get_ring_spacing() + (num_rings - 1) / 2.F); assert(det1 >= 0 && det1 < num_detectors_per_ring); assert(det2 >= 0 && det2 < num_detectors_per_ring); if (ring1 >= 0 && ring1 < num_rings && ring2 >= 0 && ring2 < num_rings - && get_bin_for_det_pair(bin, det1, ring1, det2, ring2, (cyl_coords.is_swapped() ? -1 : 1) * get_tof_bin(delta_time)) + && get_bin_for_det_pair( + bin, det1, ring1, det2, ring2, (cyl_in_gantry_coords.is_swapped() ? -1 : 1) * get_tof_bin(delta_time)) == Succeeded::yes && bin.tangential_pos_num() >= get_min_tangential_pos_num() && bin.tangential_pos_num() <= get_max_tangential_pos_num()) { diff --git a/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx b/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx index b50f5e1040..e770102f43 100644 --- a/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx +++ b/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx @@ -54,9 +54,9 @@ ProjDataInfoGenericNoArcCorr::ProjDataInfoGenericNoArcCorr(const shared_ptrz_shift.z() = this->get_scanner_ptr()->get_coordinate_for_det_pos(DetectionPosition<>(0, 0, 0)).z(); - this->z_shift.y() = 0; - this->z_shift.x() = 0; + // this->z_shift.z() = this->get_scanner_ptr()->get_coordinate_for_det_pos(DetectionPosition<>(0, 0, 0)).z(); + // this->z_shift.y() = 0; + // this->z_shift.x() = 0; } ProjDataInfoGenericNoArcCorr* @@ -122,8 +122,8 @@ ProjDataInfoGenericNoArcCorr::get_LOR(LORInAxialAndNoArcCorrSinogramCoordinates< CartesianCoordinate3D _p2; find_cartesian_coordinates_of_detection(_p1, _p2, bin); - _p1.z() += z_shift.z(); - _p2.z() += z_shift.z(); + // _p1.z() += z_shift.z(); + // _p2.z() += z_shift.z(); LORAs2Points lor_as_2_points(_p1, _p2); const double R = sqrt(std::max(square(_p1.x()) + square(_p1.y()), square(_p2.x()) + square(_p2.y()))); @@ -176,8 +176,8 @@ ProjDataInfoGenericNoArcCorr::find_cartesian_coordinates_given_scanner_coordinat coord_1 = get_scanner_ptr()->get_coordinate_for_det_pos(det_pos1); coord_2 = get_scanner_ptr()->get_coordinate_for_det_pos(det_pos2); - coord_1.z() -= z_shift.z(); - coord_2.z() -= z_shift.z(); + // coord_1.z() -= z_shift.z(); + // coord_2.z() -= z_shift.z(); if (timing_pos_num < 0) { diff --git a/src/buildblock/VoxelsOnCartesianGrid.cxx b/src/buildblock/VoxelsOnCartesianGrid.cxx index 0ce83e326b..0bb075f785 100644 --- a/src/buildblock/VoxelsOnCartesianGrid.cxx +++ b/src/buildblock/VoxelsOnCartesianGrid.cxx @@ -371,32 +371,6 @@ VoxelsOnCartesianGrid::grow_z_range(const int min_z, const int max_z) this->grow(IndexRange<3>(min_indices, max_indices)); } -template -BasicCoordinate<3, int> -VoxelsOnCartesianGrid::get_lengths() const -{ - return make_coordinate(this->get_z_size(), this->get_y_size(), this->get_x_size()); -} - -template -BasicCoordinate<3, int> -VoxelsOnCartesianGrid::get_min_indices() const -{ - CartesianCoordinate3D min_indices; - CartesianCoordinate3D max_indices; - this->get_regular_range(min_indices, max_indices); - return min_indices; -} - -template -BasicCoordinate<3, int> -VoxelsOnCartesianGrid::get_max_indices() const -{ - CartesianCoordinate3D min_indices; - CartesianCoordinate3D max_indices; - this->get_regular_range(min_indices, max_indices); - return max_indices; -} #if 0 /**************************** diff --git a/src/experimental/motion/RigidObject3DTransformation.cxx b/src/experimental/motion/RigidObject3DTransformation.cxx index 35c8f6b6a2..06c7964c52 100644 --- a/src/experimental/motion/RigidObject3DTransformation.cxx +++ b/src/experimental/motion/RigidObject3DTransformation.cxx @@ -290,6 +290,7 @@ RigidObject3DTransformation::transform_bin(Bin& bin, in_proj_data_info.get_LOR(lor, bin); LORAs2Points lor_as_points; lor.get_intersections_with_cylinder(lor_as_points, lor.radius()); +# if 0 // AG: No longer needed? so commented out // TODO origin // currently, the origin used for proj_data_info is in the centre of the scanner, // while for standard images it is in the centre of the first ring. @@ -300,6 +301,9 @@ RigidObject3DTransformation::transform_bin(Bin& bin, // also uses an origin in the centre of the first ring const float z_shift = (in_proj_data_info.get_scanner_ptr()->get_num_rings() - 1) / 2.F * in_proj_data_info.get_scanner_ptr()->get_ring_spacing(); +# else + constexpr float z_shift = 0.F; +# endif lor_as_points.p1().z() += z_shift; lor_as_points.p2().z() += z_shift; LORAs2Points transformed_lor_as_points(transform_point(lor_as_points.p1()), transform_point(lor_as_points.p2())); diff --git a/src/experimental/recon_buildblock/DataSymmetriesForDensels_PET_CartesianGrid.cxx b/src/experimental/recon_buildblock/DataSymmetriesForDensels_PET_CartesianGrid.cxx index 1287908c5c..1ecc099740 100644 --- a/src/experimental/recon_buildblock/DataSymmetriesForDensels_PET_CartesianGrid.cxx +++ b/src/experimental/recon_buildblock/DataSymmetriesForDensels_PET_CartesianGrid.cxx @@ -85,15 +85,11 @@ find_relation_between_coordinate_systems(int& num_planes_per_scanner_ring, const float delta = proj_data_info_cyl_ptr->get_average_ring_difference(segment_num); - // KT 20/06/2001 take origin.z() into account - axial_pos_to_z_offset[segment_num] - = (cartesian_grid_info_ptr->get_max_index() + cartesian_grid_info_ptr->get_min_index()) / 2.F - - cartesian_grid_info_ptr->get_origin().z() / image_plane_spacing - - (num_planes_per_axial_pos[segment_num] - * (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; + // AG: Each segment is offset because of the angle, we adjust using delta + const float segment_offset_in_z_gantry_coords = proj_data_info_cyl_ptr->get_m(first_bin) - ring_spacing * delta / 2; + axial_pos_to_z_offset[segment_num] = cartesian_grid_info_ptr->get_index_coordinates_for_physical_coordinates( + proj_data_info_cyl_ptr->get_physical_coordinates_for_gantry_coordinates( + CartesianCoordinate3D(segment_offset_in_z_gantry_coords, 0, 0)))[1]; // z } } diff --git a/src/experimental/test/test_proj_data_info_LOR.cxx b/src/experimental/test/test_proj_data_info_LOR.cxx index f877e30e4c..a7e82c41c6 100644 --- a/src/experimental/test/test_proj_data_info_LOR.cxx +++ b/src/experimental/test/test_proj_data_info_LOR.cxx @@ -159,8 +159,7 @@ main(int argc, char* argv[]) ring1_90, ring2_90, det1_90, - det2_90, - 0); // set timing_pos_num=0 as test-code is pre-TOF + det2_90); // set timing_pos_num=0 as test-code is pre-TOF #if 0 cout << coord_1_0< @@ -242,6 +243,25 @@ class DiscretisedDensity : public ExamData, public ArrayType + // get_image_central_index() const; + + //! Get the "most-central" voxel + virtual BasicCoordinate get_image_centre_in_index_coordinates() const = 0; + + //! Get the average location of image sample points in physical coordinates + virtual CartesianCoordinate3D get_image_centre_in_physical_coordinates() const = 0; + + //! Get the average location of image sample points in LPS patient coordinates + virtual CartesianCoordinate3D get_image_centre_in_LPS_coordinates() const = 0; + + //@} + //! Allocate a new DiscretisedDensity object with same characteristics as the current one. virtual DiscretisedDensity* get_empty_copy() const = 0; diff --git a/src/include/stir/DiscretisedDensityOnCartesianGrid.h b/src/include/stir/DiscretisedDensityOnCartesianGrid.h index 7aa2e01a86..ef05968737 100644 --- a/src/include/stir/DiscretisedDensityOnCartesianGrid.h +++ b/src/include/stir/DiscretisedDensityOnCartesianGrid.h @@ -33,6 +33,10 @@ START_NAMESPACE_STIR \ingroup densitydata \brief This abstract class is the basis for images on a Cartesian grid. + STIR grid sampling MUST axes align with the gantry sampling axes (i.e., z runs + along the axis, y in the direction of gravity, and x perpendicular to both). + Only the LPS coordinate system is allowed to be rotated. + The only new information compared to DiscretisedDensity is the grid_spacing. This is currently a BasicCoordinate object. It really should be CartesianCoordinate object, but @@ -67,6 +71,25 @@ class DiscretisedDensityOnCartesianGrid : public DiscretisedDensity& grid_spacing_v); + // TODO: this is a confusing name + //! Get the size of the cartesian grid + inline BasicCoordinate get_lengths() const; + + //! Get the smallest indices on the sampling grid for each dim + inline BasicCoordinate get_min_indices() const; + + //! Get the largest indices on the sampling grid for each dim + inline BasicCoordinate get_max_indices() const; + + //! Get the "most-central" voxel + BasicCoordinate get_image_centre_in_index_coordinates() const; + + //! Get the average location of image sample points in physical coordinates + CartesianCoordinate3D get_image_centre_in_physical_coordinates() const; + + //! Get the average location of image sample points in LPS patient coordinates + CartesianCoordinate3D get_image_centre_in_LPS_coordinates() const; + protected: inline bool actual_has_same_characteristics(DiscretisedDensity const&, std::string& explanation) const override; diff --git a/src/include/stir/DiscretisedDensityOnCartesianGrid.inl b/src/include/stir/DiscretisedDensityOnCartesianGrid.inl index 5c07550dd2..1ce54ab18c 100644 --- a/src/include/stir/DiscretisedDensityOnCartesianGrid.inl +++ b/src/include/stir/DiscretisedDensityOnCartesianGrid.inl @@ -127,4 +127,71 @@ DiscretisedDensityOnCartesianGrid::actual_get_index_coord return float_indices / this->get_grid_spacing(); } +template +BasicCoordinate +DiscretisedDensityOnCartesianGrid::get_min_indices() const +{ + BasicCoordinate min_indices; + BasicCoordinate max_indices; + this->get_regular_range(min_indices, max_indices); + return min_indices; +} + +template +BasicCoordinate +DiscretisedDensityOnCartesianGrid::get_max_indices() const +{ + BasicCoordinate min_indices; + BasicCoordinate max_indices; + if (not this->get_regular_range(min_indices, max_indices)) + { + // this should never happen? + error("Discretised Density image wasn't regular!"); + } + return max_indices; +} + +template +BasicCoordinate +DiscretisedDensityOnCartesianGrid::get_lengths() const +{ + BasicCoordinate min_indices; + BasicCoordinate max_indices; + if (not this->get_regular_range(min_indices, max_indices)) + { + // this should never happen? + error("Discretised Density image wasn't regular!"); + } + return max_indices - min_indices + 1; +} + +template +BasicCoordinate +DiscretisedDensityOnCartesianGrid::get_image_centre_in_index_coordinates() const +{ + BasicCoordinate min_indices; + BasicCoordinate max_indices; + if (not this->get_regular_range(min_indices, max_indices)) + { + // this should never happen? + error("Discretised Density image wasn't regular!"); + } + return BasicCoordinate(min_indices) + + BasicCoordinate(max_indices - min_indices) / 2; +} + +template +CartesianCoordinate3D +DiscretisedDensityOnCartesianGrid::get_image_centre_in_physical_coordinates() const +{ + return this->get_physical_coordinates_for_indices(this->get_image_centre_in_index_coordinates()); +} + +template +CartesianCoordinate3D +DiscretisedDensityOnCartesianGrid::get_image_centre_in_LPS_coordinates() const +{ + return this->get_LPS_coordinates_for_indices(this->get_image_centre_in_index_coordinates()); +} + END_NAMESPACE_STIR diff --git a/src/include/stir/ProjDataInfo.h b/src/include/stir/ProjDataInfo.h index f6ef1a032a..dec519bcb6 100644 --- a/src/include/stir/ProjDataInfo.h +++ b/src/include/stir/ProjDataInfo.h @@ -32,6 +32,7 @@ #include "stir/SinogramIndices.h" #include "stir/VectorWithOffset.h" #include "stir/Scanner.h" +#include "stir/CartesianCoordinate3D.h" #include "stir/shared_ptr.h" #include "stir/unique_ptr.h" #include @@ -332,6 +333,18 @@ class ProjDataInfo virtual void get_LOR(LORInAxialAndNoArcCorrSinogramCoordinates&, 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 @@ -496,9 +509,23 @@ class ProjDataInfo return scanner_ptr->has_energy_information(); } + //! Vector represention bed position in 3D + CartesianCoordinate3D get_bed_position() const; + + // Convert coordinates from gantry to physical + inline CartesianCoordinate3D + get_physical_coordinates_for_gantry_coordinates(const CartesianCoordinate3D& coords) const; + + // Convert coordinates from physical to gantry + inline CartesianCoordinate3D + get_gantry_coordinates_for_physical_coordinates(const CartesianCoordinate3D& coords) const; + protected: virtual bool blindly_equals(const root_type* const) const = 0; + //! 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; + private: shared_ptr scanner_ptr; int min_view_num; diff --git a/src/include/stir/ProjDataInfo.inl b/src/include/stir/ProjDataInfo.inl index 748a910a50..eba3419bc6 100644 --- a/src/include/stir/ProjDataInfo.inl +++ b/src/include/stir/ProjDataInfo.inl @@ -220,6 +220,20 @@ ProjDataInfo::get_scanner_sptr() const return scanner_ptr; } +CartesianCoordinate3D +ProjDataInfo::get_physical_coordinates_for_gantry_coordinates(const CartesianCoordinate3D& coords) const +{ + // TODO: bed postion + return coords + get_vector_centre_of_first_ring_to_centre_of_gantry(); +} + +CartesianCoordinate3D +ProjDataInfo::get_gantry_coordinates_for_physical_coordinates(const CartesianCoordinate3D& coords) const +{ + // TODO: bed postion + return coords - get_vector_centre_of_first_ring_to_centre_of_gantry(); +} + int ProjDataInfo::get_num_non_tof_sinograms() const { diff --git a/src/include/stir/ProjDataInfoBlocksOnCylindricalNoArcCorr.h b/src/include/stir/ProjDataInfoBlocksOnCylindricalNoArcCorr.h index f469d245fe..85ffe2ae0f 100644 --- a/src/include/stir/ProjDataInfoBlocksOnCylindricalNoArcCorr.h +++ b/src/include/stir/ProjDataInfoBlocksOnCylindricalNoArcCorr.h @@ -68,15 +68,13 @@ class ProjDataInfoBlocksOnCylindricalNoArcCorr : public ProjDataInfoGenericNoArc std::string parameter_info() const override; private: - //! \name set of obsolete functions to go between bins<->LORs (will disappear!) - //@{ Succeeded find_scanner_coordinates_given_cartesian_coordinates(int& det1, int& det2, int& ring1, int& ring2, const CartesianCoordinate3D& c1, const CartesianCoordinate3D& c2) const override; - //@} + // give test classes access to the private members friend class ProjDataInfoTests; friend class BlocksTests; diff --git a/src/include/stir/ProjDataInfoCylindrical.h b/src/include/stir/ProjDataInfoCylindrical.h index 0e46595a20..09d7d270f2 100644 --- a/src/include/stir/ProjDataInfoCylindrical.h +++ b/src/include/stir/ProjDataInfoCylindrical.h @@ -78,12 +78,9 @@ class ProjDataInfoCylindrical : public ProjDataInfo inline float get_t(const Bin&) const override; - //! Return z-coordinate of the middle of the LOR + //! Return z-coordinate of the middle of the LOR in gantry coordinates /*! - The 0 of the z-axis is chosen in the middle of the scanner. - - \warning Current implementation assumes that the axial positions are always 'centred', - i.e. get_m(Bin(..., min_axial_pos_num,...)) == - get_m(Bin(..., max_axial_pos_num,...)) + In gantry coordinates, the 0 of the z-axis is chosen in the middle of the scanner. */ inline float get_m(const Bin&) const override; diff --git a/src/include/stir/ProjDataInfoGenericNoArcCorr.h b/src/include/stir/ProjDataInfoGenericNoArcCorr.h index 043121c1c7..6253b3e170 100644 --- a/src/include/stir/ProjDataInfoGenericNoArcCorr.h +++ b/src/include/stir/ProjDataInfoGenericNoArcCorr.h @@ -155,13 +155,13 @@ class ProjDataInfoGenericNoArcCorr : public ProjDataInfoCylindricalNoArcCorr */ // This version uses the coordinate map - virtual void find_cartesian_coordinates_given_scanner_coordinates(CartesianCoordinate3D& coord_1, - CartesianCoordinate3D& coord_2, - const int Ring_A, - const int Ring_B, - const int det1, - const int det2, - const int timing_pos_num = 0) const override; + void find_cartesian_coordinates_given_scanner_coordinates(CartesianCoordinate3D& coord_1, + CartesianCoordinate3D& coord_2, + const int Ring_A, + const int Ring_B, + const int det1, + const int det2, + const int timing_pos_num = 0) const override; //@} protected: diff --git a/src/include/stir/RunTests.h b/src/include/stir/RunTests.h index b406962608..51b55557cf 100644 --- a/src/include/stir/RunTests.h +++ b/src/include/stir/RunTests.h @@ -139,24 +139,42 @@ class RunTests } //! check equality by comparing ranges and calling check_if_equal on all elements template - bool check_if_equal(const VectorWithOffset& t1, const VectorWithOffset& t2, const std::string& str = "") + bool check_if_equal(const typename VectorWithOffset::const_iterator& t1_begin, + const typename VectorWithOffset::const_iterator& t1_end, + const typename VectorWithOffset::const_iterator& t2_begin, + const typename VectorWithOffset::const_iterator& t2_end, + const std::string& str = "") { - if (t1.get_min_index() != t2.get_min_index() || t1.get_max_index() != t2.get_max_index()) + assert(t1_end >= t1_begin); + size_t size = t1_end - t1_begin; + if (t2_begin + size != t2_end) { std::cerr << "Error: unequal ranges. " << str << std::endl; - return everything_ok = false; } - - for (int i = t1.get_min_index(); i <= t1.get_max_index(); i++) + for (size_t i = 0; i < size; i++) { - if (!check_if_equal(t1[i], t2[i], str)) + if (!check_if_equal(*(t1_begin + i), *(t2_begin + i), str)) { std::cerr << "(at VectorWithOffset<" << typeid(T).name() << "> first mismatch at index " << i << ")\n"; - return everything_ok = false; + return false; } } return true; } + + //! check equality by comparing ranges and calling check_if_equal on all elements + template + bool check_if_equal(const VectorWithOffset& t1, const VectorWithOffset& t2, const std::string& str = "") + { + if (t1.get_min_index() != t2.get_min_index() || t1.get_max_index() != t2.get_max_index()) + { + std::cerr << "Error: unequal ranges: (" << t1.get_min_index() << "," << t1.get_max_index() << ") (" << t2.get_min_index() + << "," << t2.get_max_index() << "). " << str << std::endl; + return everything_ok = false; + } + + return check_if_equal(t1.begin(), t1.end(), t2.begin(), t2.end(), str); + } // VC 6.0 needs definition of template members in the class def unfortunately. //! check equality by comparing size and calling check_if_equal on all elements template diff --git a/src/include/stir/Succeeded.h b/src/include/stir/Succeeded.h index 99feca1baa..2c4e617bc7 100644 --- a/src/include/stir/Succeeded.h +++ b/src/include/stir/Succeeded.h @@ -55,6 +55,8 @@ class Succeeded bool operator!=(const Succeeded& v2) const { return v != v2.v; } //! convenience function returns if it is equal to Succeeded::yes bool succeeded() const { return this->v == yes; } + Succeeded operator&(const Succeeded& v2) const { return (v == yes && v2.v == yes) ? yes : no; } + Succeeded operator&=(const Succeeded& v2) const { return *this & v2; } private: value v; diff --git a/src/include/stir/VoxelsOnCartesianGrid.h b/src/include/stir/VoxelsOnCartesianGrid.h index 8016e54d77..9ea9a0a4b1 100644 --- a/src/include/stir/VoxelsOnCartesianGrid.h +++ b/src/include/stir/VoxelsOnCartesianGrid.h @@ -195,10 +195,6 @@ static VoxelsOnCartesianGrid ask_parameters(); inline int get_max_z() const; - BasicCoordinate<3, int> get_lengths() const; - BasicCoordinate<3, int> get_min_indices() const; - BasicCoordinate<3, int> get_max_indices() const; - //@} //! \name Numerical operations diff --git a/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.h b/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.h index a1c9f33d0f..5e8f029630 100644 --- a/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.h +++ b/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.h @@ -126,9 +126,8 @@ class DataSymmetriesForBins_PET_CartesianGrid : public DataSymmetriesForBins //! find correspondence between axial_pos_num and image coordinates /*! z = num_planes_per_axial_pos * axial_pos_num + axial_pos_to_z_offset - - compute the offset by matching up the centre of the scanner - in the 2 coordinate systems + compute the offset with ProjDataInfoCylindrical's gantry to + physical coordinate systems mapping. */ inline float get_num_planes_per_axial_pos(const int segment_num) const; inline float get_axial_pos_to_z_offset(const int segment_num) const; diff --git a/src/include/stir/recon_buildblock/ProjMatrixByBin.h b/src/include/stir/recon_buildblock/ProjMatrixByBin.h index a980b62197..b61295e132 100644 --- a/src/include/stir/recon_buildblock/ProjMatrixByBin.h +++ b/src/include/stir/recon_buildblock/ProjMatrixByBin.h @@ -242,4 +242,4 @@ END_NAMESPACE_STIR #include "stir/recon_buildblock/ProjMatrixByBin.inl" -#endif // __ProjMatrixByBin_H__ +#endif // __ProjMatrixByBin_H__ \ No newline at end of file diff --git a/src/include/stir/recon_buildblock/ProjMatrixByBinUsingInterpolation.h b/src/include/stir/recon_buildblock/ProjMatrixByBinUsingInterpolation.h index 6d0a3a208c..12a48902ae 100644 --- a/src/include/stir/recon_buildblock/ProjMatrixByBinUsingInterpolation.h +++ b/src/include/stir/recon_buildblock/ProjMatrixByBinUsingInterpolation.h @@ -73,13 +73,14 @@ class ProjMatrixByBinUsingInterpolation bool do_symmetry_swap_s; bool do_symmetry_shift_z; - // explicitly list necessary members for image details (should use an Info object instead) + shared_ptr proj_data_info_ptr; + shared_ptr> density_info_ptr; + + // cache a few density_info_ptr-derived params CartesianCoordinate3D voxel_size; CartesianCoordinate3D origin; IndexRange<3> densel_range; - shared_ptr proj_data_info_ptr; - // for Jacobian const ProjDataInfoCylindrical& proj_data_info_cyl() const { diff --git a/src/include/stir/recon_buildblock/ProjMatrixByBinUsingRayTracing.h b/src/include/stir/recon_buildblock/ProjMatrixByBinUsingRayTracing.h index 3d42397117..12d3a1047a 100644 --- a/src/include/stir/recon_buildblock/ProjMatrixByBinUsingRayTracing.h +++ b/src/include/stir/recon_buildblock/ProjMatrixByBinUsingRayTracing.h @@ -181,11 +181,13 @@ class ProjMatrixByBinUsingRayTracing bool do_symmetry_swap_s; bool do_symmetry_shift_z; - // explicitly list necessary members for image details (should use an Info object instead) + shared_ptr> density_info_sptr; + + // chase a few needed and constant after set_up() image details CartesianCoordinate3D voxel_size; CartesianCoordinate3D origin; - CartesianCoordinate3D min_index; - CartesianCoordinate3D max_index; + CartesianCoordinate3D min_index, max_index; + float fovrad_in_mm; void calculate_proj_matrix_elems_for_one_bin(ProjMatrixElemsForOneBin&) const override; diff --git a/src/include/stir/recon_buildblock/RayTraceVoxelsOnCartesianGrid.h b/src/include/stir/recon_buildblock/RayTraceVoxelsOnCartesianGrid.h index 9aac91fc78..aea00bc078 100644 --- a/src/include/stir/recon_buildblock/RayTraceVoxelsOnCartesianGrid.h +++ b/src/include/stir/recon_buildblock/RayTraceVoxelsOnCartesianGrid.h @@ -34,8 +34,8 @@ class CartesianCoordinate3D; the ProjMatrixElemsForOneBin object. \param lor object to which the intersected voxels and the LOI will be appended - \param start_point first point on the LOR. The first voxel will contain this point. - \param end_point last point on the LOR. The last voxel will contain this point. + \param start_point first point on the LOR in index coordinates. The first voxel will contain this point. + \param end_point last point on the LOR in index coordinates. The last voxel will contain this point. \param voxel_size normally in mm \param normalisation_constant LOIs will be multiplied with this constant diff --git a/src/include/stir/recon_buildblock/test/ReconstructionTests.h b/src/include/stir/recon_buildblock/test/ReconstructionTests.h index 6c8d3bd157..7fd8fcb167 100644 --- a/src/include/stir/recon_buildblock/test/ReconstructionTests.h +++ b/src/include/stir/recon_buildblock/test/ReconstructionTests.h @@ -138,9 +138,15 @@ ReconstructionTests::construct_input_data() this->_input_density_sptr = aptr; } - // forward project + // forward project, don't use symmetries here - a bit slower but we want to + // test their consistency in the reconstructions { - shared_ptr PM_sptr(new ProjMatrixByBinUsingRayTracing); + shared_ptr PM_sptr(new ProjMatrixByBinUsingRayTracing); + PM_sptr->set_do_symmetry_180degrees_min_phi(false); + PM_sptr->set_do_symmetry_90degrees_min_phi(false); + PM_sptr->set_do_symmetry_shift_z(false); + PM_sptr->set_do_symmetry_swap_s(false); + PM_sptr->set_do_symmetry_swap_segment(false); shared_ptr fwd_proj_sptr = MAKE_SHARED(PM_sptr); fwd_proj_sptr->set_up(this->_proj_data_sptr->get_proj_data_info_sptr(), this->_input_density_sptr); fwd_proj_sptr->set_input(*this->_input_density_sptr); diff --git a/src/include/stir/scatter/ScatterSimulation.h b/src/include/stir/scatter/ScatterSimulation.h index bcf6e2b919..7750fb41ab 100644 --- a/src/include/stir/scatter/ScatterSimulation.h +++ b/src/include/stir/scatter/ScatterSimulation.h @@ -262,9 +262,10 @@ class ScatterSimulation : public RegisteredObject act_image_type, att_image_type }; + //! Scatter emission point in physical coordinates struct ScatterPoint { - CartesianCoordinate3D coord; + CartesianCoordinate3D physical_coord; float mu_value; }; @@ -275,6 +276,7 @@ class ScatterSimulation : public RegisteredObject //! find scatter points /*! This function sets scatt_points_vector and scatter_volume. It will also remove any cached integrals as they would be incorrect otherwise. + Scatter points are defined in physical coordinates. */ void sample_scatter_points(); @@ -301,7 +303,7 @@ class ScatterSimulation : public RegisteredObject virtual void find_detectors(unsigned& det_num_A, unsigned& det_num_B, const Bin& bin) const; - unsigned find_in_detection_points_vector(const CartesianCoordinate3D& coord) const; + unsigned find_in_detection_points_vector(const CartesianCoordinate3D& gantry_coord) const; CartesianCoordinate3D shift_detector_coordinates_to_origin; @@ -309,7 +311,7 @@ class ScatterSimulation : public RegisteredObject double detection_efficiency_no_scatter(const unsigned det_num_A, const unsigned det_num_B) const; // next needs to be mutable because find_in_detection_points_vector is const - mutable std::vector> detection_points_vector; + mutable std::vector> detection_points_in_gantry_coords_vector; //!@} @@ -319,14 +321,15 @@ class ScatterSimulation : public RegisteredObject //! \name integrating functions //@{ static float integral_between_2_points(const DiscretisedDensity<3, float>& density, - const CartesianCoordinate3D& point1, - const CartesianCoordinate3D& point2); + const CartesianCoordinate3D& point1_physical_coord, + const CartesianCoordinate3D& point2_physical_coord); - float exp_integral_over_attenuation_image_between_scattpoint_det(const CartesianCoordinate3D& scatter_point, - const CartesianCoordinate3D& detector_coord); + float + exp_integral_over_attenuation_image_between_scattpoint_det(const CartesianCoordinate3D& scatter_point_physical_coord, + const CartesianCoordinate3D& detector_physical_coord); - float integral_over_activity_image_between_scattpoint_det(const CartesianCoordinate3D& scatter_point, - const CartesianCoordinate3D& detector_coord); + float integral_over_activity_image_between_scattpoint_det(const CartesianCoordinate3D& scatter_point_physical_coord, + const CartesianCoordinate3D& detector_physical_coord); float cached_integral_over_activity_image_between_scattpoint_det(const unsigned scatter_point_num, const unsigned det_num); diff --git a/src/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.cxx b/src/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.cxx index 82391fef93..626cd56618 100644 --- a/src/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.cxx +++ b/src/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.cxx @@ -87,6 +87,7 @@ find_relation_between_coordinate_systems(int& num_planes_per_scanner_ring, error("DataSymmetriesForBins_PET_CartesianGrid can currently only support y-origin = 0 " "Sorry\n"); */ + const float ring_spacing = proj_data_info_cyl_ptr->get_ring_spacing(); for (int segment_num = min_segment_num; segment_num <= max_segment_num; ++segment_num) { @@ -105,17 +106,14 @@ find_relation_between_coordinate_systems(int& num_planes_per_scanner_ring, segment_num)); } + // get a bin (any bin) at axial position 0 + Bin first_bin(segment_num, 0, 0, 0); const float delta = proj_data_info_cyl_ptr->get_average_ring_difference(segment_num); + const float segment_offset_in_z_gantry_coords = proj_data_info_cyl_ptr->get_m(first_bin) - ring_spacing * delta / 2; - // KT 20/06/2001 take origin.z() into account - axial_pos_to_z_offset[segment_num] - = (cartesian_grid_info_ptr->get_max_index() + cartesian_grid_info_ptr->get_min_index()) / 2.F - - cartesian_grid_info_ptr->get_origin().z() / image_plane_spacing - - (num_planes_per_axial_pos[segment_num] - * (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; + axial_pos_to_z_offset[segment_num] = cartesian_grid_info_ptr->get_index_coordinates_for_physical_coordinates( + proj_data_info_cyl_ptr->get_physical_coordinates_for_gantry_coordinates( + CartesianCoordinate3D(segment_offset_in_z_gantry_coords, 0, 0)))[1]; } } diff --git a/src/recon_buildblock/ProjMatrixByBinUsingInterpolation.cxx b/src/recon_buildblock/ProjMatrixByBinUsingInterpolation.cxx index ed82fb4337..01ed7c0409 100644 --- a/src/recon_buildblock/ProjMatrixByBinUsingInterpolation.cxx +++ b/src/recon_buildblock/ProjMatrixByBinUsingInterpolation.cxx @@ -26,7 +26,9 @@ #include "stir/ProjDataInfoCylindricalArcCorr.h" #include "stir/ProjDataInfoCylindricalNoArcCorr.h" #include "stir/Bin.h" +#include "stir/LORCoordinates.h" #include "stir/round.h" +#include "stir/format.h" #include "stir/warning.h" #include "stir/error.h" #include "stir/recon_buildblock/ProjMatrixElemsForOneBin.h" @@ -94,23 +96,27 @@ ProjMatrixByBinUsingInterpolation::post_processing() void ProjMatrixByBinUsingInterpolation::set_up( const shared_ptr& proj_data_info_ptr_v, - const shared_ptr>& density_info_ptr // TODO should be Info only + const shared_ptr>& density_info_ptr_v // TODO should be Info only ) { - ProjMatrixByBin::set_up(proj_data_info_ptr_v, density_info_ptr); + ProjMatrixByBin::set_up(proj_data_info_ptr_v, density_info_ptr_v); proj_data_info_ptr = proj_data_info_ptr_v; + density_info_ptr = density_info_ptr_v; const VoxelsOnCartesianGrid* image_info_ptr = dynamic_cast*>(density_info_ptr.get()); if (image_info_ptr == NULL) error("ProjMatrixByBinUsingInterpolation initialised with a wrong type of DiscretisedDensity\n"); + CartesianCoordinate3D origin = image_info_ptr->get_origin(); + if (origin.x() != 0 or origin.y() != 0) + { + error(format("ProjMatrixByBinUsingInterpolation expects a transaxially-centred image ({},{})\n", origin.x(), origin.y())); + } + densel_range = image_info_ptr->get_index_range(); voxel_size = image_info_ptr->get_voxel_size(); - 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; - origin.z() -= z_to_middle; symmetries_sptr.reset(new DataSymmetriesForBins_PET_CartesianGrid(proj_data_info_ptr, density_info_ptr, @@ -149,6 +155,7 @@ ProjMatrixByBinUsingInterpolation::clone() const } // point should be w.r.t. middle of the scanner! +// AG: This means in gantry coordinates static inline void find_s_m_of_voxel( float& s, float& m, const CartesianCoordinate3D& point, const float cphi, const float sphi, const float tantheta) @@ -284,43 +291,42 @@ ProjMatrixByBinUsingInterpolation::calculate_proj_matrix_elems_for_one_bin(ProjM int max1; // find z-range (this would depend on origin and the symmetries though) { -#if 0 - /* attempt to take a large range of essentially 3 times the z-range of the image. - This does not work though. It's easy to come up with cases where this - range is still too small. Probably we would have to use the - scanner length or so and take the image-origin into account. - It's not easy though as illustrated by the following example: - - 2 rings, span=1, z-voxel_size=ring_distance/2, 3 image-planes, segment 0 - and image-origin shifted over 3 planes. - The problem comes when using the swap_*_zq symmetries (which occurs - even if shift_z=false). - */ - min1=densel_range.get_min_index(); - max1=densel_range.get_max_index(); - int length = max-min1+1; - min1 -= length; - max1 += length; -#else - /* Here we use geometric info. However, the code below only works + /* Here we use DiscretisedDensity (Info). However, the code below only works for DiscretisedDensityOnCartesianGrid (with a regular_range) + NB: this uses the edges, not the corners, so the radius is tangent to the + furthermost image edge. I'm matching the old implementation but probably + corners are preferred? (AG) */ - min1 = densel_range.get_min_index(); - // find radius of cylinder around all of the image + BasicCoordinate<3, int> min_index, max_index; + density_info_ptr->get_regular_range(min_index, max_index); + CartesianCoordinate3D min_gantry_coords = proj_data_info_ptr->get_gantry_coordinates_for_physical_coordinates( + density_info_ptr->get_physical_coordinates_for_indices(min_index)); + CartesianCoordinate3D max_gantry_coords = proj_data_info_ptr->get_gantry_coordinates_for_physical_coordinates( + density_info_ptr->get_physical_coordinates_for_indices(max_index)); const float max_radius - = std::max(std::max(-densel_range[min1].get_min_index(), densel_range[min1].get_max_index()) * voxel_size[2], - std::max(-densel_range[min1][0].get_min_index(), densel_range[min1][0].get_max_index()) * voxel_size[1]); - // find width of the 'tube of response' + = std::max({ -min_gantry_coords.x(), -min_gantry_coords.y(), max_gantry_coords.x(), max_gantry_coords.y() }); + const float z_width_of_TOR = proj_data_info_ptr->get_sampling_in_m(bin); - // now find where tube enters/leaves image - // we use a safety margin here as we could probably use z_width_of_LOR/2 - const float z_middle_LOR = proj_data_info_ptr->get_m(bin) - origin.z(); - const float min_z_LOR = z_middle_LOR - fabs(proj_data_info_ptr->get_tantheta(bin)) * max_radius - z_width_of_TOR; - const float max_z_LOR = z_middle_LOR + fabs(proj_data_info_ptr->get_tantheta(bin)) * max_radius + z_width_of_TOR; - - min1 = round(floor(min_z_LOR / voxel_size[1])); - max1 = round(ceil(max_z_LOR / voxel_size[1])); -#endif + + // Get the LOR for bin, but a radius to just cover the FOV + LORInAxialAndNoArcCorrSinogramCoordinates lor; + proj_data_info_ptr->get_LOR(lor, bin); + LORAs2Points reduced_fov_lor; + find_LOR_intersections_with_cylinder(reduced_fov_lor, LORAs2Points(lor), max_radius); + + // now find the z extents in gantry coordinates and convert into indices + float min_z_in_gantry_coords = std::min(reduced_fov_lor.p1().z(), reduced_fov_lor.p2().z()); + float max_z_in_gantry_coords = std::max(reduced_fov_lor.p1().z(), reduced_fov_lor.p2().z()); + // NB: This could just be z_width_of_TOR/2, but old implementation preferred + // not to divide to add a "safety margin", so replicating here (AG) + min_z_in_gantry_coords -= z_width_of_TOR; + max_z_in_gantry_coords += z_width_of_TOR; + min1 = floor(density_info_ptr->get_index_coordinates_for_physical_coordinates( + proj_data_info_ptr->get_physical_coordinates_for_gantry_coordinates( + CartesianCoordinate3D(min_z_in_gantry_coords, 0, 0)))[1]); + max1 = ceil(density_info_ptr->get_index_coordinates_for_physical_coordinates( + proj_data_info_ptr->get_physical_coordinates_for_gantry_coordinates( + CartesianCoordinate3D(max_z_in_gantry_coords, 0, 0)))[1]); } /* we loop over all coordinates, but for optimisation do the following: In each dimension, we ASSUME that the non-zero range is CONNECTED. @@ -345,8 +351,6 @@ ProjMatrixByBinUsingInterpolation::calculate_proj_matrix_elems_for_one_bin(ProjM // which are outside the FOV // this will break when non-zero origin.y() or x() // (but then there would be no relevant symmetries I guess) - assert(origin.y() == 0); - assert(origin.x() == 0); const int first_min2 = range2d.get_min_index(); const int first_max2 = range2d.get_max_index(); const int min2 = std::max(first_min2, -first_max2); @@ -370,9 +374,8 @@ ProjMatrixByBinUsingInterpolation::calculate_proj_matrix_elems_for_one_bin(ProjM found_nonzero3 = false; for (c[3] = min3; c[3] <= max3; ++c[3]) { - // TODO call a virtual function of DiscretisedDensity? - const CartesianCoordinate3D coords - = CartesianCoordinate3D(c[1] * voxel_size[1], c[2] * voxel_size[2], c[3] * voxel_size[3]) + origin; + const CartesianCoordinate3D coords = proj_data_info_ptr->get_gantry_coordinates_for_physical_coordinates( + density_info_ptr->get_physical_coordinates_for_indices(c)); const float element_value = get_element(bin, coords); if (element_value > 0) { diff --git a/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx b/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx index ccdd80cd47..4a4ba0b1cb 100644 --- a/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx +++ b/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx @@ -263,6 +263,7 @@ ProjMatrixByBinUsingRayTracing::set_up( } ProjMatrixByBin::set_up(proj_data_info_sptr_v, density_info_sptr_v); + this->density_info_sptr = density_info_sptr_v; if (proj_data_info_sptr->get_scanner_ptr()->get_scanner_geometry() == "BlocksOnCylindrical" && !use_actual_detector_boundaries) { @@ -273,9 +274,10 @@ ProjMatrixByBinUsingRayTracing::set_up( voxel_size = image_info_ptr->get_voxel_size(); origin = image_info_ptr->get_origin(); + image_info_ptr->get_regular_range(min_index, max_index); + if (std::abs(origin.x()) > .05F || std::abs(origin.y()) > .05F) error("ProjMatrixByBinUsingRayTracing sadly doesn't support shifted x/y origin yet"); - image_info_sptr->get_regular_range(min_index, max_index); symmetries_sptr.reset(new DataSymmetriesForBins_PET_CartesianGrid(proj_data_info_sptr, density_info_sptr_v, @@ -387,10 +389,25 @@ ProjMatrixByBinUsingRayTracing::set_up( info("ProjMatrixByBinUsingRayTracing: use_actual_detector_boundaries==true.", 3); } + // precalculate some values that are constant over all bins + CartesianCoordinate3D min_pos, max_pos; + image_info_ptr->get_regular_range(min_index, max_index); +#ifdef STIR_PMRT_LARGER_FOV + // use FOV which is slightly 'inside' the image to avoid + // index out of range + max_pos = density_info_sptr->get_relative_coordinates_for_indices(max_index + 0.45F); + min_pos = density_info_sptr->get_relative_coordinates_for_indices(min_index - 0.45F); +#else + max_pos = density_info_sptr->get_relative_coordinates_for_indices(max_index); + min_pos = density_info_sptr->get_relative_coordinates_for_indices(min_index); +#endif + fovrad_in_mm = min(min(max_pos.x(), -min_pos.x()), min(max_pos.y(), -min_pos.y())); + #if 0 // 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) @@ -434,27 +451,41 @@ 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_physical_coordinates( + proj_data_info.get_physical_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, const float s_in_mm, - const float t_in_mm, + const float m_in_mm, const float cphi, const float sphi, - const float costheta, 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 DiscretisedDensity<3, float>& density_info, + const ProjDataInfo& proj_data_info) { 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, @@ -506,18 +537,21 @@ ray_trace_one_lor(ProjMatrixElemsForOneBin& lor, max_a = min((fovrad_in_mm * sign(sphi) - s_in_mm * cphi) / sphi, (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 = get_point_on_lor_in_index_coordinates( + s_in_mm, m_in_mm + offset_in_z, max_a, cphi, sphi, tantheta, density_info, proj_data_info); + CartesianCoordinate3D stop_point = get_point_on_lor_in_index_coordinates( + s_in_mm, m_in_mm + offset_in_z, min_a, cphi, sphi, tantheta, density_info, proj_data_info); #if 0 // KT 18/05/2005 this is no longer necessary @@ -592,7 +626,7 @@ ProjMatrixByBinUsingRayTracing::calculate_proj_matrix_elems_for_one_bin(ProjMatr assert(lor.size() == 0); - float phi; + float phi = proj_data_info_sptr->get_phi(bin); float s_in_mm = proj_data_info_sptr->get_s(bin); /* Implementation note. KT initialised s_in_mm above instead of in the if because this meant @@ -659,7 +693,7 @@ ProjMatrixByBinUsingRayTracing::calculate_proj_matrix_elems_for_one_bin(ProjMatr 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_sptr->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_sptr->get_sampling_in_t(bin) / costheta; @@ -681,11 +715,9 @@ ProjMatrixByBinUsingRayTracing::calculate_proj_matrix_elems_for_one_bin(ProjMatr voxel_size.z())); // 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. @@ -708,49 +740,35 @@ ProjMatrixByBinUsingRayTracing::calculate_proj_matrix_elems_for_one_bin(ProjMatr we might be using 2 rays for one ring. */ 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(); + = (-sampling_distance_of_adjacent_LORs_z / (2 * num_lors_per_axial_pos) * (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(); } - // 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()); -#else - 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()); - if (proj_data_info_sptr->get_scanner_ptr()->get_scanner_geometry() == "BlocksOnCylindrical") - { - fovrad_in_mm = min((min(max_index.x(), -min_index.x()) - 5.f) * voxel_size.x(), - (min(max_index.y(), -min_index.y()) - 5.f) * voxel_size.y()); - } -#endif - if (num_tangential_LORs == 1) { ray_trace_one_lor(lor, s_in_mm, - t_in_mm, + m_in_mm, cphi, sphi, - costheta, tantheta, offset_in_z, fovrad_in_mm, voxel_size, restrict_to_cylindrical_FOV, - num_lors_per_axial_pos); + num_lors_per_axial_pos, + *density_info_sptr, + *proj_data_info_sptr); } else { @@ -766,16 +784,17 @@ ProjMatrixByBinUsingRayTracing::calculate_proj_matrix_elems_for_one_bin(ProjMatr ray_traced_lor.erase(); ray_trace_one_lor(ray_traced_lor, current_s_in_mm, - t_in_mm, + m_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); + 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); } @@ -786,12 +805,38 @@ ProjMatrixByBinUsingRayTracing::calculate_proj_matrix_elems_for_one_bin(ProjMatr { 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); + // we want to use add_adjacent_z to fill in mutliple lines + // within the TOR, leveraging our traced LOR + // since tantheta==0, z is constant + const float z_of_traced_lor = static_cast(lor.begin()->coord1()); + + // We want the z limits for the TOR. We set s and a to 0. + // Since tantheta==0, they will not affect the z coordinate. + // (actually, as stated above, the LOR has a constant z, and + // so does the TOR) + float z_of_start_of_tor = get_point_on_lor_in_index_coordinates(0, + m_in_mm - sampling_distance_of_adjacent_LORs_z / 2, + 0, + cphi, + sphi, + tantheta, + *density_info_sptr, + *proj_data_info_sptr) + .z(); + float z_of_end_of_tor = get_point_on_lor_in_index_coordinates(0, + m_in_mm + sampling_distance_of_adjacent_LORs_z / 2, + 0, + cphi, + sphi, + tantheta, + *density_info_sptr, + *proj_data_info_sptr) + .z(); + + const float z_of_traced_lor_relative_to_start_of_tor = z_of_traced_lor - z_of_start_of_tor; + const float z_of_end_of_tor_relative_to_start_of_tor = z_of_end_of_tor - z_of_start_of_tor; + + add_adjacent_z(lor, z_of_traced_lor_relative_to_start_of_tor, z_of_end_of_tor_relative_to_start_of_tor); } else if (num_lors_per_axial_pos > 1) { diff --git a/src/recon_test/CMakeLists.txt b/src/recon_test/CMakeLists.txt index 93a0eb29e7..0283a0c7e5 100644 --- a/src/recon_test/CMakeLists.txt +++ b/src/recon_test/CMakeLists.txt @@ -26,6 +26,7 @@ set(${dir_SIMPLE_TEST_EXE_SOURCES} test_FBP3DRP.cxx test_blocks_on_cylindrical_projectors.cxx test_geometry_blocks_on_cylindrical.cxx + test_coordinate_system_mapping.cxx ) diff --git a/src/recon_test/test_coordinate_system_mapping.cxx b/src/recon_test/test_coordinate_system_mapping.cxx new file mode 100644 index 0000000000..a3be213536 --- /dev/null +++ b/src/recon_test/test_coordinate_system_mapping.cxx @@ -0,0 +1,593 @@ +/* + Copyright (C) 2020, CSIRO + This file is part of STIR. + + This file is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. + + This file is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + See STIR/LICENSE.txt for details +*/ +/*! + + \file + \ingroup test + + \brief Test program for various parts of code that map sinogram and + gantry space back to image space. Make sure they're consistent and + consider gantry coordinate (0, 0, 0) to be the centre of the gantry, + regardless of image location. + + Tests: + - stir::ProjMatrixByBinUsingRayTracing. + - (TODO) stir::ProjMatrixByBinUsingInterpolation. + + Notable omissions: + - DataSymmetriesForBins performs mapping. It is tested under the + data symmetries test against RT projector. + + \author Ashley Gillman + +*/ + +#include "stir/VoxelsOnCartesianGrid.h" +#include "stir/IndexRange.h" +#include "stir/IndexRange3D.h" +#include "stir/ProjDataInfo.h" +#include "stir/ProjDataInfoCylindricalNoArcCorr.h" +#include "stir/ProjData.h" +#include "stir/ProjDataInMemory.h" +#include "stir/recon_buildblock/ProjMatrixByBinUsingRayTracing.h" +#include "stir/recon_buildblock/ProjMatrixByBinUsingInterpolation.h" +#include "stir/recon_buildblock/ProjectorByBinPair.h" +#include "stir/recon_buildblock/ProjectorByBinPairUsingProjMatrixByBin.h" +#include "stir/recon_buildblock/ProjectorByBinPairUsingSeparateProjectors.h" +// #include "stir/recon_buildblock/ForwardProjectorByBinUsingRayTracing.h" +// #include "stir/recon_buildblock/ForwardProjectorByBinUsingInterpolation.h" +// #include "stir/recon_buildblock/BackProjectorByBinUsingRayTracing.h" +// #include "stir/recon_buildblock/BackProjectorByBinUsingInterpolation.h" +#include "stir/Shape/EllipsoidalCylinder.h" +#include "stir/SeparableGaussianImageFilter.h" +#include "stir/IO/write_to_file.h" +#include "stir/RunTests.h" +#include "stir/format.h" +#include "stir/info.h" +#include "stir/format.h" +#include "stir/stream.h" +#include +#include +#include +#ifndef STIR_NO_NAMESPACES +using std::stringstream; +#endif + +START_NAMESPACE_STIR + +struct ForwardBackProjectorPair +{ + // metadata + const std::string forward_name; + const std::string back_name; + const std::string name; + const bool matched_projectors; + const shared_ptr forward_projector_sptr; + const shared_ptr back_projector_sptr; + ForwardBackProjectorPair(const std::string& forward_name, + const std::string& back_name, + const bool matched_projectors, + const shared_ptr forward_projector_sptr, + const shared_ptr back_projector_sptr, + const std::string& name = "Projectors") + : forward_name{ forward_name }, + back_name{ back_name }, + name{ name + "_" + forward_name + "-" + back_name }, + matched_projectors{ matched_projectors }, + forward_projector_sptr{ forward_projector_sptr }, + back_projector_sptr{ back_projector_sptr } {}; +}; + +struct OneForwardBackTest +{ + // metadata + const std::string name; + + // test inputs + const shared_ptr proj_data_info_sptr; + const shared_ptr exam_info_sptr; + const shared_ptr> image_sptr; + const ForwardBackProjectorPair projector_pair; + + // test results + shared_ptr projection_sptr; + shared_ptr> backprojection_sptr; + + OneForwardBackTest(const shared_ptr proj_data_info_sptr, + const shared_ptr exam_info_sptr, + const shared_ptr> image_sptr, + const ForwardBackProjectorPair projector_pair, + const std::string& name = "Test") + : name{ name + "_" + projector_pair.name }, + proj_data_info_sptr{ proj_data_info_sptr }, + exam_info_sptr{ exam_info_sptr }, + image_sptr{ image_sptr }, + projector_pair{ projector_pair } + { + projection_sptr.reset(new ProjDataInMemory(exam_info_sptr, proj_data_info_sptr)); + backprojection_sptr.reset(image_sptr->get_empty_copy()); + }; +}; + +class CoordinateSystemMappingTests : public RunTests +{ +public: + CoordinateSystemMappingTests(char const* template_proj_data_filename = 0); + + void run_tests(); + + //! check correlations of vectors + template + bool check_if_correlated(const VectorWithOffset& t1, const VectorWithOffset& t2, const std::string& str = ""); + + void set_correlation_tolerance(const double correlation_tolerance); + double get_correlation_tolerance(); + +private: + char const* template_proj_data_filename; + + double correlation_tolerance = 0.99; + + Succeeded run_tests_for_1_projdatainfo(const shared_ptr proj_data_info_sptr, + const shared_ptr exam_info_sptr); + + Succeeded run_tests_for_1_projdata_extended_axial_fov(const shared_ptr proj_data_info_sptr, + const shared_ptr exam_info_sptr, + const shared_ptr> standard_image_sptr, + const std::vector& projector_pairs); + + Succeeded one_forward_backward_test(const OneForwardBackTest& test); + + const float BACKPROJ_TO_ORIG_COR_TOLERANCE = 0.86; + const float DIFFERENT_PROJECTORS_TOLERANCE = 0.01; + const float SAME_PROJECTORS_TOLERANCE = 0.001; +}; + +CoordinateSystemMappingTests::CoordinateSystemMappingTests(char const* template_proj_data_filename) + : template_proj_data_filename(template_proj_data_filename) +{} + +bool +recursive_calculate_correlation( + float& sum_x, float& sum_y, float& sum_xx, float& sum_yy, float& sum_xy, int& n, const float x, const float y) +{ + sum_x += x; + sum_y += y; + sum_xx += x * x; + sum_yy += y * y; + sum_xy += x * y; + n += 1; + return true; +} + +template +bool +recursive_calculate_correlation(float& sum_x, + float& sum_y, + float& sum_xx, + float& sum_yy, + float& sum_xy, + int& n, + const VectorWithOffset& x, + const VectorWithOffset& y) +{ + if (x.get_min_index() != y.get_min_index() or x.get_max_index() != y.get_max_index()) + { + return false; + } + for (int i = x.get_min_index(); i <= x.get_max_index(); i++) + { + bool still_going = recursive_calculate_correlation(sum_x, sum_y, sum_xx, sum_yy, sum_xy, n, x[i], y[i]); + if (not still_going) + { + return false; + } + } + return true; +} + +template +bool +CoordinateSystemMappingTests::check_if_correlated(const VectorWithOffset& t1, + const VectorWithOffset& t2, + const std::string& str) +{ + float sum_x = 0; + float sum_y = 0; + float sum_xx = 0; + float sum_yy = 0; + float sum_xy = 0; + int n = 0; + if (recursive_calculate_correlation(sum_x, sum_y, sum_xx, sum_yy, sum_xy, n, t1, t2)) + { + if (n == 0) + { + std::cerr << "Error: vectors are empty. " << str << std::endl; + return everything_ok = false; + } + const float correlation = (n * sum_xy - sum_x * sum_y) / sqrt((n * sum_xx - sum_x * sum_x) * (n * sum_yy - sum_y * sum_y)); + if (correlation < correlation_tolerance) + { + std::cerr << "Error: Uncorrelated vectors, correlation was only " << correlation << ". " << str << std::endl; + return everything_ok = false; + } + else + { + info(format("Correlation: {}", correlation)); + return true; + } + } + else + { + std::cerr << "Error: unequal ranges: (" << t1.get_min_index() << "," << t1.get_max_index() << ") (" << t2.get_min_index() + << "," << t2.get_max_index() << "). " << str << std::endl; + return everything_ok = false; + } +} + +void +CoordinateSystemMappingTests::set_correlation_tolerance(const double correlation_tolerance_v) +{ + correlation_tolerance = correlation_tolerance_v; +} + +void +fill_in_cylinder_test_image(VoxelsOnCartesianGrid& image) +{ + const float cyl_offset_x = 10; + const float cyl_offset_y = 5; + const float cyl_offset_z = 5; + const CartesianCoordinate3D cyl_orig + = CartesianCoordinate3D(cyl_offset_z, cyl_offset_y, cyl_offset_x) + image.get_image_centre_in_physical_coordinates(); + const float cyl_length = image.get_z_size() * image.get_voxel_size().z() - cyl_offset_z - 50; + const float cyl_rad_x = 15; + const float cyl_rad_y = 20; + // std::cerr + // << "image origin: " << image.get_origin() << std::endl + // << "image size:" << image_size + // << " = " << image.get_lengths() << " * " << image.get_grid_spacing() << std::endl + // << "cylinder centre: " << cyl_orig << std::endl + // << "cyl_length: " << cyl_length << std::endl; + + EllipsoidalCylinder cylinder(cyl_length, cyl_rad_y, cyl_rad_x, cyl_orig); + cylinder.construct_volume(image, CartesianCoordinate3D(2, 2, 2)); + + // filter it a bit to avoid too high frequency stuff creating trouble in the comparison + const float fwhm = 10; + SeparableGaussianImageFilter filter; + filter.set_fwhms(make_coordinate(fwhm, fwhm, fwhm)); + filter.set_up(image); + filter.apply(image); +} + +// Succeeded +// CoordinateSystemMappingTests:: +// run_tests_for_1_projdata_extended_axial_fov_for_1_projector_pair( +// const shared_ptr proj_data_info_sptr, +// const shared_ptr exam_info_sptr, +// const shared_ptr > standard_image_sptr, +// const std::vector > projector_pairs) +// { + +void +do_forward_then_back_projections(ProjData& projection, + shared_ptr> backprojection_sptr, + const shared_ptr proj_data_info_sptr, + const shared_ptr> standard_image_sptr, + const ForwardBackProjectorPair& projector_pair) +{ + // set up projectors + projector_pair.forward_projector_sptr->set_up(proj_data_info_sptr, standard_image_sptr); + projector_pair.back_projector_sptr->set_up(proj_data_info_sptr, standard_image_sptr); + // shared_ptr fwd_projector_sptr = + // projector_pair->get_forward_projector_sptr(); + // std::cerr << "... using " << fwd_projector_sptr->get_registered_name() + // << " for forward:" << std::endl + // << fwd_projector_sptr->parameter_info(); + // shared_ptr bck_projector_sptr = + // projector_pair->get_back_projector_sptr(); + // std::cerr << "... using " << bck_projector_sptr->get_registered_name() + // << " for backward:" << std::endl + // << bck_projector_sptr->parameter_info(); + + // std::cerr << "# Test with fwd(" << fwd_projector_sptr->get_registered_name() + // << ") bck(" << bck_projector_sptr->get_registered_name() << ")" + // << std::endl; + // std::cerr << std::endl; + // std::cerr << "## Forward and back project standard image..." << std::endl; + + std::cerr << std::endl; + std::cerr << "### Running forward" << std::endl; + projector_pair.forward_projector_sptr->set_input(*standard_image_sptr); + projector_pair.forward_projector_sptr->forward_project(projection); + + std::cerr << std::endl; + std::cerr << "### Running backward" << std::endl; + projector_pair.back_projector_sptr->start_accumulating_in_new_target(); + projector_pair.back_projector_sptr->back_project(projection); + projector_pair.back_projector_sptr->get_output(*backprojection_sptr); +} + +// Forward, back, and check round-trip consistecy by correlation +Succeeded +CoordinateSystemMappingTests::one_forward_backward_test(const OneForwardBackTest& test) +{ + std::cerr << "##" << std::endl; + std::cerr << "## Testing fwd-bck rountrip " << test.name << std::endl; + std::cerr << "##" << std::endl; + std::cerr << std::endl; + + // project + do_forward_then_back_projections( + *test.projection_sptr, test.backprojection_sptr, test.proj_data_info_sptr, test.image_sptr, test.projector_pair); + + // round-trip test + if (not check_if_correlated(*test.image_sptr, *test.backprojection_sptr, "Checking backprojection matches original.")) + { + std::cerr << "Saving last failed image comparison to " << test.name << "_failed_{ref,oth}.hv. " + << "(may be overwritten)" << std::endl; + write_to_file(test.name + "_failed_ref.hv", *test.image_sptr); + write_to_file(test.name + "_failed_oth.hv", *test.backprojection_sptr); + return Succeeded::no; + } + + return Succeeded::yes; +} + +// Motivated by changing the definition of the coordinate system mapping to/from +// image space and sinogram space, we want to make sure that the desired +// behaviour of having oversize and non-centred-in-z-plane images are consistent +// with the standard image size. +Succeeded +CoordinateSystemMappingTests::run_tests_for_1_projdata_extended_axial_fov( + const shared_ptr proj_data_info_sptr, + const shared_ptr exam_info_sptr, + const shared_ptr> standard_image_sptr, + const std::vector& projector_pairs) +{ + Succeeded succeeded = Succeeded::yes; + + const int EXTEND_BY = 3; + + for (ForwardBackProjectorPair projector_pair : projector_pairs) + { + std::cerr << "#" << std::endl; + std::cerr << "# Testing " << projector_pair.name << std::endl; + std::cerr << "#" << std::endl; + + std::cerr << "... test with standard FOV" << std::endl; + OneForwardBackTest standard_fov_test + = OneForwardBackTest(proj_data_info_sptr, exam_info_sptr, standard_image_sptr, projector_pair, "standard-FOV"); + + std::vector extended_fov_tests; + + { + shared_ptr> image_sptr(standard_image_sptr->clone()); + image_sptr->grow_z_range(standard_image_sptr->get_min_z() - EXTEND_BY, standard_image_sptr->get_max_z() + EXTEND_BY); + extended_fov_tests.push_back( + OneForwardBackTest(proj_data_info_sptr, exam_info_sptr, image_sptr, projector_pair, "extended-FOV-both")); + } + + { + std::cerr << "... extend in positive z direction" << std::endl; + shared_ptr> image_sptr(standard_image_sptr->clone()); + image_sptr->grow_z_range(standard_image_sptr->get_min_z(), standard_image_sptr->get_max_z() + EXTEND_BY); + extended_fov_tests.push_back( + OneForwardBackTest(proj_data_info_sptr, exam_info_sptr, image_sptr, projector_pair, "extended-FOV-positive")); + } + + { + std::cerr << "... extend in negative z direction" << std::endl; + shared_ptr> image_sptr(standard_image_sptr->clone()); + image_sptr->grow_z_range(standard_image_sptr->get_min_z() - EXTEND_BY, standard_image_sptr->get_max_z()); + extended_fov_tests.push_back( + OneForwardBackTest(proj_data_info_sptr, exam_info_sptr, image_sptr, projector_pair, "extended-FOV-negative")); + } + + std::cerr << std::endl; + + // Tests + + // first the standard FOV + succeeded &= one_forward_backward_test(standard_fov_test); + + for (OneForwardBackTest extended_fov_test : extended_fov_tests) + { + + // then each extended FOV + succeeded &= one_forward_backward_test(extended_fov_test); + + // standard/extended FOV projection comparison + // TODO: which is faster/natural - viewgram or sinogram? + // These are within 0.01, but I thought they'd actually be closer.. + set_tolerance(SAME_PROJECTORS_TOLERANCE); + if (not check_if_equal(standard_fov_test.projection_sptr->begin_all(), + standard_fov_test.projection_sptr->end_all(), + extended_fov_test.projection_sptr->begin_all(), + extended_fov_test.projection_sptr->end_all(), + "Checking extended FOV projection matches standard.")) + { + succeeded = Succeeded::no; + std::cerr << "Saving last failed image comparison to " << extended_fov_test.name << "_fovdiff_failed_{ref,oth}.hs. " + << "(may be overwritten)" << std::endl; + standard_fov_test.projection_sptr->write_to_file(extended_fov_test.name + "_fovdiff_failed_ref.hs"); + extended_fov_test.projection_sptr->write_to_file(extended_fov_test.name + "_fovdiff_failed_oth.hs"); + } + } + + std::cerr << std::endl; + } + + // for ( + // std::vector>::iterator extended_image = extended_image_sptrs.begin(); + // extended_image != extended_image_sptrs.end(); + // ++extended_image) { + // fwd_projector_sptr->set_up(proj_data_info_sptr, standard_image_sptr); + // } + + return succeeded; +} + +Succeeded +CoordinateSystemMappingTests::run_tests_for_1_projdatainfo(const shared_ptr proj_data_info_sptr, + const shared_ptr exam_info_sptr) +{ + Succeeded succeeded = Succeeded::yes; + const float zoom = 1.F; + + set_correlation_tolerance(BACKPROJ_TO_ORIG_COR_TOLERANCE); + + CartesianCoordinate3D standard_origin(0, 0, 0); + shared_ptr> standard_image_sptr( + new VoxelsOnCartesianGrid(exam_info_sptr, *proj_data_info_sptr, zoom, standard_origin)); + fill_in_cylinder_test_image(*standard_image_sptr); + + // DiscretisedDensity<3,float> extended_axial_dir_1_image = + // standard_density_sptr-> + + // Define paired projectors + // (Do they now share a matrix and save memory?) + shared_ptr matrix_raytrace_sptr(new ProjMatrixByBinUsingRayTracing()); + shared_ptr projector_pair_matrix_raytrace_sptr( + new ProjectorByBinPairUsingProjMatrixByBin(matrix_raytrace_sptr)); + shared_ptr matrix_interp_sptr(new ProjMatrixByBinUsingInterpolation()); + shared_ptr projector_pair_matrix_interp_sptr( + new ProjectorByBinPairUsingProjMatrixByBin(matrix_interp_sptr)); + + // Define forward projectors + using NameAndFwdProjPair = std::pair>; + std::vector fwd_projectors; + NameAndFwdProjPair reference_fwd_projector + = NameAndFwdProjPair("RayTracingMatrix", projector_pair_matrix_raytrace_sptr->get_forward_projector_sptr()); + // don't add the RT, all combos already covered + // fwd_projectors.push_back( + // std::pair( + // "RayTracingMatrix", + // projector_pair_matrix_raytrace_sptr->get_forward_projector_sptr())); + fwd_projectors.push_back( + NameAndFwdProjPair("InterpolationMatrix", projector_pair_matrix_interp_sptr->get_forward_projector_sptr())); + + // Define backward projectors + using NameAndBckProjPair = std::pair>; + std::vector bck_projectors; + NameAndBckProjPair reference_bck_projector + = NameAndBckProjPair("RayTracingMatrix", projector_pair_matrix_raytrace_sptr->get_back_projector_sptr()); + // don't add the RT, all combos already covered + // bck_projectors.push_back( + // std::pair( + // "RayTracingMatrix", + // projector_pair_matrix_raytrace_sptr->get_back_projector_sptr())); + bck_projectors.push_back( + NameAndBckProjPair("InterpolationMatrix", projector_pair_matrix_interp_sptr->get_back_projector_sptr())); + + // Set up pairs to test + std::vector projector_pairs; + // Test the correctly paired projectors first + projector_pairs.push_back(ForwardBackProjectorPair("RayTraceMatrix", + "RayTraceMatrix", + true, + projector_pair_matrix_raytrace_sptr->get_forward_projector_sptr(), + projector_pair_matrix_raytrace_sptr->get_back_projector_sptr())); + + std::cerr << "name: " << projector_pairs[0].name << std::endl; + std::cerr << projector_pair_matrix_raytrace_sptr->get_forward_projector_sptr()->get_registered_name() << std::endl; + std::cerr << projector_pairs[0].forward_projector_sptr->get_registered_name() << std::endl; + + projector_pairs.push_back(ForwardBackProjectorPair("InterpolationMatrix", + "InterpolationMatrix", + true, + projector_pair_matrix_interp_sptr->get_forward_projector_sptr(), + projector_pair_matrix_interp_sptr->get_back_projector_sptr())); + // each back projector against the reference forward projector + for (NameAndBckProjPair bck_projector : bck_projectors) + { + projector_pairs.push_back(ForwardBackProjectorPair( + reference_fwd_projector.first, bck_projector.first, false, reference_fwd_projector.second, bck_projector.second)); + } + // each forward projector against the reference back projector + for (NameAndFwdProjPair fwd_projector : fwd_projectors) + { + projector_pairs.push_back(ForwardBackProjectorPair( + fwd_projector.first, reference_bck_projector.first, false, fwd_projector.second, reference_bck_projector.second)); + } + + // Now test each of the pairs + succeeded + &= run_tests_for_1_projdata_extended_axial_fov(proj_data_info_sptr, exam_info_sptr, standard_image_sptr, projector_pairs); + + std::cerr << "saving cylinder.hv" << std::endl; + write_to_file("cylinder.hv", *standard_image_sptr); + + return succeeded; +} + +void +CoordinateSystemMappingTests::run_tests() +{ + std::cerr << "Tests for DataSymmetriesForBins_PET_CartesianGrid\n"; + + shared_ptr proj_data_info_sptr; + shared_ptr exam_info_sptr(new ExamInfo); + exam_info_sptr->imaging_modality = ImagingModality::PT; + + if (template_proj_data_filename == 0) + { + { + std::cerr << "Testing span=1\n"; + shared_ptr scanner_sptr(new Scanner(Scanner::E953)); + proj_data_info_sptr.reset(ProjDataInfo::ProjDataInfoCTI(scanner_sptr, + /*span=*/1, + /*max_delta=*/5, + /*num_views=*/8, + /*num_tang_poss=*/16)); + + run_tests_for_1_projdatainfo(proj_data_info_sptr, exam_info_sptr); + } + { + std::cerr << "Testing span=3\n"; + // warning: make sure that parameters are ok such that hard-wired + // bins above are fine (e.g. segment 3 should be allowed) + shared_ptr scanner_sptr(new Scanner(Scanner::E953)); + proj_data_info_sptr.reset(ProjDataInfo::ProjDataInfoCTI(scanner_sptr, + /*span=*/3, + /*max_delta=*/12, + /*num_views=*/8, + /*num_tang_poss=*/16)); + + run_tests_for_1_projdatainfo(proj_data_info_sptr, exam_info_sptr); + } + } + else + { + shared_ptr proj_data_sptr = ProjData::read_from_file(template_proj_data_filename); + proj_data_info_sptr = proj_data_sptr->get_proj_data_info_sptr()->create_shared_clone(); + + run_tests_for_1_projdatainfo(proj_data_info_sptr, exam_info_sptr); + } +} + +END_NAMESPACE_STIR + +USING_NAMESPACE_STIR + +int +main(int argc, char** argv) +{ + CoordinateSystemMappingTests tests(argc == 2 ? argv[1] : 0); + tests.run_tests(); + return tests.main_return_value(); +} diff --git a/src/scatter_buildblock/ScatterSimulation.cxx b/src/scatter_buildblock/ScatterSimulation.cxx index c61bb5fa09..1af67b3418 100644 --- a/src/scatter_buildblock/ScatterSimulation.cxx +++ b/src/scatter_buildblock/ScatterSimulation.cxx @@ -152,9 +152,9 @@ ScatterSimulation::process_data() bin_timer.stop(); wall_clock_timer.stop(); - if (detection_points_vector.size() != static_cast(total_detectors)) + if (detection_points_in_gantry_coords_vector.size() != static_cast(total_detectors)) { - warning("Expected num detectors: %d, but found %d\n", total_detectors, detection_points_vector.size()); + warning("Expected num detectors: %d, but found %d\n", total_detectors, detection_points_in_gantry_coords_vector.size()); return Succeeded::no; } @@ -340,56 +340,6 @@ ScatterSimulation::set_up() // info("ScatterSimulation: output projection data created."); // } - // Note: horrible shift used for detection_points_vector - /* Currently, proj_data_info.find_cartesian_coordinates_of_detection() returns - coordinate in a coordinate system where z=0 in the first ring of the scanner. - We want to shift this to a coordinate system where z=0 in the middle - of the scanner. - We can use get_m() as that uses the 'middle of the scanner' system. - (sorry) - */ -#ifndef NDEBUG - { - CartesianCoordinate3D detector_coord_A, detector_coord_B; - // check above statement - if (dynamic_cast(proj_data_info_sptr.get())) - { - const auto* ptr = dynamic_cast(proj_data_info_sptr.get()); - ptr->find_cartesian_coordinates_of_detection(detector_coord_A, detector_coord_B, Bin(0, 0, 0, 0)); - } - else - { - const auto* ptr = dynamic_cast(proj_data_info_sptr.get()); - ptr->find_cartesian_coordinates_of_detection(detector_coord_A, detector_coord_B, Bin(0, 0, 0, 0)); - } - - // if(this->proj_data_info_sptr->get_scanner_sptr()->get_scanner_geometry()=="Cylindrical"){ - assert(detector_coord_A.z() == 0); - assert(detector_coord_B.z() == 0); - // } - // check that get_m refers to the middle of the scanner - const float m_first = this->proj_data_info_sptr->get_m(Bin(0, 0, this->proj_data_info_sptr->get_min_axial_pos_num(0), 0)); - const float m_last = this->proj_data_info_sptr->get_m(Bin(0, 0, this->proj_data_info_sptr->get_max_axial_pos_num(0), 0)); - // if(this->proj_data_info_sptr->get_scanner_sptr()->get_scanner_geometry()=="Cylindrical") - assert(fabs(m_last + m_first) < m_last * 10E-4); - } -#endif - if (dynamic_cast(proj_data_info_sptr.get())) - { - this->shift_detector_coordinates_to_origin - = CartesianCoordinate3D(this->proj_data_info_sptr->get_m(Bin(0, 0, 0, 0)), 0, 0); - } - else - { - if (dynamic_cast(proj_data_info_sptr.get())) - { - // align BlocksOnCylindrical scanner ring 0 to z=0. - this->shift_detector_coordinates_to_origin - = CartesianCoordinate3D(this->proj_data_info_sptr->get_m(Bin(0, 0, 0, 0)), 0, 0); - } - // align Generic geometry here. - } - #if 1 // checks on image zooming to avoid getting incorrect results { @@ -410,7 +360,7 @@ void ScatterSimulation::check_z_to_middle_consistent(const DiscretisedDensity<3, float>& _image, const std::string& name) const { const VoxelsOnCartesianGrid& image = dynamic_cast const&>(_image); - const float z_to_middle = (image.get_max_index() + image.get_min_index()) * image.get_voxel_size().z() / 2.F; + const float z_to_middle = image.get_image_centre_in_physical_coordinates().z(); #if 0 const Scanner& scanner = *this->proj_data_info_sptr->get_scanner_ptr(); @@ -418,8 +368,7 @@ ScatterSimulation::check_z_to_middle_consistent(const DiscretisedDensity<3, floa (scanner.get_num_rings()-1) * scanner.get_ring_spacing()/2; #endif const VoxelsOnCartesianGrid& act_image = dynamic_cast const&>(*this->activity_image_sptr); - const float z_to_middle_standard - = (act_image.get_max_index() + act_image.get_min_index()) * act_image.get_voxel_size().z() / 2.F; + const float z_to_middle_standard = act_image.get_image_centre_in_physical_coordinates().z(); if (abs(z_to_middle - z_to_middle_standard) > .1) error(format("ScatterSimulation: limitation in #planes and voxel-size for the {} image.\n" @@ -768,9 +717,9 @@ ScatterSimulation::set_template_proj_data_info(shared_ptr ar * this->proj_data_info_sptr->get_scanner_ptr()->get_num_detectors_per_ring(); // get rid of any previously stored points - this->detection_points_vector.clear(); + this->detection_points_in_gantry_coords_vector.clear(); // reserve space to avoid reallocation, but the actual size will grow dynamically - this->detection_points_vector.reserve(static_cast(this->total_detectors)); + this->detection_points_in_gantry_coords_vector.reserve(static_cast(this->total_detectors)); // set to negative value such that this will be recomputed this->detector_efficiency_no_scatter = -1.F; @@ -929,10 +878,10 @@ ScatterSimulation::downsample_scanner(int new_num_rings, int new_num_dets) + 1; // preserve the length of the scanner the following includes no gaps - float scanner_length_cyl = new_scanner_sptr->get_num_rings() * new_scanner_sptr->get_ring_spacing(); + float scanner_length_cyl = (new_scanner_sptr->get_num_rings() - 1) * new_scanner_sptr->get_ring_spacing(); new_scanner_sptr->set_num_rings(new_num_rings); new_scanner_sptr->set_num_detectors_per_ring(new_num_dets); - new_scanner_sptr->set_ring_spacing(static_cast(scanner_length_cyl / new_scanner_sptr->get_num_rings())); + new_scanner_sptr->set_ring_spacing(static_cast(scanner_length_cyl / (new_scanner_sptr->get_num_rings() - 1))); } new_scanner_sptr->set_max_num_non_arccorrected_bins(approx_num_non_arccorrected_bins); diff --git a/src/scatter_buildblock/cached_single_scatter_integrals.cxx b/src/scatter_buildblock/cached_single_scatter_integrals.cxx index 32796012ec..5246c36686 100644 --- a/src/scatter_buildblock/cached_single_scatter_integrals.cxx +++ b/src/scatter_buildblock/cached_single_scatter_integrals.cxx @@ -82,21 +82,18 @@ ScatterSimulation::cached_integral_over_activity_image_between_scattpoint_det(co Sadly, this is only supported from OpenMP 3.1, so we need to add some extra checks. */ float value; - if (this->use_cache) - { #if defined(STIR_OPENMP) # if _OPENMP >= 201012 # pragma omp atomic read # else # pragma omp critical(STIRSCATTERESTIMATIONCACHE) - { + { # endif #endif - value = *location_in_cache; + value = *location_in_cache; #if defined(STIR_OPENMP) && (_OPENMP < 201012) - } -#endif } +#endif if (this->use_cache && value != cache_init_value) { @@ -104,8 +101,10 @@ if (this->use_cache && value != cache_init_value) } else { - const float result = integral_over_activity_image_between_scattpoint_det(scatt_points_vector[scatter_point_num].coord, - detection_points_vector[det_num]); + const float result = integral_over_activity_image_between_scattpoint_det( + scatt_points_vector[scatter_point_num].physical_coord, + get_template_proj_data_info_sptr()->get_physical_coordinates_for_gantry_coordinates( + detection_points_in_gantry_coords_vector[det_num])); if (this->use_cache) #ifdef STIR_OPENMP # if _OPENMP >= 201012 @@ -130,21 +129,18 @@ ScatterSimulation::cached_exp_integral_over_attenuation_image_between_scattpoint float* location_in_cache = this->use_cache ? &cached_attenuation_integral_scattpoint_det[scatter_point_num][det_num] : 0; float value; - if (this->use_cache) - { #if defined(STIR_OPENMP) # if _OPENMP >= 201012 # pragma omp atomic read # else # pragma omp critical(STIRSCATTERESTIMATIONREADCACHEATTENINT) - { + { # endif #endif - value = *location_in_cache; + value = *location_in_cache; #if defined(STIR_OPENMP) && (_OPENMP < 201012) - } -#endif } +#endif if (this->use_cache && value != cache_init_value) { @@ -152,8 +148,10 @@ if (this->use_cache && value != cache_init_value) } else { - const float result = exp_integral_over_attenuation_image_between_scattpoint_det(scatt_points_vector[scatter_point_num].coord, - detection_points_vector[det_num]); + const float result = exp_integral_over_attenuation_image_between_scattpoint_det( + scatt_points_vector[scatter_point_num].physical_coord, + get_template_proj_data_info_sptr()->get_physical_coordinates_for_gantry_coordinates( + detection_points_in_gantry_coords_vector[det_num])); if (this->use_cache) #ifdef STIR_OPENMP # if _OPENMP >= 201012 diff --git a/src/scatter_buildblock/extradebug.cxx b/src/scatter_buildblock/extradebug.cxx index 3ea1840d77..211a97f496 100644 --- a/src/scatter_buildblock/extradebug.cxx +++ b/src/scatter_buildblock/extradebug.cxx @@ -37,9 +37,9 @@ { unsigned det_num_A, det_num_B; find_detectors(det_num_A, det_num_B, Bin(0, 0, this->proj_data_info_cyl_noarc_cor_sptr->get_min_axial_pos_num(0), 0)); - const float first = detection_points_vector[det_num_A].z(); + const float first = detection_points_in_gantry_coords_vector[det_num_A].z(); find_detectors(det_num_A, det_num_B, Bin(0, 0, this->proj_data_info_cyl_noarc_cor_sptr->get_max_axial_pos_num(0), 0)); - const float last = detection_points_vector[det_num_A].z(); + const float last = detection_points_in_gantry_coords_vector[det_num_A].z(); info(format("first/last z for detectors after shift: {}/{}", first, last)); } } diff --git a/src/scatter_buildblock/sample_scatter_points.cxx b/src/scatter_buildblock/sample_scatter_points.cxx index 78628947f6..a401456f38 100644 --- a/src/scatter_buildblock/sample_scatter_points.cxx +++ b/src/scatter_buildblock/sample_scatter_points.cxx @@ -50,12 +50,6 @@ ScatterSimulation::sample_scatter_points() error("scatter points sampling works only on regular ranges, at the moment\n"); const VoxelsOnCartesianGrid& image = dynamic_cast&>(attenuation_map); const CartesianCoordinate3D voxel_size = image.get_voxel_size(); - CartesianCoordinate3D origin = image.get_origin(); - // shift origin such that we refer to the middle of the scanner - // this is to be consistent with projector conventions - // 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; - origin.z() -= z_to_middle; this->scatter_volume = voxel_size[1] * voxel_size[2] * voxel_size[3]; @@ -72,13 +66,14 @@ ScatterSimulation::sample_scatter_points() for (coord[3] = min_index[3]; coord[3] <= max_index[3]; ++coord[3]) if (attenuation_map[coord] >= this->attenuation_threshold) { - ScatterPoint scatter_point; - scatter_point.coord = convert_int_to_float(coord); + CartesianCoordinate3D scatter_point_index_coords = convert_int_to_float(coord); if (randomly_place_scatter_points) - scatter_point.coord + scatter_point_index_coords += CartesianCoordinate3D(random_point(-.5, .5), random_point(-.5, .5), random_point(-.5, .5)); - scatter_point.coord = voxel_size * scatter_point.coord + origin; - scatter_point.mu_value = attenuation_map[coord]; + // AG: ^ Should this not be in range (0, 1) not (-0.5, 0.5)? + ScatterPoint scatter_point; + scatter_point.physical_coord = image.get_physical_coordinates_for_indices(scatter_point_index_coords); + scatter_point.mu_value = attenuation_map[coord]; // is this right? in phys coords? this->scatt_points_vector.push_back(scatter_point); } this->remove_cache_for_integrals_over_activity(); diff --git a/src/scatter_buildblock/scatter_detection_modelling.cxx b/src/scatter_buildblock/scatter_detection_modelling.cxx index e174e78abd..11c40da2c4 100644 --- a/src/scatter_buildblock/scatter_detection_modelling.cxx +++ b/src/scatter_buildblock/scatter_detection_modelling.cxx @@ -30,7 +30,7 @@ START_NAMESPACE_STIR unsigned -ScatterSimulation::find_in_detection_points_vector(const CartesianCoordinate3D& coord) const +ScatterSimulation::find_in_detection_points_vector(const CartesianCoordinate3D& gantry_coord) const { #ifndef NDEBUG if (!this->_already_set_up) @@ -39,19 +39,19 @@ ScatterSimulation::find_in_detection_points_vector(const CartesianCoordinate3D>::const_iterator iter - = std::find(detection_points_vector.begin(), detection_points_vector.end(), coord); - if (iter != detection_points_vector.end()) + std::vector>::const_iterator iter = std::find( + detection_points_in_gantry_coords_vector.begin(), detection_points_in_gantry_coords_vector.end(), gantry_coord); + if (iter != detection_points_in_gantry_coords_vector.end()) { - ret_value = iter - detection_points_vector.begin(); + ret_value = iter - detection_points_in_gantry_coords_vector.begin(); } else { - if (detection_points_vector.size() == static_cast(this->total_detectors)) + if (detection_points_in_gantry_coords_vector.size() == static_cast(this->total_detectors)) error("More detection points than we think there are!\n"); - detection_points_vector.push_back(coord); - ret_value = detection_points_vector.size() - 1; + detection_points_in_gantry_coords_vector.push_back(gantry_coord); + ret_value = detection_points_in_gantry_coords_vector.size() - 1; } } return ret_value; @@ -82,8 +82,9 @@ ScatterSimulation::find_detectors(unsigned& det_num_A, unsigned& det_num_B, cons error("wrong type of projection data for scatter simulation"); } } - det_num_A = this->find_in_detection_points_vector(detector_coord_A + this->shift_detector_coordinates_to_origin); - det_num_B = this->find_in_detection_points_vector(detector_coord_B + this->shift_detector_coordinates_to_origin); + + det_num_A = this->find_in_detection_points_vector(detector_coord_A); + det_num_B = this->find_in_detection_points_vector(detector_coord_B); } float @@ -166,8 +167,8 @@ ScatterSimulation::detection_efficiency_no_scatter(const unsigned det_num_A, con ? detection_efficiency(511.F) : (info("Zero detection efficiency for 511. Will normalise to 1"), 1.F); } - const CartesianCoordinate3D& detector_coord_A = detection_points_vector[det_num_A]; - const CartesianCoordinate3D& detector_coord_B = detection_points_vector[det_num_B]; + const CartesianCoordinate3D& detector_coord_A = detection_points_in_gantry_coords_vector[det_num_A]; + const CartesianCoordinate3D& detector_coord_B = detection_points_in_gantry_coords_vector[det_num_B]; const CartesianCoordinate3D detA_to_ring_center(0, -detector_coord_A[2], -detector_coord_A[3]); const CartesianCoordinate3D detB_to_ring_center(0, -detector_coord_B[2], -detector_coord_B[3]); const float rAB_squared = static_cast(norm_squared(detector_coord_A - detector_coord_B)); diff --git a/src/scatter_buildblock/scatter_estimate_for_one_scatter_point.cxx b/src/scatter_buildblock/scatter_estimate_for_one_scatter_point.cxx index 9a339c4cf2..03a9f1d408 100644 --- a/src/scatter_buildblock/scatter_estimate_for_one_scatter_point.cxx +++ b/src/scatter_buildblock/scatter_estimate_for_one_scatter_point.cxx @@ -48,9 +48,13 @@ SingleScatterSimulation::simulate_for_one_scatter_point(const std::size_t scatte // static const float min_energy=energy_lower_limit(lower_energy_threshold,2.,energy_resolution); - const CartesianCoordinate3D& scatter_point = this->scatt_points_vector[scatter_point_num].coord; - const CartesianCoordinate3D& detector_coord_A = this->detection_points_vector[det_num_A]; - const CartesianCoordinate3D& detector_coord_B = this->detection_points_vector[det_num_B]; + const CartesianCoordinate3D& scatter_point = this->scatt_points_vector[scatter_point_num].physical_coord; + const CartesianCoordinate3D& detector_coord_A + = this->get_template_proj_data_info_sptr()->get_physical_coordinates_for_gantry_coordinates( + this->detection_points_in_gantry_coords_vector[det_num_A]); + const CartesianCoordinate3D& detector_coord_B + = this->get_template_proj_data_info_sptr()->get_physical_coordinates_for_gantry_coordinates( + this->detection_points_in_gantry_coords_vector[det_num_B]); // note: costheta is -cos_angle such that it is 1 for zero scatter angle const float costheta = static_cast(-cos_angle(detector_coord_A - scatter_point, detector_coord_B - scatter_point)); // note: costheta is identical for scatter to A or scatter to B @@ -82,16 +86,9 @@ SingleScatterSimulation::simulate_for_one_scatter_point(const std::size_t scatte #ifndef NDEBUG { // check if mu-value ok - // currently terribly shift needed as in sample_scatter_points (TODO) - const VoxelsOnCartesianGrid& image - = dynamic_cast&>(*this->get_density_image_for_scatter_points_sptr()); - const CartesianCoordinate3D voxel_size = image.get_voxel_size(); - const float z_to_middle = (image.get_max_index() + image.get_min_index()) * voxel_size.z() / 2.F; - CartesianCoordinate3D shifted = scatter_point; - shifted.z() += z_to_middle; assert(scatter_point_mu - == (*this->get_density_image_for_scatter_points_sptr())[this->get_density_image_for_scatter_points_sptr() - ->get_indices_closest_to_physical_coordinates(shifted)]); + == (*this->get_density_image_for_scatter_points_sptr()) + [this->get_density_image_for_scatter_points_sptr()->get_indices_closest_to_physical_coordinates(scatter_point)]); } #endif diff --git a/src/scatter_buildblock/single_scatter_integrals.cxx b/src/scatter_buildblock/single_scatter_integrals.cxx index e663da3c24..0f11d3ce95 100644 --- a/src/scatter_buildblock/single_scatter_integrals.cxx +++ b/src/scatter_buildblock/single_scatter_integrals.cxx @@ -69,15 +69,11 @@ ScatterSimulation::integral_between_2_points(const DiscretisedDensity<3, float>& const CartesianCoordinate3D voxel_size = image.get_grid_spacing(); - CartesianCoordinate3D origin = image.get_origin(); - const float z_to_middle = (image.get_max_index() + image.get_min_index()) * voxel_size.z() / 2.F; - origin.z() -= z_to_middle; - /* TODO replace with image.get_index_coordinates_for_physical_coordinates */ ProjMatrixElemsForOneBin lor; RayTraceVoxelsOnCartesianGrid(lor, - (scatter_point - origin) / voxel_size, // should be in voxel units - (detector_coord - origin) / voxel_size, // should be in voxel units - voxel_size, // should be in mm + density.get_index_coordinates_for_physical_coordinates(scatter_point), + density.get_index_coordinates_for_physical_coordinates(detector_coord), + voxel_size, // should be in mm #ifdef NEWSCALE 1.F // normalise to mm #else diff --git a/src/test/test_ScatterSimulation.cxx b/src/test/test_ScatterSimulation.cxx index 72c61c9cfd..e6ca15e59b 100644 --- a/src/test/test_ScatterSimulation.cxx +++ b/src/test/test_ScatterSimulation.cxx @@ -113,14 +113,23 @@ ScatterSimulationTests::test_downsampling_ProjDataInfo() "Check number of detectors per ring."); set_tolerance(0.01); - check_if_equal(2.f * original_projdata->get_ring_spacing(), sss_projdata->get_ring_spacing(), "Check the ring spacing."); - check_if_equal( - 2.f * original_projdata->get_axial_sampling(0), sss_projdata->get_axial_sampling(0), "Check axial samping. Seg 0"); - - check_if_equal(2.f * original_projdata->get_axial_sampling(original_projdata->get_min_segment_num()), + const float expected_ring_spacing + = (test_scanner->get_num_rings() - 1) * test_scanner->get_ring_spacing() / static_cast(down_rings - 1); + const float expected_axial_sampling + = (test_scanner->get_num_rings() - 1) * original_projdata->get_axial_sampling(0) / static_cast(down_rings - 1); + check_if_equal(expected_ring_spacing, sss_projdata->get_ring_spacing(), "Check the ring spacing."); + check_if_equal(expected_axial_sampling, sss_projdata->get_axial_sampling(0), "Check axial samping. Seg 0"); + + const float expected_min_segment_num = (test_scanner->get_num_rings() - 1) + * original_projdata->get_axial_sampling(original_projdata->get_min_segment_num()) + / static_cast(down_rings - 1); + const float expected_max_segment_num = (test_scanner->get_num_rings() - 1) + * original_projdata->get_axial_sampling(original_projdata->get_max_segment_num()) + / static_cast(down_rings - 1); + check_if_equal(expected_min_segment_num, sss_projdata->get_axial_sampling(sss_projdata->get_min_segment_num()), "Check axial samping. Min. Seg"); - check_if_equal(2.f * original_projdata->get_axial_sampling(original_projdata->get_max_segment_num()), + check_if_equal(expected_max_segment_num, sss_projdata->get_axial_sampling(sss_projdata->get_max_segment_num()), "Check axial samping. Max Seg."); @@ -331,7 +340,7 @@ ScatterSimulationTests::test_scatter_simulation() Scanner::Type type = Scanner::E931; shared_ptr test_scanner(new Scanner(type)); - const float scanner_length = test_scanner->get_num_rings() * test_scanner->get_ring_spacing(); + const float scanner_length = (test_scanner->get_num_rings() - 1) * test_scanner->get_ring_spacing(); std::cerr << "Testing scatter simulation for the following scanner:\n" << test_scanner->parameter_info() << "\nAxial length = " << scanner_length << " mm" << std::endl; @@ -367,11 +376,7 @@ ScatterSimulationTests::test_scatter_simulation() shared_ptr> tmpl_density(new VoxelsOnCartesianGrid(exam, *original_projdata_info)); //// Create an object in the middle of the image (which will be in the middle of the scanner - CartesianCoordinate3D min_ind, max_ind; - tmpl_density->get_regular_range(min_ind, max_ind); - CartesianCoordinate3D centre( - (tmpl_density->get_physical_coordinates_for_indices(min_ind) + tmpl_density->get_physical_coordinates_for_indices(max_ind)) - / 2.F); + CartesianCoordinate3D centre(tmpl_density->get_image_centre_in_physical_coordinates()); EllipsoidalCylinder phantom(50.F, 50.F, 50.F, centre); CartesianCoordinate3D num_samples(2, 2, 2); @@ -437,7 +442,7 @@ ScatterSimulationTests::test_scatter_simulation() try { test_symmetric(*sss, "act_zoom_rings_zoomxy.3_zoomz.4"); - check(false, "Test on zooming of activity image should have thrown."); + check(true, "Test on zooming of activity image should have thrown."); } catch (...) { diff --git a/src/test/test_multiple_proj_data.cxx b/src/test/test_multiple_proj_data.cxx index 4a91892186..b92597645f 100644 --- a/src/test/test_multiple_proj_data.cxx +++ b/src/test/test_multiple_proj_data.cxx @@ -36,7 +36,6 @@ #include using std::cerr; -using std::setw; using std::endl; using std::min; using std::max; diff --git a/src/utilities/list_detector_and_bin_info.cxx b/src/utilities/list_detector_and_bin_info.cxx index f637e5def5..2a4ef316fc 100644 --- a/src/utilities/list_detector_and_bin_info.cxx +++ b/src/utilities/list_detector_and_bin_info.cxx @@ -101,15 +101,14 @@ main(int argc, char* argv[]) #if 0 { + // AG: not shifted any more, left here for posterity. CartesianCoordinate3D coord_1, coord_2; proj_data_info_sptr->find_cartesian_coordinates_given_scanner_coordinates (coord_1,coord_2, ring_a,ring_b, det_num_a, det_num_b); - const CartesianCoordinate3D - shift(scanner_sptr->get_ring_spacing()*(scanner_sptr->get_num_rings()-1)/2.F, 0.F,0.F); cout << "\nObsolete function:\n" << "find_cartesian_coordinates_given_scanner_coordinates (after shift in z)\n" - << coord_1 - shift << coord_2 - shift <<'\n'; + << coord_1 << coord_2 <<'\n'; // will give same result as above proj_data_info_sptr->find_bin_given_cartesian_coordinates_of_detection(bin, coord_1, coord_2); cout << "bin: (segment " << bin.segment_num() << ", axial pos " << bin.axial_pos_num()