#### NOTE
New to DeepInverse? Get started with the basics with the
[5 minute quickstart tutorial](https://deepinv.org/auto_examples/basics/demo_quickstart.html.md#sphx-glr-auto-examples-basics-demo-quickstart-py)..
# Positron emission tomography (PET) in 3D
This demo shows how to define a non time-of-flight PET scanner, simulate measurements
and reconstruct a volume from them.
The (unnormalized) PET forward model is defined as
$$
y \sim \gamma \mathcal{P}(c \circ H(g*x) + b)
$$
where $H \in \mathbb{R}_{+}^{m \times n}$ is the projection operator,
$g \in \mathbb{R}_{+}^{n}$ is a Gaussian blur kernel, $x\in\mathbb{R}_{+}^{n}$
is the emission image, $b \in \mathbb{R}_{+}^{m}$ is the (expected) background,
$\mathcal{P}$ denotes Poisson noise,
$c=\exp(-H\mu)\in \mathbb{R}_{+}^{m}$ is an (optional) attenuation term
with $\mu \in \mathbb{R}_{+}^{n}$ an attenuation map (typically obtained through an auxiliary CT scan).
#### NOTE
In this example, we consider the unnormalized case, which allows to obtain quantitative reconstructions (i.e., $x$ has real
physical units). The operator also can be used in a normalized setting (forcing $\|A\|_2=1$ and normalizing counts to be between 0 and 1).
See also the [normalized 2D PET example](https://deepinv.org/auto_examples/physics/demo_pet2d.html.md#sphx-glr-auto-examples-physics-demo-pet2d-py).
When using deep learning-based reconstruction methods, it is often easier to consider the normalized case, but a special attention is required
to denormalize the reconstructions and obtain physical units.
#### TIP
If you prefer to get started with PET on a simpler 2D problem, please check out [the 2D PET demo](https://deepinv.org/auto_examples/physics/demo_pet2d.html.md#sphx-glr-auto-examples-physics-demo-pet2d-py).
#### NOTE
This operator requires the `parallelproj` package to be installed.
This in turn requires [installing deepinv via pixi or conda](https://deepinv.org/index.html.md#install),
but not pypi/uv (as `parallelproj` is not currently available on pypi).
If you are working on a conda environment, you can install `parallelproj` as
```default
conda install -c conda-forge parallelproj
```
If you are working on a pixi installation, simply do
```default
pixi install -e full
```
which installs all optional dependencies.
Check the `parallelproj` documentation for more details: [https://parallelproj.readthedocs.io/en/stable/](https://parallelproj.readthedocs.io/en/stable/).
```Python
import time
import matplotlib.pyplot as plt
import deepinv as dinv
from deepinv.physics import PET
from deepinv.utils.phantoms import generate_pet_phantom
import torch
import parallelproj
from array_api_compat import torch as torch_compat
```
## Setup a minimal non-TOF PET projector
Here we define each voxel to have size $3\times 3\times 3$ mm
such that the total volume to reconstruct is of size $38.4\times 38.4\times 7.2$ cm
which fits approximately a portion of a human chest.
The maximum achievable resolution (in high count settings) is typically proportional to the full-width at half
maximum (FWHM) of the Gaussian blur kernel, which here is set to 4 mm.
We use a PET scanner with 8 rings of detectors, each ring being a polygon of
32 sides, and each side containing 16 detectors. This gives us a total of 32\*16=512 detectors per ring.
#### TIP
You can play with different geometries and voxel sizes to get a good grasp of
the scanner geometry.
```Python
device = "cuda" if torch.cuda.is_available() else "cpu"
# DeepInverse uses the (D, H, W) convention, as opposed to parallelproj which uses (H, W, D)
img_size = (24, 128, 128)
voxel_size = (3, 3, 3)
# number of sides of the polygon approximating a circle
num_sides = 32
# number of detectors per polygon side
num_lor_endpoints_per_side = 16
# number of rings of detectors on the depth axes
num_rings = 8
scanner = parallelproj.pet_scanners.DemoPETScannerGeometry(
torch_compat,
dev=device,
num_rings=num_rings,
num_sides=num_sides,
num_lor_endpoints_per_side=num_lor_endpoints_per_side,
)
# FWHM of the Gaussian blur kernel in mm
fwhm_data_mm = 4
# gain of the device:
# higher gains are associated to lower dose and/or shorter acquisition times,
# while lower gains are associated to higher dose and/or longer acquisition times.
# larger gain -> more poisson noise -> harder reconstruction
gain = 0.001
physics = PET(
device=device,
voxel_size=voxel_size,
scanner=scanner,
fwhm_data_mm=fwhm_data_mm,
img_size=img_size,
normalize_counts=True,
normalize=True,
gain=gain,
)
physics.plot_geometry()
```
## Define a phantom and attenuation map
We define a 3D phantom and attenuation map, whose shape is the same as the phantom.
In practice, the attenuation is typically obtained with an auxiliary CT scan of the patient.
```Python
x, attenuation = generate_pet_phantom(img_size, device=device)
mid_slice = img_size[0] // 2
dinv.utils.plot(
[x[:, :, mid_slice], attenuation[:, :, mid_slice]],
titles=["Emission image", "Attenuation image"],
)
```
## Simulating measurements
The shape of measurements is approximately `(B, 1, N, N/2, R^2)` where
`N=num_lor_endpoints_per_side*num_sides` is the number of detectors per ring
and `R` is the number of rings.
This provides one measurement for every possible Line of Response (LOR), or in other words ‘rays’, connecting
two detectors in the scanner, which are arranged in a sinogram format, with the first axis
corresponding to the angle of the ray, the second axis corresponding to the distance of the ray to the center of the field of view
and the last axis corresponding to the depth of the ray (i.e., which rings of detectors are connected by the ray)
#### TIP
The size of measurements is independent of the chosen `img_size`
```Python
y = physics(x)
sinogram_mid_slice = y.shape[-1] // 2
print(
f"Measurements shape={tuple(y.shape)}, range=({y.min().item():.2f},{y.max().item():.2f})"
)
```
```none
Measurements shape=(1, 1, 507, 256, 64), range=(0.00,0.29)
```
## Setting up background and attenuation
The attenuation term reduces the amount of signal measured in rays that
go through highly attenuating regions, such as bones. This makes the reconstruction more challenging, but also more realistic.
In PET, we generally have access to a realization of the background,
i.e., $\tilde{s} \sim \mathcal{P}(s)$, which is a Poisson random variable with mean $s$.
Both attenuation and background are stored as “physics parameters” which are patient dependent
and can be updated via [`physics.update(...)`](https://deepinv.org/api/stubs/deepinv.physics.Physics.html.md#deepinv.physics.Physics.update) or by passing them as kwargs in
[`physics(x, ...)`](https://deepinv.org/api/stubs/deepinv.physics.Physics.html.md#deepinv.physics.Physics.forward), [`physics.A(x, ...)`](https://deepinv.org/api/stubs/deepinv.physics.Physics.html.md#deepinv.physics.Physics.A) or
[`physics.A_adjoint(y, ...)`](https://deepinv.org/api/stubs/deepinv.physics.LinearPhysics.html.md#deepinv.physics.LinearPhysics.A_adjoint).
#### NOTE
The attenuation is stored in the physics in sinogram space as $\exp(-\mu)$ to speed up computations,
but it can be provided either in image space, i.e., $\mu$, to the physics, or in sinogram space, i.e., $\exp(-\mu)$.
The class figures out the attenuation space by comparing it to `img_size`.
```Python
expected_background = torch.ones_like(y) * x.max() * 0.05
background = physics.generate_background(expected_background)
physics.update(attenuation=attenuation, background=background)
y = physics(x)
y2 = y - background
dinv.utils.plot(
[
physics.attenuation[..., sinogram_mid_slice],
y[..., sinogram_mid_slice],
y2[..., sinogram_mid_slice],
],
["sino. atten.", "meas.", "corrected meas."],
figsize=(6, 6),
)
```
```none
/local/jtachell/deepinv/deepinv/deepinv/physics/functional/matrix.py:42: UserWarning: Power iteration: convergence not reached
warnings.warn("Power iteration: convergence not reached")
```
## Backprojection and sensitivities
We backproject the data to visualize the sensitivity map of the scanner.
The sensitivity map is defined as the back-projection of a sinogram of ones $s = A^\top \mathbf{1}$, which corresponds to the number of rays intersecting each voxel.
Here we also obtain a simple linear least-squares reconstruction by using
[`A_dagger`](https://deepinv.org/api/stubs/deepinv.physics.LinearPhysics.html.md#deepinv.physics.LinearPhysics.A_dagger).
```Python
with torch.no_grad():
x_dag = physics.A_dagger(y - background)
sensitivities = physics.A_adjoint(torch.ones_like(y))
print(f"Norm operator: {physics.compute_norm(x):.2f}")
dinv.utils.plot(sensitivities[:, :, mid_slice], ["sensitivities"])
```
```none
/local/jtachell/deepinv/deepinv/deepinv/physics/forward.py:639: DeprecationWarning: Using `compute_norm(squared=True)` is deprecated. Use `compute_sqnorm()` instead to compute the squared spectral norm (||A^T A||_2). In a future version, `compute_norm()` will compute the non-squared spectral norm (||A||_2) by default.
warnings.warn(
Power iteration converged at iteration 8, ||A^T A||_2=1.00
Norm operator: 1.00
```
## MLEM reconstruction
We run the standard MLEM reconstruction algorithm [1](#footcite-sheppmaximumlikelihoodreconstruction1982)
to obtain a reconstructed emission volume.
The algorithm can be seen as a preconditioned gradient descent on the negative log-likelihood of the Poisson model:
$$
x^{(k+1)} = x^{(k)} - P \nabla f(Ax^{(k)}+b,y)
$$
where $f$ is the Poisson data-fidelity term, $P=\mathrm{diag}(\frac{x}{A^T\mathbf{1}})$ is a preconditioner
and $b$ is the background.
```Python
gain = physics.noise_model.gain
mlem_iter = 40
def _sync():
if torch.device(device).type == "cuda":
torch.cuda.synchronize()
nrmse = dinv.metric.NRMSE()
# With ``denormalize=True``, the likelihood is evaluated in the count domain.
# The background produced by PET is in the normalized measurement domain, so it
# must be divided by the gain as well.
data_fidelity = dinv.optim.PoissonLikelihood(
gain=gain,
bkg=background / gain,
denormalize=True,
)
model_mlem = dinv.optim.MLEM(
data_fidelity=data_fidelity,
prior=None,
max_iter=mlem_iter,
custom_metrics={
"nrmse": lambda _values, _x_prev, x_cur: nrmse(x_cur.unsqueeze(0), x).item()
},
)
with torch.no_grad():
_sync()
start = time.perf_counter()
x_mlem, metrics_mlem = model_mlem(
y, physics, init=torch.ones_like(x), compute_metrics=True
)
_sync()
mlem_time = time.perf_counter() - start
print(f"MLEM runtime: {mlem_time:.2f} s for {mlem_iter} iterations")
psnr = dinv.metric.PSNR(max_pixel=x.max().item())
psnr_mlem = psnr(x_mlem, x)
psnr_dag = psnr(x_dag, x)
nrmse_mlem = nrmse(x_mlem, x)
nrmse_dag = nrmse(x_dag, x)
dinv.utils.plot(
[x[:, :, mid_slice], x_mlem[:, :, mid_slice], x_dag[:, :, mid_slice]],
["Ground truth", f"MLEM ({mlem_iter} it.)", "L2 pseudoinv."],
subtitles=[
"Reference",
f"PSNR: {psnr_mlem.item():.2f} dB\n" f"NRMSE: {100 * nrmse_mlem.item():.2f}%",
f"PSNR: {psnr_dag.item():.2f} dB\n" f"NRMSE: {100 * nrmse_dag.item():.2f}%",
],
rescale_mode="clip",
vmin=0,
vmax=x.max().item(),
figsize=(8, 4),
cbar=True,
)
```
```none
MLEM runtime: 0.57 s for 40 iterations
/local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:550: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator. Otherwise, ticks may be mislabeled.
colbar.set_ticklabels(true_labels)
/local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:550: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator. Otherwise, ticks may be mislabeled.
colbar.set_ticklabels(true_labels)
/local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:550: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator. Otherwise, ticks may be mislabeled.
colbar.set_ticklabels(true_labels)
/local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:570: UserWarning: This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust.
plt.subplots_adjust(hspace=0.2, wspace=0.2)
```
## Accelerating MLEM with ordered subsets
One MLEM iteration evaluates the forward and adjoint operators using the complete
PET acquisition. This can be slow when the scanner geometry and its measurements
are large, especially in 3D.
Ordered-Subsets Expectation-Maximization (OSEM) [2](#footcite-hudsonacceleratedimagereconstruction1994)
partitions the measurements and the forward operator into $L$ matching angular subsets indexed by
$l=1,\ldots,L$:
$$
y = (y_1,\ldots,y_L), \qquad
A = (A_1^\top,\ldots,A_L^\top)^\top.
$$
OSEM applies an EM update successively with each pair $(y_l,A_l)$. Each subset
update is cheaper than a full MLEM iteration and updates the complete image, so OSEM
generally reaches a useful reconstruction in fewer passes over the complete acquisition.
## OSEM reconstruction
OSEM accepts the full measurements and physics and splits them internally. Alternatively,
pre-split inputs can be created with [`deepinv.physics.split_measurements()`](https://deepinv.org/api/stubs/deepinv.physics.split_measurements.html.md#deepinv.physics.split_measurements) and
[`deepinv.physics.split_physics()`](https://deepinv.org/api/stubs/deepinv.physics.split_physics.html.md#deepinv.physics.split_physics) and passed directly to [`deepinv.optim.OSEM`](https://deepinv.org/api/stubs/deepinv.optim.OSEM.html.md#deepinv.optim.OSEM).
```Python
osem_iter = 3
num_subsets = 16
model_osem = dinv.optim.OSEM(
data_fidelity=data_fidelity,
prior=None,
max_iter=osem_iter,
num_subsets=num_subsets,
custom_metrics={
"nrmse": lambda _values, _x_prev, x_cur: nrmse(x_cur.unsqueeze(0), x).item()
},
)
with torch.no_grad():
_sync()
start = time.perf_counter()
x_osem, metrics_osem = model_osem(
y,
physics,
init=torch.ones_like(x),
compute_metrics=True,
)
_sync()
osem_time = time.perf_counter() - start
print(
f"OSEM runtime: {osem_time:.2f} s "
f"for {osem_iter} iterations "
f"({num_subsets} subsets)"
)
print(f"Reconstruction speedup MLEM/OSEM: {mlem_time / osem_time:.2f}x")
psnr_osem = psnr(x_osem, x)
nrmse_osem = nrmse(x_osem, x)
dinv.utils.plot(
[
x[:, :, mid_slice],
x_mlem[:, :, mid_slice],
x_osem[:, :, mid_slice],
],
[
"Ground truth",
f"MLEM ({mlem_iter} it.)",
f"OSEM ({osem_iter} it.)",
],
subtitles=[
"Reference",
f"PSNR: {psnr_mlem.item():.2f} dB\n" f"NRMSE: {100 * nrmse_mlem.item():.2f}%",
f"PSNR: {psnr_osem.item():.2f} dB\n" f"NRMSE: {100 * nrmse_osem.item():.2f}%",
],
rescale_mode="clip",
vmin=0,
vmax=x.max().item(),
figsize=(8, 4),
cbar=True,
)
# We also compare the Poisson objective and NRMSE after every iteration.
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
axes[0].plot(metrics_mlem["cost"][0], label="MLEM")
axes[0].plot(metrics_osem["cost"][0], label="OSEM")
axes[0].set_xlabel("Iteration")
axes[0].set_ylabel("Poisson NLL")
axes[0].legend()
axes[1].plot([100 * v for v in metrics_mlem["nrmse"][0]], label="MLEM")
axes[1].plot([100 * v for v in metrics_osem["nrmse"][0]], label="OSEM")
axes[1].set_xlabel("Iteration")
axes[1].set_ylabel("NRMSE (%)")
axes[1].yaxis.set_major_formatter("{x:.0f}%")
axes[1].legend()
fig.tight_layout()
```
```none
OSEM runtime: 2.69 s for 3 iterations (16 subsets)
Reconstruction speedup MLEM/OSEM: 0.21x
/local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:550: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator. Otherwise, ticks may be mislabeled.
colbar.set_ticklabels(true_labels)
/local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:550: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator. Otherwise, ticks may be mislabeled.
colbar.set_ticklabels(true_labels)
/local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:550: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator. Otherwise, ticks may be mislabeled.
colbar.set_ticklabels(true_labels)
/local/jtachell/deepinv/deepinv/deepinv/utils/plotting.py:570: UserWarning: This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust.
plt.subplots_adjust(hspace=0.2, wspace=0.2)
```
## What next?
Now that you master the basics of PET, you can go further by
- Reconstructing PET with learning-based methods ([PnP](https://deepinv.org/user_guide/reconstruction/iterative.html.md#iterative), [diffusion](https://deepinv.org/user_guide/reconstruction/sampling.html.md#sampling), [unrolled](https://deepinv.org/user_guide/reconstruction/unfolded.html.md#unfolded), etc.)
- Playing with the scanner setup: changing number of detectors, voxel size, etc.
* **References:**
* Lawrence A Shepp and Yehuda Vardi. Maximum likelihood reconstruction for emission tomography. *IEEE Transactions on Medical Imaging*, 1(2):113–122, 1982.
* H.M. Hudson and R.S. Larkin. Accelerated image reconstruction using ordered subsets of projection data. *IEEE Transactions on Medical Imaging*, 13(4):601–609, December 1994. [doi:10.1109/42.363108](https://doi.org/10.1109/42.363108).
**Total running time of the script:** (0 minutes 6.268 seconds)
🏷 Tags: Tomography, PET/SPECT, Low-photon imaging
[`Download Jupyter notebook: demo_pet3d.ipynb`](https://deepinv.org/auto_examples/physics/_downloads/2775468dc95494f8de216087fccc1343/demo_pet3d.ipynb)
[`Download Python source code: demo_pet3d.py`](https://deepinv.org/auto_examples/physics/_downloads/6e7db66cf381ecc0de5f82bdee647998/demo_pet3d.py)
[`Download zipped: demo_pet3d.zip`](https://deepinv.org/auto_examples/physics/_downloads/9b4465fb87591e11f2ba6914b75fc93b/demo_pet3d.zip)
[Gallery generated by Sphinx-Gallery](https://sphinx-gallery.github.io)