Low-level TOF listmode projection example

Note

parallelproj aims to be compatible with the python array API and supports different python array backends (numpy, cupy, pytorch). You can change the array backend in this example by commenting / uncommenting the import lines.

Note

The example below shows how to do a simple TOF forward projection along a set of known lines of response with known start and end points in listmode.

  1"""minimal example that shows how to use the joseph3d TOF forward and back projector in listmode"""
  2
  3# parallelproj supports the numpy, cupy and pytorch array API and different devices
  4# choose your preferred array API uncommenting the corresponding line
  5import array_api_compat.numpy as xp
  6#import array_api_compat.cupy as xp
  7#import array_api_compat.torch as xp
  8
  9import parallelproj
 10from array_api_compat import to_device, device
 11
 12# choose a device (CPU or CUDA GPU)
 13if 'numpy' in xp.__name__:
 14    # using numpy, device must be cpu
 15    dev = 'cpu'
 16elif 'cupy' in xp.__name__:
 17    # using cupy, only cuda devices are possible
 18    dev = xp.cuda.Device(0)
 19elif 'torch' in xp.__name__:
 20    # using torch valid choices are 'cpu' or 'cuda'
 21    dev = 'cuda'
 22
 23#---------------------------------------------------------------
 24#--- setup a simple test image ---------------------------------
 25#---------------------------------------------------------------
 26
 27# setup the image dimensions
 28n0, n1, n2 = (7, 7, 7)
 29img_dim = (n0, n1, n2)
 30
 31# define the voxel sizes (in physical units)
 32voxel_size = to_device(xp.asarray([2., 2., 2.], dtype=xp.float32), dev)
 33# define the origin of the image (location of voxel (0,0,0) in physical units)
 34img_origin = (
 35    (-to_device(xp.asarray(img_dim, dtype=xp.float32), dev) / 2 + 0.5) *
 36    voxel_size)
 37
 38# create a simple test image
 39img = to_device(xp.zeros((n0, n1, n2), dtype=xp.float32), dev)
 40img[n0 // 2, n1 // 2, n2 // 2] = 1
 41
 42#---------------------------------------------------------------
 43#--- setup the LOR start and end points ------------------------
 44#---------------------------------------------------------------
 45
 46# Every line of response (LOR) along which we want to project is
 47# defined by its start point (3 element array) and end point (3 element array).
 48# Here we define 2 LORs and group all start and end points in two
 49# 2D arrays of shape (2,3).
 50
 51# We first define the LORs start/end points in voxel coordinates (for convenience)
 52# and convert them later to physical units (as required for the projectors)
 53
 54# define start/end points in voxel coordinates
 55vstart = to_device(
 56    xp.asarray(
 57        [
 58            [n0 // 2, -1, n2 // 2],  # 
 59            [n0 // 2, n1 // 2, -1],  # 
 60        ],
 61        dtype=xp.float32),
 62    dev)
 63
 64vend = to_device(
 65    xp.asarray(
 66        [
 67            [n0 // 2, n1, n2 // 2],  #           
 68            [n0 // 2, n1 // 2, n2],  # 
 69        ],
 70        dtype=xp.float32),
 71    dev)
 72
 73# convert the LOR coordinates to world coordinates (physical units)
 74xstart = vstart * voxel_size + img_origin
 75xend = vend * voxel_size + img_origin
 76
 77#---------------------------------------------------------------
 78#--- setup the TOF related parameters --------------------------
 79#---------------------------------------------------------------
 80
 81# the width of the TOF bins in spatial physical units
 82# same unit as voxel size
 83tofbin_width = 1.5
 84
 85# the number of TOF bins
 86num_tof_bins = 17
 87
 88# number of sigmas after which TOF kernel is truncated
 89nsigmas = 3.
 90
 91# FWHM of the Gaussian TOF kernel in physical units
 92fwhm_tof = 6.
 93
 94# sigma of the Gaussian TOF kernel in physical units
 95# if this is an array of length 1, the same sigma is used
 96# for all LORs
 97sigma_tof = to_device(xp.asarray([fwhm_tof / 2.35], dtype=xp.float32), dev)
 98
 99# TOF center offset for the central TOF bin in physical units
100# if this is an array of length 1, the same offset is used
101# for all LORs
102tofcenter_offset = to_device(xp.asarray([0], dtype=xp.float32), dev)
103
104# setup an array containing the TOF bin of each event
105tof_bin = to_device(xp.zeros(xstart.shape[0], dtype=xp.int16), dev)
106tof_bin[-1] = 1
107
108#---------------------------------------------------------------
109#--- call the forward projector --------------------------------
110#---------------------------------------------------------------
111
112img_fwd = parallelproj.joseph3d_fwd_tof_lm(xstart, xend, img, img_origin,
113                                           voxel_size, tofbin_width, sigma_tof,
114                                           tofcenter_offset, nsigmas, tof_bin)
115
116print(img_fwd)
117print(type(img_fwd))
118print(device(img_fwd))
119print('')
120
121#---------------------------------------------------------------
122#--- call the adjoint of the forward projector -----------------
123#---------------------------------------------------------------
124
125# setup a list of ones to be back projected
126lst = to_device(xp.ones(img_fwd.shape, dtype=xp.float32), dev)
127
128back_img = parallelproj.joseph3d_back_tof_lm(xstart, xend, img_dim, img_origin,
129                                             voxel_size, lst, tofbin_width,
130                                             sigma_tof, tofcenter_offset,
131                                             nsigmas, tof_bin)
132
133print(back_img[:, :, 3])
134print(type(back_img))
135print(device(back_img))