Low-level projector interface `parallelproj.backend`

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

LOR descriptors and sinogram definition

Parameters:

xp (ModuleType) – array module to use for storing the LOR endpoints
dev (str) – device to use for storing the LOR endpoints
radius (float) – radius of the scanner
num_sides (int) – number of sides (faces) of each regular polygon
num_lor_endpoints_per_side (int) – number of LOR endpoints in each side (face) of each polygon
lor_spacing (float) – spacing between the LOR endpoints in each side
ring_positions (Array) – 1D array with the coordinate of the rings along the ring axis
symmetry_axis (int) – the ring axis (0,1,2)
phis (None | Array, optional) – angle of each side, by default None means that the sides are equally spaced around a circle

property all_lor_endpoints: Array: the world coordinates of all LOR endpoints

property all_lor_endpoints_index_in_ring: Array: the index within the ring (regular polygon) of all LOR endpoints

property all_lor_endpoints_index_offset: Array: the offset in the linear (flattend) index for all LOR endpoints

property all_lor_endpoints_module_number: Array: the module number of all LOR endpoints

property all_lor_endpoints_ring_number: Array: the ring (regular polygon) number of all LOR endpoints

property dev: str: device to use for storing the LOR endpoints

get_lor_endpoints(module: Array, index_in_module: Array) → Array

get the coordinates for LOR endpoints defined by module and index in module

Parameters:

module (Array) – the module number of the LOR endpoints
index_in_module (Array) – the index in module number of the LOR endpoints

Returns:

the 3 world coordinates of the LOR endpoints

Return type:

Array

linear_lor_endpoint_index(module: Array, index_in_module: Array) → Array

transform the module + index_in_modules indices into a flattened / linear LOR endpoint index

Parameters:

module (Array) – containing module numbers
index_in_module (Array) – containing index in modules

Returns:

the flattened LOR endpoint index

Return type:

Array

property lor_spacing: float: the spacing between the LOR endpoints in every side (face) of each polygon

property modules: tuple[PETScannerModule]: tuple of modules defining the scanner

property num_lor_endpoints: int: the total number of LOR endpoints in the scanner

property num_lor_endpoints_per_module: Array: numpy array showing how many LOR endpoints are in every module

property num_lor_endpoints_per_ring: int: the number of LOR endpoints per ring (regular polygon)

property num_lor_endpoints_per_side: int: number of LOR endpoints per side (face) in each polygon

property num_modules: int: the number of modules defining the scanner

property num_rings: int: number of rings (regular polygons)

property num_sides: int: number of sides (faces) of each polygon

property radius: float: radius of the scanner

property ring_positions: Array: the ring (regular polygon) positions

setup_all_lor_endpoints() → None

calculate the position of all lor endpoints by iterating over the modules and calculating the transformed coordinates of all module endpoints

Return type:: None

show_lor_endpoints(ax: Axes, show_linear_index: bool = True, **kwargs) → None

show all LOR endpoints in a 3D plot

Parameters:

ax (plt.Axes) – a 3D matplotlib axes
show_linear_index (bool, optional) – annotate the LOR endpoints with the linear LOR endpoint index
**kwargs (keyword arguments) – passed to show_lor_endpoints() of the scanner module

Return type:

None

property symmetry_axis: int: The symmetry axis. Also called axial (or ring) direction.

property xp: ModuleType: array module to use for storing the LOR endpoints

class parallelproj.pet_scanners.RegularPolygonPETScannerModule(xp: ModuleType, dev: str, radius: float, num_sides: int, num_lor_endpoints_per_side: int, lor_spacing: float, ax0: int = 2, ax1: int = 1, affine_transformation_matrix: Array | None = None, phis: None | Array = None)[source]

Bases: PETScannerModule

Regular polygon PET scanner module (detectors on a regular polygon)

Parameters:

