Rename filtered_backprojection -> backproject_filtered

roflmaostc · Feb 9, 2024 · fcec9cc · fcec9cc · roflmaostc · Feb 9, 2024
1 parent 91cadf2
commit fcec9cc
Show file tree

Hide file tree

Showing 14 changed files with 51 additions and 36 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "RadonKA"
 uuid = "86de8297-835b-47df-b249-c04e8db91db5"
 authors = ["Felix Wechsler (roflmaostc) <fxw+git@mailbox.org>"]
-version = "0.4.0"
+version = "0.5.0"
 
 [deps]
 Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458"

diff --git a/README.md b/README.md
@@ -1,5 +1,5 @@
 # RadonKA.jl
-A simple yet sufficiently fast Radon and inverse Radon (backproject) transform implementation using [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl).
+A simple yet sufficiently fast Radon and adjoint Radon (backproject) transform implementation using [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl).
 It offers multithreading and CUDA support and outperforms any existing Julia Radon transforms (at least the ones we are aware of). 
 On CUDA it is faster much than Matlab and it offers the same or faster speed than ASTRA.
 

diff --git a/docs/src/functions.md b/docs/src/functions.md
@@ -6,7 +6,7 @@ radon
 ## IRadon
 ```@docs
 backproject
-filtered_backprojection
+backproject_filtered
 ```
 
 ## Specifying Geometries

diff --git a/docs/src/index.md b/docs/src/index.md
@@ -1,5 +1,5 @@
 # RadonKA.jl
-A simple yet sufficiently fast Radon and inverse Radon (backproject) transform implementation using [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl).
+A simple yet sufficiently fast Radon and adjoint Radon (backproject) transform implementation using [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl).
 
 ```@raw html
 <a  href="assets/RadonKA_logo.png"><img src="assets/RadonKA_logo.png"  width="200"></a>

diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md
@@ -24,7 +24,7 @@ angles = range(0f0, 2f0π, 500)[begin:end-1]
 ```
 ![](../assets/sinogram.png)
 
-# inverse Radon (backproject) transform
+# adjoint Radon (backproject) transform
 ```julia
 # 0.268649 seconds (147 allocations: 1.015 MiB)
 @time backproject = RadonKA.backproject(sinogram, angles);
@@ -42,7 +42,7 @@ Left is the original sample and right the simple backprojection.
 In the absence of noise, the filtered backprojection works well:
 ```julia
 #   0.252915 seconds (171 allocations: 13.664 MiB)
-@time filtered_backproject = RadonKA.filtered_backprojection(sinogram, angles);
+@time filtered_backproject = RadonKA.backproject_filtered(sinogram, angles);
 ```
 ![](../assets/filtered.png)
 

diff --git a/examples/0_example_radon_iradon.jl b/examples/0_example_radon_iradon.jl
@@ -124,7 +124,7 @@ simshow(Array(backproject_cu[:, :, i_z3]))
 md"# Filtered Backprojection"
 
 # ╔═╡ 32b15077-5e09-4693-8f12-3b2029fe63cc
-@mytime filtered_bproj = RadonKA.filtered_backprojection(sinogram_c, togoc(angles));
+@mytime filtered_bproj = RadonKA.backproject_filtered(sinogram_c, togoc(angles));
 
 # ╔═╡ bc6e2d40-fcd1-4d7b-8f96-3d4d9e4336de
 @bind i_z4 Slider(1:size(sinogram, 3), show_value=true)

diff --git a/examples/2_CT_with_optimizer.jl b/examples/2_CT_with_optimizer.jl
@@ -70,7 +70,7 @@ simshow(measurement)
 md"""# Simple Backprojection"""
 
 # ╔═╡ 464f022d-0773-43dd-81a7-ca2f6fc91634
-img_bp = filtered_backprojection(measurement, angles);
+img_bp = backproject_filtered(measurement, angles);
 
 # ╔═╡ 49e59001-0e8d-4872-ae50-47c38486b3fd
 img_backproject = backproject(measurement, angles);

diff --git a/examples/99_docs_example.jl b/examples/99_docs_example.jl
@@ -42,7 +42,7 @@ simshow(sinogram)
 [simshow(img) simshow(backproject)]
 
 # ╔═╡ ebd1170c-84ea-420d-a2e1-9afe0acb637b
-@time filtered_backproject = RadonKA.filtered_backprojection(sinogram, angles);
+@time filtered_backproject = RadonKA.backproject_filtered(sinogram, angles);
 
 # ╔═╡ 00756bf0-e309-4991-84a8-4afed5cdfa93
 simshow(filtered_backproject)

diff --git a/src/RadonKA.jl b/src/RadonKA.jl
@@ -11,7 +11,7 @@ include("geometry.jl")
 include("utils.jl")
 include("radon.jl")
 include("backproject.jl")
-include("filtered_backprojection.jl")
+include("backproject_filtered.jl")
 
 
 # PrecompileTools
@@ -22,10 +22,10 @@ include("filtered_backprojection.jl")
     @compile_workload begin
         r = radon(Float32.(img), Float32.(angles)) 
         backproject(r, Float32.(angles)) 
-        RadonKA.filtered_backprojection(r, angles)
+        RadonKA.backproject_filtered(r, angles)
         r = radon(img, angles) 
         backproject(r, angles) 
-        RadonKA.filtered_backprojection(r, angles)
+        RadonKA.backproject_filtered(r, angles)
 
         r = radon(Float32.(img), Float32.(angles), μ=0.1f0) 
         backproject(r, Float32.(angles), μ=0.1f0) 

diff --git a/src/backproject.jl b/src/backproject.jl
@@ -13,6 +13,8 @@ end
 Conceptually the adjoint operation of [`radon`](@ref). Intuitively, the `backproject` smears rays back into the space.
 See also [`radon`](@ref).
 
+For filtered backprojection see [`backproject_filtered`](@ref).
+
 # Example
 ```jldoctest
 julia> arr = zeros((5,2)); arr[2,:] .= 1; arr[4, :] .= 1

