pct_tools.matrix

System-matrix construction from proton trajectories.

Each trajectory is a polyline of sampled points (x, y, z). Its chord length through every voxel is the length of the polyline inside that voxel, and these chord lengths are the entries of one row of the system matrix.

Conventions shared by all functions in this module:

  • Lengths and coordinates are in the same unit (millimetres in practice); the chord lengths come out in that unit too.

  • img_shape is (n_i, n_j, n_k) and coordinate_origin is (x0, y0, z0). A point falls into the voxel i = floor((y - y0) / slice_thickness), j = floor((x0 - x) / pixel_width), k = floor((z - z0) / pixel_width). Slices are therefore stacked along y, x decreases with increasing j and z increases with k. The flat voxel index used for matrix columns is (i * n_j + j) * n_k + k, i.e. C order of img_shape.

  • Each segment between consecutive samples is cut at the voxel boundaries it crosses, and every piece contributes its length to the voxel containing it. A segment may cross at most one boundary along each axis, so consecutive samples must be less than one voxel apart along every axis; a step of half the smallest voxel dimension is a safe choice. Segments that cross more raise ValueError. Within that limit the chord lengths are exact for the polyline, so the sampling only needs to be fine enough for the polyline to follow the true path.

  • Pieces that lie outside the image contribute nothing.

  • The coordinate arrays must be float32; other dtypes raise TypeError instead of being copied into float32 silently (see the package docstring).

pct_tools.matrix.compute_matrix_elements(x, y, z, img_shape, coordinate_origin, pixel_width, slice_thickness, *, verbose=False, n_threads=1)[source]

Compute the dense array of chord lengths for a single proton trajectory.

Mainly useful for testing and visualisation; use compute_matrix_triplets() or compute_and_store_sparse_matrix() to build system matrices. See the module docstring for the voxel and sampling conventions.

Parameters:

  • x (NDArray[np.float32]) – The x-coordinates of the sampled trajectory, of shape (n_samples,).

  • y (NDArray[np.float32]) – The y-coordinates of the sampled trajectory, of shape (n_samples,).

  • z (NDArray[np.float32]) – The z-coordinates of the sampled trajectory, of shape (n_samples,).

  • img_shape (tuple[int, int, int]) – The number of voxels (n_i, n_j, n_k) along y, x and z.

  • coordinate_origin (tuple[float, float, float]) – The coordinates (x0, y0, z0) of the voxel (0, 0, 0) corner.

  • pixel_width (float) – The voxel size along x and z.

  • slice_thickness (float) – The voxel size along y.

  • verbose (bool) – Print debugging information for every segment of the trajectory.

  • n_threads (int) – The number of OpenMP threads to use.

Returns:

Array of shape img_shape holding the chord length of the trajectory through each voxel, as 32-bit floats.

Raises:

TypeError – If x, y or z is not a float32 array.

Return type:

NDArray[np.float32]

pct_tools.matrix.compute_matrix_triplets(x, y, z, img_shape, coordinate_origin, pixel_width, slice_thickness, *, verbosity_level=0, n_threads=1)[source]

Compute the non-zero system-matrix entries for a block of proton trajectories.

See the module docstring for the voxel and sampling conventions.

Parameters:

  • x (NDArray[np.float32]) – The x-coordinates of the trajectories, of shape (n_trajectories, n_samples) with one row per trajectory. All trajectories in a block must have the same number of samples.

  • y (NDArray[np.float32]) – The y-coordinates of the trajectories, of shape (n_trajectories, n_samples).

  • z (NDArray[np.float32]) – The z-coordinates of the trajectories, of shape (n_trajectories, n_samples).

  • img_shape (tuple[int, int, int]) – The number of voxels (n_i, n_j, n_k) along y, x and z.

  • coordinate_origin (tuple[float, float, float]) – The coordinates (x0, y0, z0) of the voxel (0, 0, 0) corner.

  • pixel_width (float) – The voxel size along x and z.

  • slice_thickness (float) – The voxel size along y.

  • verbosity_level (int) – 0 is silent; 1 prints a summary, 2 and 3 print per-block and per-segment details.

  • n_threads (int) – The number of OpenMP threads to use; trajectories are processed in parallel.

Returns:

A list of (row, column, chord_length) triplets, where row is the trajectory index within the block and column the flat voxel index. A voxel crossed by several segments of the same trajectory appears in several triplets; summing duplicates is left to the caller (scipy.sparse does so by default).

Raises:

TypeError – If x, y or z is not a float32 array.

Return type:

list[tuple[int, int, float]]

pct_tools.matrix.compute_and_store_sparse_matrix(filename, x, y, z, img_shape, coordinate_origin, pixel_width, slice_thickness, verbosity_level=0, n_threads=8, compression_level=9, shuffle_method=ShuffleMethod.BLOSC_SHUFFLE, codec='lz4hc')[source]

Compute the sparse system matrix for a block of proton trajectories and store it.

Calls compute_matrix_triplets(), assembles the triplets into a sparse matrix of shape (n_trajectories, n_i * n_j * n_k) with duplicate entries summed, and stores it in the compressed format described in pct_tools.io. Read it back with pct_tools.io.read_compressed_matrix().

Parameters:

  • filename (str) – The filename of the output file.

  • x (NDArray[np.float32]) – The x-coordinates of the trajectories, of shape (n_trajectories, n_samples) with one row per trajectory.

  • y (NDArray[np.float32]) – The y-coordinates of the trajectories, of shape (n_trajectories, n_samples).

  • z (NDArray[np.float32]) – The z-coordinates of the trajectories, of shape (n_trajectories, n_samples).

  • img_shape (tuple[int, int, int]) – The number of voxels (n_i, n_j, n_k) along y, x and z.

  • coordinate_origin (tuple[float, float, float]) – The coordinates (x0, y0, z0) of the voxel (0, 0, 0) corner.

  • pixel_width (float) – The voxel size along x and z.

  • slice_thickness (float) – The voxel size along y.

  • verbosity_level (int) – 0 is silent; higher values print progressively more detail.

  • n_threads (int) – The number of OpenMP threads to use.

  • compression_level (int) – The Blosc compression level 0-9. Default is 9.

  • shuffle_method (ShuffleMethod) – The Blosc shuffle method. Default is byte shuffling.

  • codec (Codec) – The Blosc codec to use. Default is LZ4HC.

Raises:

TypeError – If x, y or z is not a float32 array.