xp (ModuleType) – array module to use for storing the LOR endpoints
device (str) – device to use for storing the LOR endpoints
radius (float) – inner radius of the regular polygon
num_sides (int) – number of sides of the regular polygon
num_lor_endpoints_per_sides (int) – number of LOR endpoints per side
lor_spacing (float) – spacing between the LOR endpoints in the polygon direction
ax0 (int, optional) – axis number for the first direction, by default 2
ax1 (int, optional) – axis number for the second direction, by default 1
affine_transformation_matrix (Array | None, optional) – 4x4 affine transformation matrix applied to the LOR endpoint coordinates, default None if None, the 4x4 identity matrix is used
phis (None | Array, optional) – angle of each side, by default None means that the sides are equally spaced around a circle
dev (str)
num_lor_endpoints_per_side (int)

property affine_transformation_matrix: Array

4x4 affine transformation matrix

Return type:: Array

property ax0: int

axis number for the first module direction

Return type:: int

property ax1: int

axis number for the second module direction

Return type:: int

property dev: str: device to use for storing the LOR endpoints

get_lor_endpoints(inds: Array | None = None) → Array

mapping from LOR endpoint indices within module to an array of “transformed” world coordinates

Parameters:: inds (Array | None, optional) – an non-negative integer array of indices, default None if None means all possible indices [0, … , num_lor_endpoints - 1]
Returns:: a 3 x len(inds) float array with the world coordinates of the LOR endpoints including an affine transformation
Return type:: Array

get_raw_lor_endpoints(inds: Array | None = None) → Array[source]

mapping from LOR endpoint indices within module to an array of “raw” world coordinates

Parameters:: inds (Array | None, optional) – an non-negative integer array of indices, default None if None means all possible indices [0, … , num_lor_endpoints - 1]
Returns:: a 3 x len(inds) float array with the world coordinates of the LOR endpoints
Return type:: Array

property lor_endpoint_numbers: Array

array enumerating all the LOR endpoints in the module

Return type:: Array

property lor_spacing: float

spacing between the LOR endpoints in a module along the polygon

Return type:: float

property num_lor_endpoints: int

total number of LOR endpoints in the module

Return type:: int

property num_lor_endpoints_per_side: int

number of LOR endpoints per side

Return type:: int

property num_sides: int

number of sides of the regular polygon

Return type:: int

property phis: Array

azimuthal angle of each side

Return type:: Array

property radius: float

inner radius of the regular polygon

Return type:: float

show_lor_endpoints(ax: Axes, annotation_fontsize: float = 0, annotation_prefix: str = '', annotation_offset: int = 0, transformed: bool = True, **kwargs) → None

show the LOR coordinates in a 3D scatter plot

Parameters:

ax (plt.Axes) – 3D matplotlib axes
annotation_fontsize (float, optional) – fontsize of LOR endpoint number annotation, by default 0
annotation_prefix (str, optional) – prefix for annotation, by default ‘’
annotation_offset (int, optional) – number to add to crystal number, by default 0
transformed (bool, optional) – use transformed instead of raw coordinates, by default True

Return type:

None

property xp: ModuleType: array module to use for storing the LOR endpoints

PET LOR / sinogram descriptors `parallelproj.pet_lors`

description of PET LORs (and sinograms bins) consisting of two detector endpoints

class parallelproj.pet_lors.EqualBlockPETLORDescriptor(scanner: ModularizedPETScannerGeometry, all_block_pairs: Array)[source]

Bases: PETLORDescriptor

LOR descriptor for scanner consisting of block modules where each block module has the same number of LOR endpoints

Parameters:

scanner (ModularizedPETScannerGeometry) – A modularized PET scanner consisting of block modules with the same number of LOR endpoints.
all_block_pairs (Array) – An array containing pairs of integer numbers encoding which block pairs are in coincidence and form valid LORs.

Return type:

None

property all_block_pairs: Array: all block pairs in coincidence

property dev: str: device to use for storing the LOR endpoints

get_lor_coordinates(block_pair_nums: None | Array = None) → tuple[Array, Array][source]

Get the coordinates of LORs for the given block pair numbers.

Parameters:

block_pair_nums (None or Array, optional) – The block pair numbers for which to retrieve the LOR coordinates. If None, all block pair numbers will be used.

Returns:

tuple[Array, Array]
A tuple containing two arrays –
- the start coordinates of the LORs, with shape (N, 3), where N is the total number of LORs.
- the end coordinates of the LORs, with shape (N, 3)

Return type:

tuple[Array, Array]

