Note
Go to the end to download the full example code.
Low-level TOF sinogram 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
125 img_fwd = parallelproj.joseph3d_fwd_tof_sino(
126 xstart,
127 xend,
128 img,
129 img_origin,
130 voxel_size,
131 tofbin_width,
132 sigma_tof,
133 tofcenter_offset,
134 nsigmas,
135 num_tof_bins,
136 )
137
138 print(img_fwd)
139 print(type(img_fwd))
140 print(device(img_fwd))
141 print("")
[[0. 0. 0. 0.00698569 0.03164449 0.10241873
0.23692152 0.3918243 0.46335313 0.3918243 0.23692152 0.10241873
0.03164449 0.00698569 0.00110137 0. 0. ]
[0. 0. 0. 0.00698569 0.03164449 0.10241873
0.23692152 0.3918243 0.46335313 0.3918243 0.23692152 0.10241873
0.03164449 0.00698569 0.00110137 0. 0. ]]
<class 'numpy.ndarray'>
cpu
Call the adjoint of the forward projector
147 # setup a "TOF sinogram"
148 sino = to_device(xp.zeros(img_fwd.shape, dtype=xp.float32), dev)
149 sino[:, num_tof_bins // 2] = 1
150
151 back_img = parallelproj.joseph3d_back_tof_sino(
152 xstart,
153 xend,
154 img_dim,
155 img_origin,
156 voxel_size,
157 sino,
158 tofbin_width,
159 sigma_tof,
160 tofcenter_offset,
161 nsigmas,
162 num_tof_bins,
163 )
164
165 print(back_img[:, :, 3])
166 print(type(back_img))
167 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.92670625 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.003 seconds)
Related examples
