Low-level non-TOF 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 non-TOF forward projection along a set of known lines of response with known start and end points.

  1"""minimal example that shows how to use the joseph3d non-TOF forward and back projector"""
  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 = (2, 3, 4)
 29img_dim = (n0, n1, n2)
 30
 31# define the voxel sizes (in physical units)
 32voxel_size = to_device(xp.asarray([4., 3., 2.], dtype=xp.float32), dev)
 33# define the origin of the image (location of voxel (0,0,0) in physical units)
 34img_origin = (-to_device(xp.asarray(img_dim, dtype=xp.float32), dev) / 2 +
 35              0.5) * voxel_size
 36
 37# create a simple test image
 38img = to_device(
 39    xp.reshape(xp.arange(n0 * n1 * n2, dtype=xp.float32), (n0, n1, n2)), dev)
 40
 41#---------------------------------------------------------------
 42#--- setup the LOR start and end points ------------------------
 43#---------------------------------------------------------------
 44
 45# Every line of response (LOR) along which we want to project is
 46# defined by its start point (3 element array) and end point (3 element array).
 47# Here we define 10 LORs and group all start and end points in two
 48# 2D arrays of shape (10,3).
 49
 50# We first define the LORs start/end points in voxel coordinates (for convenience)
 51# and convert them later to physical units (as required for the projectors)
 52
 53# define start/end points in voxel coordinates
 54vstart = to_device(
 55    xp.asarray(
 56        [
 57            [0, -1, 0],  # 
 58            [0, -1, 0],  #
 59            [0, -1, 1],  #
 60            [0, -1, 0.5],  #
 61            [0, 0, -1],  #
 62            [-1, 0, 0],  #
 63            [n0 - 1, -1, 0],  # 
 64            [n0 - 1, -1, n2 - 1],  #
 65            [n0 - 1, 0, -1],  #
 66            [n0 - 1, n1 - 1, -1]
 67        ],
 68        dtype=xp.float32),
 69    dev)
 70
 71vend = to_device(
 72    xp.asarray(
 73        [
 74            [0, n1, 0],  #           
 75            [0, n1, 0],  #           
 76            [0, n1, 1],  #          
 77            [0, n1, 0.5],  #         
 78            [0, 0, n2],  #          
 79            [n0, 0, 0],  #          
 80            [n0 - 1, n1, 0],  #      
 81            [n0 - 1, n1, n2 - 1],  # 
 82            [n0 - 1, 0, n2],  #     
 83            [n0 - 1, n1 - 1, n2]
 84        ],
 85        dtype=xp.float32),
 86    dev)
 87
 88# convert the LOR coordinates to world coordinates (physical units)
 89xstart = vstart * voxel_size + img_origin
 90xend = vend * voxel_size + img_origin
 91
 92#---------------------------------------------------------------
 93#--- call the forward projector --------------------------------
 94#---------------------------------------------------------------
 95
 96# allocate memory for the forward projection array
 97# call the forward projector
 98img_fwd = parallelproj.joseph3d_fwd(xstart, xend, img, img_origin, voxel_size)
 99
100print(img_fwd)
101print(type(img_fwd))
102print(device(img_fwd))
103print('')
104
105#---------------------------------------------------------------
106#--- call the adjoint of the forward projector -----------------
107#---------------------------------------------------------------
108
109# setup a "sinogram" full of ones
110sino = to_device(xp.ones(img_fwd.shape, dtype=xp.float32), dev)
111
112# call the back projector
113back_img = parallelproj.joseph3d_back(xstart, xend, img_dim, img_origin,
114                                      voxel_size, sino)
115
116print(back_img)
117print(type(back_img))
118print(device(back_img))