Merge pull request #22 from KrisThielemans/GeometricDescription

scanner geometry

.gitignore

@@ -10,3 +10,6 @@ python/prd/
 ~$*
 *~
 *.bak
+\#*
+*#
+.#*

cpp/prd_analysis.cpp

@@ -1,11 +1,20 @@
 /*
   Copyright (C) 2022-2023 Microsoft Corporation
-  Copyright (C) 2023 University College London
+  Copyright (C) 2023-2024 University College London
 
   SPDX-License-Identifier: Apache-2.0
 */
 
-#include "generated/hdf5/protocols.h"
+// (un)comment if you want HDF5 or binary output
+#define USE_HDF5
+
+#ifdef USE_HDF5
+#  include "generated/hdf5/protocols.h"
+using prd::hdf5::PrdExperimentReader;
+#else
+#  include "generated/binary/protocols.h"
+using prd::binary::PrdExperimentReader;
+#endif
 #include <xtensor/xview.hpp>
 #include <xtensor/xio.hpp>
 #include <iostream>
@@ -21,16 +30,22 @@ main(int argc, char* argv[])
     }
 
   // Open the file
-  prd::hdf5::PrdExperimentReader reader(argv[1]);
+  PrdExperimentReader reader(argv[1]);
   prd::Header header;
   reader.ReadHeader(header);
 
   std::cout << "Processing file: " << argv[1] << std::endl;
   if (header.exam) // only do this if present
     std::cout << "Subject ID: " << header.exam->subject.id << std::endl;
-  std::cout << "Number of detectors: " << header.scanner.NumberOfDetectors() << std::endl;
+  std::cout << "Types of modules: " << header.scanner.scanner_geometry.replicated_modules.size() << std::endl;
+  std::cout << "Number of modules of first type: " << header.scanner.scanner_geometry.replicated_modules[0].transforms.size()
+            << std::endl;
+  std::cout << "Number of types of detecting elements in modules of first type: "
+            << 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 << "Number of TOF bins: " << header.scanner.NumberOfTOFBins() << std::endl;
-  std::cout << "Number of energy bins: " <<header.scanner.NumberOfEnergyBins() << std::endl;
+  std::cout << "Number of energy bins: " << header.scanner.NumberOfEnergyBins() << std::endl;
 
   const auto& tof_bin_edges = header.scanner.tof_bin_edges;
   std::cout << "TOF bin edges: " << tof_bin_edges << std::endl;

cpp/prd_generator.cpp

@@ -1,62 +1,138 @@
 /*
   Copyright (C) 2022-2023 Microsoft Corporation
-  Copyright (C) 2023 University College London
+  Copyright (C) 2023-2024 University College London
 
   SPDX-License-Identifier: Apache-2.0
 */
 
 #include <iostream>
 #include <cmath>
 #include <random>
-#include "generated/hdf5/protocols.h"
+
+// (un)comment if you want HDF5 or binary output
+#define USE_HDF5
+
+#ifdef USE_HDF5
+#  include "generated/hdf5/protocols.h"
+using prd::hdf5::PrdExperimentWriter;
+#else
+#  include "generated/binary/protocols.h"
+using prd::binary::PrdExperimentWriter;
+#endif
 
 // these are constants for now
-const uint32_t NUMBER_OF_ENERGY_BINS = 3;
-const uint32_t NUMBER_OF_TOF_BINS = 300;
-const float RADIUS = 400.F;
-const uint32_t NUMBER_OF_TIME_BLOCKS = 6;
-const float COUNT_RATE = 500.F;
+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 uint32_t NUMBER_OF_TIME_BLOCKS = 6;
+constexpr float COUNT_RATE = 500.F;
+
+//! return a cuboid volume
+prd::BoxSolidVolume
+get_crystal()
+{
+  using prd::Coordinate;
+  prd::BoxShape crystal_shape{ Coordinate{ { 0, 0, 0 } },
+                               Coordinate{ { 0, 0, CRYSTAL_LENGTH[2] } },
+                               Coordinate{ { 0, CRYSTAL_LENGTH[1], CRYSTAL_LENGTH[2] } },
+                               Coordinate{ { 0, CRYSTAL_LENGTH[1], 0 } },
+                               Coordinate{ { CRYSTAL_LENGTH[0], 0, 0 } },
+                               Coordinate{ { CRYSTAL_LENGTH[0], 0, CRYSTAL_LENGTH[2] } },
+                               Coordinate{ { CRYSTAL_LENGTH[0], CRYSTAL_LENGTH[1], CRYSTAL_LENGTH[2] } },
+                               Coordinate{ { CRYSTAL_LENGTH[0], CRYSTAL_LENGTH[1], 0 } } };
+
+  prd::BoxSolidVolume crystal{ crystal_shape, /* material_id */ 1 };
+  return crystal;
+}
+
+//! return a module of NUM_CRYSTALS_PER_MODULE cuboids
+prd::DetectorModule
+get_detector_module()
+{
+  prd::ReplicatedBoxSolidVolume rep_volume;
+  {
+    rep_volume.object = get_crystal();
+    constexpr auto N0 = NUM_CRYSTALS_PER_MODULE[0];
+    constexpr auto N1 = NUM_CRYSTALS_PER_MODULE[1];
+    constexpr auto N2 = NUM_CRYSTALS_PER_MODULE[2];
+    for (int rep0 = 0; rep0 < N0; ++rep0)
+      for (int rep1 = 0; rep1 < N1; ++rep1)
+        for (int rep2 = 0; rep2 < N2; ++rep2)
+          {
+            prd::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] } } };
+            rep_volume.transforms.push_back(transform);
+            rep_volume.ids.push_back(rep0 + N0 * (rep1 + N1 * rep2));
+          }
+  }
+
+  prd::DetectorModule detector_module;
+  detector_module.detecting_elements.push_back(rep_volume);
+  detector_module.detecting_element_ids.push_back(0);
+
+  return detector_module;
+}
+
+//! return scanner build by rotating a module around the (0,0,1) axis
+prd::ScannerGeometry
+get_scanner_geometry()
+{
+  prd::ReplicatedDetectorModule rep_module;
+  {
+    rep_module.object = get_detector_module();
+    int module_id = 0;
+    std::vector<float> angles;
+    for (int i = 0; i < 10; ++i)
+      {
+        angles.push_back(static_cast<float>(2 * M_PI * i / 10));
+      }
+    for (auto angle : angles)
+      {
+        prd::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);
+      }
+  }
+  prd::ScannerGeometry scanner_geometry;
+  scanner_geometry.replicated_modules.push_back(rep_module);
+  scanner_geometry.ids.push_back(0);
+  return scanner_geometry;
+}
 
