PET projectors `parallelproj.projectors`¶

Convergence comparison: SGD vs SVRG with regularization (sinogram and listmode)

PDHG and LM-SPDHG to optimize the Poisson logL and total variation

Parameters:

event_start_coordinates (Array) – float world coordinates of event LOR start points, shape (num_events, 3)
event_end_coordinates (Array) – float world coordinates of event LOR end points, shape (num_events, 3)
img_shape (tuple[int, int, int]) – shape of the image to be projected
voxel_size (tuple[float, float, float]) – the voxel size of the image to be projected
img_origin (None | Array, optional) – the origin of the image to be projected, by default None means that the center of the image is at world coordinate (0,0,0)

Examples

PET listmode projector (non-TOF and TOF)

Convergence comparison: SGD vs SVRG with regularization (sinogram and listmode)

PDHG and LM-SPDHG to optimize the Poisson logL and total variation

property H: AdjointLinearOperator¶: adjoint operator \(A^H\)

__call__(x: Array) → Array¶

alias to apply(x)

Examples

Parameters:: x (Array)
Return type:: Array

adjoint(y: Array) → Array¶

(scaled) adjoint step \(x = \overline{\alpha} A^H y\)

Parameters:: y (Array)
Return type:: Array

Examples

PET listmode projector (non-TOF and TOF)

adjointness_test(xp: ModuleType | None = None, dev: str | None = None, verbose: bool = False, iscomplex: bool = False, dtype: type | None = None, **kwargs) → bool¶

test whether the adjoint is correctly implemented

Parameters:

xp (ModuleType, optional) – array module to use. Defaults to the operator’s own namespace (self.xp) when it exposes one; otherwise must be given.
dev (str, optional) – device (cpu or cuda). Defaults to the operator’s own device (self.dev) when it exposes one, else the default device.
verbose (bool, optional) – verbose output
iscomplex (bool, optional) – use complex arrays
dtype (type | None, optional) – data type of the arrays
**kwargs (dict) – passed to np.isclose

Returns:

whether the adjoint is correctly implemented

Return type:

bool

Examples

apply(x: Array) → Array¶

(scaled) forward step \(y = \alpha A x\)

Parameters:: x (Array)
Return type:: Array

Examples

property event_end_coordinates: Array¶: coordinates of LOR end points

property event_start_coordinates: Array¶: coordinates of LOR start points

property event_tofbins: None | Array¶

Integer TOF bin index for each event, or None for non-TOF mode.

Assigning an array enables per-event TOF binning; its length must match num_events. Assigning None clears the TOF bins and disables TOF.

property in_shape: tuple[int, int, int]¶: Image shape (n0, n1, n2).

norm(xp: ModuleType | None = None, dev: str | None = None, num_iter: int = 30, iscomplex: bool = False, verbose: bool = False) → float¶

estimate norm of the linear operator using power iterations

Parameters:

xp (ModuleType, optional) – array module to use. Defaults to the operator’s own namespace (self.xp) when it exposes one; otherwise must be given.
dev (str, optional) – device (cpu or cuda). Defaults to the operator’s own device (self.dev) when it exposes one, else the default device.
num_iter (int, optional) – number of power iterations
iscomplex (bool, optional) – use complex arrays
verbose (bool, optional) – verbose output

Returns:

the norm of the linear operator

Return type:

float

Examples

property num_events: int¶: number of events

property out_shape: tuple[int, ...]¶: (num_events,) — one value per detected event.

property scale: float | complex¶: scalar factor applied to the linear operator

property tof: bool¶

Enable or disable TOF mode.

Must set tof_parameters and event_tofbins before setting to True; raises ValueError otherwise. Setting event_tofbins to None automatically disables TOF.

property tof_parameters: TOFParameters | None¶

TOF kernel parameters, or None for non-TOF mode.

Assigning None disables TOF.

property voxel_size: Array¶: voxel size

property xp: ModuleType¶: array module

class parallelproj.projectors.ParallelViewProjector2D(image_shape: tuple[int, int], radial_positions: Array, view_angles: Array, radius: float, image_origin: tuple[float, float], voxel_size: tuple[float, float])[source]¶

Bases: LinearOperator

2D non-TOF parallel view projector

Examples

2D non-TOF filtered back projection (FBP) of Poisson data

Set up a 2D parallel-beam projector using Joseph’s ray-tracing method.

