High-level 2D 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 non-TOF forward projection using a predefined 2D parallelview projector.

 1from __future__ import annotations
 2
 3# parallelproj supports the numpy, cupy and pytorch array API
 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
10import matplotlib.pyplot as plt
11from array_api_compat import to_device
12
13# choose our device depending on the array API and the availability of CUDA
14if 'numpy' in xp.__name__:
15    dev = 'cpu'
16elif 'cupy' in xp.__name__:
17    dev = xp.cuda.Device(0)
18elif 'torch' in xp.__name__:
19    # using torch valid choises are cpu or cuda
20    dev = 'cuda'
21
22print(f'running on {dev} device using {xp.__name__}')
23
24# shape of the 2D image
25image_shape = (128, 128)
26# voxel size
27voxel_size = (2., 2.)
28# world coordinates of the [0, 0] pixel
29image_origin = (-127., -127.)
30
31# radial positions of the projection lines
32radial_positions = to_device(xp.linspace(-128, 128, 200), dev)
33# projection angles
34view_angles = to_device(xp.linspace(0, xp.pi, 180, endpoint=False), dev)
35# distance between the image center and the start / end of the center line
36radius = 200.
37
38proj2d = parallelproj.ParallelViewProjector2D(image_shape, radial_positions,
39                                              view_angles, radius,
40                                              image_origin, voxel_size)
41
42img = to_device(xp.zeros(proj2d.in_shape, dtype=xp.float32), dev)
43img[32:96, 32:64] = 1.
44
45img_fwd = proj2d(img)
46img_fwd_back = proj2d.adjoint(img_fwd)
47
48#----------------------------------------------------------------------------
49# show the geometry of the projector and the projections
50
51fig = proj2d.show_views(image=img)
52fig.show()
53
54fig2, ax2 = plt.subplots(1, 3, figsize=(12, 4))
55ax2[0].imshow(to_device(img, 'cpu'))
56ax2[1].imshow(to_device(img_fwd, 'cpu'))
57ax2[2].imshow(to_device(img_fwd_back, 'cpu'))
58ax2[0].set_title('x')
59ax2[1].set_title('A x')
60ax2[2].set_title('A^H A x')
61fig2.tight_layout()
62fig2.show()