-// single ring as example
 prd::ScannerInformation
 get_scanner_info()
 {
-  float radius = RADIUS;
-  std::vector<float> angles;
-  for (int i = 0; i < 10; ++i)
-    {
-      angles.push_back(static_cast<float>(2 * M_PI * i / 10));
-    }
-  std::vector<prd::Detector> detectors;
-  int detector_id = 0;
-  for (auto angle : angles)
-    {
-      // Create a new detector
-      prd::Detector d;
-      d.x = radius * std::sin(angle);
-      d.y = radius * std::cos(angle);
-      d.z = 0.;
-      d.id = detector_id++;
-      detectors.push_back(d);
-    }
-
-  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(tof_bin_edges_shape);
-  for (std::size_t i = 0; i < tof_bin_edges.size(); ++i)
-    tof_bin_edges[i] = (i - NUMBER_OF_TOF_BINS / 2.F) / NUMBER_OF_TOF_BINS * 2 * radius;
-  FArray1D::shape_type energy_bin_edges_shape = { NUMBER_OF_ENERGY_BINS + 1 };
-  FArray1D energy_bin_edges(energy_bin_edges_shape);
-  for (std::size_t i = 0; i < energy_bin_edges.size(); ++i)
-    energy_bin_edges[i] = 430.F + i * (650.F - 430.F) / NUMBER_OF_ENERGY_BINS;
   prd::ScannerInformation scanner_info;
-  scanner_info.detectors = detectors;
-  scanner_info.tof_bin_edges = tof_bin_edges;
-  scanner_info.tof_resolution = 9.4F; // in mm
-  scanner_info.energy_bin_edges = energy_bin_edges;
-  scanner_info.energy_resolution_at_511 = .11F;    // as fraction of 511
-  scanner_info.listmode_time_block_duration = 1.F; // ms
+  scanner_info.model_name = "PETSIRD_TEST";
+
+  scanner_info.scanner_geometry = get_scanner_geometry();
+
+  // TODO scanner_info.bulk_materials
+
+  // TOF and energy information
+  {
+    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(tof_bin_edges_shape);
+    for (std::size_t i = 0; i < tof_bin_edges.size(); ++i)
+      tof_bin_edges[i] = (i - NUMBER_OF_TOF_BINS / 2.F) / NUMBER_OF_TOF_BINS * 2 * RADIUS;
+    FArray1D::shape_type energy_bin_edges_shape = { NUMBER_OF_ENERGY_BINS + 1 };
+    FArray1D energy_bin_edges(energy_bin_edges_shape);
+    for (std::size_t i = 0; i < energy_bin_edges.size(); ++i)
+      energy_bin_edges[i] = 430.F + i * (650.F - 430.F) / NUMBER_OF_ENERGY_BINS;
+    scanner_info.tof_bin_edges = tof_bin_edges;
+    scanner_info.tof_resolution = 9.4F; // in mm
+    scanner_info.energy_bin_edges = energy_bin_edges;
+    scanner_info.energy_resolution_at_511 = .11F;    // as fraction of 511
+    scanner_info.listmode_time_block_duration = 1.F; // ms
+  }
+
   return scanner_info;
 }
 
@@ -99,13 +175,13 @@ get_random_tof_value()
 }
 
 std::vector<prd::CoincidenceEvent>