LOR start and end points are computed from the scanner radius and the supplied radial positions and view angles. The image is treated as a 2D slice; internally a unit axial dimension is prepended so the 3D parallelproj_core kernels can be used directly.

Parameters:

image_shape (tuple[int, int]) – Shape of the 2D input image (n1, n2).
radial_positions (Array) – Radial positions of the detector bins in world coordinates (mm).
view_angles (Array) – View angles in radians, one per projection view.
radius (float) – Scanner radius in mm (distance from centre to detector).
image_origin (tuple[float, float]) – World coordinates of the [0, 0] voxel centre in mm (o1, o2). Internally promoted to (0, o1, o2) so the 3D backend can be used.
voxel_size (tuple[float, float]) – Voxel size (d1, d2) in mm. Internally promoted to (1, d1, d2) with a unit axial dimension.

Examples

2D non-TOF filtered back projection (FBP) of Poisson data

property H: AdjointLinearOperator¶: adjoint operator \(A^H\)

__call__(x: Array) → Array¶

alias to apply(x)

Examples

Parameters:: x (Array)
Return type:: Array

adjoint(y: Array) → Array¶

(scaled) adjoint step \(x = \overline{\alpha} A^H y\)

Parameters:: y (Array)
Return type:: Array

Examples

2D non-TOF filtered back projection (FBP) of Poisson data

adjointness_test(xp: ModuleType | None = None, dev: str | None = None, verbose: bool = False, iscomplex: bool = False, dtype: type | None = None, **kwargs) → bool¶

test whether the adjoint is correctly implemented

Parameters:

xp (ModuleType, optional) – array module to use. Defaults to the operator’s own namespace (self.xp) when it exposes one; otherwise must be given.
dev (str, optional) – device (cpu or cuda). Defaults to the operator’s own device (self.dev) when it exposes one, else the default device.
verbose (bool, optional) – verbose output
iscomplex (bool, optional) – use complex arrays
dtype (type | None, optional) – data type of the arrays
**kwargs (dict) – passed to np.isclose

Returns:

whether the adjoint is correctly implemented

Return type:

bool

Examples

apply(x: Array) → Array¶

(scaled) forward step \(y = \alpha A x\)

Parameters:: x (Array)
Return type:: Array

Examples

property dev: str¶: device used for storage of LOR endpoints

property image_origin: Array¶

World coordinates of the [0, 0] voxel, shape (3,).

The first element is always 0 (the prepended axial dimension used internally by the 3D backend); elements 1 and 2 correspond to the image_origin passed at construction.

property image_shape: tuple[int, int]¶: image shape

property in_shape: tuple[int, ...]¶: Image shape (n1, n2).

norm(xp: ModuleType | None = None, dev: str | None = None, num_iter: int = 30, iscomplex: bool = False, verbose: bool = False) → float¶

estimate norm of the linear operator using power iterations

Parameters:

xp (ModuleType, optional) – array module to use. Defaults to the operator’s own namespace (self.xp) when it exposes one; otherwise must be given.
dev (str, optional) – device (cpu or cuda). Defaults to the operator’s own device (self.dev) when it exposes one, else the default device.
num_iter (int, optional) – number of power iterations
iscomplex (bool, optional) – use complex arrays
verbose (bool, optional) – verbose output

Returns:

the norm of the linear operator

Return type:

float

Examples

property num_rad: int¶: number of radial elements

property num_views: int¶: number of views

property out_shape: tuple[int, ...]¶: Projection shape (num_rad, num_views).

property scale: float | complex¶: scalar factor applied to the linear operator

show_views(views_to_show: None | Array = None, image: None | Array = None, **kwargs: Any) → Figure[source]¶

Visualize the geometry of selected projection views.

Parameters:

views_to_show (None | Array) – view numbers to show
image (None | Array) – show an image inside the projector geometry
**kwargs (dict) – passed to matplotlib.pyplot.imshow

Returns:

the matplotlib figure

Return type:

plt.Figure

Examples

property voxel_size: Array¶

Voxel size, shape (3,).

The first element is always 1 (the unit axial dimension used internally by the 3D backend); elements 1 and 2 correspond to the voxel_size passed at construction.

property xend: Array¶: coordinates of LOR end points

property xp: ModuleType¶: array module

property xstart: Array¶: coordinates of LOR start points

class parallelproj.projectors.ParallelViewProjector3D(image_shape: tuple[int, int, int], radial_positions: Array, view_angles: Array, radius: float, image_origin: tuple[float, float, float], voxel_size: tuple[float, float, float], ring_positions: Array, michelogram: Michelogram)[source]¶

