From 0d543ad2ca7b5147703897e0e6b542c53bd2338e Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Sun, 27 Oct 2024 12:12:53 +0000 Subject: [PATCH] Update Python example with DetectorEfficiencies --- python/petsird_analysis.py | 20 ++++++++- python/petsird_generator.py | 80 +++++++++++++++++++++++++++++++--- python/petsird_helpers.py | 86 +++++++++++++++++++++++++++++++++++++ 3 files changed, 180 insertions(+), 6 deletions(-) create mode 100644 python/petsird_helpers.py diff --git a/python/petsird_analysis.py b/python/petsird_analysis.py index 64ca5f5..88b1fc7 100644 --- a/python/petsird_analysis.py +++ b/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}") diff --git a/python/petsird_generator.py b/python/petsird_generator.py index 8860072..757542a 100644 --- a/python/petsird_generator.py +++ b/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), diff --git a/python/petsird_helpers.py b/python/petsird_helpers.py new file mode 100644 index 0000000..69183ab --- /dev/null +++ b/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