Low-level TOF sinogram 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 sinogram mode (using a set of fixed TOF bins along each LOR).

  1"""minimal example that shows how to use the joseph3d TOF forward and back projector in sinogram mode"""
  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#---------------------------------------------------------------
105#--- call the forward projector --------------------------------
106#---------------------------------------------------------------
107
108img_fwd = parallelproj.joseph3d_fwd_tof_sino(xstart, xend, img, img_origin,
109                                             voxel_size, tofbin_width,
110                                             sigma_tof, tofcenter_offset,
111                                             nsigmas, num_tof_bins)
112
113print(img_fwd)
114print(type(img_fwd))
115print(device(img_fwd))
116print('')
117
118#---------------------------------------------------------------
119#--- call the adjoint of the forward projector -----------------
120#---------------------------------------------------------------
121
122# setup a "TOF sinogram"
123sino = to_device(xp.zeros(img_fwd.shape, dtype=xp.float32), dev)
124sino[:, num_tof_bins // 2] = 1
125
126back_img = parallelproj.joseph3d_back_tof_sino(xstart, xend, img_dim,
127                                               img_origin, voxel_size, sino,
128                                               tofbin_width, sigma_tof,
129                                               tofcenter_offset, nsigmas,
130                                               num_tof_bins)
131
132print(back_img[:, :, 3])
133print(type(back_img))
134print(device(back_img))