property num_block_pairs: int: number of block pairs in coincidence

property num_lorendpoints_per_block: int: number of LOR endpoints per block

property num_lors_per_block_pair: int: number of LORs per block pair

property scanner: ModularizedPETScannerGeometry: the scanner for which coincidences are described

show_block_pair_lors(ax: Axes, block_pair_nums: Array, lw: float = 0.2, **kwargs) → None[source]

show all LORs connecting all endpoints between blocks forming a block pairs

Parameters:

ax (plt.Axes) – a 3D matplotlib axes
block_pair_nums (int) – the block pair numbers to show
lw (float, optional) – the line width, by default 0.2

Return type:

None

property xp: ModuleType: array module to use for storing the LOR endpoints

class parallelproj.pet_lors.PETLORDescriptor(scanner: ModularizedPETScannerGeometry)[source]

Bases: ABC

abstract base class to describe which modules / indices in modules of a modularized PET scanner are in coincidence; defining geometrical LORs

Parameters:: scanner (ModularizedPETScannerGeometry) – a modularized PET scanner

property dev: str: device to use for storing the LOR endpoints

abstract get_lor_coordinates() → tuple[Array, Array][source]

return the start and end coordinates of all (or a subset of) LORs

Return type:: tuple[Array, Array]

property scanner: ModularizedPETScannerGeometry: the scanner for which coincidences are described

property xp: ModuleType: array module to use for storing the LOR endpoints

class parallelproj.pet_lors.RegularPolygonPETLORDescriptor(scanner: RegularPolygonPETScannerGeometry, radial_trim: int = 3, max_ring_difference: int | None = None, sinogram_order: SinogramSpatialAxisOrder = SinogramSpatialAxisOrder.RVP)[source]

Bases: PETLORDescriptor

LOR descriptor for a regular polygon PET scanner where we have coincidences within and between “rings (polygons of modules)” The geometrical LORs can be sorted into a sinogram having a “plane”, “view” and “radial” axis.

Examples

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

LOR descriptors and sinogram definition

Parameters:

scanner (RegularPolygonPETScannerGeometry) – a regular polygon PET scanner
radial_trim (int, optional) – number of geometrial LORs to disregard in the radial direction, by default 3
max_ring_difference (int | None, optional) – maximim ring difference to consider for coincidences, by default None means all ring differences are included
sinogram_order (SinogramSpatialAxisOrder, optional) – the order of the sinogram axes, by default SinogramSpatialAxisOrder.RVP

property dev: str: device to use for storing the LOR endpoints

property end_in_ring_index: Array: end index within ring for all views - shape (num_view, num_rad)

property end_plane_index: Array: end plane for all planes

get_distributed_views_and_slices(num_subsets: int, num_dim: int) → tuple[list[Array], list[tuple[slice, ...]]][source]

distribute sinogram views numbers into subsets

Parameters:

num_subsets (int) – number of subsets
num_dim (int) – number of dimensions of the sinogram to setup the subset slices (e.g. 3 for non-TOF, 4 for TOF)

Returns:

subset views numbers and subset slices

Return type:

tuple[list[Array], list[tuple[slice, …]]]

get_lor_coordinates(views: None | Array = None) → tuple[Array, Array][source]

return the start and end coordinates of all LORs / or a subset of views

Parameters:: views (None | Array, optional) – the views to consider, by default None means all views
Returns:: xstart, xend – 2 dimensional floating point arrays containing the start and end coordinates of all LORs
Return type:: Array

property max_ring_difference: int: the maximum ring difference

property num_planes: int: number of planes in the sinogram

property num_rad: int: number of radial elements in the sinogram

property num_views: int: number of views in the sinogram

property plane_axis_num: int: the axis number of the plane axis

property radial_axis_num: int: the axis number of the radial axis

property radial_trim: int: number of geometrial LORs to disregard in the radial direction

property scanner: ModularizedPETScannerGeometry: the scanner for which coincidences are described

show_views(ax: Axes, views: Array, planes: Array, lw: float = 0.2, **kwargs) → None[source]

show all LORs of a single view in a given plane

Parameters:

