From 016df4d266c3d6bff1c7caba172424b4e4b4c7c1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Coussat?= Date: Wed, 12 Nov 2025 17:56:23 +0100 Subject: [PATCH 1/3] WIP: New PCT format --- documentation/docs/pct_format.md | 39 ++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/documentation/docs/pct_format.md b/documentation/docs/pct_format.md index 86df4cb..f5263b7 100644 --- a/documentation/docs/pct_format.md +++ b/documentation/docs/pct_format.md @@ -4,6 +4,45 @@ PCT uses its own data format to store list-mode proton CT data as [MetaImage](ht Some PCT executables produce such data, e.g. `pctpairprotons` from [ROOT](https://root.cern/) data produced by a [GATE](https://github.com/OpenGATE/opengate) simulation, whereas other ones take such data as input, e.g. `pctbinning`. PCT uses an [ITK](https://itk.org/) image internally (`itk::Image`) to read from / write to these [MetaImage](https://itk.org/Wiki/ITK/MetaIO/Documentation) files. +## Data format + +A PCT list-mode dataset is a two-dimensional $N\times M$ MetaImage where +* $N$ is the number of proton pairs, i.e. each row corresponds to one proton pair, and +* $M$ is the number of columns in the dataset. + +More precisely, the header contains a series of key/value pairs that specify at which column indexes a given information can be found. For instance, if the header specifies that `UpstreamPositionU = 2`, then the third element of each proton pair corresponds to the $u$ coordinate of the proton as detected by the upstream detector. + +**Important note:** the data is stored as an ITK image file to facilitate its processing using ITK, but it does **not** represent an image! For this reason, some fields that typically appear in usual header files are meaningless here, such as `ElementSpacing` or `Offset`. Their values are not interpreted in any way. + +Below is a complete description of the fields that PCT can handle. For backward-compatibility, some fields have a default value which corresponds to the legacy format (detailed below). + +| Header key | Description | Default value | +|------------------------|--------------------------------------------------------------|---------------| +| `UpstreamPositionU` | $u$ coordinate of the proton at the upstream detector | 0 | +| `UpstreamPositionV` | $v$ coordinate of the proton at the upstream detector | 1 | +| `UpstreamPositionW` | $w$ coordinate of the proton at the upstream detector | 2 | +| `DownstreamPositionU` | $u$ coordinate of the proton at the downstream detector | 3 | +| `DownstreamPositionV` | $v$ coordinate of the proton at the downstream detector | 4 | +| `DownstreamPositionW` | $w$ coordinate of the proton at the downstream detector | 5 | +| `UpstreamDirectionU` | Direction along $u$ of the proton at the upstream detector | 6 | +| `UpstreamDirectionV` | Direction along $v$ of the proton at the upstream detector | 7 | +| `UpstreamDirectionW` | Direction along $w$ of the proton at the upstream detector | 8 | +| `DownstreamDirectionU` | Direction along $u$ of the proton at the downstream detector | 9 | +| `DownstreamDirectionV` | Direction along $v$ of the proton at the downstream detector | 10 | +| `DownstreamDirectionW` | Direction along $w$ of the proton at the downstream detector | 11 | +| `UpstreamEnergy` | Kinetic energy of the proton at the upstream detector | 12 | +| `DownstreamEnergy` | Kinetic energy of the proton at the downstream detector | 13 | +| `TrackID` | `TrackID` of the proton | 14 | +| `WEPL` | Water-equivalent path length | ‒ | +| `CreatorProcess` | Process which created the exit particle | ‒ | +| `NuclearProcess` | Whether the particle encountered a nuclear interaction | ‒ | +| `Order` | Number of particle interactions | ‒ | +| `TOF` | Time-of-flight between the two detectors | ‒ | + +## Legacy data format + +PCT used to use a different data format, which was discontinued due to its rigidity. The new data format should be backward-compatible with files that use the old data format. The description below provides details about this legacy data format. + Proton pairs are stored in 2D images of 3D float vectors, in which the first dimension is a series of 5 or 6 vectors of 3 floats each, one per proton pair, and the second dimension corresponds to the number of proton pairs. Each vector stores the following data: 1. $(u_{\text{in}},v_{\text{in}},w_{\text{in}})$ position of the proton at the entrance detector; From 752b1e103a26101292e8e1508b9b250acc4a54ec Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Coussat?= Date: Thu, 13 Nov 2025 11:40:27 +0100 Subject: [PATCH 2/3] WIP: First Python POC that uses the new format --- poc_new_format.py | 55 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) create mode 100644 poc_new_format.py diff --git a/poc_new_format.py b/poc_new_format.py new file mode 100644 index 0000000..faa1e5c --- /dev/null +++ b/poc_new_format.py @@ -0,0 +1,55 @@ +# Dummy example of usage of the new PCT format, +# which defines a dataset that only contains upstream position of the protons. +import itk +import numpy as np + +TEST_FILE = "/tmp/test.mhd" + +# Create a list-mode dataset that uses the new format +data = itk.image_from_array(np.random.rand(10, 3)) +data["UpstreamPositionU"] = 0 +data["UpstreamPositionV"] = 1 +data["UpstreamPositionW"] = 2 + +itk.imwrite(data, TEST_FILE) + + +# Read from a file that uses the new format +class ProtonPairs: + + pct_fields = [ + "UpstreamPositionU", + "UpstreamPositionV", + "UpstreamPositionW", + "DownstreamPositionU", + "DownstreamPositionV", + "DownstreamPositionW", + ] + + def __init__(self, file_path): + im = itk.imread(file_path) + + self.data = itk.array_from_image(im) + self.metadata = dict(im) + + def __getitem__(self, idx): + # __getitem__ is enough to implement the iterator protocol + data_idx = self.data[idx] + item = {} + for pct_field in ProtonPairs.pct_fields: + try: + item[pct_field] = data_idx[int(self.metadata[pct_field])] + except KeyError: + pass + return item + + +proton_pairs = ProtonPairs(TEST_FILE) + +for i, p in enumerate(proton_pairs): + print("Proton pair", i, ": ", p) + +p_iter = iter(proton_pairs) +print(next(p_iter)) +print(next(p_iter)) +print(next(p_iter)) From 33036764121e235f93b6cc9e0844d3b13bf3e7e9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Coussat?= Date: Tue, 25 Nov 2025 15:46:15 +0100 Subject: [PATCH 3/3] WIP: First C++ POC that uses the new format --- applications/CMakeLists.txt | 1 + applications/pctpocnewformat/CMakeLists.txt | 30 ++++++++++ .../pctpocnewformat/pctpocnewformat.cxx | 59 +++++++++++++++++++ .../pctpocnewformat/pctpocnewformat.ggo | 4 ++ include/pctProtonPairsIterator.h | 39 ++++++++++++ include/pctProtonPairsIterator.hxx | 44 ++++++++++++++ 6 files changed, 177 insertions(+) create mode 100644 applications/pctpocnewformat/CMakeLists.txt create mode 100644 applications/pctpocnewformat/pctpocnewformat.cxx create mode 100644 applications/pctpocnewformat/pctpocnewformat.ggo create mode 100644 include/pctProtonPairsIterator.h create mode 100644 include/pctProtonPairsIterator.hxx diff --git a/applications/CMakeLists.txt b/applications/CMakeLists.txt index 53f75e4..d44bb52 100644 --- a/applications/CMakeLists.txt +++ b/applications/CMakeLists.txt @@ -20,6 +20,7 @@ include(${ITK_USE_FILE}) add_subdirectory(pctbackprojectionbinning) add_subdirectory(pctbinning) add_subdirectory(pctfdk) +add_subdirectory(pctpocnewformat) add_subdirectory(pctpaircuts) add_subdirectory(pctpairgeometry) add_subdirectory(pctprojections) diff --git a/applications/pctpocnewformat/CMakeLists.txt b/applications/pctpocnewformat/CMakeLists.txt new file mode 100644 index 0000000..4274c6d --- /dev/null +++ b/applications/pctpocnewformat/CMakeLists.txt @@ -0,0 +1,30 @@ +wrap_ggo(pctpocnewformat_GGO_C pctpocnewformat.ggo) +add_executable( + pctpocnewformat + pctpocnewformat.cxx + ${pctpocnewformat_GGO_C} +) +target_link_libraries(pctpocnewformat PCT) + +# Ensure binary is placed in the central PCT binary dir +set_target_properties( + pctpocnewformat + PROPERTIES + RUNTIME_OUTPUT_DIRECTORY + "${PCT_BINARY_DIR}/bin" +) + +# Installation code +install( + TARGETS + pctpocnewformat + RUNTIME + DESTINATION ${PCT_INSTALL_RUNTIME_DIR} + COMPONENT Runtime + LIBRARY + DESTINATION ${PCT_INSTALL_LIB_DIR} + COMPONENT RuntimeLibraries + ARCHIVE + DESTINATION ${PCT_INSTALL_ARCHIVE_DIR} + COMPONENT Development +) diff --git a/applications/pctpocnewformat/pctpocnewformat.cxx b/applications/pctpocnewformat/pctpocnewformat.cxx new file mode 100644 index 0000000..a39678d --- /dev/null +++ b/applications/pctpocnewformat/pctpocnewformat.cxx @@ -0,0 +1,59 @@ +#include "pctpocnewformat_ggo.h" +#include "pctProtonPairsIterator.h" + +#include +#include + +const std::string TEST_FILE = "/tmp/test.mhd"; + +int +main(int argc, char * argv[]) +{ + GGO(pctpocnewformat, args_info); // RTK macro parsing options from .ggo file (rtkMacro.h) + + using VectorType = itk::Vector; + using ImageType = itk::Image; + + // Create a list-mode dataset that uses the new format + + auto img = ImageType::New(); + const ImageType::IndexType start = { { 0 } }; + const ImageType::SizeType size = { { 10 } }; + ImageType::RegionType region; + region.SetSize(size); + region.SetIndex(start); + img->SetRegions(region); + img->Allocate(); + + itk::ImageRegionIterator it(img, img->GetLargestPossibleRegion()); + float p = 0.; + for (it.GoToBegin(); !it.IsAtEnd(); ++it) + { // Fill the data with dummy data + float x = p + 1.; + float y = p + 2.; + float z = p + 3.; + it.Set(itk::MakeVector(x, y, z)); + p++; + } + + auto metaDataDictionary = img->GetMetaDataDictionary(); + itk::EncapsulateMetaData(metaDataDictionary, "UpstreamPositionU", 0); + itk::EncapsulateMetaData(metaDataDictionary, "UpstreamPositionV", 1); + itk::EncapsulateMetaData(metaDataDictionary, "UpstreamPositionW", 2); + img->SetMetaDataDictionary(metaDataDictionary); + + itk::WriteImage(img, TEST_FILE); + + // Read from a file that uses the new format + + auto img2 = itk::ReadImage(TEST_FILE); + + pct::ProtonPairsIterator it2(img2.GetPointer(), img2->GetLargestPossibleRegion()); + for (it2.GoToBegin(); !it2.IsAtEnd(); ++it2) + { + auto pos = it2.GetUpstreamPosition(); + std::cout << pos[0] << " " << pos[1] << " " << pos[2] << std::endl; + } + + return EXIT_SUCCESS; +} diff --git a/applications/pctpocnewformat/pctpocnewformat.ggo b/applications/pctpocnewformat/pctpocnewformat.ggo new file mode 100644 index 0000000..926ee26 --- /dev/null +++ b/applications/pctpocnewformat/pctpocnewformat.ggo @@ -0,0 +1,4 @@ +package "pct" +version "Test new format" + +option "verbose" v "Verbose execution" flag off diff --git a/include/pctProtonPairsIterator.h b/include/pctProtonPairsIterator.h new file mode 100644 index 0000000..0801184 --- /dev/null +++ b/include/pctProtonPairsIterator.h @@ -0,0 +1,39 @@ +#ifndef __pctProtonPairsIterator_h +#define __pctProtonPairsIterator_h + +#include "itkImageRegionConstIterator.h" + +namespace pct +{ + +template +class ITK_TEMPLATE_EXPORT ProtonPairsIterator : public itk::ImageRegionConstIterator> +{ +public: + using ValueType = typename TVectorType::ValueType; + using ImageType = itk::Image; + using Superclass = itk::ImageRegionConstIterator; + using RegionType = typename Superclass::RegionType; + + ProtonPairsIterator(const ImageType * img, const RegionType & region); + + itk::Vector + GetUpstreamPosition(); + +protected: + int + GetFieldIndex(const itk::MetaDataDictionary & metaDataDictionary, const std::string key); + + int m_upsteamPositionU; + int m_upsteamPositionV; + int m_upsteamPositionW; + +}; // end of class + +} // end namespace pct + +#ifndef ITK_MANUAL_INSTANTIATION +# include "pctProtonPairsIterator.hxx" +#endif + +#endif diff --git a/include/pctProtonPairsIterator.hxx b/include/pctProtonPairsIterator.hxx new file mode 100644 index 0000000..9e9d54a --- /dev/null +++ b/include/pctProtonPairsIterator.hxx @@ -0,0 +1,44 @@ +#include "pctProtonPairsIterator.h" + +#include "itkMetaDataObject.h" + +namespace pct +{ + +template +ProtonPairsIterator::ProtonPairsIterator(const ImageType * img, const RegionType & region) + : itk::ImageRegionConstIterator(img, region) +{ + auto metaDataDictionary = img->GetMetaDataDictionary(); + this->m_upsteamPositionU = GetFieldIndex(metaDataDictionary, "UpstreamPositionU"); + this->m_upsteamPositionV = GetFieldIndex(metaDataDictionary, "UpstreamPositionV"); + this->m_upsteamPositionW = GetFieldIndex(metaDataDictionary, "UpstreamPositionW"); +} + +template +int +ProtonPairsIterator::GetFieldIndex(const itk::MetaDataDictionary & metaDataDictionary, + const std::string key) +{ + std::string metaData; + if (itk::ExposeMetaData(metaDataDictionary, key, metaData)) + { + // will throw a std::invalid_argument if number cannot be parsed + return std::stoi(metaData); + } + return -1; // indicates that the key is not defined in the metadata dictionary +} + +template +itk::Vector +ProtonPairsIterator::GetUpstreamPosition() +{ + auto currentVec = this->Get(); + itk::Vector vec; + vec.SetNthComponent(0, currentVec[m_upsteamPositionU]); + vec.SetNthComponent(1, currentVec[m_upsteamPositionV]); + vec.SetNthComponent(2, currentVec[m_upsteamPositionW]); + return vec; +} + +} // namespace pct