PET LOR / sinogram descriptors parallelproj.pet_lors - parallelproj 2.0.0dev0 documentation

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

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

Examples

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

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]

Examples

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

Examples

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

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

class parallelproj.pet_lors.Michelogram(num_rings: int, max_ring_difference: int, span: int = 1)[source]¶

Bases: object

Axial plane layout for a cylindrical PET scanner under odd span.

Encapsulates the segment / axial-position combinatorics that map every valid ring pair \((s, e)\) onto a sinogram plane index under span conventions.

For span \(S\) (odd) and a maximum ring difference \(D\), each ring pair with \(|e - s| \le D\) is assigned a segment via ring_diff_to_segment(). Ring pairs sharing the same \((\text{segment},\; s + e)\) collapse into a single plane. Planes are ordered by \((|\text{seg}|,\; -\text{seg},\; s + e)\) \([0, +1, -1, +2, -2, \ldots]\) with axial bins increasing in \(s + e\) (equivalently in z for equispaced rings).

The class knows nothing about ring z-positions, scanner radius, or sinogram axis ordering — it operates on pure integer indices. Consumers (e.g. RegularPolygonPETLORDescriptor, SinogramAxialCompressionOperator) combine it with the geometry- and array-API-specific information they need.

For span = 1 the layout reduces to the unspanned Michelogram (each ring pair is its own plane with max_multiplicity == 1); the ordering rd = 0, +1, -1, +2, -2, ... with each ring difference sorted by ring sum.

Parameters:

num_rings (int) – Number of detector rings (\(R \ge 1\)).
max_ring_difference (int) – Maximum ring difference \(|e - s|\) considered (\(\ge 0\)). Values larger than num_rings - 1 have no extra effect.
span (int, optional) – Axial compression factor — must be odd and \(\ge 1\). Default 1 (no compression).

Examples

>>> m = Michelogram(num_rings=3, max_ring_difference=2, span=3)
>>> int(m.num_planes)
7
>>> int(m.max_multiplicity)
2
>>> int(m.ring_diff_to_segment(0)), int(m.ring_diff_to_segment(2)),         int(m.ring_diff_to_segment(-2))
(0, 1, -1)

Examples

Convergence comparison: SGD vs SVRG with quadratic regularization

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

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

average_z_per_plane(ring_positions) → tuple[ndarray, ndarray][source]¶

Mean ring z-coordinate per plane, separately for start and end rings.

Equivalent to averaging ring_positions over the contributing ring pairs of each plane. For span =1 planes this is trivially the single contributing ring’s z; for span > 1 planes it produces the averaged-LOR z-position used by the spanned setup of RegularPolygonPETLORDescriptor and by show_segment_lors().

Parameters:

ring_positions (array-like, shape (num_rings,)) – z-coordinate of each ring (any backend; converted via np.asarray).

Returns:

start_z (np.ndarray, shape (num_planes,), dtype float32)
end_z (np.ndarray, shape (num_planes,), dtype float32)

Return type:

tuple[ndarray, ndarray]

Examples

compression_index_maps(target_span: int) → tuple[ndarray, ndarray, ndarray, ndarray][source]¶

Build gather/scatter index maps to a higher-span Michelogram.

Convenience wrapper around compression_index_maps_to() that builds the target Michelogram internally as Michelogram(self.num_rings, self.max_ring_difference, span=target_span).

