Note
Go to the end to download the full example code.
Low-level non-TOF projection example
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.
Tip
parallelproj is python array API compatible meaning it supports different
array backends (e.g. numpy, cupy, torch, …) and devices (CPU or GPU).
Choose your preferred array API xp and device dev below.
18 import array_api_compat.numpy as xp
19
20 # import array_api_compat.cupy as xp
21 # import array_api_compat.torch as xp
22
23 import parallelproj
24 from array_api_compat import to_device, device
25
26 # choose a device (CPU or CUDA GPU)
27 if "numpy" in xp.__name__:
28 # using numpy, device must be cpu
29 dev = "cpu"
30 elif "cupy" in xp.__name__:
31 # using cupy, only cuda devices are possible
32 dev = xp.cuda.Device(0)
33 elif "torch" in xp.__name__:
34 # using torch valid choices are 'cpu' or 'cuda'
35 dev = "cuda"
Setup a simple test image
41 # setup the image dimensions
42 n0, n1, n2 = (2, 3, 4)
43 img_dim = (n0, n1, n2)
44
45 # define the voxel sizes (in physical units)
46 voxel_size = to_device(xp.asarray([4.0, 3.0, 2.0], dtype=xp.float32), dev)
47 # define the origin of the image (location of voxel (0,0,0) in physical units)
48 img_origin = (
49 -to_device(xp.asarray(img_dim, dtype=xp.float32), dev) / 2 + 0.5
50 ) * voxel_size
51
52 # create a simple test image
53 img = to_device(
54 xp.reshape(xp.arange(n0 * n1 * n2, dtype=xp.float32), (n0, n1, n2)), dev
55 )
Setup the LOR start and end points
Every line of response (LOR) along which we want to project is defined by its start point (3 element array) and end point (3 element array). Here we define 10 LORs and group all start and end points in two 2D arrays of shape (10,3).
We first define the LORs start/end points in voxel coordinates (for convenience) and convert them later to physical units (as required for the projectors)
69 # define start/end points in voxel coordinates
70 vstart = to_device(
71 xp.asarray(
72 [
73 [0, -1, 0], #
74 [0, -1, 0], #
75 [0, -1, 1], #
76 [0, -1, 0.5], #
77 [0, 0, -1], #
78 [-1, 0, 0], #
79 [n0 - 1, -1, 0], #
80 [n0 - 1, -1, n2 - 1], #
81 [n0 - 1, 0, -1], #
82 [n0 - 1, n1 - 1, -1],
83 ],
84 dtype=xp.float32,
85 ),
86 dev,
87 )
88
89 vend = to_device(
90 xp.asarray(
91 [
92 [0, n1, 0],
93 [0, n1, 0],
94 [0, n1, 1],
95 [0, n1, 0.5],
96 [0, 0, n2],
97 [n0, 0, 0],
98 [n0 - 1, n1, 0],
99 [n0 - 1, n1, n2 - 1], #
100 [n0 - 1, 0, n2],
101 [n0 - 1, n1 - 1, n2],
102 ],
103 dtype=xp.float32,
104 ),
105 dev,
106 )
107
108 # convert the LOR coordinates to world coordinates (physical units)
109 xstart = vstart * voxel_size + img_origin
110 xend = vend * voxel_size + img_origin
Call the forward projector
116 # allocate memory for the forward projection array
117 # call the forward projector
118 img_fwd = parallelproj.joseph3d_fwd(xstart, xend, img, img_origin, voxel_size)
119
120 print(img_fwd)
121 print(type(img_fwd))
122 print(device(img_fwd))
123 print("")
[ 36. 36. 45. 40.5 12. 48. 144. 171. 108. 172. ]
<class 'numpy.ndarray'>
cpu
Call the adjoint of the forward projector
129 # setup a "sinogram" full of ones
130 sino = to_device(xp.ones(img_fwd.shape, dtype=xp.float32), dev)
131
132 # call the back projector
133 back_img = parallelproj.joseph3d_back(
134 xstart, xend, img_dim, img_origin, voxel_size, sino
135 )
136
137 print(back_img)
138 print(type(back_img))
139 print(device(back_img))
[[[13.5 6.5 2. 2. ]
[ 7.5 4.5 0. 0. ]
[ 7.5 4.5 0. 0. ]]
[[ 9. 2. 2. 5. ]
[ 3. 0. 0. 3. ]
[ 5. 2. 2. 5. ]]]
<class 'numpy.ndarray'>
cpu
Total running time of the script: (0 minutes 0.005 seconds)
Related examples
Non-TOF and TOF projections using a modularized (block) PET scanner geometry