Inverse radon (#2)

pixel-driven inverse radon transform

Author: kczimm

Project.toml

@@ -2,6 +2,9 @@ name = "ImageReconstruction"
 uuid = "d5d0ac58-a118-45a5-925e-339d09fa0c5f"
 version = "0.1.0"
 
+[deps]
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+
 [compat]
 julia = "1"
 

docs/make.jl

@@ -1,17 +1,12 @@
 using Documenter, ImageReconstruction
 
 makedocs(;
-    modules=[ImageReconstruction],
-    format=Documenter.HTML(),
-    pages=[
-        "Home" => "index.md",
-    ],
-    repo="https://github.com/kczimm/ImageReconstruction.jl/blob/{commit}{path}#L{line}",
-    sitename="ImageReconstruction.jl",
-    authors="Kevin C. Zimmerman",
-    assets=String[],
+    modules = [ImageReconstruction],
+    format = Documenter.HTML(),
+    pages = ["Home" => "index.md"],
+    repo = "https://github.com/JuliaImages/ImageReconstruction.jl/blob/{commit}{path}#L{line}",
+    sitename = "ImageReconstruction.jl",
+    assets = String[],
 )
 
-deploydocs(;
-    repo="github.com/kczimm/ImageReconstruction.jl",
-)
+deploydocs(; repo = "github.com/JuliaImages/ImageReconstruction.jl")

src/ImageReconstruction.jl

@@ -1,20 +1,22 @@
 module ImageReconstruction
 
+using FFTW
+
 export radon, iradon
 
 """
 radon transform
 
 https://en.wikipedia.org/wiki/Radon_transform
 """
-radon
+function radon end
 
 """
 inverse radon transform
 
 https://en.wikipedia.org/wiki/Radon_transform
 """
-iradon
+function iradon end
 
 
 """
@@ -24,24 +26,75 @@ function radon(image::AbstractMatrix, θ::AbstractRange, t::AbstractRange)
     P = zeros(eltype(image), length(t), length(θ))
     nr, nc = size(image)
 
-    for i = 1:nr, j = 1:nc
+    for i in 1:nr, j in 1:nc
         x = j - nc / 2 + 0.5
         y = i - nr / 2 + 0.5
         @inbounds for (k, θₖ) in enumerate(θ)
             t′ = x * cos(θₖ) + y * sin(θₖ)
 
-            a = convert(Int, round((t′ - t.start) / step(t) + 1))
+            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]
 
             (a > length(t) + 1) && continue
-            P[a+1, k] += α * image[i, j]
+            P[a + 1, k] += α * image[i, j]
         end
     end
 
     P
 end
 
+function _ramp_spatial(N::Int, τ)
+    h = zeros(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
+end
+
+_zero_pad(p::AbstractVector, N::Int) = vcat(p, zeros(N))
+
+"""
+Inverse radon transform using a ramp filter and a pixel-driven algorithm.
+"""
+function iradon(sinogram::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
+    j = i + N - 1
+
+    image = zeros(eltype(sinogram), pixels, pixels)
+    Q = Vector{eltype(sinogram)}(undef, N)
+    for (k, θₖ) in enumerate(θ)
+        # filter projection
+        Q[:] .= τ .* real.(ifft(fft(_zero_pad(sinogram[:, k], Npad - N)) .* ramp)[i:j])
+
+        # backproject
+        for c in CartesianIndices(image)
+            x = c.I[2] - pixels ÷ 2 + 0.5
+            y = c.I[1] - pixels ÷ 2 + 0.5
+            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]
+        end
+    end
+    @. image * π / K
+end
+
 end # module