Note
Go to the end to download the full example code.
Low-level TOF listmode projection example
A minimal example that shows how to use the joseph3d TOF forward and back projector in sinogram mode.
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.
17 import array_api_compat.numpy as xp
18
19 # import array_api_compat.cupy as xp
20 # import array_api_compat.torch as xp
21
22 import parallelproj
23 from array_api_compat import to_device, device
24
25 # choose a device (CPU or CUDA GPU)
26 if "numpy" in xp.__name__:
27 # using numpy, device must be cpu
28 dev = "cpu"
29 elif "cupy" in xp.__name__:
30 # using cupy, only cuda devices are possible
31 dev = xp.cuda.Device(0)
32 elif "torch" in xp.__name__:
33 # using torch valid choices are 'cpu' or 'cuda'
34 dev = "cuda"
Setup a simple test image
40 # setup the image dimensions
41 n0, n1, n2 = (7, 7, 7)
42 img_dim = (n0, n1, n2)
43
44 # define the voxel sizes (in physical units)
45 voxel_size = to_device(xp.asarray([2.0, 2.0, 2.0], dtype=xp.float32), dev)
46 # define the origin of the image (location of voxel (0,0,0) in physical units)
47 img_origin = (
48 -to_device(xp.asarray(img_dim, dtype=xp.float32), dev) / 2 + 0.5
49 ) * voxel_size
50
51 # create a simple test image
52 img = to_device(xp.zeros((n0, n1, n2), dtype=xp.float32), dev)
53 img[n0 // 2, n1 // 2, n2 // 2] = 1
Setup the LOR start and end points
59 # Every line of response (LOR) along which we want to project is
60 # defined by its start point (3 element array) and end point (3 element array).
61 # Here we define 2 LORs and group all start and end points in two
62 # 2D arrays of shape (2,3).
63
64 # We first define the LORs start/end points in voxel coordinates (for convenience)
65 # and convert them later to physical units (as required for the projectors)
66
67 # define start/end points in voxel coordinates
68 vstart = to_device(
69 xp.asarray(
70 [
71 [n0 // 2, -1, n2 // 2], #
72 [n0 // 2, n1 // 2, -1], #
73 ],
74 dtype=xp.float32,
75 ),
76 dev,
77 )
78
79 vend = to_device(
80 xp.asarray(
81 [
82 [n0 // 2, n1, n2 // 2],
83 [n0 // 2, n1 // 2, n2], #
84 ],
85 dtype=xp.float32,
86 ),
87 dev,
88 )
89
90 # convert the LOR coordinates to world coordinates (physical units)
91 xstart = vstart * voxel_size + img_origin
92 xend = vend * voxel_size + img_origin
Call the forward projector
129 img_fwd = parallelproj.joseph3d_fwd_tof_lm(
130 xstart,
131 xend,
132 img,
133 img_origin,
134 voxel_size,
135 tofbin_width,
136 sigma_tof,
137 tofcenter_offset,
138 nsigmas,
139 tof_bin,
140 )
141
142 print(img_fwd)
143 print(type(img_fwd))
144 print(device(img_fwd))
145 print("")
[0.46335313 0.3918243 ]
<class 'numpy.ndarray'>
cpu
Call the adjoint of the forward projector
151 # setup a list of ones to be back projected
152 lst = to_device(xp.ones(img_fwd.shape, dtype=xp.float32), dev)
153
154 back_img = parallelproj.joseph3d_back_tof_lm(
155 xstart,
156 xend,
157 img_dim,
158 img_origin,
159 voxel_size,
160 lst,
161 tofbin_width,
162 sigma_tof,
163 tofcenter_offset,
164 nsigmas,
165 tof_bin,
166 )
167
168 print(back_img[:, :, 3])
169 print(type(back_img))
170 print(device(back_img))
[[0. 0. 0. 0. 0. 0.
0. ]
[0. 0. 0. 0. 0. 0.
0. ]
[0. 0. 0. 0. 0. 0.
0. ]
[0.03164449 0.14060059 0.34391424 0.8551774 0.34391424 0.14060059
0.03164449]
[0. 0. 0. 0. 0. 0.
0. ]
[0. 0. 0. 0. 0. 0.
0. ]
[0. 0. 0. 0. 0. 0.
0. ]]
<class 'numpy.ndarray'>
cpu
Total running time of the script: (0 minutes 0.004 seconds)
Related examples
