diff --git a/Project.toml b/Project.toml index 20cd372..9279ccc 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "RadonKA" uuid = "86de8297-835b-47df-b249-c04e8db91db5" authors = ["Felix Wechsler (roflmaostc) "] -version = "0.4.0" +version = "0.5.0" [deps] Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458" diff --git a/README.md b/README.md index 93fc141..98cb9c1 100644 --- 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 index 59d52e8..aa827f3 100644 --- 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 index c691d41..a41f67c 100644 --- 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 diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md index 86aea1a..9a3b87a 100644 --- 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 index 648e6e6..3d4f9ad 100644 --- 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 index 56c8fd0..3870b92 100644 --- 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 index be73b1f..0ca7f1e 100644 --- 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 index 9f5bf57..d14e7c1 100644 --- 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 index ed50903..e9eb68b 100644 --- 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 new file mode 100644 index 0000000..97c2429 --- /dev/null +++ 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 deleted file mode 100644 index 9ecabb9..0000000 --- a/src/filtered_backprojection.jl +++ /dev/null @@ -1,21 +0,0 @@ -export filtered_backprojection - -""" - filtered_backprojection(sinogram, θs; - geometry, μ) - -Calculates the simple Filtered Backprojection in CT with applying a ramp filter -in Fourier space. - -See [`radon`](@ref) for the explanation of the keyword arguments -""" -function filtered_backprojection(sinogram::AbstractArray{T}, θs::AbstractVector; - geometry=RadonParallelCircle(size(sinogram,1) + 1,-(size(sinogram,1))÷2:(size(sinogram,1))÷2), μ=nothing) 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 diff --git a/test/notebook.jl b/test/notebook.jl index f7810dd..ad2372c 100644 --- 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 index 99d7a76..3b341bf 100644 --- 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