diff --git a/documentation/release_5.2.htm b/documentation/release_5.2.htm index a11ca8f5d6..e3d005fb89 100644 --- a/documentation/release_5.2.htm +++ b/documentation/release_5.2.htm @@ -136,6 +136,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

    diff --git a/src/include/stir/LORCoordinates.inl b/src/include/stir/LORCoordinates.inl index f7754c3981..e17685a0a0 100644 --- a/src/include/stir/LORCoordinates.inl +++ b/src/include/stir/LORCoordinates.inl @@ -342,12 +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 - intersection_coords.p1() = d*l1 + c1; - intersection_coords.p2() = 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); - CartesianCoordinate3D p2_minus_p1 = _p2 - _p1; - return p2_minus_p1.z() / (sqrt(square(p2_minus_p1.x())+square(p2_minus_p1.y()))); + LORInAxialAndNoArcCorrSinogramCoordinates lor; + get_LOR(lor, bin); + + return (lor.z2() - lor.z1()) / (2 * lor.radius()); } 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)