Non-TOF 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 non-TOF forward and back projector"""
 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 = (2, 3, 4)
17img_dim = xp.array([n0, n1, n2])
18
19# define the voxel sizes (in physical units)
20voxel_size = xp.array([4., 3., 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.arange(n0 * n1 * n2, dtype=xp.float32).reshape((n0, n1, n2))
26
27
28#---------------------------------------------------------------
29#--- setup the LOR start and end points ------------------------
30#---------------------------------------------------------------
31
32# Every line of response (LOR) along which we want to project is
33# defined by its start point (3 element array) and end point (3 element array).
34# Here we define 10 LORs and group all start and end points in two
35# 2D arrays of shape (10,3).
36
37# We first define the LORs start/end points in voxel coordinates (for convenience)
38# and convert them later to physical units (as required for the projectors)
39
40# define start/end points in voxel coordinates
41vstart = xp.array([
42    [0, -1, 0], # 
43    [0, -1, 0], #
44    [0, -1, 1], #
45    [0, -1, 0.5], #
46    [0, 0, -1], #
47    [-1, 0, 0], #
48    [n0 - 1, -1, 0], # 
49    [n0 - 1, -1, n2 - 1], #
50    [n0 - 1, 0, -1], #
51    [n0 - 1, n1 - 1, -1]
52])
53
54vend = xp.array([
55    [0, n1, 0], #           
56    [0, n1, 0], #           
57    [0, n1, 1], #          
58    [0, n1, 0.5], #         
59    [0, 0, n2], #          
60    [n0, 0, 0], #          
61    [n0 - 1, n1, 0], #      
62    [n0 - 1, n1, n2 - 1], # 
63    [n0 - 1, 0, n2], #     
64    [n0 - 1, n1 - 1, n2]
65])
66
67# convert the LOR coordinates to world coordinates (physical units)
68xstart = (vstart * voxel_size + img_origin).astype(xp.float32)
69xend = (vend * voxel_size + img_origin).astype(xp.float32)
70
71
72#---------------------------------------------------------------
73#--- call the forward projector --------------------------------
74#---------------------------------------------------------------
75
76# allocate memory for the forward projection array
77img_fwd = xp.zeros(xstart.shape[0], dtype=xp.float32)
78
79# call the forward projector
80parallelproj.joseph3d_fwd(xstart, xend, img, img_origin, voxel_size,
81                          img_fwd)
82
83
84#---------------------------------------------------------------
85#--- call the adjoint of the forward projector -----------------
86#---------------------------------------------------------------
87
88# setup a "sinogram" full of ones
89sino = xp.ones_like(img_fwd)
90
91# allocate memory for the back projection array
92back_img = xp.zeros_like(img)
93
94# call the back projector
95parallelproj.joseph3d_back(xstart, xend, back_img, img_origin, voxel_size,
96                           sino)