Parameters:: target_span (int) – Odd integer >= self.span; additionally (target_span // self.span) must be odd.
Return type:: See compression_index_maps_to().
Raises:: ValueError – If target_span is not a positive odd integer. Further validation errors are raised by compression_index_maps_to().

Examples

compression_index_maps_to(target: Michelogram) → tuple[ndarray, ndarray, ndarray, ndarray][source]¶

Build gather/scatter index maps to a pre-built target Michelogram.

Returns the integer index structures needed to map planes of this Michelogram onto planes of target.

Both Michelograms must describe the same scanner geometry (target.num_rings == self.num_rings), and target.span must be an integer multiple of self.span. Because both spans are odd by construction, the ratio target.span / self.span is then automatically odd, which guarantees that every ring pair of any input plane shares the same target plane — so the operation is a single-valued gather.

The target’s max_ring_difference must be at least self.max_ring_difference so every input ring pair has a target plane. If it is strictly greater, the resulting maps still work but some output planes will have zero multiplicity (output bins that no input ring pair contributes to).

Parameters:

target (Michelogram) – Pre-built target Michelogram. Validation rules above.

Returns:

target_for_p1 (np.ndarray, shape (self.num_planes,), dtype int64) – For each plane of this Michelogram, the corresponding plane index in target.
idx2d (np.ndarray, shape (target.num_planes, target_max_mult), int64) – For each target plane, the indices in this Michelogram that contribute, right-padded with 0. Use mask to filter.
mask (np.ndarray, same shape as idx2d, dtype float32) – 1.0 for valid entries, 0.0 for right-padding.
target_multiplicity (np.ndarray, shape (target.num_planes,), int32) – Number of self-planes folded into each target plane.

Raises:

TypeError – If target is not a Michelogram instance.
ValueError – If target.num_rings differs from self.num_rings; if target.span < self.span; if self.span does not divide target.span; or if target.max_ring_difference < self.max_ring_difference.

Return type:

tuple[ndarray, ndarray, ndarray, ndarray]

Examples

property max_multiplicity: int¶: Largest plane multiplicity (most ring pairs in any one plane).

property max_ring_difference: int¶: Maximum ring difference \(|e - s|\).

property num_planes: int¶: Total number of sinogram planes.

property num_rings: int¶: Number of rings.

property plane_axial_midpoint_int: ndarray¶: Integer axial midpoint \(s + e\) (= twice the actual midpoint) for each plane, shape (num_planes,), dtype int32.

property plane_end_rings: ndarray¶

Contributing end ring indices per plane, right-padded with 0.

Shape (num_planes, max_multiplicity), dtype int32. Use plane_mask to identify the valid entries.

plane_for_ring_pair(s: int, e: int) → int[source]¶

Plane index for the ring pair (s, e).

Raises:

IndexError – If either s or e is outside [0, num_rings).
ValueError – If \(|e - s| > \text{max\_ring\_difference}\).

Parameters:

s (int)
e (int)

Return type:

int

Examples

property plane_for_ring_pair_table: ndarray¶: (num_rings, num_rings) lookup table whose entry [s, e] is the plane index for ring pair (s, e), or -1 if \(|e - s| > \text{max\_ring\_difference}\).

property plane_mask: ndarray¶

Validity mask for plane_start_rings / plane_end_rings.

Shape (num_planes, max_multiplicity), dtype float32. Entries are 1.0 for valid contributing ring pairs and 0.0 for right-padding.

property plane_multiplicity: ndarray¶: Number of ring pairs contributing to each plane, shape (num_planes,), dtype int32.

property plane_segment: ndarray¶: Signed segment number for each plane, shape (num_planes,), dtype int32.

property plane_start_rings: ndarray¶

Contributing start ring indices per plane, right-padded with 0.

Shape (num_planes, max_multiplicity), dtype int32. Use plane_mask to identify the valid entries.

ring_diff_to_segment(rd: int) → int[source]¶

Signed segment number for a given ring difference \(e - s\).

Returns:: 0 if \(|rd| \le \text{half\_span}\), otherwise \(\pm k\) with \(k = \lceil (|rd| - \text{half\_span}) / S \rceil\) and sign equal to that of \(rd\).
Return type:: int
Parameters:: rd (int)

Examples

show(ax: Axes, show_merge_lines: bool = True, plane_index_fontsize: float = 6, **kwargs) → None[source]¶

Draw the Michelogram scatter plot onto ax.

Each point represents a valid ring pair (s, e), colored by abs(segment). For span > 1, ring pairs that share the same (segment, s + e) and therefore collapse into the same sinogram plane are connected by a thin grey line when show_merge_lines is True.

Parameters:

ax (plt.Axes) – 2-D matplotlib axes (not 3-D).
show_merge_lines (bool, optional) – Draw lines connecting ring pairs that merge into the same plane. Defaults to True. Only has a visible effect for span > 1.
plane_index_fontsize (float, optional) – Font size of the per-plane index annotations placed at each ring-pair (or merged-group) centroid. Defaults to 6. Useful knob when the Michelogram is large (lower to avoid overlap) or small (raise for readability).
**kwargs – Forwarded to ax.scatter (e.g. s=4, cmap="RdBu_r").

Return type:

None

Examples

show_segment_lors(ring_positions, axs=None, uncompressed_lor_kwargs: dict | None = None, compressed_lor_kwargs: dict | None = None, inset_plane_index_fontsize: float = 4)[source]¶

Side-view LOR diagram per segment with a Michelogram inset.

Mirrors the descriptor’s RegularPolygonPETLORDescriptor.show_segment_lors(), but takes ring_positions explicitly so the Michelogram can be visualised standalone (e.g. with np.arange(num_rings) for a purely schematic plot, or with the user’s actual ring z-positions).

Subplots are arranged in a 2-row grid (when negative segments exist):

columns indexed by abs(segment): 0, 1, 2, …
row 0 non-negative segments (0, +1, +2, …)
row 1, col 0 Michelogram inset
row 1, col >= 1 negative segments (-1, -2, …)

Each LOR subplot shows the uncompressed (per-ring-pair) LORs as solid black lines and the compressed (axially-averaged) LORs as dashed coloured lines.

Parameters:

ring_positions (array-like, shape (num_rings,)) – z-coordinate of each ring.
axs (2-D array-like of Axes, optional) – Pre-existing axes of shape (n_rows, n_cols). If None, a new figure is created.
uncompressed_lor_kwargs (dict, optional) – Style overrides for the uncompressed LOR lines.
compressed_lor_kwargs (dict, optional) – Style overrides for the compressed LOR lines.
inset_plane_index_fontsize (float)

Return type:

matplotlib.figure.Figure

Examples

property span: int¶: Axial compression factor (odd).

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

Examples

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

Convergence comparison: SGD vs SVRG with quadratic regularization

PDHG and SPDHG for PET reconstruction with a directional TV prior

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

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

Examples

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

Convergence comparison: SGD vs SVRG with quadratic regularization

PDHG and SPDHG for PET reconstruction with a directional TV prior

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

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

abstractmethod get_lor_coordinates() → tuple[Array, Array][source]¶

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

Examples

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, michelogram: Michelogram | None = None, radial_trim: int = 3, 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

Convergence comparison: SGD vs SVRG with quadratic 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-SPHG to optimize the Poisson logL and total variation

Parameters:

scanner (RegularPolygonPETScannerGeometry) – a regular polygon PET scanner.
michelogram (Michelogram, optional) – the axial plane layout — the single source of truth for the spanning combinatorics (segments, axial midpoints, ring-pair grouping, ordering). If None (default), a span-1 layout with no constraint on the ring difference is used, i.e. Michelogram(scanner.num_rings, scanner.num_rings - 1, span=1). The Michelogram must have num_rings == scanner.num_rings.
radial_trim (int, optional) – number of geometrial LORs to disregard in the radial direction. Defaults to 3.
sinogram_order (SinogramSpatialAxisOrder, optional) – the order of the sinogram axes. Defaults to SinogramSpatialAxisOrder.RVP.

Examples

Convergence comparison: SGD vs SVRG with quadratic 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-SPHG to optimize the Poisson logL and total variation

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 ring index for all planes (only defined for span=1)

property end_plane_z: Array¶: end z-coordinate for all planes (averaged over constituent ring pairs for span > 1)

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, …]]]

