Update Python example with DetectorEfficiencies

python/petsird_analysis.py

@@ -6,7 +6,11 @@
 import sys
 import numpy
 import petsird
-
+from petsird_helpers import (
+    get_num_det_els,
+    get_module_and_element,
+    get_detection_efficiency,
+)
 
 if __name__ == "__main__":
     with petsird.BinaryPETSIRDReader(sys.stdin.buffer) as reader:
@@ -25,6 +29,9 @@
         print(
             f"Number of elements of first type in modules of first type: {len(header.scanner.scanner_geometry.replicated_modules[0].object.detecting_elements[0].transforms)}"
         )
+        print(
+            f"Total number of 'crystals' {get_num_det_els(header.scanner.scanner_geometry)}"
+        )
         print(f"Number of TOF bins: {header.scanner.number_of_tof_bins()}")
         print(f"Number of energy bins: {header.scanner.number_of_energy_bins()}")
 
@@ -43,6 +50,17 @@
                 energy_1 += energy_mid_points[event.energy_indices[0]]
                 energy_2 += energy_mid_points[event.energy_indices[1]]
 
+                print(event)
+                print(
+                    "   ",
+                    get_module_and_element(
+                        header.scanner.scanner_geometry, event.detector_ids
+                    ),
+                )
+                print(
+                    "    efficiency:", get_detection_efficiency(header.scanner, event)
+                )
+
         print(f"Last time block at {last_time} ms")
         print(f"Number of prompt events: {num_prompts}")
         print(f"Average energy_1: {energy_1 / num_prompts}")

python/petsird_generator.py

@@ -10,6 +10,7 @@
 import numpy.typing as npt
 from collections.abc import Iterator
 import petsird
+from petsird_helpers import get_num_det_els, get_module_and_element
 
 # these are constants for now
 NUMBER_OF_ENERGY_BINS = 3
@@ -99,6 +100,58 @@ def get_scanner_geometry() -> petsird.ScannerGeometry:
     return petsird.ScannerGeometry(replicated_modules=[rep_module], ids=[0])
 
 
+def get_detection_efficiencies(
+    scanner: petsird.ScannerInformation,
+) -> petsird.DetectionEfficiencies:
+    num_det_els = get_num_det_els(scanner.scanner_geometry)
+    det_el_efficiencies = numpy.ones(
+        (num_det_els, scanner.number_of_energy_bins()), dtype=numpy.float32
+    )
+
+    # only rotations of 1 type of module in the current scanner
+    assert len(scanner.scanner_geometry.replicated_modules) == 1
+    rep_module = scanner.scanner_geometry.replicated_modules[0]
+    num_modules = len(rep_module.transforms)
+    # We only need to store geometric efficiencies of the first module with all others.
+    # Assume all module-pairs are in coincidence.
+    # We can now use as SGID=abs(mod_1 - mod_2) -1
+    # This gives SGID = -1 for "self-coincidences", which are normally not recorded
+    num_SGIDs = num_modules
+    module_pair_SGID_LUT = numpy.ndarray((num_modules, num_modules), dtype="int32")
+    for mod1 in range(num_modules):
+        for mod2 in range(num_modules):
+            module_pair_SGID_LUT[mod1, mod2] = int(abs(mod1 - mod2) - 1)
+
+    module_pair_efficiencies_vector = []
+    assert len(rep_module.object.detecting_elements) == 1
+    detecting_elements = rep_module.object.detecting_elements[0]
+    mod1 = 0
+    num_det_els_in_module = len(detecting_elements.transforms)
+    for mod2 in range(1, num_modules):
+        SGID = abs(mod1 - mod2) - 1
+        module_pair_efficiencies = numpy.ones(
+            (
+                num_det_els_in_module,
+                scanner.number_of_energy_bins(),
+                num_det_els_in_module,
+                scanner.number_of_energy_bins(),
+            ),
+            dtype=numpy.float32,
+        )
+        # give some (non-physical) value
+        module_pair_efficiencies *= SGID
+        module_pair_efficiencies_vector.append(
+            petsird.ModulePairEfficiencies(values=module_pair_efficiencies, sgid=SGID)
+        )
+        assert len(module_pair_efficiencies_vector) == SGID + 1
+
+    return petsird.DetectionEfficiencies(
+        det_el_efficiencies=det_el_efficiencies,
+        module_pair_sgidlut=module_pair_SGID_LUT,
+        module_pair_efficiencies_vector=module_pair_efficiencies_vector,
+    )
+
+
 def get_scanner_info() -> petsird.ScannerInformation:
 
     scanner_geometry = get_scanner_geometry()
@@ -112,7 +165,9 @@ def get_scanner_info() -> petsird.ScannerInformation:
     energyBinEdges = numpy.linspace(
         430, 650, NUMBER_OF_ENERGY_BINS + 1, dtype="float32"
     )