-get_events(const prd::Header& header, std::size_t num_events)
+get_events(const prd::Header&, std::size_t num_events)
 {
   std::vector<prd::CoincidenceEvent> events;
   events.reserve(num_events);
   for (std::size_t i = 0; i < num_events; ++i)
     {
-      const auto detectors = get_random_pair(header.scanner.NumberOfDetectors());
+      const auto detectors = get_random_pair(1); // TODO header.scanner.NumberOfDetectors());
       prd::CoincidenceEvent e;
       e.detector_ids[0] = detectors.first;
       e.detector_ids[1] = detectors.second;
@@ -129,7 +205,7 @@ main(int argc, char* argv[])
 
   std::string outfile = argv[1];
   std::remove(outfile.c_str());
-  prd::hdf5::PrdExperimentWriter writer(outfile);
+  PrdExperimentWriter writer(outfile);
 
   const auto header = get_header();
   writer.WriteHeader(header);

model/DetectorInformation.yml

@@ -0,0 +1,42 @@
+# A shape filled with a uniform material
+SolidVolume<Shape>: !record
+  fields:
+    shape: Shape
+    # identifier referring to `ScannerInformation.bulkMaterials` list
+    materialId: uint
+
+BoxSolidVolume: SolidVolume<BoxShape>
+GenericSolidVolume: SolidVolume<GeometricShape>
+
+# A list of identical SolidVolumes<BoxShape> at different locations
+ReplicatedBoxSolidVolume: ReplicatedObject< BoxSolidVolume >
+
+# A list of identical SolidVolumes<BGeometricShape> at different locations
+ReplicatedGenericSolidVolume: ReplicatedObject< GenericSolidVolume >
+
+# Top-level detector structure, consisting of one or more lists of detecting elements (or "crystals")
+# This allows having different types of detecting elements (e.g. for phoswich detectors)
+# TODO this could be made into a hierarchical structure
+DetectorModule: !record
+  fields:
+    detectingElements: ReplicatedBoxSolidVolume*
+    # list of unique ids for every replicated solid volume
+    # constraint: size(detectingElements) == size(detectingElementsIds)
+    detectingElementIds: uint*
+    # optional list describing shielding/optical reflectors etc
+    nonDetectingElements: ReplicatedGenericSolidVolume*
+
+# A list of identical modules at different locations
+ReplicatedDetectorModule: ReplicatedObject< DetectorModule >
+
+# Full definition of the geometry of the scanner, consisting of
+# one of more types of modules replicated in space and (optional) other structures (e.g. side-shielding)
+ScannerGeometry: !record
+  fields:
+    # list of different types of replicated modules
+    replicatedModules: ReplicatedDetectorModule*
+    # list of unique ids for every replicated module
+    # constraint: size(replicated_modules) == size(ids)
+    ids: uint*
+    # shielding etc
+    nonDetectingVolumes: GenericSolidVolume*?

model/GeometryInformation.yml

@@ -0,0 +1,64 @@
+# Type definitions related to geometry and coordinates
+
+# 3D coordinates (in mm)
+Coordinate: !record
+  fields:
+    c: float[3]
+
+# 3D direction vector (normalized to 1)
+Direction: !record
+  fields:
+    c: float[3]
+
+# Orthonormal matrix
+# direction_of_first_axis = matrix * [1, 0 ,0] (as a column vector)
+DirectionMatrix: !record
+  fields:
+    matrix: float[3, 3]
+
+# Rigid transformation, encoded via homogenous transformation
+# transformed_coord = matrix * [c, 1] (where [c,1] is a column vector)
+# with `c` of type `Coordinate`
+RigidTransformation: !record
+  fields:
+    matrix: float[3, 4]
+
+# A list of identical objects at different locations
+ReplicatedObject<T>: !record
+  fields:
+    object: T
+    # list of transforms
+    # constraint: length >= 1
+    transforms: RigidTransformation*
+    # list of unique ids for every replicated solid volume
+    # constraint: size(transforms) == size(ids)
+    ids: uint*
+  computedFields:
+    numberOfObjects: size(transforms)
+
+# A box-shape specified by 8 corners (e.g. cuboid, wedge, etc.)
+# TODO need to think about a clear definition of planes
+# We do not want to have to check about intersection planes
+# Potential mechanisms:
+# - lexicographical ordering of corner coordinates?
+# - first 4 coordinates give first plane, 5th and 6th need to define plane with first 2, etc.
+BoxShape: !record
+  fields:
+    corners: Coordinate*8
+
+# Annulus of certain thickness centered at [0,0,0] and oriented along the [0,0,1] axis
+# in radians. An angle of 0 corresponds to the [1,0,0] axis, Pi/2 corresponds to the [0,1,0] axis.
+AnnulusShape: !record
+  fields:
+    # inner radius (in mm)
+    innerRadius: float
+    # outer radius (in mm)
+    outerRadius: float
+    # thickness of the annulus, i.e. length along the axis (in mm)
+    thickness: float
+    # start-stop angle (in radians)
+    angularRange: float*2
+    # center point of the cylinder defining the annulus
+
+# Union of all possible shapes
+GeometricShape: [BoxShape, AnnulusShape]