diff --git a/src/backproject_filtered.jl b/src/backproject_filtered.jl
@@ -0,0 +1,34 @@
+export backproject_filtered
+
+"""
+    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`](@ref) for the explanation of the keyword arguments
+"""
+function backproject_filtered(sinogram::AbstractArray{T}, θs::AbstractVector;
+        geometry=RadonParallelCircle(size(sinogram,1) + 1,-(size(sinogram,1))÷2:(size(sinogram,1))÷2), μ=nothing,
+        filter=nothing) where {T}
+    _backproject(sinogram, θs, geometry, μ, filter, T)
+end
+
+
+function _backproject(sinogram, θs, geometry, μ, filter::Nothing, ::Type{T}) where {T}
+    filter = similar(sinogram, (size(sinogram, 1),))
+    filter .= rr(T, (size(sinogram, 1), )) 
+
+    p = plan_fft(sinogram, (1,))
+    sinogram = real(inv(p) * (p * sinogram .* ifftshift(filter)))
+    return backproject(sinogram, θs; geometry, μ)
+end
+
+
+function _backproject(sinogram, θs, geometry, μ, filter::AbstractArray, ::Type{T}) where {T}
+    p = plan_fft(sinogram, (1,))
+    sinogram = real(inv(p) * (p * sinogram .* ifftshift(filter)))
+    return backproject(sinogram, θs; geometry, μ)
+end
diff --git a/src/filtered_backprojection.jl b/src/filtered_backprojection.jl
diff --git a/test/notebook.jl b/test/notebook.jl
@@ -237,7 +237,7 @@ simshow(I_2[:,:,1])
 simshow(sinogram2)
 
 # ╔═╡ ef6de020-ba4a-4dce-8a77-f706c1910270
-array_filtered = filtered_backprojection(sinogram2, angles2)
+array_filtered = backproject_filtered(sinogram2, angles2)
 
 # ╔═╡ ecb578ff-4dca-4260-9124-0bc17499d071
 array_backproject = backproject(sinogram2, angles2)

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -94,12 +94,12 @@ using Zygote
         angles2 = range(0, π, 200);
 
         sinogram2 = radon(array3, angles2)
-        array_filtered = filtered_backprojection(sinogram2, angles2)
+        array_filtered = backproject_filtered(sinogram2, angles2)
         @test ≈(array_filtered / sum(array_filtered) .+ 0.1, array3 / sum(array3) .+ 0.1, rtol=0.06)
 
         geometry = RadonParallelCircle(32, -15:0.1:15)
         sinogram2 = radon(array3, angles2; geometry)
-        array_filtered = filtered_backprojection(sinogram2, angles2; geometry)
+        array_filtered = backproject_filtered(sinogram2, angles2; geometry)
         @test ≈(array_filtered[5:28, 5:28] / sum(array_filtered[5:28, 5:28]) .+ 0.1, array3[5:28, 5:28] / sum(array3[5:28, 5:28]) .+ 0.1, rtol=0.05)
     end