Fix and simplify get_tantheta etc for ProjDataInfoGeneric (#1231)

- Fix get_tantheta to correctly invert theta value when obtained from LOR that is being flipped when obtaining phi and s by using get_LOR
- Simplifications of get_s()
- make sure that find_LOR_intersections_with_cylinder(LORAs2Points&, ...) returns points in some order as input (they were swapped)
- expand test_projdata_info to run more views
- Updated release notes.

Fixes #1223
---------

Co-authored-by: Markus Jehl <markus.jehl@positrigo.com>

documentation/release_5.2.htm

@@ -136,6 +136,9 @@ <h3>Minor (?) bug fixes</h3>
   <li><code>sample_function_on_regular_grid</code> did not handle offset correctly in all places <br /><a href="https://github.com/UCL/STIR/issues/1178/">issue #1178</a>.
     <br /><a href="https://github.com/UCL/STIR/pull/1172/">PR #1172</a>.
   </li>
+  <li>Ray tracing projection for BlocksOnCylindrical scanner geometries contained a bug where some bins were swapped across oblique segments <br /><a href="https://github.com/UCL/STIR/issues/1223/">issue #1223</a>. This sometimes lead to large artifacts in reconstructions.
+    <br /><a href="https://github.com/UCL/STIR/pull/1231/">PR #1231</a>.
+  </li>
 </ul>
 
 <h3>Documentation changes</h3>

src/include/stir/LORCoordinates.inl

@@ -342,12 +342,11 @@ find_LOR_intersections_with_cylinder(LORAs2Points<coordT1>& intersection_coords,
   }
   const coordT2 root = static_cast<coordT2>(sqrt(argsqrt));
 
-  auto l1 = static_cast<coordT2>((-b + root) / a);
-  auto l2 = static_cast<coordT2>((-b - root) / a);
+  auto l1 = static_cast<coordT2>((-b - root) / a);
+  auto l2 = static_cast<coordT2>((-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))
 	 <square(radius)*10.E-5);

src/include/stir/ProjDataInfoGeneric.inl

@@ -80,11 +80,10 @@ ProjDataInfoGeneric::get_t(const Bin& bin) const
 float
 ProjDataInfoGeneric::get_tantheta(const Bin& bin) const
 {
-  CartesianCoordinate3D<float> _p1;
-  CartesianCoordinate3D<float> _p2;
-  find_cartesian_coordinates_of_detection(_p1, _p2, bin);
-  CartesianCoordinate3D<float> p2_minus_p1 = _p2 - _p1;
-  return p2_minus_p1.z() / (sqrt(square(p2_minus_p1.x())+square(p2_minus_p1.y()))); 
+  LORInAxialAndNoArcCorrSinogramCoordinates<float> lor;
+	get_LOR(lor, bin);
+
+  return (lor.z2() - lor.z1()) / (2 * lor.radius());
 }
 
 float

src/include/stir/ProjDataInfoGenericNoArcCorr.inl

@@ -88,43 +88,7 @@ get_s(const Bin& bin) const
 {
   LORInAxialAndNoArcCorrSinogramCoordinates<float> 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<float> _p1;
-	CartesianCoordinate3D<float> _p2;
-	find_cartesian_coordinates_of_detection(_p1, _p2, bin);
-
-  // // Get p1-p0 and p2-p1 vector.
-  CartesianCoordinate3D<float> p2_minus_p1 = _p2 - _p1;
-  CartesianCoordinate3D<float> p0_minus_p1 = CartesianCoordinate3D<float>(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

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)