-    return petsird.ScannerInformation(
+    # need energy bin info before being able to construct the detection efficiencies
+    # so we will construct a scanner without the efficiencies first
+    scanner = petsird.ScannerInformation(
         model_name="PETSIRD_TEST",
         scanner_geometry=scanner_geometry,
         tof_bin_edges=tofBinEdges,
@@ -122,6 +177,9 @@ def get_scanner_info() -> petsird.ScannerInformation:
         listmode_time_block_duration=1,  # ms
     )
 
+    scanner.detection_efficiencies = get_detection_efficiencies(scanner)
+    return scanner
+
 
 def get_header() -> petsird.Header:
     subject = petsird.Subject(id="123456")
@@ -138,13 +196,25 @@ def get_header() -> petsird.Header:
 def get_events(
     header: petsird.Header, num_events: int
 ) -> Iterator[petsird.CoincidenceEvent]:
-    detector_count = 5  # TODO header.scanner.number_of_detectors()
+    """Generate some random events"""
+    detector_count = get_num_det_els(header.scanner.scanner_geometry)
     for _ in range(num_events):
-        yield petsird.CoincidenceEvent(
-            detector_ids=[
+        # Generate random detector_ids, where the corresponding modules are distinct
+        # This is because in the generated scanner.detection_efficiencies, we assume that
+        # "within module" coincidences are not recorded
+        while True:
+            detector_ids = [
                 random.randrange(0, detector_count),
                 random.randrange(0, detector_count),
-            ],
+            ]
+            mod_and_els = get_module_and_element(
+                header.scanner.scanner_geometry, detector_ids
+            )
+            if abs(mod_and_els[0].module - mod_and_els[1].module) > 0:
+                break
+
+        yield petsird.CoincidenceEvent(
+            detector_ids=detector_ids,
             energy_indices=[
                 random.randrange(0, NUMBER_OF_ENERGY_BINS),
                 random.randrange(0, NUMBER_OF_ENERGY_BINS),

python/petsird_helpers.py

@@ -0,0 +1,86 @@
+"""
+Preliminary helpers for PETSIRD data
+"""
+
+#  Copyright (C) 2024 University College London
+#
+#  SPDX-License-Identifier: Apache-2.0
+
+import numpy
+import petsird
+from dataclasses import dataclass
+import typing
+
+
+def get_num_det_els(scanner_geometry: petsird.ScannerGeometry) -> int:
+    """Compute total number of detecting elements in the scanner"""
+    num_det_els = 0
+    for rep_module in scanner_geometry.replicated_modules:
+        det_els = rep_module.object.detecting_elements
+        for rep_volume in det_els:
+            num_det_els += len(rep_volume.transforms) * len(rep_module.transforms)
+    return num_det_els
+
+
+@dataclass
+class ModuleAndElement:
+    """
+    Dataclass to store a module ID and element ID, where the latter runs
+    over all detecting volumes in the module
+    """
+
+    module: int
+    el: int
+
+
+def get_module_and_element(
+    scanner_geometry: petsird.ScannerGeometry, scanner_det_ids: typing.Iterable[int]
+) -> list[ModuleAndElement]:
+    """Find ModuleAndElement for a list of detector_ids"""
+    assert len(scanner_geometry.replicated_modules) == 1
+    rep_module = scanner_geometry.replicated_modules[0]
+    assert len(rep_module.object.detecting_elements) == 1
+    num_modules = len(rep_module.transforms)
+    return [
+        ModuleAndElement(module=det % num_modules, el=det // num_modules)
+        for det in scanner_det_ids
+    ]
+
+
+def get_detection_efficiency(
+    scanner: petsird.ScannerInformation, event: petsird.CoincidenceEvent
+) -> float:
+    """Compute the detection efficiency for a coincidence event"""
+    eff = 1
+
+    # per det_el efficiencies
+    det_el_efficiencies = scanner.detection_efficiencies.det_el_efficiencies
+    if det_el_efficiencies is not None:
+        eff *= (
+            det_el_efficiencies[event.detector_ids[0], event.energy_indices[0]]
+            * det_el_efficiencies[event.detector_ids[1], event.energy_indices[1]]
+        )
+
+    # per module-pair efficiencies
+    module_pair_efficiencies_vector = (
+        scanner.detection_efficiencies.module_pair_efficiencies_vector
+    )
+    if module_pair_efficiencies_vector is not None:
+        module_pair_SGID_LUT = scanner.detection_efficiencies.module_pair_sgidlut
+        assert module_pair_SGID_LUT is not None
+        mod_and_els = get_module_and_element(
+            scanner.scanner_geometry, event.detector_ids
+        )
+        assert len(scanner.scanner_geometry.replicated_modules) == 1
+        SGID = module_pair_SGID_LUT[mod_and_els[0].module, mod_and_els[1].module]
+        assert SGID >= 0
+        module_pair_efficiencies = module_pair_efficiencies_vector[SGID]
+        assert module_pair_efficiencies.sgid == SGID
+        eff *= module_pair_efficiencies.values[
+            mod_and_els[0].el,
+            event.energy_indices[0],
+            mod_and_els[1].el,
+            event.energy_indices[1],
+        ]
+
+    return eff