Merge pull request #130 from KrisThielemans/geom_helpers

- moved transform functions from plot_scanner to helpers.geometry (Python) petsird_helpers::geometry (C++)
- add equivalent C++ functions
- add get_detecting_box helpers to find the BoxShape for an ExpandedDetectionBin
- C++: add petsird_helpers/create.h to make it easier to create new ScannerInformation objects

Author: KrisThielemans

cpp/CMakeLists.txt

@@ -17,8 +17,14 @@ endif()
 
 add_executable(petsird_generator petsird_generator.cpp)
 target_link_libraries(petsird_generator petsird_generated)
+target_include_directories(petsird_generator PRIVATE "${PROJECT_SOURCE_DIR}/helpers/include")
+# needed for generated PETSIRD files
+target_include_directories(petsird_generator PRIVATE "${PROJECT_SOURCE_DIR}/")
 
 add_executable(petsird_analysis petsird_analysis.cpp)
 target_link_libraries(petsird_analysis petsird_generated)
+target_include_directories(petsird_analysis PRIVATE "${PROJECT_SOURCE_DIR}/helpers/include")
+# needed for generated PETSIRD files
+target_include_directories(petsird_analysis PRIVATE "${PROJECT_SOURCE_DIR}/")
 
 add_subdirectory(generated)

cpp/helpers/include/petsird_helpers/create.h

@@ -0,0 +1,78 @@
+/*
+  Copyright (C) 2025 University College London
+
+  SPDX-License-Identifier: Apache-2.0
+*/
+
+#include "generated/types.h"
+#include <array>
+
+namespace petsird_helpers
+{
+namespace create
+{
+
+using namespace petsird;
+
+//! Helper function to create a std::vector<T>
+/*! This function is added to have a 1D analogue construct_2D_nested_vector() */
+template <typename T>
+std::vector<T>
+construct_vector(std::size_t size)
+{
+  std::vector<T> v(size);
+  return v;
+}
+
+//! Helper function to create a std::vector<std::vector<T>> as a 2D array of size (size0, size1)
+template <typename T>
+std::vector<std::vector<T>>
+construct_2D_nested_vector(std::size_t size0, std::size_t size1)
+{
+  std::vector<std::vector<T>> v(size0);
+  for (auto& one_of_them : v)
+    {
+      one_of_them.resize(size1);
+    }
+  return v;
+}
+
+//! Set  various structures to have the correct size for the given num_types_of_modules
+/*!
+  This will set scanner.tof_bin_edges, scanner.tof_resolution,
+  scanner.event_energy_bin_edges, scanner.energy_resolution_at_511, and
+  and (optionally) scanner.detection_efficiencies.detection_bin_efficiencies,
+  scanner.detection_efficiencies.module_pair_sgidlut and
+  scanner.detection_efficiencies.module_pair_efficiencies_vectors
+  to (nested) vectors of the appropriate type and size.
+
+  Elements will be constructed via the default constructors, so you will still have
+  to fill in the actual values.
+*/
+inline void
+initialize_scanner_information_dimensions(petsird::ScannerInformation& scanner, const std::size_t num_module_types,
+                                          bool allocate_detection_bin_efficiencies, bool allocate_module_pair_efficiencies)
+{
+  scanner.tof_bin_edges = construct_2D_nested_vector<petsird::BinEdges>(num_module_types, num_module_types);
+  scanner.tof_resolution = construct_2D_nested_vector<float>(num_module_types, num_module_types);
+  scanner.event_energy_bin_edges = construct_vector<petsird::BinEdges>(num_module_types);
+  scanner.energy_resolution_at_511 = construct_vector<float>(num_module_types);
+
+  scanner.detection_efficiencies = petsird::DetectionEfficiencies();
+
+  if (allocate_detection_bin_efficiencies)
+    {
+      scanner.detection_efficiencies.detection_bin_efficiencies
+          = construct_vector<petsird::DetectionBinEfficiencies>(num_module_types);
+    }
+  if (allocate_module_pair_efficiencies)
+    {
+      scanner.detection_efficiencies.module_pair_sgidlut
+          = construct_2D_nested_vector<petsird::ModulePairSGIDLUT>(num_module_types, num_module_types);
+      scanner.detection_efficiencies.module_pair_efficiencies_vectors
+          = construct_2D_nested_vector<petsird::ModulePairEfficienciesVector>(num_module_types, num_module_types);
+    }
+}
+
+} // namespace create
+} // namespace petsird_helpers

