Skip to content

WIP: New PCT format #66

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
acoussat wants to merge 3 commits into
base: main
Choose a base branch
from
Open

WIP: New PCT format #66

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
016df4d
WIP: New PCT format
acoussat Nov 12, 2025
752b1e1
WIP: First Python POC that uses the new format
acoussat Nov 13, 2025
3303676
WIP: First C++ POC that uses the new format
acoussat Nov 25, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions applications/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
30 changes: 30 additions & 0 deletions applications/pctpocnewformat/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
)
59 changes: 59 additions & 0 deletions applications/pctpocnewformat/pctpocnewformat.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,59 @@
#include "pctpocnewformat_ggo.h"
#include "pctProtonPairsIterator.h"

#include <rtkMacro.h>
#include <rtkGgoFunctions.h>

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<float, 3>;
using ImageType = itk::Image<VectorType, 1>;

// 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<ImageType> 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<int>(metaDataDictionary, "UpstreamPositionU", 0);
itk::EncapsulateMetaData<int>(metaDataDictionary, "UpstreamPositionV", 1);
itk::EncapsulateMetaData<int>(metaDataDictionary, "UpstreamPositionW", 2);
img->SetMetaDataDictionary(metaDataDictionary);

itk::WriteImage(img, TEST_FILE);

// Read from a file that uses the new format

auto img2 = itk::ReadImage<ImageType>(TEST_FILE);

pct::ProtonPairsIterator<VectorType> 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;
}
4 changes: 4 additions & 0 deletions applications/pctpocnewformat/pctpocnewformat.ggo
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
package "pct"
version "Test new format"

option "verbose" v "Verbose execution" flag off
39 changes: 39 additions & 0 deletions documentation/docs/pct_format.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
39 changes: 39 additions & 0 deletions include/pctProtonPairsIterator.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,39 @@
#ifndef __pctProtonPairsIterator_h
#define __pctProtonPairsIterator_h

#include "itkImageRegionConstIterator.h"

namespace pct
{

template <typename TVectorType>
class ITK_TEMPLATE_EXPORT ProtonPairsIterator : public itk::ImageRegionConstIterator<itk::Image<TVectorType, 1>>
{
public:
using ValueType = typename TVectorType::ValueType;
using ImageType = itk::Image<TVectorType, 1>;
using Superclass = itk::ImageRegionConstIterator<ImageType>;
using RegionType = typename Superclass::RegionType;

ProtonPairsIterator(const ImageType * img, const RegionType & region);

itk::Vector<ValueType, 3>
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
44 changes: 44 additions & 0 deletions include/pctProtonPairsIterator.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,44 @@
#include "pctProtonPairsIterator.h"

#include "itkMetaDataObject.h"

namespace pct
{

template <typename TVectorType>
ProtonPairsIterator<TVectorType>::ProtonPairsIterator(const ImageType * img, const RegionType & region)
: itk::ImageRegionConstIterator<ImageType>(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 <typename TVectorType>
int
ProtonPairsIterator<TVectorType>::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 <typename TVectorType>
itk::Vector<typename TVectorType::ValueType, 3>
ProtonPairsIterator<TVectorType>::GetUpstreamPosition()
{
auto currentVec = this->Get();
itk::Vector<ValueType, 3> vec;
vec.SetNthComponent(0, currentVec[m_upsteamPositionU]);
vec.SetNthComponent(1, currentVec[m_upsteamPositionV]);
vec.SetNthComponent(2, currentVec[m_upsteamPositionW]);
return vec;
}

} // namespace pct
55 changes: 55 additions & 0 deletions poc_new_format.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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))
Loading