Examples

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

Examples

property max_ring_difference: int¶: the maximum ring difference

property michelogram: Michelogram¶

The Michelogram describing the axial plane layout.

This is the single source of truth for the spanning combinatorics (segments, axial midpoints, ring-pair grouping, ordering). Useful for visualization, axial compression operators, or any user code that needs access to the integer ring-pair structure.

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 plane_multiplicity: Array¶: number of ring pairs contributing to each plane (always 1 for span=1)

property plane_segment: Array¶: segment number for each plane (equals abs(rd) for span=1)

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: RegularPolygonPETScannerGeometry¶: the scanner for which coincidences are described

show_michelogram(ax: Axes, show_merge_lines: bool = True, **kwargs) → None[source]¶

Visualize the Michelogram.

Thin wrapper around Michelogram.show(); see that method for full documentation of arguments.

Examples

Parameters:

ax (Axes)
show_merge_lines (bool)

Return type:

None

show_segment_lors(axs=None, uncompressed_lor_kwargs: dict | None = None, compressed_lor_kwargs: dict | None = None)[source]¶

Side-view LOR diagram per segment with the Michelogram inset.

Thin wrapper around Michelogram.show_segment_lors(); this method supplies the scanner’s ring positions automatically. See Michelogram.show_segment_lors() for full documentation.

Examples

Parameters:

uncompressed_lor_kwargs (dict | None)
compressed_lor_kwargs (dict | None)

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

Examples

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

property span: int¶: axial compression factor (1 = no compression)