Bases: LinearOperator

3D non-TOF parallel view projector supporting any odd span.

Examples

Set up a 3D parallel-beam projector using Joseph’s ray-tracing method.

Extends ParallelViewProjector2D to 3D by adding axial rings with support for any odd sinogram span via a Michelogram. Each sinogram plane’s axial LOR position is determined by the average z-coordinate of the ring pairs contributing to that plane (exact for span=1; the standard averaged-LOR approximation for span>1).

Parameters:

image_shape (tuple[int, int, int]) – Shape of the 3D input image (n0, n1, n2) where n2 is axial.
radial_positions (Array) – Radial positions of the detector bins in world coordinates (mm).
view_angles (Array) – View angles in radians, one per projection view.
radius (float) – Scanner radius in mm.
image_origin (tuple[float, float, float]) – World coordinates of the [0, 0, 0] voxel centre in mm.
voxel_size (tuple[float, float, float]) – Voxel size (d0, d1, d_axial) in mm.
ring_positions (Array) – Axial positions of the detector rings in world coordinates (mm). Must have length equal to michelogram.num_rings.
michelogram (Michelogram) – Axial plane layout encoding the span, max ring difference, and ring-pair grouping. Use Michelogram with span=1 for uncompressed data or any odd span for compressed data.

Examples

property H: AdjointLinearOperator¶: adjoint operator \(A^H\)

__call__(x: Array) → Array¶

alias to apply(x)

Examples

Parameters:: x (Array)
Return type:: Array

adjoint(y: Array) → Array¶

(scaled) adjoint step \(x = \overline{\alpha} A^H y\)

Parameters:: y (Array)
Return type:: Array

Examples

adjointness_test(xp: ModuleType | None = None, dev: str | None = None, verbose: bool = False, iscomplex: bool = False, dtype: type | None = None, **kwargs) → bool¶

test whether the adjoint is correctly implemented

Parameters:

xp (ModuleType, optional) – array module to use. Defaults to the operator’s own namespace (self.xp) when it exposes one; otherwise must be given.
dev (str, optional) – device (cpu or cuda). Defaults to the operator’s own device (self.dev) when it exposes one, else the default device.
verbose (bool, optional) – verbose output
iscomplex (bool, optional) – use complex arrays
dtype (type | None, optional) – data type of the arrays
**kwargs (dict) – passed to np.isclose

Returns:

whether the adjoint is correctly implemented

Return type:

bool

Examples

apply(x: Array) → Array¶

(scaled) forward step \(y = \alpha A x\)

Parameters:: x (Array)
Return type:: Array

Examples

property image_origin: Array¶: image origin - world coordinates of the [0,0,0] voxel

property image_shape: tuple[int, int, int]¶: image shape

property in_shape: tuple[int, int, int]¶: Image shape (n0, n1, n2) where n2 is axial.

property max_ring_diff: int¶: maximum ring difference

property michelogram: Michelogram¶: the Michelogram defining the axial plane layout

norm(xp: ModuleType | None = None, dev: str | None = None, num_iter: int = 30, iscomplex: bool = False, verbose: bool = False) → float¶

estimate norm of the linear operator using power iterations

Parameters:

xp (ModuleType, optional) – array module to use. Defaults to the operator’s own namespace (self.xp) when it exposes one; otherwise must be given.
dev (str, optional) – device (cpu or cuda). Defaults to the operator’s own device (self.dev) when it exposes one, else the default device.
num_iter (int, optional) – number of power iterations
iscomplex (bool, optional) – use complex arrays
verbose (bool, optional) – verbose output

Returns:

the norm of the linear operator

Return type:

float

Examples

property out_shape: tuple[int, int, int]¶: Sinogram shape (num_rad, num_views, num_planes).

property scale: float | complex¶: scalar factor applied to the linear operator

property voxel_size: Array¶: the voxel size in all directions

property xend: Array¶: coordinates of LOR end points

property xp: ModuleType¶: array module

property xstart: Array¶: coordinates of LOR start points

class parallelproj.projectors.RegularPolygonPETProjector(lor_descriptor: RegularPolygonPETLORDescriptor, img_shape: tuple[int, int, int], voxel_size: tuple[float, float, float], img_origin: None | Array = None, views: None | Array = None, cache_lor_endpoints: bool = True)[source]¶

Bases: LinearOperator

geometric non-TOF and TOF sinogram projector for regular polygon PET scanners

