Merge pull request #35 from KrisThielemans/DetectionEfficiencies

add DetectionEfficiencies (aka "normalisation")

Author: KrisThielemans

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 <xtensor/xview.hpp>
 #include <xtensor/xio.hpp>
 #include <iostream>
@@ -45,6 +46,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;
 
@@ -70,11 +73,23 @@ main(int argc, char* argv[])
           auto& event_time_block = std::get<petsird::EventTimeBlock>(time_block);
           last_time = event_time_block.start;
           num_prompts += event_time_block.prompt_events.size();
+          std::cout << "=====================  Events in time block from " << last_time << " ==============\n";
 
           for (auto& event : event_time_block.prompt_events)
             {
               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";
             }
         }
     }

cpp/petsird_generator.cpp

@@ -8,6 +8,7 @@
 #include <iostream>
 #include <cmath>
 #include <random>
+#include <algorithm>
 
 // (un)comment if you want HDF5 or binary output
 #define USE_HDF5
@@ -20,12 +21,18 @@ 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;
 constexpr float RADIUS = 400.F;
-constexpr std::array<float, 3> CRYSTAL_LENGTH{ 4.F, 4.F, 20.F };
-constexpr std::array<float, 3> NUM_CRYSTALS_PER_MODULE{ 5, 6, 2 };
+constexpr std::array<float, 3> CRYSTAL_LENGTH{ 20.F, 4.F, 4.F };
+constexpr std::array<float, 3> 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;
 
@@ -62,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));
           }
@@ -85,25 +92,95 @@ get_scanner_geometry()
     rep_module.object = get_detector_module();
     int module_id = 0;
     std::vector<float> angles;
-    for (int i = 0; i < 10; ++i)
+    for (unsigned int i = 0; i < NUM_MODULES_ALONG_RING; ++i)
       {
-        angles.push_back(static_cast<float>(2 * M_PI * i / 10));
+        angles.push_back(static_cast<float>((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);
   scanner_geometry.ids.push_back(0);
   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<float>({ num_det_els, scanner.NumberOfEnergyBins() });
+
+  // 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 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<int, 2>({ 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)
+        {
+          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();
+  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();
+  for (unsigned int SGID = 0; SGID < num_SGIDs; ++SGID)
+    {
+      // 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<float, 4>(
+          { 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 +210,8 @@ get_scanner_info()
     scanner_info.event_time_block_duration = 1.F; // ms
   }
 
+  scanner_info.detection_efficiencies = get_detection_efficiencies(scanner_info);
+
   return scanner_info;
 }
 
@@ -154,12 +233,13 @@ get_header()
 }
 
 // return pair of integers between 0 and max
-std::pair<int, int>
+std::array<unsigned, 2>
 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<unsigned, 2> p{ a, b };
+  return p;
 }
 
 uint32_t
@@ -175,16 +255,27 @@ get_random_tof_value()
 }
 
 std::vector<petsird::CoincidenceEvent>
-get_events(const petsird::Header&, std::size_t num_events)
+get_events(const petsird::Header& header, std::size_t num_events)
 {
   std::vector<petsird::CoincidenceEvent> 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
+      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 (!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;
+            }
+        }
       e.energy_indices[0] = get_random_energy_value();
       e.energy_indices[1] = get_random_energy_value();
       e.tof_idx = get_random_tof_value();

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 <class T>
+inline std::vector<ModuleAndElement>
+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<ModuleAndElement> 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<unsigned>(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

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?

model/ScannerInformation.yml

@@ -59,6 +59,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