ax (plt.Axes) – a 3D matplotlib axes
view (int) – the view number
plane (int) – the plane number
lw (float, optional) – the line width, by default 0.2
views (Array)
planes (Array)

Return type:

None

property sinogram_order: SinogramSpatialAxisOrder: the order of the sinogram axes

property spatial_sinogram_shape: tuple[int, int, int]: the shape of the sinogram in spatial order

property start_in_ring_index: Array: start index within ring for all views - shape (num_view, num_rad)

property start_plane_index: Array: start plane for all planes

property view_axis_num: int: the axis number of the view axis

property xp: ModuleType: array module to use for storing the LOR endpoints

class parallelproj.pet_lors.SinogramSpatialAxisOrder(value, names=<not given>, *values, module=None, qualname=None, type=None, start=1, boundary=None)[source]

Bases: Enum

order of spatial axis in a sinogram R (radial), V (view), P (plane)

Examples

PRV = 5: [plane,radial,view]

PVR = 6: [plane,view,radial]

RPV = 2: [radial,plane,view]

RVP = 1: [radial,view,plane]

VPR = 4: [view,plane,radial]

VRP = 3: [view,radial,plane]

PET TOF parameters `parallelproj.tof`

PET time-of-flight (TOF) related classes and functions

class parallelproj.tof.TOFParameters(num_tofbins: int = 29, tofbin_width: float = 25.37139696, sigma_tof: float = 24.50529087048832, num_sigmas: float = 3.0, tofcenter_offset: float = 0)[source]

Bases: object

generic time of flight (TOF) parameters for a scanner with 385ps FWHM TOF

num_tofbins: int: number of time of flight bins
tofbin_width: float: width of the TOF bin in spatial units (mm)
sigma_tof: float: standard deviation of Gaussian TOF kernel in spatial units (mm)
num_sigmas: float: number of sigmas after which TOF kernel is truncated
tofcenter_offset: float: offset of center of central TOF bin from LOR center in spatial units (mm)

Parameters:

num_tofbins (int)
tofbin_width (float)
sigma_tof (float)
num_sigmas (float)
tofcenter_offset (float)

num_sigmas: float = 3.0

num_tofbins: int = 29

sigma_tof: float = 24.50529087048832

tofbin_width: float = 25.37139696

tofcenter_offset: float = 0

PET projectors `parallelproj.projectors`

high-level geometrical forward and back projectors

class parallelproj.projectors.EqualBlockPETProjector(lor_descriptor: EqualBlockPETLORDescriptor, img_shape: tuple[int, int, int], voxel_size: tuple[float, float, float], img_origin: None | Array = None)[source]

Bases: LinearOperator

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

Examples

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

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)

__call__(x: Array) → Array

alias to apply(x)

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

adjointness_test(xp: ModuleType, dev: str, verbose: bool = False, iscomplex: bool = False, **kwargs) → bool

test whether the adjoint is correctly implemented

Parameters:

xp (ModuleType) – array module to use
dev (str) – device (cpu or cuda)
verbose (bool, optional) – verbose output
iscomplex (bool, optional) – use complex arrays
**kwargs (dict) – passed to np.isclose

Returns:

whether the adjoint is correctly implemented

Return type:

apply(x: Array) → Array

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

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

property dev: str: device

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

property in_shape: tuple[int, int, int]: shape of the input array

property lor_descriptor: EqualBlockPETLORDescriptor: LOR descriptor

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

estimate norm of the linear operator using power iterations

Parameters:

xp (ModuleType) – array module to use
dev (str) – device (cpu or cuda)
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:

property out_shape: tuple[int, int, int]: shape of the output array

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

show_geometry(ax: Axes, 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 (plt.Axes) – 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

property tof: bool: bool indicating whether to use TOF or not

property tof_parameters: TOFParameters | None: TOF parameters

property voxel_size: Array: voxel size

property xp: ModuleType: array module

class parallelproj.projectors.ListmodePETProjector(event_start_coordinates: Array, event_end_coordinates: Array, img_shape: tuple[int, int, int], voxel_size: tuple[float, float, float], img_origin: None | Array = None)[source]

Bases: LinearOperator

non-TOF and TOF listmode projector for regular polygon PET scanners

Examples

PDHG and LM-SPHG 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)

