iradon: threading, interpolations, and tests (#3)

add threading to iradon. add tests.

Author: kczimm

Project.toml

@@ -4,12 +4,14 @@ version = "0.1.0"
 
 [deps]
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 
 [compat]
 julia = "1"
 
 [extras]
+Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["Images", "Test"]

src/ImageReconstruction.jl

@@ -1,100 +1,104 @@
 module ImageReconstruction
 
+using Base.Threads
+using Interpolations
 using FFTW
 
 export radon, iradon
 
 """
-radon transform
+    radon(I::AbstractMatrix, θ::AbstractRange, t::AbstractRange)
 
-https://en.wikipedia.org/wiki/Radon_transform
+Radon transform of a image `I` producing a sinogram from view angles `θ` in
+radians and detector sampling `t`.
 """
-function radon end
-
-"""
-inverse radon transform
-
-https://en.wikipedia.org/wiki/Radon_transform
-"""
-function iradon end
-
-
-"""
-Radon transform of an image using a pixel-driven algorithm.
-"""
-function radon(image::AbstractMatrix, θ::AbstractRange, t::AbstractRange)
-    P = zeros(eltype(image), length(t), length(θ))
-    nr, nc = size(image)
-
-    for i in 1:nr, j in 1:nc
-        x = j - nc / 2 + 0.5
-        y = i - nr / 2 + 0.5
+function radon(I::AbstractMatrix, θ::AbstractRange, t::AbstractRange)
+    P = zeros(eltype(I), length(t), length(θ))
+    ax1, ax2 = axes(I)
+
+    nax1, nax2 = length(ax1), length(ax2)
+    for j in ax2, i in ax1
+        x = j - nax2 / 2 + 0.5
+        y = i - nax1 / 2 + 0.5
         @inbounds for (k, θₖ) in enumerate(θ)
             t′ = x * cos(θₖ) + y * sin(θₖ)
 
             a = convert(Int, round((t′ - minimum(t)) / step(t) + 1))
 
             (a < 1 || a > length(t)) && continue
             α = abs(t′ - t[a])
-            P[a, k] += (1 - α) * image[i, j]
+
+            I′ = I[i, j]
+            P[a, k] += (1 - α) * I′
 
             (a > length(t) + 1) && continue
-            P[a + 1, k] += α * image[i, j]
+            P[a+1, k] += α * I′
         end
     end
 
     P
 end
 
-function _ramp_spatial(N::Int, τ)
-    h = zeros(N)
+function _ramp_spatial(N::Int, τ, Npad::Int = N)
+    @assert Npad ≥ N
     N2 = N ÷ 2
-    for i in eachindex(h)
-        n = i - N2 - 1
-        if mod(n, 2) != 0
-            h[i] = -1 / (π * n * τ)^2
-        elseif n == 0
-            h[i] = 1 / (4 * τ^2)
-        end
-    end
-    h
+    hval(n) = n == 0 ? 1 / (4*τ^2) : - mod(n, 2)/(π * n * τ)^2
+    [i ≤ N ? hval(i - N2 - 1) : 0. for i = 1:Npad]
 end
 
-_zero_pad(p::AbstractVector, N::Int) = vcat(p, zeros(N))
-
 """
-Inverse radon transform using a ramp filter and a pixel-driven algorithm.
+	iradon(P::AbstractMatrix, θ::AbstractRange, t::AbstractRange)
+
+Inverse radon transform of a sinogram `P` with view angles `θ` in radians and
+detector sampling `t` producing an image on a 128x128 matrix.
 """
-function iradon(sinogram::AbstractMatrix, θ::AbstractRange, t::AbstractRange)
+function iradon(P::AbstractMatrix, θ::AbstractRange, t::AbstractRange)
     pixels = 128
 
     N = length(t)
     K = length(θ)
     Npad = nextpow(2, 2 * N - 1)
     τ = step(t)
-    ramp = fft(_zero_pad(_ramp_spatial(N, τ), Npad - N))
-    i = div(N, 2) + 1
+    ramp = fft(_ramp_spatial(N, τ, Npad))
+    i = N ÷ 2 + 1
     j = i + N - 1
 
-    image = zeros(eltype(sinogram), pixels, pixels)
-    Q = Vector{eltype(sinogram)}(undef, N)
-    for (k, θₖ) in enumerate(θ)
+    T = eltype(P)
+    I = [zeros(T, pixels, pixels) for _ = 1:nthreads()]
+    P′ = [Vector{Complex{T}}(undef, Npad) for _ = 1:nthreads()]
+    Q = [Vector{T}(undef, N) for _ = 1:nthreads()]
+
+    l = SpinLock()
+    @inbounds @threads for (k, θₖ) in collect(enumerate(θ))
+        id = threadid()
+        Pid, Qid, Iid = P′[id], Q[id], I[id]
+
         # filter projection
-        Q[:] .= τ .* real.(ifft(fft(_zero_pad(sinogram[:, k], Npad - N)) .* ramp)[i:j])
+        Pid[1:N] .= P[:, k]
+        Pid[N+1:end] .= 0
+
+        # Need to prevent multiple thread execution during fft/ifft.
+        # double free occurs otherwise.
+        # https://github.com/JuliaMath/FFTW.jl/issues/134
+        lock(l)
+        fft!(Pid)
+        Pid .*= ramp
+        ifft!(Pid)
+        unlock(l)
+        Qid .= τ .* real.(@view Pid[i:j])
+
+        Qₖ = LinearInterpolation(t, Qid)
 
         # backproject
-        for c in CartesianIndices(image)
+        for c in CartesianIndices(first(I))
             x = c.I[2] - pixels ÷ 2 + 0.5
             y = c.I[1] - pixels ÷ 2 + 0.5
+            x^2 + y^2 ≥ pixels^2 / 4 && continue
             t′ = x * cos(θₖ) + y * sin(θₖ)
-            # linear interpolation
-            #   image pixel (x, y) at θₖ is projected on t[z] and t[z+1]
-            z = convert(Int, round((t′ - minimum(t)) / step(t) + 1))
-            α = abs(t′ - t[z])
-            image[c] += (1 - α) * Q[z] + α * Q[z + 1]
+            Iid[c] += Qₖ(t′)
         end
     end
-    @. image * π / K
+    sum(I) .* π ./ K
 end
 
 end # module

test/runtests.jl

@@ -1,7 +1,7 @@
 using ImageReconstruction
 using Test
 
-@testset "radon: pixel-driven" begin
+@testset "radon - centered impulse response" begin
     pixels = 129
     views = 200
     I = zeros(pixels, pixels)
@@ -13,3 +13,14 @@ using Test
 
     @test all(P[101, :] .== 1)
 end
+
+using Images: shepp_logan, sad
+@testset "iradon - shepp logan" begin
+    Igt = shepp_logan(128)
+    views = 200
+    θ = range(0, 2π, length=views)
+    t = -150:150
+    P = radon(Igt, θ, t)
+    I = iradon(P, θ, t)
+	@test sad(I, Igt) < 500
+end