pct-tools¶
Tools for proton CT image reconstruction with very large, sparse system matrices.
The package is a thin Python layer over a C++ extension. The system matrix of a proton CT scan has one row per proton and one column per voxel, and for realistic scans it does not fit in memory. pct-tools therefore treats it as a set of row blocks, one per block of proton trajectories. Each block is stored compressed on disk and is read back only while it is needed; the operations loop over blocks in parallel with OpenMP.
The trajectories are typically produced with protonct-sim, a Geant4 application that simulates proton CT acquisitions.
Installation¶
uv pip install pct-tools
Wheels are published for CPython 3.12 on Linux (x86_64, aarch64) and macOS (arm64). Building from source, and the development set-up, are described in the README.
Usage¶
The example below builds the system-matrix block for two protons crossing a small image, stores it together with the measured water-equivalent path lengths, and takes one gradient step towards the reconstructed image. A real reconstruction repeats the last two calls over many blocks and iterations.
import numpy as np
from pct_tools.io import read_compressed_matrix, store_vector
from pct_tools.matrix import compute_and_store_sparse_matrix
from pct_tools.recon import gradient_descent, squared_residuals
# A 2 x 3 x 3 voxel image of 1 mm voxels: 2 slices along y, 3 pixels along x and z.
# The corner of voxel (0, 0, 0) sits at x = 3, y = 0, z = 0; x decreases with the j index.
img_shape = (2, 3, 3)
origin = (3.0, 0.0, 0.0)
# Two straight trajectories along z, sampled every 0.25 mm.
z = np.arange(0.0, 3.01, 0.25, dtype=np.float32)
x = np.stack([np.full_like(z, 1.5), np.full_like(z, 2.5)])
y = np.stack([np.full_like(z, 0.5), np.full_like(z, 1.5)])
z = np.stack([z, z])
compute_and_store_sparse_matrix(
"block_000.mat", x, y, z, img_shape, origin, pixel_width=1.0, slice_thickness=1.0
)
A = read_compressed_matrix("block_000.mat") # scipy CSR matrix of shape (2, 18)
# Measured water-equivalent path lengths, one per proton.
store_vector(np.array([3.0, 3.0], dtype=np.float32), "block_000.vec")
x_image = np.zeros(A.shape[1], dtype=np.float32)
x_image = gradient_descent(x_image, 0.1, ["block_000.mat"], ["block_000.vec"])
print(squared_residuals(x_image, ["block_000.mat"], ["block_000.vec"]))
Every function that loops over blocks or samples takes an n_threads keyword. The optional verbosity_level arguments print progress from the C++ side; 0 is silent.
Array arguments must be numpy arrays of dtype float32 (int32 for the sparse indices passed to construct_matrix). Other dtypes and plain lists raise TypeError rather than being converted: the C++ side works in float32, and a silent conversion would copy the array, which for a block of trajectories or an image vector can be hundreds of megabytes. Convert once with astype(np.float32) where the data is produced.
The chord lengths are exact for the sampled polyline as long as no segment between consecutive samples crosses more than one voxel boundary along any axis; a step of half the smallest voxel dimension is a safe choice, and segments that cross more raise ValueError. examples/benchmark_chords.py times the system-matrix construction for a block of trajectories at several steps and thread counts.
API reference¶
Lengths are in millimetres and angles in radians throughout. Array arguments must be numpy arrays of dtype float32 (int32 for sparse indices); other dtypes raise TypeError rather than being converted. Submodules are not imported by the package itself; import the one you need.
Module |
Purpose |
|---|---|
Chord lengths of trajectories through a voxel grid, i.e. the system-matrix entries; builds and stores a block. |
|
Most likely path of a proton from its entry and exit position and direction. |
|
Binary matrix and vector formats, compressed with Blosc2. |
|
One steepest-descent step and the squared residual norm over all blocks. |
|
|
|
Recompression and row subsampling of stored blocks. |