__call__(x: Array) → Array

alias to apply(x)

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

adjointness_test(xp: ModuleType, dev: str, verbose: bool = False, iscomplex: bool = False, **kwargs) → bool

test whether the adjoint is correctly implemented

Parameters:

xp (ModuleType) – array module to use
dev (str) – device (cpu or cuda)
verbose (bool, optional) – verbose output
iscomplex (bool, optional) – use complex arrays
**kwargs (dict) – passed to np.isclose

Returns:

whether the adjoint is correctly implemented

Return type:

apply(x: Array) → Array

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

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

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: TOF bin of each event

property in_shape: tuple[int, int, int]: shape of the input array

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

estimate norm of the linear operator using power iterations

Parameters:

xp (ModuleType) – array module to use
dev (str) – device (cpu or cuda)
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:

property num_events: int: number of events

property out_shape: tuple[int]: shape of the output array

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

property tof: bool: bool indicating whether to use TOF projections or not

property tof_parameters: TOFParameters | None: TOF parameters

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

init method

Parameters:

image_shape (tuple[int, int]) – shape of the input image (n1, n2)
radial_positions (Array) – radial positions of the projection views in world coordinates
view_angles (Array) – angles of the projection views in radians
radius (float) – radius of the scanner
image_origin (tuple[float, float]) – world coordinates of the [0,0] voxel
voxel_size (tuple[float, float]) – the voxel size in both directions

__call__(x: Array) → Array

alias to apply(x)

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

adjointness_test(xp: ModuleType, dev: str, verbose: bool = False, iscomplex: bool = False, **kwargs) → bool

test whether the adjoint is correctly implemented

Parameters:

xp (ModuleType) – array module to use
dev (str) – device (cpu or cuda)
verbose (bool, optional) – verbose output
iscomplex (bool, optional) – use complex arrays
**kwargs (dict) – passed to np.isclose

Returns:

whether the adjoint is correctly implemented

Return type:

apply(x: Array) → Array

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

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

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

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

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

property in_shape: tuple[int, int]: shape of the input array

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

estimate norm of the linear operator using power iterations

Parameters:

xp (ModuleType) – array module to use
dev (str) – device (cpu or cuda)
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:

property num_rad: int: number of radial elements

property num_views: int: number of views

property out_shape: tuple[int, int]: shape of the output array

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

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

visualize the geometry of certrain projection views

Parameters:

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

Return type:

Figure

property voxel_size: Array: voxel size

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], ring_positions: Array, span: int = 1, max_ring_diff: int | None = None)[source]

Bases: LinearOperator

3D non-TOF parallel view projector

init method

Parameters:

image_shape (tuple[int, int, int]) – shape of the input image (n0, n1, n2) (last direction is axial)
radial_positions (Array) – radial positions of the projection views in world coordinates
view_angles (Array) – angles of the projection views in radians
radius (float) – radius of the scanner
image_origin (tuple[float, float, float]) – world coordinates of the [0,0,0] voxel
voxel_size (tuple[float, float, float]) – the voxel size in all directions (last direction is axial)
ring_positions (Array) – position of the rings in world coordinates
span (int) – span of the sinogram - default is 1
max_ring_diff (int | None) – maximum ring difference - default is None (no limit)

__call__(x: Array) → Array

alias to apply(x)

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

adjointness_test(xp: ModuleType, dev: str, verbose: bool = False, iscomplex: bool = False, **kwargs) → bool

test whether the adjoint is correctly implemented

Parameters:

xp (ModuleType) – array module to use
dev (str) – device (cpu or cuda)
verbose (bool, optional) – verbose output
iscomplex (bool, optional) – use complex arrays
**kwargs (dict) – passed to np.isclose

Returns:

whether the adjoint is correctly implemented

Return type:

apply(x: Array) → Array

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

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

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]: shape of the input array

property max_ring_diff: int: maximum ring difference

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

estimate norm of the linear operator using power iterations

Parameters:

xp (ModuleType) – array module to use
dev (str) – device (cpu or cuda)
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:

property out_shape: tuple[int, int, int]: shape of the output array

property scale: int | 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

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

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

__call__(x: Array) → Array

alias to apply(x)

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