cpp/helpers/include/petsird_helpers/geometry.h

@@ -0,0 +1,98 @@
+#include <xtensor/xarray.hpp>
+#include <xtensor/xview.hpp>
+#include <xtensor/xio.hpp>
+#include <xtensor-blas/xlinalg.hpp>
+#include <vector>
+#include "generated/types.h"
+
+namespace petsird_helpers
+{
+namespace geometry
+{
+
+//! Convert RigidTransformation to 4x4 matrix
+inline xt::xarray<float>
+transform_to_mat44(const RigidTransformation& transform)
+{
+  xt::xarray<float> bottom_row = { { 0.0f, 0.0f, 0.0f, 1.0f } };
+  return xt::concatenate(xt::xtuple(transform.matrix, bottom_row), 0);
+}
+
+//! Convert 4x4 matrix to RigidTransformation
+inline RigidTransformation
+mat44_to_transform(const xt::xarray<float>& mat)
+{
+  RigidTransformation transform;
+  transform.matrix = xt::view(mat, xt::range(0, 3), xt::all());
+  return transform;
+}
+
+//! Convert Coordinate to homogeneous vector
+inline xt::xarray<float>
+coordinate_to_homogeneous(const Coordinate& coord)
+{
+  return xt::concatenate(xt::xtuple(coord.c, xt::xarray<float>{ 1.0f }), 0);
+}
+
+//! Convert homogeneous vector to Coordinate
+inline Coordinate
+homogeneous_to_coordinate(const xt::xarray<float>& hom_coord)
+{
+  Coordinate coord;
+  coord.c = xt::view(hom_coord, xt::range(0, 3));
+  return coord;
+}
+
+//! Multiply a list of transformations
+inline RigidTransformation
+mult_transforms(const std::vector<RigidTransformation>& transforms)
+{
+  xt::xarray<float> mat = xt::eye<float>(4);
+  for (auto it = transforms.rbegin(); it != transforms.rend(); ++it)
+    {
+      mat = xt::linalg::dot(transform_to_mat44(*it), mat);
+    }
+  return mat44_to_transform(mat);
+}
+
+//! Apply transformations to a coordinate
+inline Coordinate
+mult_transforms_coord(const std::vector<RigidTransformation>& transforms, const Coordinate& coord)
+{
+  xt::xarray<float> hom = xt::linalg::dot(transform_to_mat44(mult_transforms(transforms)), coordinate_to_homogeneous(coord));
+  return homogeneous_to_coordinate(hom);
+}
+
+//! Transform a BoxShape
+inline BoxShape
+transform_BoxShape(const RigidTransformation& transform, const BoxShape& box_shape)
+{
+  BoxShape new_box;
+  for (size_t i = 0; i < box_shape.corners.size(); ++i)
+    {
+      new_box.corners[i] = mult_transforms_coord({ transform }, box_shape.corners[i]);
+    }
+  return new_box;
+}
+
+//! find the BoxShape corresponding to a ExpandedDetectionBin
+inline BoxShape
+get_detecting_box(const ScannerInformation& scanner, const TypeOfModule& type_of_module,
+                  const ExpandedDetectionBin& expanded_detection_bin)
+{
+  const auto& rep_module = scanner.scanner_geometry.replicated_modules[type_of_module];
+  const auto& det_els = rep_module.object.detecting_elements;
+  const auto& mod_transform = rep_module.transforms[expanded_detection_bin.module_index];
+  const auto& transform = det_els.transforms[expanded_detection_bin.element_index];
+  return transform_BoxShape(mult_transforms({ mod_transform, transform }), det_els.object.shape);
+}
+
+//! find the BoxShape corresponding to a DetectionBin
+inline BoxShape
+get_detecting_box(const ScannerInformation& scanner, const TypeOfModule& type_of_module, const DetectionBin& detection_bin)
+{
+  return get_detecting_box(scanner, type_of_module, expand_detection_bin(scanner, type_of_module, detection_bin));
+}
+
+} // namespace geometry
+} // namespace petsird_helpers

