rework functions for interpolating projdata for scatter, mostly for B…

…locksOnCylindrical - Fixed inverse SSRB by always linearly interpolating between the two closest slices. This makes it work for BlocksOnCylindrical as well - Rework extend_* functions: new extend_segment function that works on more cases, deprecated others This works with BlocksOnCylindrical and view-offset. Fixes #1177 - interpolate_projdata rewrite, including use linear interpolation in axial direction for BlocksOnCylindrical - Deprecate interpolate_axial_position which used nearest neighbour - Remove use_view_offset parameter from interpolate_projdata functions (always use view offset) - Fixed #1178 offset handling in sample_function_on_regular_grid, add sample_function_using_index_converter - add test function - Updated release notes. Fixes #1109 Fixes #1161
UCL · Apr 8, 2023 · 412d934 · 412d934
1 parent 7f6de79
commit 412d934
Show file tree

Hide file tree

Showing 14 changed files with 816 additions and 234 deletions.
diff --git a/documentation/release_5.2.htm b/documentation/release_5.2.htm
@@ -33,6 +33,11 @@ <h3>New functionality</h3>
 <ul>
   <li>Parallelised function <code>set_fan_data_add_gaps_help</code> across segments to reduce computation time.
   <br /><a href="https://github.com/UCL/STIR/pull/1168/">PR #1168</a>.
+  <li>Scatter estimation is now smoothed in axial direction for BlocksOnCylindrical scanners.
+    <br /><a href="https://github.com/UCL/STIR/pull/1172/">PR #1172</a>.
+  </li>
+  <li><code>InverseSSRB</code> now works for BlocksOnCylindrical after a rewrite.
+    <br /><a href="https://github.com/UCL/STIR/pull/1172/">PR #1172</a>.
   </li>
 </ul>
 
@@ -59,6 +64,9 @@ <h3>Changed functionality</h3>
   <li>Errors now throw <code>std::runtime_error</code> instead of <code>std::string</code>.
   <br /><a href="https://github.com/UCL/STIR/pull/1131/">PR #1131</a>.
   </li>
+  <li>The parameter <code>use_view_offset</code> was removed from the <code>interpolate_projdata</code> functions.
+  <br /><a href="https://github.com/UCL/STIR/pull/1172/">PR #1172</a>.
+  </li>
 </ul>
 
 <h3>Build system and dependencies</h3>
@@ -80,6 +88,12 @@ <h3>Minor (?) bug fixes</h3>
   <li>Small change in scatter simulation to how non-arccorrected bins are computed. Added a check in the construction of non-arccorrected projdata that the number of tangential bins is not larger than the maximum non-arccorrected bins.
     <br /><a href="https://github.com/UCL/STIR/pull/1152/">PR #1152</a>.
   </li>
+  <li><code>extend_segment_in_views</code> does not handle view offsets correctly and does not work for BlocksOnCylindrical scanners <br /><a href="https://github.com/UCL/STIR/issues/1177/">issue #1177</a>. A new function <code>extend_segment</code> was added that works for Cylindrical and BlocksOnCylindrical and allows extension in tangential and axial direction as well.
+    <br /><a href="https://github.com/UCL/STIR/pull/1172/">PR #1172</a>.
+  </li>
+  <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>
 </ul>
 
 <h3>Documentation changes</h3>

diff --git a/src/buildblock/extend_projdata.cxx b/src/buildblock/extend_projdata.cxx
@@ -2,6 +2,7 @@
 //
 /*
   Copyright (C) 2005- 2011, Hammersmith Imanet Ltd
+  Copyright 2023, Positrigo AG, Zurich
   This file is part of STIR.
 
   SPDX-License-Identifier: Apache-2.0
@@ -15,7 +16,7 @@
 
 \author Kris Thielemans
 \author Charalampos Tsoumpas
-  
+\author Markus Jehl
 */
 #include "stir/Array.h"
 #include "stir/SegmentBySinogram.h"
