radon(img, θs; <kwargs>)

Calculates the parallel Radon transform of the AbstractArray I. 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 equally sized. 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.

Keywords

μ=nothing - Exponential IRadon 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.

geometry=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.

Please note: the implementation is not quite optimized for cache efficiency and it is a very naive algorithm. But still, it is fast!

See also iradon.

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

iradon(sinogram, θs; <kwargs>)

See also radon.

Example

# Examples

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> iradon(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.1        0.0  0.1        0.0
 0.0  0.1  0.2        0.1  0.2        0.0232051
 0.0  0.0  0.1        0.0  0.1        0.0
 0.0  0.1  0.2        0.1  0.2        0.0232051
 0.0  0.0  0.0232051  0.0  0.0232051  0.0

filtered_backprojection(sinogram, θs; 
                        geometry, μ)

Calculates the simple Filtered Backprojection in CT with applying a ramp filter in Fourier space.

See radon for the explanation of the keyword arguments

abstract type RadonGeometry end

List of geometries:

• RadonParallelCircle
• RadonFlexibleCircle

RadonParallelCircle(N, in_height)

N is the size of the first and second dimension of the array. in_height is a vector or range indicating the positions of the detector with respect to the midpoint which is located at N ÷ 2 + 1.

So 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))

RadonFlexibleCircle(N, in_height, out_height)

N is the size of the first and second dimension of the array. 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. 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.

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 cone beam tomography.