Linear operators parallelproj.operators

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

Joint activity and attenuation reconstruction (MLAA) for TOF PET

init method

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

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)

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

Joint activity and attenuation reconstruction (MLAA) for TOF PET

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

__call__(x: Array) → Array¶

alias to apply(x)

Examples

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

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

get the i-th operator \(A_i\)

Examples

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

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)

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 in_shape: tuple[int, ...]¶: Input shape of the innermost (last) operator.

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

PDHG and SPDHG for PET reconstruction with a directional TV prior

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

property out_shape: tuple[int, ...]¶: Output shape of the outermost (first) operator.

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

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

init method

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

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)

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

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 in_shape: tuple[int, ...]¶: Shape of the diagonal values array (operator is square).

property iscomplex: bool¶: bool whether the operator is complex

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, ...]¶: Shape of the diagonal values array (operator is square).

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

Forward finite-difference gradient operator for 1-D to 4-D images.

Maps an image of shape in_shape to a gradient field of shape (ndim, *in_shape) where axis 0 enumerates the spatial directions. At each boundary the difference wraps to zero: the last difference along each axis is forced to zero in the forward pass (as if the image is padded with a copy of the border value, so the difference there vanishes). The adjoint is the negative discrete divergence, consistent with the standard TV regularisation convention. Self-adjointness is verified by adjointness_test().

Examples

Convergence comparison: SGD vs SVRG with logcosh regularization

PDHG and SPDHG for PET reconstruction with a directional TV prior

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

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

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

Joint activity and attenuation reconstruction (MLAA) for TOF PET

Parameters:: in_shape (tuple[int, ...]) – Shape of the input image. Must have 1 to 4 dimensions; raises ValueError for higher-dimensional inputs.

Examples

Convergence comparison: SGD vs SVRG with logcosh regularization

PDHG and SPDHG for PET reconstruction with a directional TV prior

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

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

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

Joint activity and attenuation reconstruction (MLAA) for TOF PET

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 in_shape: tuple[int, ...]¶: shape of the input array

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

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, ...]¶: shape of the output array

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

class parallelproj.operators.GaussianFilterOperator(in_shape: tuple[int, ...], **kwargs)[source]¶

Bases: LinearOperator

Isotropic Gaussian smoothing operator (self-adjoint).

Works with NumPy, CuPy, and PyTorch (CPU and CUDA) arrays. GPU arrays (CuPy, and PyTorch CUDA round-tripped through CuPy via DLPack) are filtered with cupyx.scipy.ndimage.gaussian_filter; CPU arrays use scipy.ndimage.gaussian_filter. Routing the GPU path through cupyx makes it independent of scipy’s array-API delegation, so it works regardless of the SCIPY_ARRAY_API env var or import order. All keyword arguments accepted by gaussian_filter (e.g. sigma, mode, truncate) are forwarded through **kwargs. Because the Gaussian kernel is symmetric, the adjoint equals the forward application.

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

Exact vs. “safe epsilon” mode of the negative Poisson log-likelihood

Exact vs. "safe epsilon" mode of the negative Poisson log-likelihood

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

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

Joint activity and attenuation reconstruction (MLAA) for TOF PET

init method

Parameters:

in_shape (tuple[int, ...]) – shape of the input array
**kwargs (dict) – passed to scipy.ndimage.gaussian_filter; most commonly sigma (standard deviation in pixels), plus optional mode, truncate, etc.

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

Exact vs. “safe epsilon” mode of the negative Poisson log-likelihood

Exact vs. "safe epsilon" mode of the negative Poisson log-likelihood

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

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

Joint activity and attenuation reconstruction (MLAA) for TOF PET

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

Exact vs. “safe epsilon” mode of the negative Poisson log-likelihood

Exact vs. "safe epsilon" mode of the negative Poisson log-likelihood

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 in_shape: tuple[int, ...]¶: Shape of the input (and output) array (operator is square).

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, ...]¶: Shape of the output (and input) array (operator is square).

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

