add listmode examples

examples/01_pet_geometry/README.txt

@@ -1,2 +1,2 @@
-PET scanner and sinogram geometries
------------------------------------
+PET scanner and sinogram geometry examples
+------------------------------------------

examples/02_pet_sinogram_projections/README.txt

@@ -1,2 +1,2 @@
-High-level PET sinogram projections
------------------------------------
+High-level PET sinogram projector examples
+------------------------------------------

examples/03_pet_listmode_projections/README.txt

@@ -0,0 +1,2 @@
+High-level PET listmode projector examples
+------------------------------------------

examples/03_pet_listmode_projections/demo_01_pet_non_tof_listmode_projector.py

@@ -0,0 +1,170 @@
+"""
+PET non-TOF listmode projector
+==============================
+
+In this example we will show how to setup and use a non-TOF 
+PET listmode projector including geometrical forward projection, 
+resolution model and correction for attenuation.
+"""
+
+# %%
+# parallelproj supports the numpy, cupy and pytorch array API and different devices
+# choose your preferred array API uncommenting the corresponding line
+
+import array_api_compat.numpy as xp
+
+# import array_api_compat.cupy as xp
+# import array_api_compat.torch as xp
+
+# %%
+import parallelproj
+from array_api_compat import to_device
+import array_api_compat.numpy as np
+import matplotlib.pyplot as plt
+
+# choose a device (CPU or CUDA GPU)
+if "numpy" in xp.__name__:
+    # using numpy, device must be cpu
+    dev = "cpu"
+elif "cupy" in xp.__name__:
+    # using cupy, only cuda devices are possible
+    dev = xp.cuda.Device(0)
+elif "torch" in xp.__name__:
+    # using torch valid choices are 'cpu' or 'cuda'
+    dev = "cuda"
+
+# %%
+# setup a small regular polygon PET scanner with 5 rings (polygons)
+# -----------------------------------------------------------------
+
+num_rings = 4
+scanner = parallelproj.RegularPolygonPETScannerGeometry(
+    xp,
+    dev,
+    radius=65.0,
+    num_sides=12,
+    num_lor_endpoints_per_side=15,
+    lor_spacing=2.3,
+    ring_positions=xp.linspace(-10, 10, num_rings),
+    symmetry_axis=1,
+)
+
+# %%
+# generate 4 arbitrary listmode events
+# ------------------------------------
+
+start_ring = xp.asarray([2, 1, 0, 3], device=dev)
+start_xtal = xp.asarray([0, 59, 143, 75], device=dev)
+
+end_ring = xp.asarray([2, 0, 1, 3], device=dev)
+end_xtal = xp.asarray([79, 140, 33, 147], device=dev)
+
+event_start_coordinates = scanner.get_lor_endpoints(start_ring, start_xtal)
+event_end_coordinates = scanner.get_lor_endpoints(end_ring, end_xtal)
+
+print(event_start_coordinates)
+print(event_end_coordinates)
+
+# %%
+# show the scanner geometry and the events
+# ----------------------------------------
+
+fig = plt.figure(figsize=(8, 8))
+ax = fig.add_subplot(111, projection="3d")
+scanner.show_lor_endpoints(ax)
+for i in range(event_start_coordinates.shape[0]):
+    ax.plot(
+        [float(event_start_coordinates[i, 0]), float(event_end_coordinates[i, 0])],
+        [float(event_start_coordinates[i, 1]), float(event_end_coordinates[i, 1])],
+        [float(event_start_coordinates[i, 2]), float(event_end_coordinates[i, 2])],
+    )
+fig.tight_layout()
+fig.show()
+
+
+# %%
+# setup a listmode projector and a test image
+# -------------------------------------------
+
+img_shape = (40, 9, 40)
+voxel_size = (2.0, 3.0, 2.0)
+
+lm_proj = parallelproj.ListmodePETProjector(
+    event_start_coordinates, event_end_coordinates, img_shape, voxel_size
+)
+
+x = xp.ones(img_shape, dtype=xp.float32, device=dev)
+
+# %% peform listmode forward and back projections
+# -----------------------------------------------
+
+x_fwd = lm_proj(x)
+print(x_fwd)
+
+# back project a list of ones
+ones_list = xp.ones(lm_proj.num_events, dtype=xp.float32, device=dev)
+y_back = lm_proj.adjoint(ones_list)
+
+# %%
+# show the backprojected list of ones (events)
+# --------------------------------------------
+
+fig2, ax2 = plt.subplots(3, 3, figsize=(8, 8))
+vmax = float(xp.max(y_back))
+for i in range(ax2.size):
+    if i < y_back.shape[1]:
+        axx = ax2.ravel()[i]
+        axx.imshow(
+            np.asarray(to_device(y_back[:, i, :].T, "cpu")),
+            cmap="Greys",
+            vmin=0,
+            vmax=vmax,
+        )
+        axx.set_title(f"img plane {i}", fontsize="medium")
+    else:
+        ax2.ravel()[i].set_axis_off()
+fig2.tight_layout()
+fig2.show()
+
+# %% combine the listmode projector with a resolution and attenuation model
+# -------------------------------------------------------------------------
+
+# setup a simple image-based resolution model with an Gaussian FWHM of 4.5mm
+res_model = parallelproj.GaussianFilterOperator(
+    lm_proj.in_shape, sigma=4.5 / (2.35 * lm_proj.voxel_size)
+)
+
+# define arbritrary attenuation factors
+att_list = xp.asarray([0.3, 0.4, 0.2, 0.6], device=dev)
+att_op = parallelproj.ElementwiseMultiplicationOperator(att_list)
+
+
+lm_proj_with_res_model_and_att = parallelproj.CompositeLinearOperator(
+    (att_op, lm_proj, res_model)
+)
+
+x_fwd2 = lm_proj_with_res_model_and_att(x)
+print(x_fwd2)
+
+y_back2 = lm_proj_with_res_model_and_att.adjoint(ones_list)
+
+# %%
+# show the backprojected list of ones (events)
+# --------------------------------------------
+
+fig3, ax3 = plt.subplots(3, 3, figsize=(8, 8))
+vmax = float(xp.max(y_back2))
+for i in range(ax3.size):
+    if i < y_back.shape[1]:
+        axx = ax3.ravel()[i]
+        axx.imshow(
+            np.asarray(to_device(y_back2[:, i, :].T, "cpu")),
+            cmap="Greys",
+            vmin=0,
+            vmax=vmax,
+        )
+        axx.set_title(f"img plane {i}", fontsize="medium")
+    else:
+        ax3.ravel()[i].set_axis_off()
+fig3.tight_layout()
+fig3.show()