cpp/petsird_analysis.cpp

@@ -16,6 +16,7 @@ using petsird::hdf5::PETSIRDReader;
 using petsird::binary::PETSIRDReader;
 #endif
 #include "petsird_helpers.h"
+#include "petsird_helpers/geometry.h"
 #include <xtensor/xview.hpp>
 #include <xtensor/xio.hpp>
 #include <iostream>
@@ -99,7 +100,16 @@ main(int argc, char const* argv[])
             << std::endl;
   std::cout << "Total number of 'crystals' in modules of first type : "
             << petsird_helpers::get_num_det_els(header.scanner, type_of_module) << std::endl;
-
+  // example printing of coordinates
+  {
+    const petsird::ExpandedDetectionBin expanded_detection_bin{ 0, 0, 0 };
+    const auto box_shape = petsird_helpers::geometry::get_detecting_box(header.scanner, type_of_module, expanded_detection_bin);
+    std::cout << "Coordinates of first detecting bin:\n";
+    for (auto& corner : box_shape.corners)
+      {
+        std::cout << "  [" << corner.c[0] << ", " << corner.c[1] << ", " << corner.c[2] << "]\n";
+      }
+  }
   const auto& tof_bin_edges = header.scanner.tof_bin_edges[type_of_module][type_of_module];
   const auto num_tof_bins = tof_bin_edges.NumberOfBins();
   std::cout << "Number of TOF bins: " << num_tof_bins << std::endl;

cpp/petsird_generator.cpp

@@ -22,6 +22,7 @@ using petsird::binary::PETSIRDWriter;
 #endif
 
 #include "petsird_helpers.h"
+#include "petsird_helpers/create.h"
 
 // these are constants for now
 constexpr uint32_t NUMBER_OF_EVENT_ENERGY_BINS = 3;
@@ -108,33 +109,30 @@ get_scanner_geometry()
   return scanner_geometry;
 }
 
-petsird::DetectionEfficiencies
-get_detection_efficiencies(const petsird::ScannerInformation& scanner)
+// set some example efficiencies in the ScannerInformation object.
+void
+set_detection_efficiencies(petsird::ScannerInformation& scanner)
 {
   const auto num_module_types = scanner.scanner_geometry.replicated_modules.size();
   // only 1 type of module in the current scanner
   assert(num_module_types == 1);
   const petsird::TypeOfModule type_of_module{ 0 };
   const auto num_detection_bins = petsird_helpers::get_num_detection_bins(scanner, type_of_module);
-  petsird::DetectionEfficiencies detection_efficiencies;
 
   const auto& event_energy_bin_edges = scanner.event_energy_bin_edges[type_of_module];
   const auto num_event_energy_bins = event_energy_bin_edges.NumberOfBins();
 
-  std::vector<petsird::DetectionBinEfficiencies> detection_bin_efficiencies(num_module_types);
-  detection_efficiencies.detection_bin_efficiencies = detection_bin_efficiencies;
-  detection_efficiencies.module_pair_sgidlut = std::vector<std::vector<petsird::ModulePairSGIDLUT>>(num_module_types);
-  for (auto& elem : *detection_efficiencies.module_pair_sgidlut)
+  // set all detection_bin_efficiencies to 1 in this example
+  if (scanner.detection_efficiencies.detection_bin_efficiencies)
     {
-      elem.resize(num_module_types);
+      auto& bin_effs = (*scanner.detection_efficiencies.detection_bin_efficiencies)[type_of_module];
+      yardl::resize(bin_effs, { num_detection_bins });
+      std::fill(begin(bin_effs), end(bin_effs), 1.F);
     }
-  detection_efficiencies.module_pair_efficiencies_vectors
-      = std::vector<std::vector<petsird::ModulePairEfficienciesVector>>(num_module_types);
-  for (auto& elem : *detection_efficiencies.module_pair_efficiencies_vectors)
-    {
-      elem.resize(num_module_types);
-    }
-  (*detection_efficiencies.detection_bin_efficiencies)[type_of_module] = xt::ones<float>({ num_detection_bins });
+
+  // check if the caller wants to have module-pair stuff. If not, return.
+  if (!scanner.detection_efficiencies.module_pair_efficiencies_vectors)
+    return;
 
   const auto& rep_module = scanner.scanner_geometry.replicated_modules[type_of_module];
   const auto num_modules = rep_module.transforms.size();
@@ -149,9 +147,9 @@ get_detection_efficiencies(const petsird::ScannerInformation& scanner)
   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)[type_of_module][type_of_module]