class parallelproj.operators.GradientFieldProjectionOperator(gradient_field: Array, eta: float = 0.0)[source]¶

Bases: LinearOperator

Gradient field projection operator (self-adjoint).

Projects a gradient field onto the subspace orthogonal to a normalised structural prior gradient \(\xi_n\):

\[P_{\xi_n}x = x - \langle \xi_n, x \rangle \xi_n, \qquad \xi_n = g_n / \| g_n \|_{\eta}\]

where \(g_n\) is the joint gradient field and \(\eta\) is a smoothing parameter for the pointwise gradient norm. The operator is self-adjoint (its own adjoint) because orthogonal projection operators are symmetric.

See Ehrhardt and Betcke, “Multicontrast MRI Reconstruction with Structure-Guided Total Variation” (doi: 10.1137/15M1047325).

Examples

PDHG and SPDHG for PET reconstruction with a directional TV prior

Parameters:

gradient_field (Array) – a real gradient field. In 3D, the shape would be [3,n0,n1,n2]. In 2D, the shape would be [2,n0,n1]. This can be e.g. the output of the FiniteForwardDifference operator applied to a structural prior image.
eta (float, optional) – smoothing parameter used in the pointwise gradient norm default 0.0

Examples

PDHG and SPDHG for PET reconstruction with a directional TV prior

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 dev: str¶: device of the operator

property eta: float¶: smoothing parameter

property in_shape: tuple[int, ...]¶: Shape of the gradient field (ndim, *spatial_shape).

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 normalized_gradient_field: Array¶: normalized gradient field

property out_shape: tuple[int, ...]¶: Shape of the gradient field (ndim, *spatial_shape) (same as in_shape).

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

property xp: ModuleType¶: array module of the operator

class parallelproj.operators.LinearOperator[source]¶

Bases: ABC

Abstract base class for array-API-compatible linear operators.

Subclasses implement _apply() (\(y = Ax\)) and _adjoint() (\(x = A^H y\)). The public apply() and adjoint() methods apply an optional scalar scale factor (\(\alpha A\) and \(\overline{\alpha} A^H\)). Utility methods adjointness_test() and norm() are provided for validation and step-size estimation.

Examples

Michelograms and axial sinogram compression

Non-TOF and TOF projections using a modularized (block) PET scanner geometry

Listmode to sinogram unlisting

Convergence comparison: SGD vs SVRG with logcosh regularization

PDHG and SPDHG for PET reconstruction with a directional TV prior

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

DePierro’s algorithm to optimize the Poisson logL with quadratic intensity prior

DePierro's algorithm to optimize the Poisson logL with quadratic intensity prior

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

Exact vs. “safe epsilon” mode of the negative Poisson log-likelihood

Exact vs. "safe epsilon" mode of the negative Poisson log-likelihood

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

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

pytorch parallelproj projection layer

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

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

alias to apply(x)

Examples

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

Examples

Michelograms and axial sinogram compression

Non-TOF and TOF projections using a modularized (block) PET scanner geometry

Convergence comparison: SGD vs SVRG with logcosh regularization

PDHG and SPDHG for PET reconstruction with a directional TV prior

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

DePierro’s algorithm to optimize the Poisson logL with quadratic intensity prior

DePierro's algorithm to optimize the Poisson logL with quadratic intensity prior

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

Exact vs. “safe epsilon” mode of the negative Poisson log-likelihood

Exact vs. "safe epsilon" mode of the negative Poisson log-likelihood

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

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[source]¶

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

Non-TOF and TOF projections using a modularized (block) PET scanner geometry

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

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

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

Examples

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

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

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

PDHG and SPDHG for PET reconstruction with a directional TV prior

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

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

Convergence comparison: SGD vs SVRG with logcosh regularization

PDHG and SPDHG for PET reconstruction with a directional TV prior

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

init method

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

Examples

Convergence comparison: SGD vs SVRG with logcosh regularization

PDHG and SPDHG for PET reconstruction with a directional TV prior