@@ -29,6 +30,7 @@ START_NAMESPACE_STIR
 
 namespace detail 
 {
+#if STIR_VERSION < 060000
   /* This function takes symmetries in the sinogram space into account
      to find data in the negative segment if necessary.
      However, it needs testing if it would work for non-direct sinograms.
@@ -144,6 +146,108 @@ extend_segment_in_views(const SegmentBySinogram<float>& sino,
     }
   return out;
 }
+#endif
+
+/* This function takes symmetries into account to extend segments in any direction.
+   However, it needs testing if it would work for non-direct sinograms.
+*/
+Array<3,float>
+extend_segment(const SegmentBySinogram<float>& segment, const int view_extension,
+               const int axial_extension, const int tangential_extension)
+{
+  Array<3,float> out(segment);
+  BasicCoordinate<3,int> min_dim, max_dim;
+  min_dim[1] = segment.get_min_index() - axial_extension;
+  min_dim[2] = segment[0].get_min_index() - view_extension;
+  min_dim[3] = segment[0][0].get_min_index() - tangential_extension;
+  max_dim[1] = segment.get_max_index() + axial_extension;
+  max_dim[2] = segment[0].get_max_index() + view_extension;
+  max_dim[3] = segment[0][0].get_max_index() + tangential_extension;
+  out.grow(IndexRange<3>(min_dim, max_dim));
+
+  // fill the axial extensions with the same entries from the border
+  for (int axial_edge = 0; axial_edge < axial_extension; axial_edge++)
+  {
+    out[min_dim[1] + axial_edge] = out[min_dim[1] + axial_extension];
+    out[max_dim[1] - axial_edge] = out[max_dim[1] - axial_extension];
+  }
+
+  // check, whether the projection data cover 180° or 360°
+  bool flip_views = false;
+  bool extend_without_wrapping = false;
+  float min_phi=_PI, max_phi=-_PI;
+  for (auto view = segment.get_proj_data_info_sptr()->get_min_view_num(); view <= segment.get_proj_data_info_sptr()->get_max_view_num(); view++)
+  {
+    auto phi = segment.get_proj_data_info_sptr()->get_phi(Bin(0, view, 0, 0));
+    if (phi < min_phi)
+      min_phi = phi;
+    if (phi > max_phi)
+      max_phi = phi;
+  }
+  const auto phi_range = max_phi - min_phi;
+  const auto average_phi_sampling = phi_range / (segment.get_proj_data_info_sptr()->get_num_views() - 1);
+  if (abs(phi_range - 2 * _PI) < 5 * average_phi_sampling)
+    flip_views = false; // if views cover 360°, we can simply wrap around
+  else if (abs(phi_range - _PI) < 5 * average_phi_sampling)
+    flip_views = true;  // if views cover 180°, the tangential positions need to be flipped
+  else
+  {
+    extend_without_wrapping = true;
+    warning("Extending ProjData by wrapping only works for view coverage of 180° or 360°. Instead, just extending with nearest neighbour.");
+  }
+
+  // fill the view extensions by wrapping around
+  for (int view_edge = 0; view_edge < view_extension; view_edge++)
+  {
+    for (int axial_pos = min_dim[1]; axial_pos <= max_dim[1]; axial_pos++)
+    {
+      if (extend_without_wrapping)
+      {
+        out[axial_pos][min_dim[2] + view_edge] = out[axial_pos][min_dim[2] + view_extension];
+        out[axial_pos][max_dim[2] - view_extension] = out[axial_pos][max_dim[2] - view_extension];
+      }
+      else if (flip_views)
+      {
+        const int sym_dim = std::min(abs(min_dim[3]), max_dim[3]);
+        for (int tang_pos = -sym_dim; tang_pos <= sym_dim; tang_pos++)
+        {
+          out[axial_pos][min_dim[2] + view_edge][tang_pos] = out[axial_pos][max_dim[2] - 2 * view_extension + view_edge + 1][-tang_pos];
+          out[axial_pos][max_dim[2] - view_extension + 1 + view_edge][tang_pos] = out[axial_pos][min_dim[2] + view_extension + view_edge][-tang_pos];
+        }
+        for (int tang_pos = min_dim[3]; tang_pos < -sym_dim; tang_pos++)
+        { // fill in asymmetric tangential positions at the end by just picking the nearest existing element
+          out[axial_pos][min_dim[2] + view_edge][tang_pos] = out[axial_pos][max_dim[2] - 2 * view_extension + view_edge + 1][sym_dim];
+          out[axial_pos][max_dim[2] - view_extension + 1 + view_edge][tang_pos] = out[axial_pos][min_dim[2] + view_extension + view_edge][sym_dim];
+        }
+        for (int tang_pos = max_dim[3]; tang_pos > sym_dim; tang_pos--)
+        { // fill in asymmetric tangential positions at the end by just picking the nearest existing element
+          out[axial_pos][min_dim[2] + view_edge][tang_pos] = out[axial_pos][max_dim[2] - 2 * view_extension + view_edge + 1][-sym_dim];
+          out[axial_pos][max_dim[2] - view_extension + 1 + view_edge][tang_pos] = out[axial_pos][min_dim[2] + view_extension + view_edge][-sym_dim];
+        }
+      }
+      else
+      {
+        out[axial_pos][min_dim[2] + view_edge] = out[axial_pos][max_dim[2] - 2 * view_extension + view_edge + 1];
+        out[axial_pos][max_dim[2] - view_extension + 1 + view_edge] = out[axial_pos][min_dim[2] + view_extension + view_edge];
+      }
+    }
+  }
+
+  // fill tangential extension with same entries than boundary
+  for (int tang_edge = 0; tang_edge < tangential_extension; tang_edge++)
+  {
+    for (int axial_pos = min_dim[1]; axial_pos <= max_dim[1]; axial_pos++)
+    {
+      for (int view = min_dim[2]; view <= max_dim[2]; view++)
+      {
+        out[axial_pos][view][min_dim[3] + tang_edge] = out[axial_pos][view][min_dim[3] + tangential_extension];
+        out[axial_pos][view][max_dim[3] - tang_edge] = out[axial_pos][view][max_dim[3] - tangential_extension];
+      }
+    }
+  }
+
+  return out;
+}
 
 Array<2,float>
 extend_sinogram_in_views(const Sinogram<float>& sino,

diff --git a/src/buildblock/interpolate_axial_position.cxx b/src/buildblock/interpolate_axial_position.cxx
@@ -40,6 +40,7 @@
 
 START_NAMESPACE_STIR
 
+#if STIR_VERSION < 060000
 Succeeded 
 interpolate_axial_position(ProjData& proj_data_out,
                      const ProjData& proj_data_in)
@@ -122,5 +123,6 @@ interpolate_axial_position(ProjData& proj_data_out,
       }
   return Succeeded::yes;
 }
+#endif
 
 END_NAMESPACE_STIR