examples/03_pet_listmode_projections/demo_02_pet_tof_listmode_projector.py

@@ -0,0 +1,183 @@
+"""
+PET TOF listmode projector
+==========================
+
+In this example we will show how to setup and use a non-TOF 
+PET listmode projector including geometrical forward projection, 
+resolution model and correction for attenuation.
+"""
+
+# %%
+# parallelproj supports the numpy, cupy and pytorch array API and different devices
+# choose your preferred array API uncommenting the corresponding line
+
+import array_api_compat.numpy as xp
+
+# import array_api_compat.cupy as xp
+# import array_api_compat.torch as xp
+
+# %%
+import parallelproj
+from array_api_compat import to_device
+import array_api_compat.numpy as np
+import matplotlib.pyplot as plt
+
+# choose a device (CPU or CUDA GPU)
+if "numpy" in xp.__name__:
+    # using numpy, device must be cpu
+    dev = "cpu"
+elif "cupy" in xp.__name__:
+    # using cupy, only cuda devices are possible
+    dev = xp.cuda.Device(0)
+elif "torch" in xp.__name__:
+    # using torch valid choices are 'cpu' or 'cuda'
+    dev = "cuda"
+
+# %%
+# setup a small regular polygon PET scanner with 5 rings (polygons)
+# -----------------------------------------------------------------
+
+num_rings = 4
+scanner = parallelproj.RegularPolygonPETScannerGeometry(
+    xp,
+    dev,
+    radius=65.0,
+    num_sides=12,
+    num_lor_endpoints_per_side=15,
+    lor_spacing=2.3,
+    ring_positions=xp.linspace(-10, 10, num_rings),
+    symmetry_axis=1,
+)
+
+# %%
+# generate 4 arbitrary listmode events
+# ------------------------------------
+
+start_ring = xp.asarray([2, 1, 0, 3], device=dev)
+start_xtal = xp.asarray([0, 59, 143, 75], device=dev)
+
+end_ring = xp.asarray([2, 0, 1, 3], device=dev)
+end_xtal = xp.asarray([79, 140, 33, 147], device=dev)
+
+event_start_coordinates = scanner.get_lor_endpoints(start_ring, start_xtal)
+event_end_coordinates = scanner.get_lor_endpoints(end_ring, end_xtal)
+
+print(event_start_coordinates)
+print(event_end_coordinates)
+
+# setup TOF parameters with a TOF resolution FWHM = 30mm (ca 200ps)
+tof_params = parallelproj.TOFParameters(
+    sigma_tof=30.0 / 2.35, num_tofbins=59, tofbin_width=12.5
+)
+event_tof_bins = xp.asarray([0, 1, 0, -1], device=dev)
+
+# %%
+# show the scanner geometry and the events
+# ----------------------------------------
+
+fig = plt.figure(figsize=(8, 8))
+ax = fig.add_subplot(111, projection="3d")
+scanner.show_lor_endpoints(ax)
+for i in range(event_start_coordinates.shape[0]):
+    ax.plot(
+        [float(event_start_coordinates[i, 0]), float(event_end_coordinates[i, 0])],
+        [float(event_start_coordinates[i, 1]), float(event_end_coordinates[i, 1])],
+        [float(event_start_coordinates[i, 2]), float(event_end_coordinates[i, 2])],
+    )
+fig.tight_layout()
+fig.show()
+
+
+# %%
+# setup a TOF listmode projector and a test image
+# -----------------------------------------------
+
+img_shape = (40, 9, 40)
+voxel_size = (2.0, 3.0, 2.0)
+
+lm_proj = parallelproj.ListmodePETProjector(
+    event_start_coordinates, event_end_coordinates, img_shape, voxel_size
+)
+
+# set the TOF parameters
+lm_proj.tof_parameters = tof_params
+# set the event TOF bins
+lm_proj.event_tofbins = event_tof_bins
+# enable TOF
+lm_proj.tof = True
+
+x = xp.ones(img_shape, dtype=xp.float32, device=dev)
+
+# %% peform listmode forward and back projections
+# -----------------------------------------------
+
+x_fwd = lm_proj(x)
+print(x_fwd)
+
+# back project a list of ones
+ones_list = xp.ones(lm_proj.num_events, dtype=xp.float32, device=dev)
+y_back = lm_proj.adjoint(ones_list)
+
+# %%
+# show the backprojected list of ones (events)
+# --------------------------------------------
+
+fig2, ax2 = plt.subplots(3, 3, figsize=(8, 8))
+vmax = float(xp.max(y_back))
+for i in range(ax2.size):
+    if i < y_back.shape[1]:
+        axx = ax2.ravel()[i]
+        axx.imshow(
+            np.asarray(to_device(y_back[:, i, :].T, "cpu")),
+            cmap="Greys",
+            vmin=0,
+            vmax=vmax,
+        )
+        axx.set_title(f"img plane {i}", fontsize="medium")
+    else:
+        ax2.ravel()[i].set_axis_off()
+fig2.tight_layout()
+fig2.show()
+
+# %% combine the listmode projector with a resolution and attenuation model
+# -------------------------------------------------------------------------
+
+# setup a simple image-based resolution model with an Gaussian FWHM of 4.5mm
+res_model = parallelproj.GaussianFilterOperator(
+    lm_proj.in_shape, sigma=4.5 / (2.35 * lm_proj.voxel_size)
+)
+
+# define arbritrary attenuation factors
+att_list = xp.asarray([0.3, 0.4, 0.2, 0.6], device=dev)
+att_op = parallelproj.ElementwiseMultiplicationOperator(att_list)
+
+
+lm_proj_with_res_model_and_att = parallelproj.CompositeLinearOperator(
+    (att_op, lm_proj, res_model)
+)
+
+x_fwd2 = lm_proj_with_res_model_and_att(x)
+print(x_fwd2)
+
+y_back2 = lm_proj_with_res_model_and_att.adjoint(ones_list)
+
+# %%
+# show the backprojected list of ones (events)
+# --------------------------------------------
+
+fig3, ax3 = plt.subplots(3, 3, figsize=(8, 8))
+vmax = float(xp.max(y_back2))
+for i in range(ax3.size):
+    if i < y_back.shape[1]:
+        axx = ax3.ravel()[i]
+        axx.imshow(
+            np.asarray(to_device(y_back2[:, i, :].T, "cpu")),
+            cmap="Greys",
+            vmin=0,
+            vmax=vmax,
+        )
+        axx.set_title(f"img plane {i}", fontsize="medium")
+    else:
+        ax3.ravel()[i].set_axis_off()
+fig3.tight_layout()
+fig3.show()

python/parallelproj/projectors.py

@@ -912,6 +912,11 @@ def event_end_coordinates(self) -> Array:
         """coordinates of LOR end points"""
         return self._xend
 
+    @property
+    def voxel_size(self) -> Array:
+        """voxel size"""
+        return self._voxel_size
+
     def _apply(self, x: Array) -> Array:
         dev = array_api_compat.device(x)