TOF sinogram 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 sinogram mode"""
  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 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 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, n1//2, -1], # 
 45])
 46
 47vend = xp.array([
 48    [n0//2, n1, n2//2], #           
 49    [n0//2, n1//2, n2], # 
 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# the number of TOF bins 
 66num_tof_bins = 17
 67
 68# number of sigmas after which TOF kernel is truncated
 69nsigmas = 3.
 70
 71# FWHM of the Gaussian TOF kernel in physical units
 72fwhm_tof = 6.
 73
 74# sigma of the Gaussian TOF kernel in physical units
 75# if this is an array of length 1, the same sigma is used
 76# for all LORs
 77sigma_tof = xp.array([fwhm_tof / (2 * xp.sqrt(2 * xp.log(2)))],
 78                     dtype=xp.float32)
 79
 80# TOF center offset for the central TOF bin in physical units
 81# if this is an array of length 1, the same offset is used
 82# for all LORs
 83tofcenter_offset = xp.array([0], dtype=xp.float32)
 84
 85#---------------------------------------------------------------
 86#--- call the forward projector --------------------------------
 87#---------------------------------------------------------------
 88
 89# allocate memory for the forward projection array 
 90# its shape is (number of LORs, number of TOF bins)
 91img_fwd = xp.zeros((xstart.shape[0], num_tof_bins), dtype=xp.float32)
 92
 93parallelproj.joseph3d_fwd_tof_sino(xstart, xend, img, img_origin,
 94                                   voxel_size, img_fwd, tofbin_width,
 95                                   sigma_tof, tofcenter_offset, nsigmas,
 96                                   num_tof_bins)
 97
 98
 99
100#---------------------------------------------------------------
101#--- call the adjoint of the forward projector -----------------
102#---------------------------------------------------------------
103
104# allocate memory for the back projection array
105back_img = xp.zeros((n0,n1,n2), dtype=xp.float32)
106
107# setup a "TOF sinogram"
108sino = xp.zeros_like(img_fwd)
109sino[:,num_tof_bins//2] = 1
110
111parallelproj.joseph3d_back_tof_sino(xstart, xend, back_img, img_origin,
112                                    voxel_size, sino, tofbin_width,
113                                    sigma_tof, tofcenter_offset, nsigmas,
114                                    num_tof_bins)