Examples

Michelograms and axial sinogram compression

Detector mashing: fewer, bigger virtual detectors

TOF-bin mashing: fewer, wider time-of-flight bins

Listmode to sinogram unlisting

Convergence comparison: SGD vs SVRG with logcosh regularization

PDHG and SPDHG for PET reconstruction with a directional TV prior

TOF vs non-TOF: variance reduction in a uniform cylinder

Convergence comparison: SGD vs SVRG with regularization (sinogram and listmode)

PDHG and LM-SPDHG to optimize the Poisson logL and total variation

Transmission reconstruction: MLTR, SPS and L-BFGS-B

Accelerating MLTR with ordered subsets (OS-MLTR) and SVRG

Penalised transmission reconstruction (MAPTR) with an edge-preserving prior

Joint activity and attenuation reconstruction (MLAA) for TOF PET

pytorch parallelproj projection layer

Parameters:

lor_descriptor (RegularPolygonPETLORDescriptor) – descriptor of the LOR start / end points
img_shape (tuple[int, int, int]) – shape of the image to be projected
voxel_size (tuple[float, float, float]) – the voxel size of the image to be projected
img_origin (None | Array, optional) – the origin of the image to be projected, by default None means that the center of the image is at world coordinate (0,0,0)
views (None | Array, optional) – sinogram views to be projected, by default None means that all views are being projected
cache_lor_endpoints (bool, optional) – whether to cache the LOR endpoints, by default True setting it to False will save memory but will slow down computations

Examples

Michelograms and axial sinogram compression

Detector mashing: fewer, bigger virtual detectors

TOF-bin mashing: fewer, wider time-of-flight bins

Listmode to sinogram unlisting

Convergence comparison: SGD vs SVRG with logcosh regularization

PDHG and SPDHG for PET reconstruction with a directional TV prior

TOF vs non-TOF: variance reduction in a uniform cylinder

Convergence comparison: SGD vs SVRG with regularization (sinogram and listmode)

PDHG and LM-SPDHG to optimize the Poisson logL and total variation

Transmission reconstruction: MLTR, SPS and L-BFGS-B

Accelerating MLTR with ordered subsets (OS-MLTR) and SVRG

Penalised transmission reconstruction (MAPTR) with an edge-preserving prior

Joint activity and attenuation reconstruction (MLAA) for TOF PET

pytorch parallelproj projection layer

property H: AdjointLinearOperator¶: adjoint operator \(A^H\)

__call__(x: Array) → Array¶

alias to apply(x)

Examples

Parameters:: x (Array)
Return type:: Array

adjoint(y: Array) → Array¶

(scaled) adjoint step \(x = \overline{\alpha} A^H y\)

Parameters:: y (Array)
Return type:: Array

Examples

Michelograms and axial sinogram compression

Penalised transmission reconstruction (MAPTR) with an edge-preserving prior

adjointness_test(xp: ModuleType | None = None, dev: str | None = None, verbose: bool = False, iscomplex: bool = False, dtype: type | None = None, **kwargs) → bool¶

test whether the adjoint is correctly implemented

Parameters:

xp (ModuleType, optional) – array module to use. Defaults to the operator’s own namespace (self.xp) when it exposes one; otherwise must be given.
dev (str, optional) – device (cpu or cuda). Defaults to the operator’s own device (self.dev) when it exposes one, else the default device.
verbose (bool, optional) – verbose output
iscomplex (bool, optional) – use complex arrays
dtype (type | None, optional) – data type of the arrays
**kwargs (dict) – passed to np.isclose

Returns:

whether the adjoint is correctly implemented

Return type:

bool

Examples

apply(x: Array) → Array¶

(scaled) forward step \(y = \alpha A x\)

Parameters:: x (Array)
Return type:: Array

Examples

clear_cached_lor_endpoints() → None[source]¶

clear cached LOR endpoints

Examples

Convergence comparison: SGD vs SVRG with logcosh regularization

PDHG and SPDHG for PET reconstruction with a directional TV prior

TOF vs non-TOF: variance reduction in a uniform cylinder

Convergence comparison: SGD vs SVRG with regularization (sinogram and listmode)

Return type:: None

convert_sinogram_to_crystal_index_events(sinogram: Array, shuffle: bool = False) → ndarray[source]¶

Convert a non-TOF or TOF span-1 sinogram to crystal-index events.

