Merge pull request UCL#1674 from KrisThielemans/TOF_improvements

TOF improvements: speed-up of Matrix projector by removing zeroes. Added extra tests.

Author: KrisThielemans

documentation/release_6.4.htm

@@ -51,11 +51,15 @@ <h4>Utilities</h4>
 
   <h3>Changed functionality</h3>
   <ul>
+    <li>When using a (ray tracing) matrix as projector, previously the matrix stored elements for the whole LOR,
+      while many elements would be zero. We now only store the non-zero ones, resulting in a speed-up for the Siemens Vision Quadra of about a factor 1.7. Memory usage is also reduced.<br>
+      <a href=https://github.com/UCL/STIR/pull/1674># 1674</a>.
+    </li>
     <li>
       OpenMP parallelisation was added to the function to create a 2d histogram of hits per crystal 
       from a proj data object, <code>make_fan_sum_data</code>.
       <br>
-      <a href=https://github.com/UCL/STIR/pull/1667>PR #1667</a>
+      <a href=https://github.com/UCL/STIR/pull/1667>PR #1667</a>.
     </li>
   </ul>
   <h4>Changes to examples</h4>
@@ -64,6 +68,11 @@ <h4>Changes to examples</h4>
 
   <h3>Bug fixes</h3>
   <ul>
+    <li>
+      An important bug fix for the "ray tracing matrix" (or actually any "matrix") projector for TOF PET. Previously, oblique
+      segments had problems, see <a href=https://github.com/UCL/STIR/issues/1537># 1537</a>.<br>
+      <a href=https://github.com/UCL/STIR/pull/1675># 1675</a>. Extra tests were introduced in <a href=https://github.com/UCL/STIR/pull/1674># 1674</a>.
+    </li>
   </ul>
 
   <h3>Build system</h3>
@@ -93,9 +102,7 @@ <h3>Build system</h3>
 </ul>
 
   <h3>Known problems</h3>
-  <p>See <a href=https://github.com/UCL/STIR/labels/bug>our issue tracker</a>.<br>
-    An important bug is for the "ray tracing matrix" projector for TOF PET. Oblique
-    segments currently have problems, see <a href=https://github.com/UCL/STIR/issues/1537># 1537</a>.
+  <p>See <a href=https://github.com/UCL/STIR/labels/bug>our issue tracker</a>.
   </p>
 
 

src/include/stir/recon_buildblock/ProjMatrixByBin.inl

@@ -133,6 +133,21 @@ ProjMatrixByBin::apply_tof_kernel(ProjMatrixElemsForOneBin& probabilities) const
   const CartesianCoordinate3D<float> diff = point2 - middle;
   const CartesianCoordinate3D<float> diff_unit_vector(diff / static_cast<float>(norm(diff)));
 
+  ProjMatrixElemsForOneBin tof_row(probabilities.get_bin());
+  // Estimate size of TOF row such that we can pre-allocate.
+  std::size_t max_num_elements;
+  {
+    const auto length = norm(point1 - point2);
+    const auto kernel_width = 8 / r_sqrt2_gauss_sigma;
+    const auto tof_bin_width = proj_data_info_sptr->tof_bin_boundaries_mm[probabilities.get_bin().timing_pos_num()].high_lim
+                               - proj_data_info_sptr->tof_bin_boundaries_mm[probabilities.get_bin().timing_pos_num()].low_lim;
+    const auto fraction = (kernel_width + tof_bin_width) / length;
+    // This seems to sometimes over-, sometimes underestimate, but it's probably not very important
+    // as std::vector will grow as necessary.
+    max_num_elements = std::size_t(probabilities.size() * std::min(fraction * 1.2, 1.001));
+  }
+  tof_row.reserve(max_num_elements);
+
   for (ProjMatrixElemsForOneBin::iterator element_ptr = probabilities.begin(); element_ptr != probabilities.end(); ++element_ptr)
     {
       Coordinate3D<int> c(element_ptr->get_coords());
@@ -150,8 +165,23 @@ ProjMatrixByBin::apply_tof_kernel(ProjMatrixElemsForOneBin& probabilities) const
       const float high_dist
           = ((proj_data_info_sptr->tof_bin_boundaries_mm[probabilities.get_bin().timing_pos_num()].high_lim - d2));
 
-      *element_ptr = ProjMatrixElemsForOneBin::value_type(c, element_ptr->get_value() * get_tof_value(low_dist, high_dist));
+      const auto tof_kernel_value = get_tof_value(low_dist, high_dist);
+      if (tof_kernel_value > 0)
+        {
+          if (auto non_tof_value = element_ptr->get_value())
+            tof_row.push_back(ProjMatrixElemsForOneBin::value_type(c, non_tof_value * tof_kernel_value));
+        }
+      else
+        {
+          // Optimisation would be to get out of the loop once we're "beyond" the TOF kernel,
+          // but it is tricky to do. It requires that the input is sorted
+          // "along" the LOR, i.e. d2 increases montonically, but that doesn't seem to be true.
+          // if (tof_row.size() > 0)
+          //   break;
+        }
     }
+  probabilities = tof_row;
+  // info("Estimate " + std::to_string(max_num_elements) + ", actual " + std::to_string(tof_row.size()));
 }
 
 float

src/include/stir/recon_buildblock/ProjMatrixElemsForOneBin.h

@@ -85,6 +85,8 @@ class ProjMatrixElemsForOneBin
   //! typedefs for iterator support
   typedef Element_vector::iterator iterator;
   typedef Element_vector::const_iterator const_iterator;
+  typedef Element_vector::reverse_iterator reverse_iterator;
+  typedef Element_vector::const_reverse_iterator const_reverse_iterator;
   typedef Element_vector::size_type size_type;
   typedef Element_vector::difference_type difference_type;
   typedef std::random_access_iterator_tag iterator_category;
@@ -117,6 +119,10 @@ class ProjMatrixElemsForOneBin
   inline const_iterator begin() const;
   inline iterator end();
   inline const_iterator end() const;
+  inline reverse_iterator rbegin();
+  inline const_reverse_iterator rbegin() const;
+  inline reverse_iterator rend();
+  inline const_reverse_iterator rend() const;
 
   //! reset lor to 0 length
   void erase();

src/include/stir/recon_buildblock/ProjMatrixElemsForOneBin.inl

@@ -73,6 +73,30 @@ ProjMatrixElemsForOneBin::end() const
   return elements.end();
 };
 
+ProjMatrixElemsForOneBin::reverse_iterator
+ProjMatrixElemsForOneBin::rbegin()
+{
+  return elements.rbegin();
+}
+
+ProjMatrixElemsForOneBin::const_reverse_iterator
+ProjMatrixElemsForOneBin::rbegin() const
+{
+  return elements.rbegin();
+};
+
+ProjMatrixElemsForOneBin::reverse_iterator
+ProjMatrixElemsForOneBin::rend()
+{
+  return elements.rend();
+};
+
+ProjMatrixElemsForOneBin::const_reverse_iterator
+ProjMatrixElemsForOneBin::rend() const
+{
+  return elements.rend();
+};
+
 ProjMatrixElemsForOneBin::iterator
 ProjMatrixElemsForOneBin::erase(iterator it)
 {