-      = yardl::NDArray<int, 2>({ num_modules, num_modules });
-  auto& module_pair_SGID_LUT = (*detection_efficiencies.module_pair_sgidlut)[type_of_module][type_of_module];
+
+  auto& module_pair_SGID_LUT = (*scanner.detection_efficiencies.module_pair_sgidlut)[type_of_module][type_of_module];
+  module_pair_SGID_LUT = yardl::NDArray<int, 2>({ num_modules, num_modules });
   for (unsigned int mod1 = 0; mod1 < num_modules; ++mod1)
     {
       for (unsigned int mod2 = 0; mod2 < num_modules; ++mod2)
@@ -174,7 +172,7 @@ get_detection_efficiencies(const petsird::ScannerInformation& scanner)
 
   // initialise module_pair_efficiencies
   auto& module_pair_efficiencies_vector
-      = (*detection_efficiencies.module_pair_efficiencies_vectors)[type_of_module][type_of_module];
+      = (*scanner.detection_efficiencies.module_pair_efficiencies_vectors)[type_of_module][type_of_module];
   // assign an empty vector first, and reserve correct size
   module_pair_efficiencies_vector = petsird::ModulePairEfficienciesVector();
   module_pair_efficiencies_vector.reserve(num_SGIDs);
@@ -194,8 +192,6 @@ get_detection_efficiencies(const petsird::ScannerInformation& scanner)
       module_pair_efficiencies_vector.emplace_back(module_pair_efficiencies);
       assert(module_pair_efficiencies_vector.size() == unsigned(SGID + 1));
     }
-
-  return detection_efficiencies;
 }
 
 petsird::ScannerInformation
@@ -206,47 +202,46 @@ get_scanner_info()
 
   scanner_info.scanner_geometry = get_scanner_geometry();
   const auto num_types_of_modules = scanner_info.scanner_geometry.replicated_modules.size();
