RadonKA.jl

A simple yet sufficiently fast Radon and adjoint Radon (backproject) transform implementation using KernelAbstractions.jl. It offers multithreading and CUDA support and outperforms any existing Julia Radon transforms (at least the ones we are aware of). On CUDA it is faster much than Matlab and it offers the same or faster speed than ASTRA.

⚠️ This package is still very young. I would be happy to receive any feedback or if we can improve anything, just open an issue! ⚠️

Quick Overview

• For 2D and 3D arrays
• parallel radon and backproject (?RadonParallelCircle)
• attenuated radon and backproject (see the parameter μ) and see this paper as reference)
• arbitrary 2D geometries where starting and endpoint of each ray can be specified (fan beam could be a special case of this) (?RadonFlexibleCircle)
• different strength weighting of rays
• based on KernelAbstractions.jl (tested on CPU() and CUDABackend())
• registered adjoint rules for both radon and backproject with ChainRulesCore.jl, hence automatic differentiation (AD) compatible.
• high performance however not ultra high performance. On par with ASTRA, on CUDA faster than Matlab.
• simple and extensible API

Installation

This toolbox runs with CUDA support on Linux, Windows and MacOS! Requires at least Julia 1.9

julia> ]add RadonKA

Simple use

using RadonKA, ImageShow, ImageIO, TestImages

img = Float32.(testimage("resolution_test_512"))
angles = range(0f0, 2f0π, 500)[begin:end-1]

# 0.085398 seconds (260 allocations: 1.006 MiB)
@time sinogram = radon(img, angles);
# 0.127043 seconds (251 allocations: 1.036 MiB)
@time backproject = RadonKA.backproject(sinogram, angles);

simshow(sinogram)
simshow(backproject)

using CUDA
img_c = CuArray(img)
# 0.003363 seconds (244 CPU allocations: 18.047 KiB) (7 GPU allocations: 1007.934 KiB, 0.96% memmgmt time)
CUDA.@time sinogram = radon(img_c, angles);
# 0.005928 seconds (218 CPU allocations: 16.109 KiB) (7 GPU allocations: 1.012 MiB, 0.49% memmgmt time)
CUDA.@time backproject = RadonKA.backproject(sinogram, angles);

Examples

See the documentation. You can also run the examples locally. Download this repository and then do the following in your REPL:

julia> cd("examples/")

julia> using Pkg; Pkg.activate("."); Pkg.instantiate()
  Activating project at `~/.julia/dev/RadonKA.jl/examples`

julia> using Pluto; Pluto.run()

A browser should open. The following examples show case the ability of this package:

• Simple radon and backproject: Pluto notebook
• Different geometries: Pluto notebook
• Reconstruction of a CT dataset with an optimizer: Pluto notebook
• How this package is used in Tomographic Volumetric Additive Manufacturing (3D printing): Pluto notebook

Citation

This package was created as part of scientific work. Please consider citing it :)

@article{Wechsler:24,
author = {Felix Wechsler and Carlo Gigli and Jorge Madrid-Wolff and Christophe Moser},
journal = {Opt. Express},
keywords = {3D printing; Computed tomography; Liquid crystal displays; Material properties; Ray tracing; Refractive index},
number = {8},
pages = {14705--14712},
publisher = {Optica Publishing Group},
title = {Wave optical model for tomographic volumetric additive manufacturing},
volume = {32},
month = {Apr},
year = {2024},
url = {https://opg.optica.org/oe/abstract.cfm?URI=oe-32-8-14705},
doi = {10.1364/OE.521322},
}

Development

File an issue on GitHub if you encounter any problems. You can also join my conference room. Give me a minute to join!

Similar packages

Python

There is TIGRE and ASTRA which both offer more functionality for classic CT problems. They also feature GPU acceleration, however we did not observe that they outperform this package. Also, they don't allow to calculate the attenuated Radon transform and don't allow for arbitrary ray geometries, as we do. The fastest implementation we found, is the unmaintained torch-radon. Its kernels are written in CUDA C code and offer a PyTorch interface. There is a torch-radon fork which allows to run it with newer versions. It offers no attenuated Radon transform.

Julia

There exists Sinograms.jl and in general the whole JuliaImageRecon organization. Again, no arbitrary geometries can be specified. And also no attenuated Radon transform is possible.

Matlab

Matlab has built-in a radon and iradon(...,'linear','none'); transform which is similar to our lightweight API. However, no CUDA acceleration, no 3D arrays and no attenuated Radon transform.