add make_detection_bin helpers

inverse of expand_detection_bin

Author: KrisThielemans

cpp/petsird_analysis.cpp

@@ -157,24 +157,29 @@ main(int argc, char const* argv[])
           // Note: just doing one module-type ATM
           for (auto& event : event_time_block.prompt_events[type_of_module_pair[0]][type_of_module_pair[1]])
             {
-              const auto expanded_det_bin0
+              const auto expanded_detection_bin0
                   = petsird_helpers::expand_detection_bin(header.scanner, type_of_module_pair[0], event.detection_bins[0]);
-              const auto expanded_det_bin1
+              const auto expanded_detection_bin1
                   = petsird_helpers::expand_detection_bin(header.scanner, type_of_module_pair[1], event.detection_bins[1]);
-              energy_1 += energy_mid_points[expanded_det_bin0.energy_index];
-              energy_2 += energy_mid_points[expanded_det_bin1.energy_index];
+
+              // TODO move this test to separate unit-tests
+              assert(event.detection_bins[0]
+                     == petsird_helpers::make_detection_bin(header.scanner, type_of_module_pair[0], expanded_detection_bin0));
+
+              energy_1 += energy_mid_points[expanded_detection_bin0.energy_index];
+              energy_2 += energy_mid_points[expanded_detection_bin1.energy_index];
 
               if (print_events)
                 {
                   std::cout << "CoincidenceEvent(detectionBins=[" << event.detection_bins[0] << ", " << event.detection_bins[1]
                             << "], tofIdx=" << event.tof_idx << "])\n";
                   std::cout << "    "
-                            << "[ExpandedDetectionBin(module=" << expanded_det_bin0.module_index << ", "
-                            << "el=" << expanded_det_bin0.element_index << ", "
-                            << "energy_index=" << expanded_det_bin0.energy_index
-                            << "), ExpandedDetectionBin(module=" << expanded_det_bin1.module_index << ", "
-                            << "el=" << expanded_det_bin1.element_index << ", "
-                            << "energy_index=" << expanded_det_bin0.energy_index << ")]\n";
+                            << "[ExpandedDetectionBin(module=" << expanded_detection_bin0.module_index << ", "
+                            << "el=" << expanded_detection_bin0.element_index << ", "
+                            << "energy_index=" << expanded_detection_bin0.energy_index
+                            << "), ExpandedDetectionBin(module=" << expanded_detection_bin1.module_index << ", "
+                            << "el=" << expanded_detection_bin1.element_index << ", "
+                            << "energy_index=" << expanded_detection_bin0.energy_index << ")]\n";
                   std::cout << "    efficiency:"
                             << petsird_helpers::get_detection_efficiency(header.scanner, type_of_module_pair, event) << "\n";
                 }

cpp/petsird_helpers.h

@@ -54,6 +54,7 @@ get_num_detection_bins(const ScannerInformation& scanner, const TypeOfModule& ty
   return det_els.transforms.size() * rep_module.transforms.size() * energy_bin_edges.NumberOfBins();
 }
 
+//! Create a vector of ExpandedDetectionBins from a list/vector of DetectionBins
 template <class T>
 inline std::vector<ExpandedDetectionBin>
 expand_detection_bins(const ScannerInformation& scanner, const TypeOfModule& type_of_module, const T& list_of_detection_bins)
@@ -74,6 +75,7 @@ expand_detection_bins(const ScannerInformation& scanner, const TypeOfModule& typ
   return result;
 }
 
+//! Expand a DetectionBin into a ExpandedDetectionBin
 inline ExpandedDetectionBin
 expand_detection_bin(const ScannerInformation& scanner, const TypeOfModule& type_of_module, const DetectionBin& detection_bin)
 {
@@ -83,6 +85,41 @@ expand_detection_bin(const ScannerInformation& scanner, const TypeOfModule& type
   return expanded_bins[0];
 }
 
+//! Create a vector of DetectionBins from a list/vector of ExpandedDetectionBins
+template <class T>
+inline std::vector<DetectionBin>
+make_detection_bins(const ScannerInformation& scanner, const TypeOfModule& type_of_module,
+                    const T& list_of_expanded_detection_bins)
+{
+  assert(type_of_module < scanner.scanner_geometry.replicated_modules.size());
+  const auto& rep_module = scanner.scanner_geometry.replicated_modules[type_of_module];
+  const auto& energy_bin_edges = scanner.event_energy_bin_edges[type_of_module];
+  const auto num_en = energy_bin_edges.NumberOfBins();
+  const auto num_el_per_module = rep_module.object.detecting_elements.transforms.size();
+
+  std::vector<DetectionBin> result;
+  // TODO call reserve()
+  for (auto bin : list_of_expanded_detection_bins)
+    {
+      // need to do static_cast to avoid compiler warning on narrowing from std::size_t to uint32_t
+      const auto detection_bin
+          = static_cast<DetectionBin>(bin.energy_index + (bin.element_index + bin.module_index * num_el_per_module) * num_en);
+      result.push_back({ detection_bin });
+    }
+  return result;
+}
+
+//! Create a DetectionBin from a ExpandedDetectionBin
+inline DetectionBin
+make_detection_bin(const ScannerInformation& scanner, const TypeOfModule& type_of_module,
+                   const ExpandedDetectionBin& expanded_detection_bin)
+{
+  // TODO very inefficient, but avoid re-implementation of code above.
+  std::vector<ExpandedDetectionBin> expanded_detection_bins{ expanded_detection_bin };
+  const auto detection_bins = make_detection_bins(scanner, type_of_module, expanded_detection_bins);
+  return detection_bins[0];
+}
+
 inline float
 get_detection_efficiency(const ScannerInformation& scanner, const TypeOfModulePair& type_of_module_pair,
                          const CoincidenceEvent& event)

matlab/toolbox/+petsird/+helpers/make_detection_bin.m

@@ -0,0 +1,18 @@
+function detection_bin = make_detection_bin(scanner_geometry, type_of_module, detection_bin)
+% Find DetectionBin for a ExpandedDetectionBin
+
+arguments
+    scanner_geometry (1,1) petsird.ScannerGeometry
+    type_of_module (1,1) petsird.TypeOfModule
+    expanded_detection_bin (1,1) petsird.ExpandedDetectionBin
+end
+
+rep_module = scanner_geometry.replicated_modules(type_of_module);
+num_el_per_module = length(rep_module.object.detecting_elements.transforms);
+energy_bin_edges = rep_module.event_energy_bin_edges(type_of_module);
+num_en = energy_bin_edges.number_of_bins();
+% use shorter name
+bin = expanded_detection_bin;
+detection_bin = petsird.DetectionBin(bin.energy_index + (bin.element_index + bin.module_index * num_el_per_module) * num_en);
+
+end

matlab/toolbox/+petsird/+helpers/make_detection_bins.m

@@ -0,0 +1,13 @@
+function detection_bins = make_detection_bins(scanner, type_of_module, expanded_detection_bins)
+% Find DetectionBin for a list of ExpandedDetectionBins
+
+arguments
+    scanner_geometry (1,1) petsird.ScannerInformation
+    type_of_module (1,1) petsird.TypeOfModule
+    expanded_detection_bins (:,1) petsird.ExpandedDetectionBin
+end
+
+% TODO slow but avoids re-implementation
+detection_bins = arrayfun(@(bin) make_detection_bin(scanner, type_of_module, bin)), expanded_detection_bins);
+
+end

python/petsird/helpers/__init__.py

@@ -56,6 +56,33 @@ def expand_detection_bin(
     return expand_detection_bins(scanner, type_of_module, [detection_bin])[0]
 
 
+def make_detection_bins(
+    scanner: petsird.ScannerInformation, type_of_module: petsird.TypeOfModule,
+    expanded_detection_bins: typing.Iterable[petsird.ExpandedDetectionBin]
+) -> list[petsird.ExpandedDetectionBin]:
+    """Find DetectionBin for a list of expanded_detection_bins"""
+    rep_module = scanner.scanner_geometry.replicated_modules[type_of_module]
+    num_el_per_module = len(rep_module.object.detecting_elements.transforms)
+    energy_bin_edges = scanner.event_energy_bin_edges[type_of_module]
+    num_en = energy_bin_edges.number_of_bins()
+
+    return [
+        (bin.energy_index +
+         (bin.element_index + bin.module_index * num_el_per_module) * num_en)
+        for bin in expanded_detection_bins
+    ]
+
+
+def make_detection_bin(
+    scanner: petsird.ScannerInformation, type_of_module: petsird.TypeOfModule,
+    expanded_detection_bin: petsird.DetectionBin
+) -> petsird.ExpandedDetectionBin:
+    """Find DetectionBin for an expanded_detection_bin"""
+    # TODO probably slow implementation, but avoids re-implementation for now
+    return make_detection_bins(scanner, type_of_module,
+                               [expanded_detection_bin])[0]
+
+
 def get_detection_efficiency(scanner: petsird.ScannerInformation,
                              type_of_module_pair: petsird.TypeOfModulePair,
                              event: petsird.CoincidenceEvent) -> float:

python/petsird/helpers/analysis.py

@@ -8,7 +8,7 @@
 
 import petsird
 from petsird.helpers import (expand_detection_bin, get_detection_efficiency,
-                             get_num_det_els)
+                             get_num_det_els, make_detection_bin)
 
 
 def parserCreator():
@@ -103,6 +103,11 @@ def parserCreator():
                         scanner, type_of_module_pair[1],
                         event.detection_bins[1])
 
+                    # TODO move this test to separate unit-tests
+                    assert event.detection_bins[0] == make_detection_bin(
+                        scanner, type_of_module_pair[0],
+                        expanded_detection_bin0)
+
                     energy_1 += energy_mid_points[
                         expanded_detection_bin0.energy_index]
                     energy_2 += energy_mid_points[