+  // Pre-allocate various structures to have the correct size for num_types_of_modules
+  // (We will still have to set descent values into each of these.)
+  petsird_helpers::create::initialize_scanner_information_dimensions(scanner_info, num_types_of_modules,
+                                                                     /* allocate_detection_bin_efficiencies = */ true,
+                                                                     /* allocate_module_pair_efficiencies = */ true);
 
   // TODO scanner_info.bulk_materials
 
   // TOF and energy information
   {
-    auto all_tof_bin_edges
-        = petsird_helpers::construct_2D_nested_vector<petsird::BinEdges>(num_types_of_modules, num_types_of_modules);
-    auto all_tof_resolutions = petsird_helpers::construct_2D_nested_vector<float>(num_types_of_modules, num_types_of_modules);
-    auto all_event_energy_bin_edges = petsird_helpers::construct_vector<petsird::BinEdges>(num_types_of_modules);
-    auto all_event_energy_resolutions = petsird_helpers::construct_vector<float>(num_types_of_modules);
+    auto& all_tof_bin_edges = scanner_info.tof_bin_edges;
+    auto& all_tof_resolutions = scanner_info.tof_resolution;
+    auto& all_event_energy_bin_edges = scanner_info.event_energy_bin_edges;
+    auto& all_event_energy_resolutions = scanner_info.energy_resolution_at_511;
 
     // only 1 type of module in the current scanner
     assert(num_types_of_modules == 1);
     const petsird::TypeOfModule type_of_module{ 0 };
 
     typedef yardl::NDArray<float, 1> FArray1D;
     // TOF info (in mm)
-    FArray1D::shape_type tof_bin_edges_shape = { NUMBER_OF_TOF_BINS + 1 };
-    FArray1D tof_bin_edges_arr(tof_bin_edges_shape);
+    FArray1D tof_bin_edges_arr;
+    yardl::resize(tof_bin_edges_arr, { NUMBER_OF_TOF_BINS + 1 });
     for (std::size_t i = 0; i < tof_bin_edges_arr.size(); ++i)
       tof_bin_edges_arr[i] = (i - NUMBER_OF_TOF_BINS / 2.F) / NUMBER_OF_TOF_BINS * 2 * RADIUS;
     const petsird::BinEdges tof_bin_edges{ tof_bin_edges_arr };
     all_tof_bin_edges[type_of_module][type_of_module] = tof_bin_edges;
 
     all_tof_resolutions[type_of_module][type_of_module] = 9.4F; // in mm
 
-    FArray1D::shape_type event_energy_bin_edges_shape = { NUMBER_OF_EVENT_ENERGY_BINS + 1 };
-    FArray1D event_energy_bin_edges_arr(event_energy_bin_edges_shape);
+    FArray1D event_energy_bin_edges_arr;
+    yardl::resize(event_energy_bin_edges_arr, { NUMBER_OF_EVENT_ENERGY_BINS + 1 });
     for (std::size_t i = 0; i < event_energy_bin_edges_arr.size(); ++i)
       event_energy_bin_edges_arr[i] = 430.F + i * (650.F - 430.F) / NUMBER_OF_EVENT_ENERGY_BINS;
     petsird::BinEdges event_energy_bin_edges{ event_energy_bin_edges_arr };
     all_event_energy_bin_edges[type_of_module] = event_energy_bin_edges;
     all_event_energy_resolutions[type_of_module] = .11F; // as fraction of 511
-
-    scanner_info.tof_bin_edges = all_tof_bin_edges;
-    scanner_info.tof_resolution = all_tof_resolutions;
-    scanner_info.event_energy_bin_edges = all_event_energy_bin_edges;
-    scanner_info.energy_resolution_at_511 = all_event_energy_resolutions;
   }
 
-  scanner_info.detection_efficiencies = get_detection_efficiencies(scanner_info);
+  set_detection_efficiencies(scanner_info);
 
   scanner_info.coincidence_policy = petsird::CoincidencePolicy::kRejectMultiples;
   scanner_info.delayed_coincidences_are_stored = false;

cpp/petsird_helpers.h

@@ -12,29 +12,6 @@ namespace petsird_helpers
 
 using namespace petsird;
 
-//! Helper function to create a std::vector<T>
-/*! This function is added to have a 1D analogue construct_2D_nested_vector() */
-template <typename T>
-std::vector<T>
-construct_vector(std::size_t size)
-{
-  std::vector<T> v(size);
-  return v;
-}
-
-//! Helper function to create a std::vector<std::vector<T>> as a 2D array of size (size0, size1)
-template <typename T>
-std::vector<std::vector<T>>
-construct_2D_nested_vector(std::size_t size0, std::size_t size1)
-{
-  std::vector<std::vector<T>> v(size0);
-  for (auto& one_of_them : v)
-    {
-      one_of_them.resize(size1);
-    }
-  return v;
-}
-
 // TODO remove?
 inline std::size_t
 get_num_det_els(const ScannerInformation& scanner, const TypeOfModule& type_of_module)