From 729dee235cdbc4dc8fed67efa64d9443aab712f6 Mon Sep 17 00:00:00 2001 From: Markus Jehl Date: Mon, 7 Aug 2023 13:50:24 +0000 Subject: [PATCH 1/3] Added functionality to invert theta value when obtained from LOR that is being flipped when obtaining phi and s. --- src/include/stir/LORCoordinates.inl | 15 +++++++- src/include/stir/ProjDataInfoGeneric.inl | 16 +++++++- .../stir/ProjDataInfoGenericNoArcCorr.inl | 38 +------------------ src/test/test_proj_data_info.cxx | 20 ++++++---- 4 files changed, 42 insertions(+), 47 deletions(-) diff --git a/src/include/stir/LORCoordinates.inl b/src/include/stir/LORCoordinates.inl index f7754c3981..e9b137c416 100644 --- a/src/include/stir/LORCoordinates.inl +++ b/src/include/stir/LORCoordinates.inl @@ -346,8 +346,19 @@ find_LOR_intersections_with_cylinder(LORAs2Points& intersection_coords, auto l2 = static_cast((-b - root) / a); // TODO won't work when coordT1!=coordT2 - intersection_coords.p1() = d*l1 + c1; - intersection_coords.p2() = d*l2 + c1; + const CartesianCoordinate3D distance1 = d * l1 + c1 - coords.p1(); + const CartesianCoordinate3D distance2 = d * l2 + c1 - coords.p1(); + if (square(distance1.x()) + square(distance1.y()) + square(distance1.z()) + < square(distance2.x()) + square(distance2.y()) + square(distance2.z())) + { + intersection_coords.p1() = d * l1 + c1; + intersection_coords.p2() = d * l2 + c1; + } + else + { + intersection_coords.p2() = d * l1 + c1; + intersection_coords.p1() = d * l2 + c1; + } assert(fabs(square(intersection_coords.p1().x())+square(intersection_coords.p1().y()) - square(radius)) _p1; CartesianCoordinate3D _p2; find_cartesian_coordinates_of_detection(_p1, _p2, bin); + + _p1.z() += z_shift.z(); + _p2.z() += z_shift.z(); + + // do the below transformation to check whether the points are swapped when phi is obtained + // if they are swapped theta needs to be inverted, otherwise the oblique segments contain flipped streaks + LORAs2Points lor_as_2_points(_p1, _p2); + const double R = sqrt(fmax(square(_p1.x()) + square(_p1.y()), square(_p2.x()) + square(_p2.y()))); + LORInAxialAndNoArcCorrSinogramCoordinates lor; + lor_as_2_points.change_representation(lor, R); + float sign = 1.0; + if (abs(_p1.z() - lor.z1()) + abs(_p2.z() - lor.z2()) > abs(_p1.z() - lor.z2()) + abs(_p2.z() - lor.z1())) + sign *= -1; + CartesianCoordinate3D p2_minus_p1 = _p2 - _p1; - return p2_minus_p1.z() / (sqrt(square(p2_minus_p1.x())+square(p2_minus_p1.y()))); + return sign * p2_minus_p1.z() / (sqrt(square(p2_minus_p1.x()) + square(p2_minus_p1.y()))); } float diff --git a/src/include/stir/ProjDataInfoGenericNoArcCorr.inl b/src/include/stir/ProjDataInfoGenericNoArcCorr.inl index 2a796e42a4..01c0ac5c88 100644 --- a/src/include/stir/ProjDataInfoGenericNoArcCorr.inl +++ b/src/include/stir/ProjDataInfoGenericNoArcCorr.inl @@ -88,43 +88,7 @@ get_s(const Bin& bin) const { LORInAxialAndNoArcCorrSinogramCoordinates lor; get_LOR(lor, bin); - // REQUIREMENT: Euclidean coordinate of 3 points, a,b and c. - // CALCULATION: // Equation of a line, in parametric form, given two point a, b: p(t) = a + (b-a)*t - // // Let a,b,c be points in 2D and let a,b form a line then shortest distance between c and line ab is: - // // || (p0-p1) - [(p0-p1)*(p2-p1)/||p2-p1||^2]*(p2-p1) || - // p1 = _p1, p2=_p2 and p0=(0,0,0) - // OUTPUT: replace the ring_radius*sin(lor.beta()) with distance from o to ab. - - // Can't access coordinate of detection from this class so we have to recalculate where it is: - - CartesianCoordinate3D _p1; - CartesianCoordinate3D _p2; - find_cartesian_coordinates_of_detection(_p1, _p2, bin); - - // // Get p1-p0 and p2-p1 vector. - CartesianCoordinate3D p2_minus_p1 = _p2 - _p1; - CartesianCoordinate3D p0_minus_p1 = CartesianCoordinate3D(0,0,0) - _p1; - float p0_minus_p1_dot_p2_minus_p1 = p0_minus_p1.x()* p2_minus_p1.x() + p0_minus_p1.y()*p2_minus_p1.y(); - float p2_minus_p1_magitude = square(p2_minus_p1.x()) + square(p2_minus_p1.y()); - float x = 0; - float y = 0; - float sign = sin(lor.beta()); - if (p2_minus_p1_magitude > 0.01) - { - x = p0_minus_p1.x() - (p0_minus_p1_dot_p2_minus_p1/p2_minus_p1_magitude)*p2_minus_p1.x(); - y = p0_minus_p1.y() - (p0_minus_p1_dot_p2_minus_p1/p2_minus_p1_magitude)*p2_minus_p1.y(); - } - else - { - error("get_s(): 2 detection points are too close to each other. This indicates an internal error."); - } - float s = sqrt(square(x) + square(y)); - - if(sign < 0.0){ - return -s; - }else{ - return s; - } + return lor.s(); } void diff --git a/src/test/test_proj_data_info.cxx b/src/test/test_proj_data_info.cxx index 0247b1c32a..242d8727ee 100644 --- a/src/test/test_proj_data_info.cxx +++ b/src/test/test_proj_data_info.cxx @@ -148,7 +148,7 @@ ProjDataInfoTests::test_generic_proj_data_info(ProjDataInfo& proj_data_info) #endif for (int segment_num = proj_data_info.get_min_segment_num(); segment_num <= proj_data_info.get_max_segment_num(); ++segment_num) { - for (int view_num = proj_data_info.get_min_view_num(); view_num <= proj_data_info.get_max_view_num(); view_num += 3) + for (int view_num = proj_data_info.get_min_view_num(); view_num <= proj_data_info.get_max_view_num(); view_num += 1) { // loop over axial_positions. Avoid using first and last positions, as // if there is axial compression, the central LOR of a bin might actually not @@ -176,10 +176,13 @@ ProjDataInfoTests::test_generic_proj_data_info(ProjDataInfo& proj_data_info) { const Bin new_bin = proj_data_info.get_bin(lor); #if 1 - const int diff_segment_num = intabs(org_bin.segment_num() - new_bin.segment_num()); - const int diff_view_num = intabs(org_bin.view_num() - new_bin.view_num()); + // the differences need to also consider wrap-around in views, which would flip tangential pos and segment + const int diff_segment_num = intabs(org_bin.view_num() - new_bin.view_num()) < proj_data_info.get_num_views() - intabs(org_bin.view_num() - new_bin.view_num()) ? + intabs(org_bin.segment_num() - new_bin.segment_num()) : intabs(org_bin.segment_num() + new_bin.segment_num()); + const int diff_view_num = min(intabs(org_bin.view_num() - new_bin.view_num()), proj_data_info.get_num_views() - intabs(org_bin.view_num() - new_bin.view_num())); const int diff_axial_pos_num = intabs(org_bin.axial_pos_num() - new_bin.axial_pos_num()); - const int diff_tangential_pos_num = intabs(org_bin.tangential_pos_num() - new_bin.tangential_pos_num()); + const int diff_tangential_pos_num = intabs(org_bin.view_num() - new_bin.view_num()) < proj_data_info.get_num_views() - intabs(org_bin.view_num() - new_bin.view_num()) ? + intabs(org_bin.tangential_pos_num() - new_bin.tangential_pos_num()) : intabs(org_bin.tangential_pos_num() + new_bin.tangential_pos_num()); if (new_bin.get_bin_value() > 0) { if (diff_segment_num > max_diff_segment_num) @@ -220,10 +223,13 @@ ProjDataInfoTests::test_generic_proj_data_info(ProjDataInfo& proj_data_info) lor.get_intersections_with_cylinder(lor_as_points, lor.radius()); const Bin new_bin = proj_data_info.get_bin(lor_as_points); #if 1 - const int diff_segment_num = intabs(org_bin.segment_num() - new_bin.segment_num()); - const int diff_view_num = intabs(org_bin.view_num() - new_bin.view_num()); + // the differences need to also consider wrap-around in views, which would flip tangential pos and segment + const int diff_segment_num = intabs(org_bin.view_num() - new_bin.view_num()) < proj_data_info.get_num_views() - intabs(org_bin.view_num() - new_bin.view_num()) ? + intabs(org_bin.segment_num() - new_bin.segment_num()) : intabs(org_bin.segment_num() + new_bin.segment_num()); + const int diff_view_num = min(intabs(org_bin.view_num() - new_bin.view_num()), proj_data_info.get_num_views() - intabs(org_bin.view_num() - new_bin.view_num())); const int diff_axial_pos_num = intabs(org_bin.axial_pos_num() - new_bin.axial_pos_num()); - const int diff_tangential_pos_num = intabs(org_bin.tangential_pos_num() - new_bin.tangential_pos_num()); + const int diff_tangential_pos_num = intabs(org_bin.view_num() - new_bin.view_num()) < proj_data_info.get_num_views() - intabs(org_bin.view_num() - new_bin.view_num()) ? + intabs(org_bin.tangential_pos_num() - new_bin.tangential_pos_num()) : intabs(org_bin.tangential_pos_num() + new_bin.tangential_pos_num()); if (new_bin.get_bin_value() > 0) { if (diff_segment_num > max_diff_segment_num) From 3ebd2739a95da4e88f8b8a9e3aff0bbfaa926b10 Mon Sep 17 00:00:00 2001 From: Markus Jehl Date: Thu, 17 Aug 2023 14:08:01 +0000 Subject: [PATCH 2/3] Simplifications of the changes. --- src/include/stir/LORCoordinates.inl | 20 ++++---------------- src/include/stir/ProjDataInfoGeneric.inl | 19 ++----------------- 2 files changed, 6 insertions(+), 33 deletions(-) diff --git a/src/include/stir/LORCoordinates.inl b/src/include/stir/LORCoordinates.inl index e9b137c416..e17685a0a0 100644 --- a/src/include/stir/LORCoordinates.inl +++ b/src/include/stir/LORCoordinates.inl @@ -342,23 +342,11 @@ find_LOR_intersections_with_cylinder(LORAs2Points& intersection_coords, } const coordT2 root = static_cast(sqrt(argsqrt)); - auto l1 = static_cast((-b + root) / a); - auto l2 = static_cast((-b - root) / a); + auto l1 = static_cast((-b - root) / a); + auto l2 = static_cast((-b + root) / a); - // TODO won't work when coordT1!=coordT2 - const CartesianCoordinate3D distance1 = d * l1 + c1 - coords.p1(); - const CartesianCoordinate3D distance2 = d * l2 + c1 - coords.p1(); - if (square(distance1.x()) + square(distance1.y()) + square(distance1.z()) - < square(distance2.x()) + square(distance2.y()) + square(distance2.z())) - { - intersection_coords.p1() = d * l1 + c1; - intersection_coords.p2() = d * l2 + c1; - } - else - { - intersection_coords.p2() = d * l1 + c1; - intersection_coords.p1() = d * l2 + c1; - } + intersection_coords.p1() = d * l1 + c1; + intersection_coords.p2() = d * l2 + c1; assert(fabs(square(intersection_coords.p1().x())+square(intersection_coords.p1().y()) - square(radius)) _p1; - CartesianCoordinate3D _p2; - find_cartesian_coordinates_of_detection(_p1, _p2, bin); - - _p1.z() += z_shift.z(); - _p2.z() += z_shift.z(); - - // do the below transformation to check whether the points are swapped when phi is obtained - // if they are swapped theta needs to be inverted, otherwise the oblique segments contain flipped streaks - LORAs2Points lor_as_2_points(_p1, _p2); - const double R = sqrt(fmax(square(_p1.x()) + square(_p1.y()), square(_p2.x()) + square(_p2.y()))); LORInAxialAndNoArcCorrSinogramCoordinates lor; - lor_as_2_points.change_representation(lor, R); - float sign = 1.0; - if (abs(_p1.z() - lor.z1()) + abs(_p2.z() - lor.z2()) > abs(_p1.z() - lor.z2()) + abs(_p2.z() - lor.z1())) - sign *= -1; + get_LOR(lor, bin); - CartesianCoordinate3D p2_minus_p1 = _p2 - _p1; - return sign * p2_minus_p1.z() / (sqrt(square(p2_minus_p1.x()) + square(p2_minus_p1.y()))); + return (lor.z2() - lor.z1()) / (2 * lor.radius()); } float From 4aa69f26587088d4f2da810c2f574108c1043d40 Mon Sep 17 00:00:00 2001 From: Markus Jehl Date: Thu, 17 Aug 2023 15:23:42 +0000 Subject: [PATCH 3/3] Updated release notes. --- documentation/release_5.2.htm | 3 +++ 1 file changed, 3 insertions(+) diff --git a/documentation/release_5.2.htm b/documentation/release_5.2.htm index 87d9a25ba7..b115b24639 100644 --- a/documentation/release_5.2.htm +++ b/documentation/release_5.2.htm @@ -125,6 +125,9 @@

Minor (?) bug fixes

  • sample_function_on_regular_grid did not handle offset correctly in all places
    issue #1178.
    PR #1172.
    • +
  • Ray tracing projection for BlocksOnCylindrical scanner geometries contained a bug where some bins were swapped across oblique segments
    issue #1223. This sometimes lead to large artifacts in reconstructions. +
    PR #1231. +

    • Documentation changes