property spatial_sinogram_shape: tuple[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 ring index for all planes (only defined for span=1)

property start_plane_z: Array¶: start z-coordinate for all planes (averaged over constituent ring pairs for span > 1)

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.SinogramAxialCompressionOperator(lor_descriptor: RegularPolygonPETLORDescriptor, target_span: int, mode: str = 'sum', num_tof_bins: int | None = None)[source]¶

Bases: LinearOperator

Linear operator that axially compresses a span-1 PET sinogram into a higher odd span.

For an input RegularPolygonPETLORDescriptor with span=1 and a target odd span \(S\), every span-1 ring pair \((s, e)\) is assigned to an output bin \((\text{segment}, \text{axial midpoint})\) where

segment is determined by the ring difference \(rd = e - s\) under target span \(S\) RegularPolygonPETLORDescriptor._ring_diff_to_segment()),
axial midpoint is \(s + e\) (an integer equal to twice the actual midpoint).

All span-1 ring pairs sharing the same \((\text{segment}, s + e)\) collapse into a single output plane.

Two reduction modes are supported:

mode="sum" (default). The output plane is the sum of the contributing input planes:

\[y_n \;=\; \sum_{p_1 \in \mathcal{G}(n)} x_{p_1} \qquad \left(G^T y\right)_{p_1} \;=\; y_{\,\tau(p_1)}\,.\]

This is the natural reduction for counts-like sinograms — emission data, measured counts, randoms, etc. — which add when ring pairs are grouped together.
mode="average". The output plane is the mean of the contributing input planes:

\[y_n \;=\; \frac{1}{m_n} \sum_{p_1 \in \mathcal{G}(n)} x_{p_1} \qquad \left(G_{\rm avg}^T y\right)_{p_1} \;=\; \frac{y_{\,\tau(p_1)}}{m_{\,\tau(p_1)}}\,.\]

This is the natural reduction for multiplicative-factor sinograms — attenuation factors, sensitivity / normalisation factors, geometric efficiency — which should average rather than sum when ring pairs are grouped together.

In both expressions, \(\mathcal{G}(n)\) is the set of input plane indices mapped to output plane \(n\), \(m_n = |\mathcal{G}(n)|\) is the plane multiplicity, and \(\tau(p_1)\) is the output plane index for input plane \(p_1\).

Output plane ordering matches that of RegularPolygonPETLORDescriptor constructed with the same scanner, radial_trim, max_ring_difference, and sinogram_order but with span=target_span. That companion descriptor is exposed as out_lor_descriptor for visualization (e.g. show_michelogram, show_segment_lors) or for composing the operator with a span-\(S\) projector.

The closed-form operator 2-norms are

\[\|G_{\rm sum}\|_2 = \sqrt{\max_n m_n}\,, \qquad \|G_{\rm avg}\|_2 = 1 / \sqrt{\min_n m_n}\,,\]

derived from \(G_{\rm sum} G_{\rm sum}^T = \operatorname{diag}(m_n)\) and \(G_{\rm avg} G_{\rm avg}^T = \operatorname{diag}(1/m_n)\). norm() returns these directly without power iteration.

Parameters:

lor_descriptor (RegularPolygonPETLORDescriptor) – A span=1 LOR descriptor whose sinogram is to be compressed.
target_span (int) – Odd integer >= 1 giving the target axial compression. 1 is accepted and yields an identity-like operator (each input plane maps to a single output plane in the same span-1 order).
mode ({"sum", "average"}, optional) – Reduction mode. "sum" (default) is appropriate for counts-like sinograms; "average" is appropriate for multiplicative-factor sinograms such as attenuation or sensitivity factors.
num_tof_bins (int or None, optional) – If None (default), the operator acts on the 3D spatial sinogram with shape RegularPolygonPETLORDescriptor.spatial_sinogram_shape. If a positive integer, the operator acts on a 4D TOF sinogram whose trailing axis (size num_tof_bins) is the TOF axis and is passed through unchanged.

Examples

>>> import array_api_compat.numpy as xp
>>> import parallelproj.pet_scanners as pps
>>> import parallelproj.pet_lors as ppl
>>> scanner = pps.RegularPolygonPETScannerGeometry(
...     xp, "cpu", radius=65.0, num_sides=12, num_lor_endpoints_per_side=4,
...     lor_spacing=4.0, ring_positions=xp.asarray([0.0, 1.0, 2.0]),
...     symmetry_axis=2,
... )
>>> lor_s1 = ppl.RegularPolygonPETLORDescriptor(
...     scanner, ppl.Michelogram(scanner.num_rings, 2, span=1), radial_trim=1,
... )
>>> comp = ppl.SinogramAxialCompressionOperator(lor_s1, target_span=3)
>>> comp.in_shape, comp.out_shape
((..., ..., 9), (..., ..., 7))
>>> comp.adjointness_test(xp, "cpu")
True

Examples