Each count in the sinogram becomes one row in the output array. Non-TOF rows are (d_red, r_red, d_blue, r_blue); TOF rows add a trailing tof_bin column in the projector convention (bin 0 = closest to the canonical xstart crystal). The output is ready for direct use with regular_polygon_events_to_sinogram().

The red crystal is the canonical xstart of each LOR as defined by the LOR descriptor; blue is the xend crystal. Unpack the returned array with events[:, 0] (d_red), events[:, 1] (r_red), etc.

Parameters:

sinogram (Array) – Integer span-1 sinogram. Non-TOF shape: lor_descriptor.spatial_sinogram_shape. TOF shape: (*lor_descriptor.spatial_sinogram_shape, num_tof_bins).
shuffle (bool, optional) – Randomly shuffle the output rows (default False). Uses numpy’s global random state; call numpy.random.seed before this method for reproducible results.

Returns:

events – Crystal-index events. Columns are (d_red, r_red, d_blue, r_blue) or (d_red, r_red, d_blue, r_blue, tof_bin).

Return type:

np.ndarray, shape (N, 4) or (N, 5), dtype int32

Raises:

TypeError – If sinogram does not have an integer dtype.
ValueError – If the LOR descriptor has span > 1; start_plane_index is only defined for span-1 descriptors.

Examples

Listmode to sinogram unlisting

convert_sinogram_to_listmode(sinogram: Array, shuffle: bool = False) → tuple[Array, Array, Array | None][source]¶

convert a non-TOF or TOF emission sinogram to listmode events

Parameters:

sinogram (Array) – an integer (TOF or non-TOF) emission sinogram
shuffle (bool, optional) – if True, randomly shuffle the order of the output events, by default False. Shuffling is implemented via numpy.random.permutation(num_events), which draws from numpy’s global random state. Use numpy.random.seed() before calling this method for reproducible results.

Returns:

event_start_coordinates, event_end_coordinates, event_tofbins in case of non-TOF, event_tofbins is None

Return type:

tuple[Array, Array, Array | None]

Examples

Convergence comparison: SGD vs SVRG with regularization (sinogram and listmode)

fov_mask() → Array[source]¶

Boolean cylindrical FOV mask for this projector’s image grid.

The cylinder radius equals the transaxial distance of the midpoint of the first LOR of the first sinogram view from the scanner isocenter. The cylinder axis is aligned with the scanner’s symmetry axis.

Returns:: True inside the cylindrical FOV, False outside.
Return type:: Array of bool, shape in_shape

Examples

Convergence comparison: SGD vs SVRG with logcosh regularization

TOF vs non-TOF: variance reduction in a uniform cylinder

Convergence comparison: SGD vs SVRG with regularization (sinogram and listmode)

Penalised transmission reconstruction (MAPTR) with an edge-preserving prior

property img_origin: Array¶: image origin - world coordinates of the [0,0,0] voxel

property in_shape: tuple[int, int, int]¶: Image shape (n0, n1, n2).

property lor_descriptor: RegularPolygonPETLORDescriptor¶: LOR descriptor

norm(xp: ModuleType | None = None, dev: str | None = None, num_iter: int = 30, iscomplex: bool = False, verbose: bool = False) → float¶

estimate norm of the linear operator using power iterations

Parameters:

xp (ModuleType, optional) – array module to use. Defaults to the operator’s own namespace (self.xp) when it exposes one; otherwise must be given.
dev (str, optional) – device (cpu or cuda). Defaults to the operator’s own device (self.dev) when it exposes one, else the default device.
num_iter (int, optional) – number of power iterations
iscomplex (bool, optional) – use complex arrays
verbose (bool, optional) – verbose output

Returns:

the norm of the linear operator

Return type:

float

Examples

property out_shape: tuple[int, ...]¶: Sinogram shape respecting sinogram_order, optionally with a trailing TOF axis.

property scale: float | complex¶: scalar factor applied to the linear operator

show_geometry(ax: Axes3D, color: tuple[float, float, float] = (1.0, 0.0, 0.0), edgecolor: str = 'grey', alpha: float = 0.1) → None[source]¶

show the geometry of the scanner and the FOV of the image

Parameters:

ax (Axes3D) – matplotlib axes object with projection = ‘3d’
color (tuple[float, float, float], optional) – color to use for the FOV cube, by default (1.,0.,0.)
edgecolor (str, optional) – edgecolor to use for the FOV cube, by default ‘grey’
alpha (float, optional) – alpha value of the FOV cube, by default 0.1

Return type:

None

Examples