adjointness_test(xp: ModuleType, dev: str, verbose: bool = False, iscomplex: bool = False, **kwargs) → bool

test whether the adjoint is correctly implemented

Parameters:

xp (ModuleType) – array module to use
dev (str) – device (cpu or cuda)
verbose (bool, optional) – verbose output
iscomplex (bool, optional) – use complex arrays
**kwargs (dict) – passed to np.isclose

Returns:

whether the adjoint is correctly implemented

Return type:

apply(x: Array) → Array

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

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

clear_cached_lor_endpoints() → None[source]

clear cached LOR endpoints

Return type:: None

convert_sinogram_to_listmode(sinogram: Array) → 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
Returns:: event_start_coordinates, event_end_coordinates, event_tofbin in case of non-TOF, event_tofbin is None
Return type:: tuple[Array, Array, Array | None]

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

property in_shape: tuple[int, int, int]: shape of the input array

property lor_descriptor: RegularPolygonPETLORDescriptor: LOR descriptor

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

estimate norm of the linear operator using power iterations

Parameters:

xp (ModuleType) – array module to use
dev (str) – device (cpu or cuda)
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:

property out_shape: tuple[int, int, int]: shape of the output array

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

show_geometry(ax: Axes, 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 (plt.Axes) – 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

property tof: bool: bool indicating whether to use TOF or not

property tof_parameters: TOFParameters | None: TOF parameters

property views: Array: view numbers to be projected

property voxel_size: Array: voxel size

property xend: Array | None: cached coordinates of LOR end points

property xp: ModuleType: array module

property xstart: Array | None: cached coordinates of LOR start points

Linear operators `parallelproj.operators`

basic linear operators

class parallelproj.operators.CompositeLinearOperator(operators: Sequence[LinearOperator])[source]

Bases: LinearOperator

Composite Linear Operator defined by a sequence of Linear Operators

Given a Sequence of operators

\[A^0, A^1, \ldots, A^{n-1}\]

the composite linear operator is defined as

\[A(x) = A^0( A^1( ... ( A^{n-1}(x) ) ) )\]

Examples

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

init method

Parameters:: operators (Sequence[LinearOperator, ...]) – Sequence of linear operators

__call__(x: Array) → Array

alias to apply(x)

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

__getitem__(i: int) → LinearOperator[source]

get the i-th operator \(A_i\)

Parameters:: i (int)
Return type:: LinearOperator

adjoint(y: Array) → Array

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

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

adjointness_test(xp: ModuleType, dev: str, verbose: bool = False, iscomplex: bool = False, **kwargs) → bool

test whether the adjoint is correctly implemented

Parameters:

xp (ModuleType) – array module to use
dev (str) – device (cpu or cuda)
verbose (bool, optional) – verbose output
iscomplex (bool, optional) – use complex arrays
**kwargs (dict) – passed to np.isclose

Returns:

whether the adjoint is correctly implemented

Return type:

apply(x: Array) → Array

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

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

property in_shape: tuple[int, ...]: shape of the input array

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

estimate norm of the linear operator using power iterations

Parameters:

xp (ModuleType) – array module to use
dev (str) – device (cpu or cuda)
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:

property operators: Sequence[LinearOperator]: tuple of linear operators

property out_shape: tuple[int, ...]: shape of the output array

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

class parallelproj.operators.ElementwiseMultiplicationOperator(values: Array)[source]

Bases: LinearOperator

Element-wise multiplication operator (multiplication with a diagonal matrix)

Examples

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

init method

Parameters:: values (Array) – values of the diagonal matrix

__call__(x: Array) → Array

alias to apply(x)

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

adjointness_test(xp: ModuleType, dev: str, verbose: bool = False, iscomplex: bool = False, **kwargs) → bool

test whether the adjoint is correctly implemented

Parameters:

xp (ModuleType) – array module to use
dev (str) – device (cpu or cuda)
verbose (bool, optional) – verbose output
iscomplex (bool, optional) – use complex arrays
**kwargs (dict) – passed to np.isclose

Returns:

whether the adjoint is correctly implemented

Return type:

apply(x: Array) → Array

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

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

property in_shape: tuple[int, ...]: shape of the input array

iscomplex() → bool[source]

bool whether the operator is complex

Return type:: bool

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

estimate norm of the linear operator using power iterations

Parameters:

xp (ModuleType) – array module to use
dev (str) – device (cpu or cuda)
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:

property out_shape: tuple[int, ...]: shape of the output array

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

property values: Array: values that get multiplied

property xp: ModuleType: array module of the operator

class parallelproj.operators.FiniteForwardDifference(in_shape: tuple[int, ...])[source]

Bases: LinearOperator

finite difference gradient operator

Parameters:: in_shape (tuple[int, ...])

__call__(x: Array) → Array

alias to apply(x)

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

adjointness_test(xp: ModuleType, dev: str, verbose: bool = False, iscomplex: bool = False, **kwargs) → bool

test whether the adjoint is correctly implemented

Parameters:

xp (ModuleType) – array module to use
dev (str) – device (cpu or cuda)
verbose (bool, optional) – verbose output
iscomplex (bool, optional) – use complex arrays
**kwargs (dict) – passed to np.isclose

Returns:

whether the adjoint is correctly implemented

Return type:

apply(x: Array) → Array

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

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

property in_shape: tuple[int, ...]: shape of the input array

property ndim: int: number of dimensions of the input array

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

estimate norm of the linear operator using power iterations

Parameters:

xp (ModuleType) – array module to use
dev (str) – device (cpu or cuda)
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:

property out_shape: tuple[int, ...]: shape of the output array

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

class parallelproj.operators.GaussianFilterOperator(in_shape: tuple[int, ...], sigma: float | Array, **kwargs)[source]

Bases: LinearOperator

Gaussian filter operator

Examples

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

init method

Parameters:

in_shape (tuple[int, ...]) – shape of the input array
sigma (float | array) – standard deviation of the gaussian filter
**kwargs (sometype) – passed to the ndimage gaussian_filter function

__call__(x: Array) → Array

alias to apply(x)

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

adjointness_test(xp: ModuleType, dev: str, verbose: bool = False, iscomplex: bool = False, **kwargs) → bool

test whether the adjoint is correctly implemented

Parameters:

xp (ModuleType) – array module to use
dev (str) – device (cpu or cuda)
verbose (bool, optional) – verbose output
iscomplex (bool, optional) – use complex arrays
**kwargs (dict) – passed to np.isclose

Returns:

whether the adjoint is correctly implemented

Return type:

apply(x: Array) → Array

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

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

property in_shape: tuple[int, ...]: shape of the input array

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

estimate norm of the linear operator using power iterations

Parameters:

xp (ModuleType) – array module to use
dev (str) – device (cpu or cuda)
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:

property out_shape: tuple[int, ...]: shape of the output array

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

class parallelproj.operators.LinearOperator[source]

Bases: ABC

abstract base class for linear operators

__call__(x: Array) → Array[source]

alias to apply(x)

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

adjoint(y: Array) → Array[source]

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

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

adjointness_test(xp: ModuleType, dev: str, verbose: bool = False, iscomplex: bool = False, **kwargs) → bool[source]

test whether the adjoint is correctly implemented

Parameters:

xp (ModuleType) – array module to use
dev (str) – device (cpu or cuda)
verbose (bool, optional) – verbose output
iscomplex (bool, optional) – use complex arrays
**kwargs (dict) – passed to np.isclose

Returns:

whether the adjoint is correctly implemented

Return type:

apply(x: Array) → Array[source]

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

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

abstract property in_shape: tuple[int, ...]: shape of the input array

norm(xp: ModuleType, dev: str, num_iter: int = 30, iscomplex: bool = False, verbose: bool = False) → float[source]

estimate norm of the linear operator using power iterations

Parameters:

xp (ModuleType) – array module to use
dev (str) – device (cpu or cuda)
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:

abstract property out_shape: tuple[int, ...]: shape of the output array

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

class parallelproj.operators.LinearOperatorSequence(operators: Sequence[LinearOperator])[source]

Bases: Sequence[LinearOperator]

Sequence of linear operators

\[A^0, A^1 \ldots, A^{n-1}\]

that can be evaluated independently.

Examples

Basic OSEM