gratopy 0.2.0a1

Gratopy - Graz accelerated tomographic projections for Python

Gratopy

The gratopy (Graz accelerated tomographic projections for Python) toolbox is a Python3 software package for the efficient and high-quality computation of Radon transforms, fanbeam transforms as well as the associated backprojections. The included operators are based on pixel-driven projection methods which were shown to possess favorable approximation properties. The toolbox offers a powerful parallel OpenCL/GPU implementation which admits high execution speed and allows for seamless integration into PyOpenCL. Gratopy can efficiently be combined with other PyOpenCL code and is well-suited for the development of iterative tomographic reconstruction approaches, in particular, for those involving optimization algorithms.

Highlights

• Easy-to-use tomographic projection toolbox.
• High-quality 2D projection operators.
• Fast projection due to custom OpenCL/GPU implementation.
• Seamless integration into PyOpenCL.
• Basic iterative reconstruction schemes included (Landweber, CG, total variation).
• Comprehensive documentation, tests and example code.

The fanbeam projection of a walnut and gratopy’s Landweber and total variation reconstructions (from left to right).

Installation

The toolbox can easily be installed using pip:

pip install gratopy

Alternatively, a release or snapshot archive file can directly be downloaded and unpacked. Calling

pip install .

inside the main folder then installs the toolbox. Gratopy uses uv as its project and dependency management tool, which allows setting up a dedicated development environment by running

uv sync

For more details we refer to the documentation.

Requirements

Most notably, gratopy is built using PyOpenCL, which might require additional care in order to correctly install and configure it: depending on the used platform and GPU, suitable drivers must be installed. We refer to PyOpenCL's documentation for detailed instructions.

The full list of requirements can be found in pyproject.toml.

Getting started

We refer to the extensive documentation, in particular to the getting started guide, the new operator syntax page, as well as to the test files for the Radon transform and fanbeam transform. The following rudimentary example is also included in the documentation.

# initial import
import numpy as np
import pyopencl as cl
import matplotlib.pyplot as plt

import gratopy

# discretization parameters
number_angles = 60
number_detectors = 300
Nx = 300
# Alternatively to number_angles one could give as angle input
# angles = np.linspace(0, np.pi, number_angles+1)[:-1]

# create pyopencl context
ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)

# create phantom as test image (a pyopencl.array.Array of dimensions (Nx, Nx))
phantom = gratopy.phantom(queue,Nx)

# create suitable projectionsettings
PS = gratopy.ProjectionSettings(queue, gratopy.RADON, phantom.shape,
                                number_angles, number_detectors)

# compute forward projection and backprojection of created sinogram
# results are pyopencl arrays
sino = gratopy.forwardprojection(phantom, PS)
backproj = gratopy.backprojection(sino, PS)

# plot results
plt.figure()
plt.title("Generated Phantom")
plt.imshow(phantom.get(), cmap="gray")

plt.figure()
plt.title("Sinogram")
plt.imshow(sino.get(), cmap="gray")

plt.figure()
plt.title("Backprojection")
plt.imshow(backproj.get(), cmap="gray")
plt.show()

Experimental operator syntax

Gratopy also provides an experimental operator-based API for Radon transforms via gratopy.operator.Radon. This interface supports operator algebra such as adjoints and compositions, for example R.T * R, and can also be paired with custom OpenCL kernels for experimentation. As this interface is still experimental, backward-incompatible changes may still occur without a full deprecation cycle.

import numpy as np
import pyopencl as cl
import gratopy

ctx = cl.create_some_context(interactive=False)
queue = cl.CommandQueue(ctx)

img = np.zeros((128, 128), dtype=np.float32)
R = gratopy.operator.Radon(image_domain=128, angles=180)

sino = R.apply_to(img, queue=queue)
backproj = R.T.apply_to(sino)

Authors

• Kristian Bredies, University of Graz, [e-mail]
• Richard Huber, University of Graz, [e-mail]
• Benjamin Hackl, University of Graz, [e-mail]

All authors are affiliated with the Department of Mathematics and Scientific Computing at the University of Graz.

Publications

If you find this tool useful, please cite the following associated publication.

• Kristian Bredies and Richard Huber. (2021). Convergence analysis of pixel-driven Radon and fanbeam transforms. SIAM Journal on Numerical Analysis 59(3), 1399–1432. https://doi.org/10.1137/20M1326635.
• Kristian Bredies and Richard Huber. (2021). Gratopy 0.1 [Software]. Zenodo. https://doi.org/10.5281/zenodo.5221442

Acknowledgements

The development of this software was supported by the following projects:

• Regularization Graphs for Variational Imaging, funded by the Austrian Science Fund (FWF), grant P-29192,

• International Research Training Group IGDK 1754 Optimization and Numerical Analysis for Partial Differential Equations with Nonsmooth Structures, funded by the German Research Council (DFG) and the Austrian Science Fund (FWF), grant W-1244.

The walnut data set included in this toolbox is licensed under CC BY 4.0 and available on Zenodo:

• Keijo Hämäläinen, Lauri Harhanen, Aki Kallonen, Antti Kujanpää, Esa Niemi and Samuli Siltanen. (2015). Tomographic X-ray data of a walnut (Version 1.0.0) [Data set]. Zenodo. https://doi.org/10.5281/zenodo.1254206

The phantom creation code is based on Phantominator, copyright by its contributors and licensed under GPLv3. See https://github.com/mckib2/phantominator.

License

This project is licensed under the GPLv3 license - see LICENSE for details.