[ANN] RadonKA.jl - A simple (exponential) Radon Transform in Julia

[roflmaostc]

Hi!

I created a package called RadonKA.jl which provides a parallel Radon and inverse Radon pair.

What can it do?

Further, it provides an exponential Radon transform which does not only calculate the line integral but also weights each value with an exponential decreasing factor depending on the distance to a circular boundary.

It’s all based on tracing rays through the pixel grid and weighting the absorption according to the length of the intersection with a cell.

Performance

The implementation is rather naive and totally ignores cache efficiency. But, based on KernelAbstractions.jl we utilize CPU threading and we also tested the CUDA backend.
I’m surprised that the performance is very similiar to Matlab’s Radon transform and the Astra toolbox. See details of the benchmark here.

Simple example

It’s not registered yet (in 3 days), but try it out with:

julia> ]add https://github.com/roflmaostc/RadonKA.jl

using RadonKA, ImageShow, ImageIO, TestImages

img = Float32.(testimage("resolution_test_512"))

angles = range(0f0, 2f0π, 500)[begin:end-1]

# 0.196049 seconds (145 allocations: 1009.938 KiB)
@time sinogram = radon(img, angles);

# 0.268649 seconds (147 allocations: 1.015 MiB)
@time backproject = RadonKA.iradon(sinogram, angles);
simshow(sinogram)

What’s next?

This package is relatively young and experimental, so please try it out and report bugs (also the reason why I post it here )! I’m willing to improve it. I provide some examples.
Further, I will work on an interface where rays don’t need to be parallel but ray angles can be specified separately. This include special cases such as Fan Beam Tomography.
Cone Beam would be more tricky, since I don’t trace 3D paths currently, only 2D slices in a 3D volume.

So yes I’m excited since it really helps already with my work regarding Volumetric Additive Manufacturing. Especially, the exponential Radon transform was important for my work, and I didn’t find any efficient implementation in any framework.

Best,

Felix

[roflmaostc]

I just tagged a new release of RadonKA.jl. The new features are documented here

There is some performance improvements but most notably there is now a way to specify the geometry of your ray propagation.
For that you describe the entrance and exit point of each ray inside the RadonFlexibleCircle constructor.

For example:

using RadonKA, ImageShow

angles = [0]

# output image size
N = 200

sinogram = zeros((N - 1, length(angles)))
sinogram[1:5:end] .= 1
geometry_cone = RadonFlexibleCircle(N, 
                          -(N-1)÷2:(N-1)÷2,                  # entrance position
                          range(-(N-1)÷4, (N-1)÷4, N-1))    # exit position


projected_cone = iradon(sinogram, angles; geometry=geometry_cone);

simshow(projected_cone, γ=0.01)

And also the attenuated (or exponential) Radon transform got some documentation.
Rays are getting attenuated with \exp(-x \cdot \mu) depending on the distance to the circle.

projected_exp = iradon(sinogram, angles; geometry=geometry_extreme, μ=0.04);

simshow(projected_exp)

Of course registered reverse rules for automatic differentation are available for both radon and iradon.

[roflmaostc]

Just a bit of an update:

Rename of iradon

iradon was renamed to backproject to indicate that no filtering is done.
There is also backproject_filtered to do the classical filtered backprojection.

Spatially Varying Absorption

Also I’m going to add soon the option that the absorption can be spatially varying. In this example we place a small box which absorbs a lot of light intensity.

μ = 0.5 .* box((N, N), (2, 50));

projected_1 = backproject(sinogram, angles, μ=μ);
projected_2 = backproject(sinogram, angles .+ π / 4, μ=μ);

[simshow(projected_1) simshow(projected_2)]

Benchmarks

I’m a bit disappointed that torch-radon beats us by a factor of 5-10 on both GPUs and CPUs.
From reading the source code I don’t see really why.
Maybe some motivated people spot performance penalties in my critical lines.

[sepehr]

Looks cool!