TOF listmode projection example

Note

This example can be run using numpy or cupy arrays (if a CUDA GPU is available and the cupy package is installed). To have numpy/cupy agnostic code, we use the import ... as xp lines at the top. In case you want to use numpy (even when cupy is available), simply force import numpy as xp.

  1"""minimal example that shows how to use the joseph3d TOF forward and back projector in listmode"""
  2import parallelproj
  3
  4# parallelproj tells us whether cupy is available and supported
  5# if it is, we use cupy, otherwise numpy as array module (xp)
  6if parallelproj.cupy_enabled:
  7    import cupy as xp
  8else:
  9    import numpy as xp
 10
 11#---------------------------------------------------------------
 12#--- setup a simple test image ---------------------------------
 13#---------------------------------------------------------------
 14
 15# setup the image dimensions
 16n0, n1, n2 = (7, 7, 7)
 17img_dim = xp.array([n0, n1, n2])
 18
 19# define the voxel sizes (in physical units)
 20voxel_size = xp.array([2., 2., 2.], dtype=xp.float32)
 21# define the origin of the image (location of voxel (0,0,0) in physical units)
 22img_origin = ((-img_dim / 2 + 0.5) * voxel_size).astype(xp.float32)
 23
 24# create a simple test image
 25img = xp.zeros((n0,n1,n2), dtype=xp.float32)
 26img[n0//2,n1//2,n2//2] = 1
 27
 28
 29#---------------------------------------------------------------
 30#--- setup the LOR start and end points ------------------------
 31#---------------------------------------------------------------
 32
 33# Every event line of response (LOR) along which we want to project is
 34# defined by its start point (3 element array) and end point (3 element array).
 35# Here we define 2 event LORs and group all start and end points in two
 36# 2D arrays of shape (2,3).
 37
 38# We first define the LORs start/end points in voxel coordinates (for convenience)
 39# and convert them later to physical units (as required for the projectors)
 40
 41# define start/end points in voxel coordinates
 42vstart = xp.array([
 43    [n0//2, -1, n2//2], # 
 44    [n0//2, -1, n2//2], # 
 45])
 46
 47vend = xp.array([
 48    [n0//2, n1, n2//2], #           
 49    [n0//2, n1, n2//2], #           
 50])
 51
 52# convert the LOR coordinates to world coordinates (physical units)
 53xstart = (vstart * voxel_size + img_origin).astype(xp.float32)
 54xend = (vend * voxel_size + img_origin).astype(xp.float32)
 55
 56
 57#---------------------------------------------------------------
 58#--- setup the TOF related parameters --------------------------
 59#---------------------------------------------------------------
 60
 61# the width of the TOF bins in spatial physical units
 62# same unit as voxel size
 63tofbin_width = 1.5
 64
 65# number of sigmas after which TOF kernel is truncated
 66nsigmas = 3.
 67
 68# FWHM of the Gaussian TOF kernel in physical units
 69fwhm_tof = 6.
 70
 71# sigma of the Gaussian TOF kernel in physical units
 72# if this is an array of length 1, the same sigma is used
 73# for all LORs
 74sigma_tof = xp.array([fwhm_tof / (2 * xp.sqrt(2 * xp.log(2)))],
 75                     dtype=xp.float32)
 76
 77# TOF center offset for the central TOF bin in physical units
 78# if this is an array of length 1, the same offset is used
 79# for all LORs
 80tofcenter_offset = xp.array([0], dtype=xp.float32)
 81
 82# setup an array containing the TOF bin of each event
 83tof_bin = xp.zeros(xstart.shape[0], dtype=xp.int16)
 84tof_bin[-1] = 1
 85
 86
 87#---------------------------------------------------------------
 88#--- call the forward projector --------------------------------
 89#---------------------------------------------------------------
 90
 91# allocate memory for the forward projection array 
 92# its shape is (number of LORs, number of TOF bins)
 93img_fwd = xp.zeros(xstart.shape[0], dtype=xp.float32)
 94
 95parallelproj.joseph3d_fwd_tof_lm(xstart, xend, img, img_origin,
 96                                   voxel_size, img_fwd, tofbin_width,
 97                                   sigma_tof, tofcenter_offset, nsigmas,
 98                                   tof_bin)
 99
100
101
102#---------------------------------------------------------------
103#--- call the adjoint of the forward projector -----------------
104#---------------------------------------------------------------
105
106# allocate memory for the back projection array
107back_img = xp.zeros((n0,n1,n2), dtype=xp.float32)
108
109# setup an array of ones
110sino = xp.ones_like(img_fwd)
111
112parallelproj.joseph3d_back_tof_lm(xstart, xend, back_img, img_origin,
113                                    voxel_size, sino, tofbin_width,
114                                    sigma_tof, tofcenter_offset, nsigmas,
115                                    tof_bin)