Note
Go to the endto download the full example code.
This example demonstrates the use of the MLEM algorithm to minimize the negative Poisson log-likelihood functionusing “sinogram” data from an open PET geometry.
\[f(x) = \sum_{i=1}^m \bar{y}_i (x) - y_i \log(\bar{y}_i (x))\]
subject to
\[x \geq 0\]
using the linear forward model
\[\bar{y}(x) = A x + s\]
parallelproj is python array API compatible meaning it supports differentarray backends (e.g. numpy, cupy, torch, …) and devices (CPU or GPU).Choose your preferred array API xp and device dev below.
31 from __future__ import annotations32 from array_api_strict._array_object import Array3334 import array_api_compat.numpy as xp3536 # import array_api_compat.cupy as xp37 # import array_api_compat.torch as xp3839 import parallelproj40 from array_api_compat import to_device41 import array_api_compat.numpy as np42 import matplotlib.pyplot as plt43 import matplotlib.animation as animation4445 # choose a device (CPU or CUDA GPU)46 if "numpy" in xp.__name__:47 # using numpy, device must be cpu48 dev = "cpu"49 elif "cupy" in xp.__name__:50 # using cupy, only cuda devices are possible51 dev = xp.cuda.Device(0)52 elif "torch" in xp.__name__:53 # using torch valid choices are 'cpu' or 'cuda'54 if parallelproj.cuda_present:55 dev = "cuda"56 else:57 dev = "cpu"
Setup of the forward model \(\bar{y}(x) = A x + s\)
We setup a linear forward operator \(A\) consisting of animage-based resolution model, a non-TOF PET projector and an attenuation model
Here we create an open geometry with 6 sides and 5 rings corresponding toa full geometry using 12 sides where 6 sides were removed.
70 num_rings = 171 scanner = parallelproj.RegularPolygonPETScannerGeometry(72 xp,73 dev,74 radius=65.0,75 num_sides=6,76 num_lor_endpoints_per_side=15,77 lor_spacing=2.3,78 ring_positions=xp.asarray([0.0], device=dev),79 symmetry_axis=2,80 phis=(2 * xp.pi / 12) * xp.asarray([-1, 0, 1, 5, 6, 7]),81 )
setup the LOR descriptor that defines the sinogram
86 img_shape = (40, 40, 1) 87 voxel_size = (2.0, 2.0, 2.0) 88 89 lor_desc = parallelproj.RegularPolygonPETLORDescriptor( 90 scanner, 91 radial_trim=1, 92 sinogram_order=parallelproj.SinogramSpatialAxisOrder.RVP, 93 ) 94 95 proj = parallelproj.RegularPolygonPETProjector( 96 lor_desc, img_shape=img_shape, voxel_size=voxel_size 97 ) 98 99 # setup a simple test image containing a few "hot rods"100 x_true = xp.ones(proj.in_shape, device=dev, dtype=xp.float32)101 c0 = proj.in_shape[0] // 2102 c1 = proj.in_shape[1] // 2103104 x_true[4, c1, :] = 5.0105 x_true[8, c1, :] = 5.0106 x_true[12, c1, :] = 5.0107 x_true[16, c1, :] = 5.0108109 x_true[c0, 4, :] = 5.0110 x_true[c0, 8, :] = 5.0111 x_true[c0, 12, :] = 5.0112 x_true[c0, 16, :] = 5.0113114 x_true[:2, :, :] = 0115 x_true[-2:, :, :] = 0116 x_true[:, :2, :] = 0117 x_true[:, -2:, :] = 0
Attenuation image and sinogram setup
123 # setup an attenuation image124 x_att = 0.01 * xp.astype(x_true > 0, xp.float32)125 # calculate the attenuation sinogram126 att_sino = xp.exp(-proj(x_att))
Complete PET forward model setup
We combine an image-based resolution model,a non-TOF or TOF PET projector and an attenuation modelinto a single linear operator.
137 ## enable TOF - comment if you want to run non-TOF138 # proj.tof_parameters = parallelproj.TOFParameters(139 # num_tofbins=13 * 5,140 # tofbin_width=12.0 / 5,141 # sigma_tof=12.0 / 5,142 # )143144 # setup the attenuation multiplication operator which is different145 # for TOF and non-TOF since the attenuation sinogram is always non-TOF146 if proj.tof:147 att_op = parallelproj.TOFNonTOFElementwiseMultiplicationOperator(148 proj.out_shape, att_sino149 )150 else:151 att_op = parallelproj.ElementwiseMultiplicationOperator(att_sino)152153 res_model = parallelproj.GaussianFilterOperator(154 proj.in_shape, sigma=4.5 / (2.35 * proj.voxel_size)155 )156157 # compose all 3 operators into a single linear operator158 pet_lin_op = parallelproj.CompositeLinearOperator((att_op, proj, res_model))
Visualization of the geometry
165 fig = plt.figure(figsize=(16, 8), tight_layout=True)166 ax1 = fig.add_subplot(121, projection="3d")167 ax2 = fig.add_subplot(122, projection="3d")168 proj.show_geometry(ax1)169 proj.show_geometry(ax2)170 proj.lor_descriptor.show_views(171 ax1,172 views=xp.asarray([0], device=dev),173 planes=xp.asarray([num_rings // 2], device=dev),174 lw=0.5,175 color="k",176 )177 ax1.set_title(f"view 0, plane {num_rings // 2}")178 proj.lor_descriptor.show_views(179 ax2,180 views=xp.asarray([proj.lor_descriptor.num_views // 2], device=dev),181 planes=xp.asarray([num_rings // 2], device=dev),182 lw=0.5,183 color="k",184 )185 ax2.set_title(f"view {proj.lor_descriptor.num_views // 2}, plane {num_rings // 2}")186 fig.tight_layout()187 fig.show()
Simulation of projection data
We setup an arbitrary ground truth \(x_{true}\) and simulatenoise-free and noisy data \(y\) by adding Poisson noise.
196 # simulated noise-free data197 noise_free_data = pet_lin_op(x_true)198199 # generate a contant contamination sinogram200 contamination = xp.full(201 noise_free_data.shape,202 0.5 * float(xp.mean(noise_free_data)),203 device=dev,204 dtype=xp.float32,205 )206207 noise_free_data += contamination208209 # add Poisson noise210 # np.random.seed(1)211 # y = xp.asarray(212 # np.random.poisson(parallelproj.to_numpy_array(noise_free_data)),213 # device=dev,214 # dtype=xp.float64,215 # )216217 y = noise_free_data
EM update to minimize \(f(x)\)
The EM update that can be used in MLEM or OSEM is given by cite:p:Dempster1977 [SV82] [LC84] [HL94]
\[x^+ = \frac{x}{(A^k)^H 1} (A^k)^H \frac{y^k}{A^k x + s^k}\]
to calculate the minimizer of \(f(x)\) iteratively.
To monitor the convergence we calculate the relative cost
\[\frac{f(x) - f(x^*)}{|f(x^*)|}\]
and the distance to the optimal point
\[\frac{\|x - x^*\|}{\|x^*\|}.\]
We setup a function that calculates a single MLEM/OSEMupdate given the current solution, a linear forward operator,data, contamination and the adjoint of ones.
246 def em_update(247 x_cur: Array,248 data: Array,249 op: parallelproj.LinearOperator,250 s: Array,251 adjoint_ones: Array,252 ) -> Array:253 """EM update254255 Parameters256 ----------257 x_cur : Array258 current solution259 data : Array260 data261 op : parallelproj.LinearOperator262 linear forward operator263 s : Array264 contamination265 adjoint_ones : Array266 adjoint of ones267268 Returns269 -------270 Array271 _description_272 """273 ybar = op(x_cur) + s274 return x_cur * op.adjoint(data / ybar) / adjoint_ones
Run the MLEM iterations
281 # number of MLEM iterations282 num_iter = 100283284 # initialize x285 x = xp.ones(pet_lin_op.in_shape, dtype=xp.float32, device=dev)286 # calculate A^H 1287 adjoint_ones = pet_lin_op.adjoint(288 xp.ones(pet_lin_op.out_shape, dtype=xp.float32, device=dev)289 )290291 for i in range(num_iter):292 print(f"MLEM iteration {(i + 1):03} / {num_iter:03}", end="\r")293 x = em_update(x, y, pet_lin_op, contamination, adjoint_ones)
MLEM iteration 001 / 100MLEM iteration 002 / 100MLEM iteration 003 / 100MLEM iteration 004 / 100MLEM iteration 005 / 100MLEM iteration 006 / 100MLEM iteration 007 / 100MLEM iteration 008 / 100MLEM iteration 009 / 100MLEM iteration 010 / 100MLEM iteration 011 / 100MLEM iteration 012 / 100MLEM iteration 013 / 100MLEM iteration 014 / 100MLEM iteration 015 / 100MLEM iteration 016 / 100MLEM iteration 017 / 100MLEM iteration 018 / 100MLEM iteration 019 / 100MLEM iteration 020 / 100MLEM iteration 021 / 100MLEM iteration 022 / 100MLEM iteration 023 / 100MLEM iteration 024 / 100MLEM iteration 025 / 100MLEM iteration 026 / 100MLEM iteration 027 / 100MLEM iteration 028 / 100MLEM iteration 029 / 100MLEM iteration 030 / 100MLEM iteration 031 / 100MLEM iteration 032 / 100MLEM iteration 033 / 100MLEM iteration 034 / 100MLEM iteration 035 / 100MLEM iteration 036 / 100MLEM iteration 037 / 100MLEM iteration 038 / 100MLEM iteration 039 / 100MLEM iteration 040 / 100MLEM iteration 041 / 100MLEM iteration 042 / 100MLEM iteration 043 / 100MLEM iteration 044 / 100MLEM iteration 045 / 100MLEM iteration 046 / 100MLEM iteration 047 / 100MLEM iteration 048 / 100MLEM iteration 049 / 100MLEM iteration 050 / 100MLEM iteration 051 / 100MLEM iteration 052 / 100MLEM iteration 053 / 100MLEM iteration 054 / 100MLEM iteration 055 / 100MLEM iteration 056 / 100MLEM iteration 057 / 100MLEM iteration 058 / 100MLEM iteration 059 / 100MLEM iteration 060 / 100MLEM iteration 061 / 100MLEM iteration 062 / 100MLEM iteration 063 / 100MLEM iteration 064 / 100MLEM iteration 065 / 100MLEM iteration 066 / 100MLEM iteration 067 / 100MLEM iteration 068 / 100MLEM iteration 069 / 100MLEM iteration 070 / 100MLEM iteration 071 / 100MLEM iteration 072 / 100MLEM iteration 073 / 100MLEM iteration 074 / 100MLEM iteration 075 / 100MLEM iteration 076 / 100MLEM iteration 077 / 100MLEM iteration 078 / 100MLEM iteration 079 / 100MLEM iteration 080 / 100MLEM iteration 081 / 100MLEM iteration 082 / 100MLEM iteration 083 / 100MLEM iteration 084 / 100MLEM iteration 085 / 100MLEM iteration 086 / 100MLEM iteration 087 / 100MLEM iteration 088 / 100MLEM iteration 089 / 100MLEM iteration 090 / 100MLEM iteration 091 / 100MLEM iteration 092 / 100MLEM iteration 093 / 100MLEM iteration 094 / 100MLEM iteration 095 / 100MLEM iteration 096 / 100MLEM iteration 097 / 100MLEM iteration 098 / 100MLEM iteration 099 / 100MLEM iteration 100 / 100
Calculation of the negative Poisson log-likelihood function of the reconstruction
300 # calculate the negative Poisson log-likelihood function of the reconstruction301 exp = pet_lin_op(x) + contamination302 # calculate the relative cost and distance to the optimal point303 cost = float(xp.sum(exp - y * xp.log(exp)))304 print(f"\nMLEM cost {cost:.6E} after {num_iter:03} iterations")
MLEM cost -2.586407E+05 after 100 iterations
Visualize the results
311 def _update_img(i):312 img0.set_data(x_true_np[:, :, i])313 img1.set_data(x_np[:, :, i])314 ax[0].set_title(f"true image - plane {i:02}")315 ax[1].set_title(f"MLEM iteration {num_iter} - plane {i:02}")316 return (img0, img1)317318319 x_true_np = parallelproj.to_numpy_array(x_true)320 x_np = parallelproj.to_numpy_array(x)321322 fig, ax = plt.subplots(1, 2, figsize=(10, 5))323 vmax = x_np.max()324 img0 = ax[0].imshow(x_true_np[:, :, 0], cmap="Greys", vmin=0, vmax=vmax)325 img1 = ax[1].imshow(x_np[:, :, 0], cmap="Greys", vmin=0, vmax=vmax)326 ax[0].set_title(f"true image - plane {0:02}")327 ax[1].set_title(f"MLEM iteration {num_iter} - plane {0:02}")328 fig.tight_layout()329 ani = animation.FuncAnimation(fig, _update_img, x_np.shape[2], interval=200, blit=False)330 fig.show()
Total running time of the script: (0 minutes 1.815 seconds)
Download Jupyter notebook: 07_run_mlem_open_geometry.ipynb
Download Python source code: 07_run_mlem_open_geometry.py
Related examples
TOF-MLEM with projection data
TOF-MLEM with projection data
TOF OSEM with projection data
TOF OSEM with projection data
TOF listmode MLEM with projection data
TOF listmode MLEM with projection data
TOF listmode OSEM with projection data
TOF listmode OSEM with projection data
