radon(I, θs; <kwargs>)

Calculates the parallel Radon transform of the AbstractArray I. Intuitively, the radon sums array entries of I along ray paths.

The first two dimensions are y and x. The third dimension is z, the rotational axis. size(I, 1) and size(I, 2) have to be equal. The Radon transform is rotated around the pixel size(I, 1) ÷ 2 + 1, so there is always an integer center pixel! Works either with a AbstractArray{T, 3} or AbstractArray{T, 2}.

θs is a vector or range storing the angles in radians.

In principle, all backends of KernelAbstractions.jl should work but are not tested. CUDA and CPU arrays are actively tested. Both radon and backproject are differentiable with respect to I.

Keywords

μ=nothing - Attenuated Radon Transform

If μ != nothing, then the rays are attenuated with exp(-μ * dist) where dist is the distance to the circular boundary of the field of view. μ is in units of pixel length. So μ=1 corresponds to an attenuation of exp(-1) if propagated through one pixel. If isnothing(μ), then the rays are not attenuated.

μ can be also an array of the same size as I which allows for spatially varying attenuation. For example, μ=0.01 and μ = ones(size(I)) * 0.01 are equivalent. In practice there is a slight difference because we calculate the attenuation differently.

geometry = RadonParallelCircle(-(size(img,1)-1)÷2:(size(img,1)-1)÷2)

This corresponds to a parallel Radon transform. See ?RadonGeometries for a full list of geometries. There is also the very flexible RadonFlexibleCircle.

See also backproject.

Example

The reason the sinogram has the value 1.41421 for the diagonal ray π/4 is, that such a diagonal travels a longer distance through the pixel.

julia> arr = zeros((4,4)); arr[3,3] = 1;

julia> radon(arr, [0, π/4, π/2])
3×3 view(::Array{Float64, 3}, :, :, 1) with eltype Float64:
 0.0  0.0      0.0
 1.0  1.41421  1.0
 0.0  0.0      0.0

Choose different detector

julia> array = ones((6,6))
6×6 Matrix{Float64}:
 1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0

julia> radon(array, [0])
5×1 view(::Array{Float64, 3}, :, :, 1) with eltype Float64:
 1.0
 3.7320508075688767
 5.0
 3.7320508075688767
 1.0

julia> radon(array, [0], geometry=RadonParallelCircle(6, -2:2))
5×1 view(::Array{Float64, 3}, :, :, 1) with eltype Float64:
 1.0
 3.7320508075688767
 5.0
 3.7320508075688767
 1.0

julia> radon(array, [0], geometry=RadonParallelCircle(6, -2:2:2))
3×1 view(::Array{Float64, 3}, :, :, 1) with eltype Float64:
 1.0
 5.0
 1.0

Apply some weights on the rays

julia> array = ones((6,6));

julia> radon(array, [0], geometry=RadonParallelCircle(6, -2:2, [2,1,0,1,2]))
5×1 view(::Array{Float64, 3}, :, :, 1) with eltype Float64:
 2.0
 3.7320508075688767
 0.0
 3.7320508075688767
 2.0

backproject(sinogram, θs; <kwargs>)

Conceptually the adjoint operation of radon. Intuitively, the backproject smears rays back into the space. See also radon.

For filtered backprojection see backproject_filtered.

Example

julia> arr = zeros((5,2)); arr[2,:] .= 1; arr[4, :] .= 1
2-element view(::Matrix{Float64}, 4, :) with eltype Float64:
 1.0
 1.0

julia> backproject(arr, [0, π/2])
6×6 view(::Array{Float64, 3}, :, :, 1) with eltype Float64:
 0.0  0.0  0.0       0.0  0.0       0.0
 0.0  0.0  0.0       0.0  0.0       0.0
 0.0  0.0  2.0       1.0  2.0       0.732051
 0.0  0.0  1.0       0.0  1.0       0.0
 0.0  0.0  2.0       1.0  2.0       0.732051
 0.0  0.0  0.732051  0.0  0.732051  0.0

julia> arr = ones((2,1)); 

julia> backproject(arr, [0], geometry=RadonFlexibleCircle(10, [-3, 3], [0,0]))
10×10 view(::Array{Float64, 3}, :, :, 1) with eltype Float64:
 0.0  0.0  0.0       0.0       0.0        0.0      0.0        0.0       0.0       0.0
 0.0  0.0  0.0       0.0       0.0        0.0      0.0        0.0       0.0       0.0
 0.0  0.0  0.0       0.0       0.335172   1.49876  0.335172   0.0       0.0       0.0
 0.0  0.0  0.0       0.0       1.08455    0.0      1.08455    0.0       0.0       0.0
 0.0  0.0  0.0       0.0       1.08455    0.0      1.08455    0.0       0.0       0.0
 0.0  0.0  0.0       1.00552   0.0790376  0.0      0.0790376  1.00552   0.0       0.0
 0.0  0.0  0.0       1.08455   0.0        0.0      0.0        1.08455   0.0       0.0
 0.0  0.0  0.591307  0.493247  0.0        0.0      0.0        0.493247  0.591307  0.0
 0.0  0.0  0.700352  0.0       0.0        0.0      0.0        0.0       0.700352  0.0
 0.0  0.0  0.0       0.0       0.0        0.0      0.0        0.0       0.0       0.0

backproject_filtered(sinogram, θs; 
                        geometry, μ, filter)

Reconstruct the image from the sinogram using the filtered backprojection algorithm.

• filter=nothing: The filter to be applied in Fourier space. If nothing, a ramp filter is used. filter should be a 1D array of the same length as the sinogram.

See radon for the explanation of the keyword arguments

abstract type RadonGeometry end

Abstract supertype for all geometries which are supported by radon and backproject.

List of geometries:

• RadonParallelCircle
• RadonFlexibleCircle

See radon and backproject how to apply.

RadonParallelCircle(N, in_height, weights)

• N is the size of the first and second dimension of the input array for radon.

• in_height is a vector or a range indicating the positions of the detector with respect to the midpoint which is located at N ÷ 2 + 1. The rays travel along straight parallel paths through the array. Maximum and minimum values are (N+1) ÷ 2 - 1 and -(N+1) ÷ 2 + 1 respectively.

For example, for an array of size N=10 the default definition is: RadonParallelCircle(10, -4:4) So the resulting sinogram has the shape: (9, length(angles), size(array, 3)).

• weights can weight each of the rays with different strength. Per default weights = 0 .* in_height .+ 1

See radon and backproject how to apply.

RadonFlexibleCircle(N, in_height, out_height, weights)

• N is the size of the first and second dimension of the input for radon.
• in_height is a vector or range indicating the vertical positions of the rays entering the circle with respect to the midpoint which is located at N ÷ 2 + 1. Maximum and minimum values are (N+1) ÷ 2 - 1 and -(N+1) ÷ 2 + 1 respectively.
• out_height is a vector or range indicating the vertical positions of the rays exiting the circle with respect to the midpoint which is located at N ÷ 2 + 1. Maximum and minimum values are (N+1) ÷ 2 - 1 and -(N+1) ÷ 2 + 1 respectively.

One definition could be: RadonFlexibleCircle(10, -4:4, zeros((9,))) It would describe rays which enter the circle at positions -4:4 but all of them would focus at the position 0 when leaving the circle. This is an extreme form of fan beam tomography.

• weights can weight each of the rays with different strength. Per default weights = 0 .* in_height .+ 1

See radon and backproject how to apply.