From 646c2da8b423919e9979f7baddb014f9c96990dc Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Sun, 27 Oct 2024 08:19:04 +0000 Subject: [PATCH 1/6] add DetectionEfficiencies to model --- model/DetectionEfficiencies.yml | 59 +++++++++++++++++++++++++++++++++ model/ScannerInformation.yml | 3 ++ 2 files changed, 62 insertions(+) create mode 100644 model/DetectionEfficiencies.yml diff --git a/model/DetectionEfficiencies.yml b/model/DetectionEfficiencies.yml new file mode 100644 index 0000000..3de31ed --- /dev/null +++ b/model/DetectionEfficiencies.yml @@ -0,0 +1,59 @@ +# Detection efficiencies for every detecting element and energy window +# If size along energyBinIdx is 1, the effiencies are assumed to be the same for every energy bin. +# Constraint: size(DetElEfficiencies, "energyBinIdx") == scannerInformation.numberOfEnergyBins or 1 +DetElEfficiencies: float[detElIdx, energyBinIdx] + +# Efficiency for two detecting elements (det_els) in a pair of modules. +# This is one component (often called "geometric") of the detection efficiency model. +# If size along energyBinIdx is 1, the effiencies are assumed to be the same for every energy bin. +ModulePairEfficiencies: !record + fields: + # Detection efficiency for a pair of detecting elements + # detElIdx1 and detElIdx2 run from 0 to the number of det_els in each module + values: float[detElIdx1, energyBinIdx1, detElIdx2, energyBinIdx2] + # Symmetry Group Identifier (SGID) + # This should be a number between 0 and numberOfSGIDs-1 + sgid: uint + + +# Lookup table for SGIDs +# For every module pair, give the SGID. If -1, the module-pair is not in coincidence. +# Values run from -1 ... (numberOfSGIDs-1) +ModulePairSGIDLUT: int[moduleIdx1, moduleIdx2] + +ModulePairEfficienciesVector: ModulePairEfficiencies* + +# Component-based information on detection efficiencies +# This encodes a simple model for the detection efficiency of (true) coincidences +# consisting of the product of the efficiency of the two detecting elements (det_els) +# and a (geometric) component determined by their location in the two modules. +# The former are stored in detElEfficiencies, and the latter in modulePairEfficienciesVector +# (a list of ModulePairEfficiencies, each entry corresponding to a module pair). +# +# Most PET scanners have some kind of geometric symmetry, e.g. rotation over a full +# module, or translation along the axis of the scanner. Module-pairs that are related +# by such a symmetry often have the same geometric detection efficiencies. PETSIRD +# calls this a "symmetry group" (SG). Each SG had an identifier (SGID). To save memory, +# the modulePairEfficienciesVector needs to contain only one of element for each SGID. +# The SGID for a module-pair can be found in modulePairSGIDLUT. +# +# Finding the total detection efficiency therefore follows these steps: +# (the following pseudo-code ignores energy bins for simplicity) +# 1. find modules for each det_el +# 2. find det_el indices "inside" each module +# 3. SGID = modulePairSGIDLUT[mod1, mod1] +# 4. if (SGID < 0) return 0 +# 5. module_pair_efficiencies = modulePairEfficienciesVector[SGID] +# 6. return detElEfficiencies[det_el1] * detElEfficiencies[det_el2] * module_pair_efficiencies[det_el_in_mod1, det_el_in_mod2] +# +# If either of the components is not present, its value is considered to be 1. +DetectionEfficiencies: !record + fields: + # Detection efficiencies for every detecting element + detElEfficiencies: DetElEfficiencies? + # Lookup table for SGIDs. + # Also indicates if coincidences between a module-pair are recorded. + modulePairSGIDLUT: ModulePairSGIDLUT? + # Vector of all modulePairEfficiencies (one for each SGID) + # Constraint: size(modulePairEfficienciesVector) == max(modulePairSGIDLUT) + 1 + modulePairEfficienciesVector: ModulePairEfficienciesVector? \ No newline at end of file diff --git a/model/ScannerInformation.yml b/model/ScannerInformation.yml index 48515f0..a045659 100644 --- a/model/ScannerInformation.yml +++ b/model/ScannerInformation.yml @@ -58,6 +58,9 @@ ScannerInformation: !record # Encode how the scanner handles multiple coincidences coincidencePolicy: CoincidencePolicy + # coincidence detection efficiencies + detectionEfficiencies: DetectionEfficiencies + computedFields: numberOfTOFBins: size(tofBinEdges)-1 numberOfEnergyBins: size(energyBinEdges)-1 From 0d543ad2ca7b5147703897e0e6b542c53bd2338e Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Sun, 27 Oct 2024 12:12:53 +0000 Subject: [PATCH 2/6] 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 From 98340d74f911019bac9199595ce5e4307ed9aebb Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Sun, 27 Oct 2024 22:37:39 +0000 Subject: [PATCH 3/6] Update cpp example with DetectorEfficiencies --- cpp/petsird_analysis.cpp | 14 +++++++ cpp/petsird_generator.cpp | 83 +++++++++++++++++++++++++++++++++++---- cpp/petsird_helpers.h | 81 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 170 insertions(+), 8 deletions(-) create mode 100644 cpp/petsird_helpers.h diff --git a/cpp/petsird_analysis.cpp b/cpp/petsird_analysis.cpp index ac35254..489c9a7 100644 --- a/cpp/petsird_analysis.cpp +++ b/cpp/petsird_analysis.cpp @@ -15,6 +15,7 @@ using petsird::hdf5::PETSIRDReader; # include "generated/binary/protocols.h" using petsird::binary::PETSIRDReader; #endif +#include "petsird_helpers.h" #include #include #include @@ -44,6 +45,8 @@ main(int argc, char* argv[]) << header.scanner.scanner_geometry.replicated_modules[0].object.detecting_elements.size() << std::endl; std::cout << "Number of elements of first type in modules of first type: " << header.scanner.scanner_geometry.replicated_modules[0].object.detecting_elements[0].transforms.size() << std::endl; + std::cout << "Total number of 'crystals': " << petsird_helpers::get_num_det_els(header.scanner.scanner_geometry) << std::endl; + std::cout << "Number of TOF bins: " << header.scanner.NumberOfTOFBins() << std::endl; std::cout << "Number of energy bins: " << header.scanner.NumberOfEnergyBins() << std::endl; @@ -71,6 +74,17 @@ main(int argc, char* argv[]) { energy_1 += energy_mid_points[event.energy_indices[0]]; energy_2 += energy_mid_points[event.energy_indices[1]]; + + std::cout << "CoincidenceEvent(detectorIds=[" << event.detector_ids[0] << ", " << event.detector_ids[1] + << "], tofIdx=" << event.tof_idx << ", energyIndices=[" << event.energy_indices[0] << ", " + << event.energy_indices[1] << "])\n"; + const auto module_and_elems + = petsird_helpers::get_module_and_element(header.scanner.scanner_geometry, event.detector_ids); + std::cout << " " + << "[ModuleAndElement(module=" << module_and_elems[0].module << ", " + << "el=" << module_and_elems[0].el << "), ModuleAndElement(module=" << module_and_elems[0].module << ", " + << "el=" << module_and_elems[0].el << ")]\n"; + std::cout << " efficiency:" << petsird_helpers::get_detection_efficiency(header.scanner, event) << "\n"; } } diff --git a/cpp/petsird_generator.cpp b/cpp/petsird_generator.cpp index 329aca5..9723e6f 100644 --- a/cpp/petsird_generator.cpp +++ b/cpp/petsird_generator.cpp @@ -20,6 +20,8 @@ using petsird::hdf5::PETSIRDWriter; using petsird::binary::PETSIRDWriter; #endif +#include "petsird_helpers.h" + // these are constants for now constexpr uint32_t NUMBER_OF_ENERGY_BINS = 3; constexpr uint32_t NUMBER_OF_TOF_BINS = 300; @@ -104,6 +106,58 @@ get_scanner_geometry() return scanner_geometry; } +petsird::DetectionEfficiencies +get_detection_efficiencies(const petsird::ScannerInformation& scanner) +{ + const auto num_det_els = petsird_helpers::get_num_det_els(scanner.scanner_geometry); + petsird::DetectionEfficiencies detection_efficiencies; + + detection_efficiencies.det_el_efficiencies = xt::ones({ num_det_els, scanner.NumberOfEnergyBins() }); + + // only rotations of 1 type of module in the current scanner + assert(scanner.scanner_geometry.replicated_modules.size() == 1); + const auto& rep_module = scanner.scanner_geometry.replicated_modules[0]; + const auto num_modules = rep_module.transforms.size(); + + // 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 + const auto num_SGIDs = num_modules - 1; + detection_efficiencies.module_pair_sgidlut = yardl::NDArray({ num_modules, num_modules }); + auto& module_pair_SGID_LUT = *detection_efficiencies.module_pair_sgidlut; + for (unsigned int mod1 = 0; mod1 < num_modules; ++mod1) + { + for (unsigned int mod2 = 0; mod2 < num_modules; ++mod2) + { + module_pair_SGID_LUT(mod1, mod2) = std::abs(int(mod1) - int(mod2)) - 1; + } + } + + // assign an empty vector first, and reserve correct size + detection_efficiencies.module_pair_efficiencies_vector = petsird::ModulePairEfficienciesVector(); + detection_efficiencies.module_pair_efficiencies_vector->reserve(num_SGIDs); + + assert(rep_module.object.detecting_elements.size() == 1); + const auto& detecting_elements = rep_module.object.detecting_elements[0]; + const auto num_det_els_in_module = detecting_elements.transforms.size(); + const unsigned int mod1 = 0; + for (unsigned int mod2 = 1; mod2 < num_modules; ++mod2) + { + const int SGID = std::abs(int(mod1) - int(mod2)) - 1; + petsird::ModulePairEfficiencies module_pair_efficiencies; + module_pair_efficiencies.values = yardl::NDArray( + { num_det_els_in_module, scanner.NumberOfEnergyBins(), num_det_els_in_module, scanner.NumberOfEnergyBins() }); + // give some (non-physical) value + module_pair_efficiencies.values.fill(SGID); + module_pair_efficiencies.sgid = SGID; + detection_efficiencies.module_pair_efficiencies_vector->emplace_back(module_pair_efficiencies); + assert(detection_efficiencies.module_pair_efficiencies_vector->size() == unsigned(SGID + 1)); + } + + return detection_efficiencies; +} + petsird::ScannerInformation get_scanner_info() { @@ -133,6 +187,8 @@ get_scanner_info() scanner_info.listmode_time_block_duration = 1.F; // ms } + scanner_info.detection_efficiencies = get_detection_efficiencies(scanner_info); + return scanner_info; } @@ -154,12 +210,13 @@ get_header() } // return pair of integers between 0 and max -std::pair +std::array get_random_pair(int max) { - int a = rand() % max; - int b = rand() % max; - return std::make_pair(a, b); + unsigned a = rand() % max; + unsigned b = rand() % max; + std::array p{ a, b }; + return p; } uint32_t @@ -175,16 +232,26 @@ get_random_tof_value() } std::vector -get_events(const petsird::Header&, std::size_t num_events) +get_events(const petsird::Header& header, std::size_t num_events) { std::vector events; events.reserve(num_events); + const auto num_det_els = petsird_helpers::get_num_det_els(header.scanner.scanner_geometry); for (std::size_t i = 0; i < num_events; ++i) { - const auto detectors = get_random_pair(1); // TODO header.scanner.NumberOfDetectors()); petsird::CoincidenceEvent e; - e.detector_ids[0] = detectors.first; - e.detector_ids[1] = detectors.second; + // 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) + { + e.detector_ids = get_random_pair(num_det_els); + const auto mod_and_els = petsird_helpers::get_module_and_element(header.scanner.scanner_geometry, e.detector_ids); + if (std::abs(int(mod_and_els[0].module) - mod_and_els[1].module) > 0) + { + break; + } + } e.energy_indices[0] = get_random_energy_value(); e.energy_indices[1] = get_random_energy_value(); e.tof_idx = get_random_tof_value(); diff --git a/cpp/petsird_helpers.h b/cpp/petsird_helpers.h new file mode 100644 index 0000000..b563a92 --- /dev/null +++ b/cpp/petsird_helpers.h @@ -0,0 +1,81 @@ +/* + Copyright (C) 2024 University College London + + SPDX-License-Identifier: Apache-2.0 +*/ + +#include "generated/types.h" + +namespace petsird_helpers +{ + +using namespace petsird; + +inline std::size_t +get_num_det_els(const ScannerGeometry& scanner_geometry) +{ + std::size_t num_det_els = 0; + for (const auto& rep_module : scanner_geometry.replicated_modules) + { + const auto& det_els = rep_module.object.detecting_elements; + for (const auto& rep_volume : det_els) + { + num_det_els += rep_volume.transforms.size() * rep_module.transforms.size(); + } + } + return num_det_els; +} + +struct ModuleAndElement +{ + int module; + int el; +}; + +template +inline std::vector +get_module_and_element(const ScannerGeometry& scanner_geometry, const T& scanner_det_ids) +{ + assert(scanner_geometry.replicated_modules.size() == 1); + const auto& rep_module = scanner_geometry.replicated_modules[0]; + assert(rep_module.object.detecting_elements.size() == 1); + int num_modules = rep_module.transforms.size(); + + std::vector result; + for (int det : scanner_det_ids) + { + result.push_back({ det % num_modules, det / num_modules }); + } + return result; +} + +inline float +get_detection_efficiency(const ScannerInformation& scanner, const CoincidenceEvent& event) +{ + float eff = 1.0F; + const auto& det_el_efficiencies = scanner.detection_efficiencies.det_el_efficiencies; + if (!det_el_efficiencies) + { + eff *= ((*det_el_efficiencies)(event.detector_ids[0], event.energy_indices[0]) + * (*det_el_efficiencies)(event.detector_ids[1], event.energy_indices[1])); + } + const auto& module_pair_efficiencies_vector = scanner.detection_efficiencies.module_pair_efficiencies_vector; + if (!module_pair_efficiencies_vector) + { + const auto& module_pair_SGID_LUT = scanner.detection_efficiencies.module_pair_sgidlut; + assert(!module_pair_SGID_LUT); + + const auto mod_and_els = get_module_and_element(scanner.scanner_geometry, event.detector_ids); + assert(scanner.scanner_geometry.replicated_modules.size() == 1); + const int SGID = (*module_pair_SGID_LUT)(mod_and_els[0].module, mod_and_els[1].module); + assert(SGID >= 0); + + const auto& module_pair_efficiencies = (*module_pair_efficiencies_vector)[SGID]; + assert(module_pair_efficiencies.sgid == static_cast(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; +} + +} // namespace petsird_helpers From 3cebe141f75b4153c9885a6275f8115659cb61d8 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Wed, 30 Oct 2024 00:50:08 +0000 Subject: [PATCH 4/6] Plot detecting volumes as boxes, not lines --- python/petsird_plot_scanner.py | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/python/petsird_plot_scanner.py b/python/petsird_plot_scanner.py index 73da3e7..4fed83a 100644 --- a/python/petsird_plot_scanner.py +++ b/python/petsird_plot_scanner.py @@ -11,9 +11,7 @@ import petsird import matplotlib.pyplot as plt - -# from mpl_toolkits.mplot3d import Axes3D -# import mpl_toolkits +from mpl_toolkits.mplot3d.art3d import Poly3DCollection def transform_to_mat44( @@ -71,9 +69,15 @@ def transform_BoxShape( def draw_BoxShape(ax, box: petsird.BoxShape) -> None: - coords = numpy.array([c.c for c in box.corners]) - # mpl_toolkits.mplot3d.art3d.Line3D(coords[:, 0], coords[:, 1], coords[:, 2]) - ax.plot3D(coords[:, 0], coords[:, 1], coords[:, 2]) + vertices = numpy.array([c.c for c in box.corners]) + edges = [[vertices[j] for j in [0, 1, 2, 3]], + [vertices[j] for j in [4, 5, 6, 7]], + [vertices[j] for j in [0, 1, 5, 4]], + [vertices[j] for j in [2, 3, 7, 6]], + [vertices[j] for j in [1, 2, 6, 5]], + [vertices[j] for j in [4, 7, 3, 0]]] + box = Poly3DCollection(edges, alpha=.25, linewidths=1, edgecolors='r') + ax.add_collection3d(box) if __name__ == "__main__": From 06aa18fef1754c3a7daf660662be38e76f183b54 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Wed, 30 Oct 2024 00:51:08 +0000 Subject: [PATCH 5/6] Generalise generated scanner to multiple "rings" of modules Still using only rotations for symmetries --- python/petsird_analysis.py | 2 +- python/petsird_generator.py | 83 ++++++++++++++++++++++--------------- 2 files changed, 51 insertions(+), 34 deletions(-) diff --git a/python/petsird_analysis.py b/python/petsird_analysis.py index 88b1fc7..f7e993a 100644 --- a/python/petsird_analysis.py +++ b/python/petsird_analysis.py @@ -39,7 +39,7 @@ print(f"Energy bin edges: {energy_bin_edges}") energy_mid_points = (energy_bin_edges[:-1] + energy_bin_edges[1:]) / 2 print(f"Energy mid points: {energy_mid_points}") - + print("SGID LUT:\n", header.scanner.detection_efficiencies.module_pair_sgidlut) energy_1, energy_2 = 0.0, 0.0 num_prompts = 0 last_time = 0 diff --git a/python/petsird_generator.py b/python/petsird_generator.py index 757542a..56b6140 100644 --- a/python/petsird_generator.py +++ b/python/petsird_generator.py @@ -16,10 +16,12 @@ NUMBER_OF_ENERGY_BINS = 3 NUMBER_OF_TOF_BINS = 300 RADIUS = 400 -CRYSTAL_LENGTH = (4, 4, 20) +CRYSTAL_LENGTH = (20, 4, 4) # num crystals in a module -NUM_CRYSTALS_PER_MODULE = (5, 6, 2) -NUM_MODULES = 20 +NUM_CRYSTALS_PER_MODULE = (2, 4,5) +NUM_MODULES_ALONG_RING = 20 +NUM_MODULES_ALONG_AXIS = 2 +MODULE_AXIS_SPACING = (NUM_CRYSTALS_PER_MODULE[2] + 4) * CRYSTAL_LENGTH[2] NUMBER_OF_TIME_BLOCKS = 6 NUMBER_OF_EVENTS = 1000 COUNT_RATE = 500 @@ -60,8 +62,8 @@ def get_detector_module() -> petsird.DetectorModule: matrix=numpy.array( ( (1.0, 0.0, 0.0, RADIUS + rep0 * CRYSTAL_LENGTH[0]), - (0.0, 1.0, 0.0, rep1 * CRYSTAL_LENGTH[1]), - (0.0, 0.0, 1.0, rep2 * CRYSTAL_LENGTH[2]), + (0.0, 1.0, 0.0, (rep1 - N1/2) * CRYSTAL_LENGTH[1]), + (0.0, 0.0, 1.0, (rep2 - N2/2) * CRYSTAL_LENGTH[2]), ), dtype="float32", ) @@ -78,24 +80,25 @@ def get_scanner_geometry() -> petsird.ScannerGeometry: """return a scanner build by rotating a module around the (0,0,1) axis""" detector_module = get_detector_module() radius = RADIUS - angles = [2 * math.pi * i / NUM_MODULES for i in range(NUM_MODULES)] + angles = [2 * math.pi * i / NUM_MODULES_ALONG_RING for i in range(NUM_MODULES_ALONG_RING)] rep_module = petsird.ReplicatedDetectorModule(object=detector_module) module_id = 0 for angle in angles: - transform = petsird.RigidTransformation( - matrix=numpy.array( - ( - (math.cos(angle), math.sin(angle), 0.0, 0.0), - (-math.sin(angle), math.cos(angle), 0.0, 0.0), - (0.0, 0.0, 1.0, 0.0), - ), - dtype="float32", + for ax_mod in range(NUM_MODULES_ALONG_AXIS): + transform = petsird.RigidTransformation( + matrix=numpy.array( + ( + (math.cos(angle), math.sin(angle), 0.0, 0.0), + (-math.sin(angle), math.cos(angle), 0.0, 0.0), + (0.0, 0.0, 1.0, MODULE_AXIS_SPACING*ax_mod), + ), + dtype="float32", + ) ) - ) - rep_module.ids.append(module_id) - module_id += 1 - rep_module.transforms.append(transform) + rep_module.ids.append(module_id) + module_id += 1 + rep_module.transforms.append(transform) return petsird.ScannerGeometry(replicated_modules=[rep_module], ids=[0]) @@ -108,27 +111,42 @@ def get_detection_efficiencies( (num_det_els, scanner.number_of_energy_bins()), dtype=numpy.float32 ) - # only rotations of 1 type of module in the current scanner + # only 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 + # We will only use rotational symmetries (no translation along the axis yet) + # We also assume all module-pairs are in coincidence, except those with the same angle. + # Writing a module number as (z-position, angle): + # eff((z1,a1), (z2, a2)) == eff((z1,0), (z2, abs(a2-a1))) + # or in linear indices + # eff(z1 + NZ * a1, z2 + NZ * a2) == eff(z1, z2 + NZ * abs(a2 - a1)) + # (coincident) SGIDs need to start from 0, so ignoring self-coincident angles + num_SGIDs = NUM_MODULES_ALONG_AXIS * NUM_MODULES_ALONG_AXIS * (NUM_MODULES_ALONG_RING-1) + # SGID = z1 + NZ * (z2 + NZ * abs(a2 - a1) - 1) + NZ = NUM_MODULES_ALONG_AXIS 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) - + z1 = mod1 % NZ + a1 = mod1 // NZ + z2 = mod2 % NZ + a2 = mod2 // NZ + if a1 == a2: + module_pair_SGID_LUT[mod1, mod2]= -1 + else: + module_pair_SGID_LUT[mod1, mod2] = z1 + NZ * (z2 + NZ * (abs(a2 - a1) - 1)) + + # print("SGID LUT:\n", module_pair_SGID_LUT, file=sys.stderr) + assert numpy.max(module_pair_SGID_LUT) == num_SGIDs - 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 + for SGID in range(num_SGIDs): + # extract first module_pair for this SGID. However, as this currently unused, it is commented out + # module_pair = numpy.argwhere(module_pair_SGID_LUT == SGID)[0] + # print(module_pair, file=sys.stderr) module_pair_efficiencies = numpy.ones( ( num_det_els_in_module, @@ -199,9 +217,7 @@ def get_events( """Generate some random events""" detector_count = get_num_det_els(header.scanner.scanner_geometry) for _ in range(num_events): - # 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 + # Generate random detector_ids, where the corresponding modules are in coincidence while True: detector_ids = [ random.randrange(0, detector_count), @@ -210,7 +226,8 @@ def get_events( 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: + if header.scanner.detection_efficiencies.module_pair_sgidlut[mod_and_els[0].module, mod_and_els[1].module] >= 0: + # in coincidence, we can get out of the loop break yield petsird.CoincidenceEvent( From 970e6122f80904ff630339eca54f2ee8056530cc Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Fri, 1 Nov 2024 14:16:36 +0000 Subject: [PATCH 6/6] Generalise generated scanner to multiple "rings" of modules (C++) Examples should now be in sync --- cpp/petsird_analysis.cpp | 2 +- cpp/petsird_generator.cpp | 76 +++++++++++++++++++++++++-------------- 2 files changed, 51 insertions(+), 27 deletions(-) diff --git a/cpp/petsird_analysis.cpp b/cpp/petsird_analysis.cpp index 0f06cb3..885e959 100644 --- a/cpp/petsird_analysis.cpp +++ b/cpp/petsird_analysis.cpp @@ -84,7 +84,7 @@ main(int argc, char* argv[]) << "], tofIdx=" << event.tof_idx << ", energyIndices=[" << event.energy_indices[0] << ", " << event.energy_indices[1] << "])\n"; const auto module_and_elems - = petsird_helpers::get_module_and_element(header.scanner.scanner_geometry, event.detector_ids); + = petsird_helpers::get_module_and_element(header.scanner.scanner_geometry, event.detector_ids); std::cout << " " << "[ModuleAndElement(module=" << module_and_elems[0].module << ", " << "el=" << module_and_elems[0].el << "), ModuleAndElement(module=" << module_and_elems[0].module << ", " diff --git a/cpp/petsird_generator.cpp b/cpp/petsird_generator.cpp index 703cbf0..bbcdb87 100644 --- a/cpp/petsird_generator.cpp +++ b/cpp/petsird_generator.cpp @@ -8,6 +8,7 @@ #include #include #include +#include // (un)comment if you want HDF5 or binary output #define USE_HDF5 @@ -26,8 +27,12 @@ using petsird::binary::PETSIRDWriter; constexpr uint32_t NUMBER_OF_ENERGY_BINS = 3; constexpr uint32_t NUMBER_OF_TOF_BINS = 300; constexpr float RADIUS = 400.F; -constexpr std::array CRYSTAL_LENGTH{ 4.F, 4.F, 20.F }; -constexpr std::array NUM_CRYSTALS_PER_MODULE{ 5, 6, 2 }; +constexpr std::array CRYSTAL_LENGTH{ 20.F, 4.F, 4.F }; +constexpr std::array NUM_CRYSTALS_PER_MODULE{ 2, 4, 5 }; +constexpr uint32_t NUM_MODULES_ALONG_RING{ 20 }; +constexpr uint32_t NUM_MODULES_ALONG_AXIS{ 2 }; +constexpr float MODULE_AXIS_SPACING{ (NUM_CRYSTALS_PER_MODULE[2] + 4) * CRYSTAL_LENGTH[2] }; + constexpr uint32_t NUMBER_OF_TIME_BLOCKS = 6; constexpr float COUNT_RATE = 500.F; @@ -64,8 +69,8 @@ get_detector_module() for (int rep2 = 0; rep2 < N2; ++rep2) { petsird::RigidTransformation transform{ { { 1.0, 0.0, 0.0, RADIUS + rep0 * CRYSTAL_LENGTH[0] }, - { 0.0, 1.0, 0.0, rep1 * CRYSTAL_LENGTH[1] }, - { 0.0, 0.0, 1.0, rep2 * CRYSTAL_LENGTH[2] } } }; + { 0.0, 1.0, 0.0, (rep1 - N1 / 2) * CRYSTAL_LENGTH[1] }, + { 0.0, 0.0, 1.0, (rep2 - N2 / 2) * CRYSTAL_LENGTH[2] } } }; rep_volume.transforms.push_back(transform); rep_volume.ids.push_back(rep0 + N0 * (rep1 + N1 * rep2)); } @@ -87,18 +92,19 @@ get_scanner_geometry() rep_module.object = get_detector_module(); int module_id = 0; std::vector angles; - for (int i = 0; i < 10; ++i) + for (unsigned int i = 0; i < NUM_MODULES_ALONG_RING; ++i) { - angles.push_back(static_cast(2 * M_PI * i / 10)); + angles.push_back(static_cast((2 * M_PI * i) / NUM_MODULES_ALONG_RING)); } for (auto angle : angles) - { - petsird::RigidTransformation transform{ { { std::cos(angle), std::sin(angle), 0.F, 0.F }, - { -std::sin(angle), std::cos(angle), 0.F, 0.F }, - { 0.F, 0.F, 1.F, 0.F } } }; - rep_module.ids.push_back(module_id++); - rep_module.transforms.push_back(transform); - } + for (unsigned ax_mod = 0; ax_mod < NUM_MODULES_ALONG_AXIS; ++ax_mod) + { + petsird::RigidTransformation transform{ { { std::cos(angle), std::sin(angle), 0.F, 0.F }, + { -std::sin(angle), std::cos(angle), 0.F, 0.F }, + { 0.F, 0.F, 1.F, MODULE_AXIS_SPACING * ax_mod } } }; + rep_module.ids.push_back(module_id++); + rep_module.transforms.push_back(transform); + } } petsird::ScannerGeometry scanner_geometry; scanner_geometry.replicated_modules.push_back(rep_module); @@ -114,25 +120,42 @@ get_detection_efficiencies(const petsird::ScannerInformation& scanner) detection_efficiencies.det_el_efficiencies = xt::ones({ num_det_els, scanner.NumberOfEnergyBins() }); - // only rotations of 1 type of module in the current scanner + // only 1 type of module in the current scanner assert(scanner.scanner_geometry.replicated_modules.size() == 1); const auto& rep_module = scanner.scanner_geometry.replicated_modules[0]; const auto num_modules = rep_module.transforms.size(); - // 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 - const auto num_SGIDs = num_modules - 1; + // We will only use rotational symmetries (no translation along the axis yet) + // We also assume all module-pairs are in coincidence, except those with the same angle. + // Writing a module number as (z-position, angle): + // eff((z1,a1), (z2, a2)) == eff((z1,0), (z2, abs(a2-a1))) + // or in linear indices + // eff(z1 + NZ * a1, z2 + NZ * a2) == eff(z1, z2 + NZ * abs(a2 - a1)) + // (coincident) SGIDs need to start from 0, so ignoring self-coincident angles + constexpr auto num_SGIDs = NUM_MODULES_ALONG_AXIS * NUM_MODULES_ALONG_AXIS * (NUM_MODULES_ALONG_RING - 1); + // SGID = z1 + NZ * (z2 + NZ * abs(a2 - a1) - 1) + constexpr auto NZ = NUM_MODULES_ALONG_AXIS; detection_efficiencies.module_pair_sgidlut = yardl::NDArray({ num_modules, num_modules }); auto& module_pair_SGID_LUT = *detection_efficiencies.module_pair_sgidlut; for (unsigned int mod1 = 0; mod1 < num_modules; ++mod1) { for (unsigned int mod2 = 0; mod2 < num_modules; ++mod2) { - module_pair_SGID_LUT(mod1, mod2) = std::abs(int(mod1) - int(mod2)) - 1; + const auto z1 = mod1 % NZ; + const auto a1 = mod1 / NZ; + const auto z2 = mod2 % NZ; + const auto a2 = mod2 / NZ; + if (a1 == a2) + { + module_pair_SGID_LUT(mod1, mod2) = -1; + } + else + { + module_pair_SGID_LUT(mod1, mod2) = z1 + NZ * (z2 + NZ * (std::abs(int(a2) - int(a1)) - 1)); + } } } + // assert(module_pair_SGID_LUT).max() == num_SGIDs - 1); // assign an empty vector first, and reserve correct size detection_efficiencies.module_pair_efficiencies_vector = petsird::ModulePairEfficienciesVector(); @@ -141,10 +164,10 @@ get_detection_efficiencies(const petsird::ScannerInformation& scanner) assert(rep_module.object.detecting_elements.size() == 1); const auto& detecting_elements = rep_module.object.detecting_elements[0]; const auto num_det_els_in_module = detecting_elements.transforms.size(); - const unsigned int mod1 = 0; - for (unsigned int mod2 = 1; mod2 < num_modules; ++mod2) + for (unsigned int SGID = 0; SGID < num_SGIDs; ++SGID) { - const int SGID = std::abs(int(mod1) - int(mod2)) - 1; + // extract first module_pair for this SGID. However, as this currently unused, it is commented out + // const auto& module_pair = *std::find(module_pair_SGID_LUT.begin(), module_pair_SGID_LUT.end(), SGID); petsird::ModulePairEfficiencies module_pair_efficiencies; module_pair_efficiencies.values = yardl::NDArray( { num_det_els_in_module, scanner.NumberOfEnergyBins(), num_det_els_in_module, scanner.NumberOfEnergyBins() }); @@ -241,14 +264,15 @@ get_events(const petsird::Header& header, std::size_t num_events) { petsird::CoincidenceEvent e; // 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) { e.detector_ids = get_random_pair(num_det_els); const auto mod_and_els = petsird_helpers::get_module_and_element(header.scanner.scanner_geometry, e.detector_ids); - if (std::abs(int(mod_and_els[0].module) - mod_and_els[1].module) > 0) + if (!header.scanner.detection_efficiencies.module_pair_sgidlut /* there is no LUT */ + || ((*header.scanner.detection_efficiencies.module_pair_sgidlut)(mod_and_els[0].module, mod_and_els[1].module) + >= 0)) { + // in coincidence, we can get out of the loop break; } }