diff --git a/Project.toml b/Project.toml index 572a561..c7820e1 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.3.2" +version = "0.3.3" [deps] Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458" diff --git a/README.md b/README.md index 81b9db3..7b2c6da 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,5 @@ # RadonKA.jl -A simple yet sufficiently fast Radon and inverse Radon (iradon) transform implementation using [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl). +A simple yet sufficiently fast Radon and inverse 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. @@ -14,13 +14,13 @@ On CUDA it is faster much than Matlab and it offers the same or faster speed tha # Quick Overview * [x] For 2D and 3D arrays -* [x] parallel `radon` and `iradon` (`?RadonParallelCircle`) -* [x] attenuated `radon` and `iradon` (see the parameter `μ`) and see this [paper](https://iopscience.iop.org/article/10.1088/0266-5611/17/1/309/meta) as reference) +* [x] parallel `radon` and `backproject` (`?RadonParallelCircle`) +* [x] attenuated `radon` and `backproject` (see the parameter `μ`) and see this [paper](https://iopscience.iop.org/article/10.1088/0266-5611/17/1/309/meta) as reference) * [x] arbitrary 2D geometries where starting and endpoint of each ray can be specified (fan beam could be a special case if this) (`?RadonFlexibleCircle`) * [x] different strength weighting of rays * [x] based on [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl) * [x] tested on `CPU()` and `CUDABackend()` -* [x] registered adjoint rules for both `radon` and `iradon` +* [x] registered adjoint rules for both `radon` and `backproject` * [x] high performance however not ultra high performance. On par with ASTRA, on CUDA faster than Matlab. * [x] simple and extensible API @@ -40,7 +40,7 @@ angles = range(0f0, 2f0π, 500)[begin:end-1] # 0.085398 seconds (260 allocations: 1.006 MiB) @time sinogram = radon(img, angles); # 0.127043 seconds (251 allocations: 1.036 MiB) -@time backproject = RadonKA.iradon(sinogram, angles); +@time backproject = RadonKA.backproject(sinogram, angles); simshow(sinogram) simshow(backproject) @@ -50,10 +50,10 @@ img_c = CuArray(img) # 0.003363 seconds (244 CPU allocations: 18.047 KiB) (7 GPU allocations: 1007.934 KiB, 0.96% memmgmt time) CUDA.@time sinogram = radon(img_c, angles); # 0.005928 seconds (218 CPU allocations: 16.109 KiB) (7 GPU allocations: 1.012 MiB, 0.49% memmgmt time) -CUDA.@time backproject = RadonKA.iradon(sinogram, angles); +CUDA.@time backproject = RadonKA.backproject(sinogram, angles); ``` - + # Examples See the [documentation](https://roflmaostc.github.io/RadonKA.jl/dev/tutorial). @@ -69,8 +69,8 @@ julia> using Pluto; Pluto.run() ``` A browser should open. The following examples show case the ability of this package: -* Simple `radon` and `iradon`: [Pluto notebook](examples/0_example_radon_iradon.jl) -* Different geometries: [Pluto notebook](examples/0_example_radon_iradon.jl) +* Simple `radon` and `backproject`: [Pluto notebook](examples/0_example_radon_backproject.jl) +* Different geometries: [Pluto notebook](examples/0_example_radon_backproject.jl) * Reconstruction of a CT dataset with an optimizer: [Pluto notebook](examples/2_CT_with_optimizer.jl) * How this package is used in Tomographic Volumetric Additive Manufacturing (3D printing): [Pluto notebook](examples/4_Tomographic_Volumetric_Additive_Manufacturing_with_Refraction.jl) @@ -93,4 +93,4 @@ There exists [Sinograms.jl](https://github.com/JuliaImageRecon/Sinograms.jl) and Again, no arbitrary geometries can be specified. And also no attenuated Radon transform is possible. ## Matlab -Matlab has built-in a `radon` and `iradon` transform which is similar to our lightweight API. However, no CUDA acceleration, no 3D arrays and no attenuated Radon transform. +Matlab has built-in a `radon` and `backproject` transform which is similar to our lightweight API. However, no CUDA acceleration, no 3D arrays and no attenuated Radon transform. diff --git a/docs/src/assets/astra_iradon.png b/docs/src/assets/astra_backproject.png similarity index 100% rename from docs/src/assets/astra_iradon.png rename to docs/src/assets/astra_backproject.png diff --git a/docs/src/assets/matlab_iradon.png b/docs/src/assets/matlab_backproject.png similarity index 100% rename from docs/src/assets/matlab_iradon.png rename to docs/src/assets/matlab_backproject.png diff --git a/docs/src/assets/radonka_iradon.png b/docs/src/assets/radonka_backproject.png similarity index 100% rename from docs/src/assets/radonka_iradon.png rename to docs/src/assets/radonka_backproject.png diff --git a/docs/src/assets/sample_iradon.png b/docs/src/assets/sample_backproject.png similarity index 100% rename from docs/src/assets/sample_iradon.png rename to docs/src/assets/sample_backproject.png diff --git a/docs/src/benchmark.md b/docs/src/benchmark.md index a9c3973..d48944e 100644 --- a/docs/src/benchmark.md +++ b/docs/src/benchmark.md @@ -30,11 +30,11 @@ Tested on a AMD Ryzen 9 5900X 12-Core Processor with 24 Threads and a NVIDIA GeF sinogram = radon(sample_2D, angles); @btime sinogram = radon($sample_2D, $angles); simshow(sinogram) - @btime sample_iradon = iradon($sinogram, $angles); + @btime sample_backproject = backproject($sinogram, $angles); @btime CUDA.@sync sinogram_c = radon($sample_2D_c, $angles_c); sinogram_c = radon(sample_2D_c, angles_c); - @btime CUDA.@sync sample_iradon_c = iradon($sinogram_c, $angles_c); + @btime CUDA.@sync sample_backproject_c = backproject($sinogram_c, $angles_c); sample_3D = randn(Float32, (512, 512, 100)); @@ -42,14 +42,14 @@ Tested on a AMD Ryzen 9 5900X 12-Core Processor with 24 Threads and a NVIDIA GeF sinogram_3D = radon(sample_3D, angles); @btime radon($sample_3D, $angles) - @btime iradon($sinogram_3D, $angles) + @btime backproject($sinogram_3D, $angles) sinogram_3D_c = radon(sample_3D_c, angles_c) @btime CUDA.@sync radon($sample_3D_c, $angles_c) - @btime CUDA.@sync iradon($sinogram_3D_c, $angles_c) + @btime CUDA.@sync backproject($sinogram_3D_c, $angles_c) ``` ![](../assets/radonka_sinogram.png) -![](../assets/radonka_iradon.png) +![](../assets/radonka_backproject.png) ## Matlab (R2023a) @@ -65,10 +65,10 @@ R = R / max(R(:)); imwrite(R, "/tmp/matlab_sinogram.png"); tic; -iR = iradon(R, theta, "linear", "none"); +iR = backproject(R, theta, "linear", "none"); toc; iR = iR / max(iR(:)); -imwrite(iR, "/tmp/matlab_iradon.png"); +imwrite(iR, "/tmp/matlab_backproject.png"); @@ -87,13 +87,13 @@ toc; tic; for i = 1:size(y, 1) - ix(i, :, :) = iradon(squeeze(y(i, :, :)), theta); + ix(i, :, :) = backproject(squeeze(y(i, :, :)), theta); end toc; ``` ![](../assets/matlab_sinogram.png) -![](../assets/matlab_iradon.png) +![](../assets/matlab_backproject.png) ## Astra Some of the benchmarks did not run properly or were providing non-sense. @@ -174,5 +174,5 @@ np.copy(rec) ``` ![](../assets/astra_sinogram.png) -![](../assets/astra_iradon.png) +![](../assets/astra_backproject.png) diff --git a/docs/src/functions.md b/docs/src/functions.md index 1749ac7..59d52e8 100644 --- a/docs/src/functions.md +++ b/docs/src/functions.md @@ -5,7 +5,7 @@ radon ## IRadon ```@docs -iradon +backproject filtered_backprojection ``` diff --git a/docs/src/geometries.md b/docs/src/geometries.md index 5a39b1c..79c0fe4 100644 --- a/docs/src/geometries.md +++ b/docs/src/geometries.md @@ -23,7 +23,7 @@ sinogram[1:5:end] .= 1 geometry_parallel = RadonParallelCircle(N, -(N-1)÷2:(N-1)÷2) -projection_parallel = iradon(sinogram, angles; geometry=geometry_parallel); +projection_parallel = backproject(sinogram, angles; geometry=geometry_parallel); simshow(projection_parallel) ``` @@ -37,7 +37,7 @@ sinogram_small[1:3:end] .= 1 geometry_small = RadonParallelCircle(200, -49:49) -projection_small = iradon(sinogram_small, angles; geometry=geometry_small); +projection_small = backproject(sinogram_small, angles; geometry=geometry_small); simshow(projection_small) ``` @@ -54,7 +54,7 @@ The second range indicates the position upon exit of the circle. ```julia geometry_fan = RadonFlexibleCircle(N, -(N-1)÷2:(N-1)÷2, range(-(N-1)÷4, (N-1)÷4, N-1)) -projected_fan = iradon(sinogram, angles; geometry=geometry_fan); +projected_fan = backproject(sinogram, angles; geometry=geometry_fan); simshow(projected_fan, γ=0.01) ``` @@ -63,7 +63,7 @@ simshow(projected_fan, γ=0.01) ```julia geometry_extreme = RadonFlexibleCircle(N, -(N-1)÷2:(N-1)÷2, zeros((199,))) -projected_extreme = iradon(sinogram, angles; geometry=geometry_extreme); +projected_extreme = backproject(sinogram, angles; geometry=geometry_extreme); simshow(projected_extreme, γ=0.01) ``` @@ -74,7 +74,7 @@ For example, if in your application some rays are stronger than others you can i ```julia geometry_weight = RadonParallelCircle(N, -(N-1)÷2:(N-1)÷2, abs.(-(N-1)÷2:(N-1)÷2)) -projection_weight = iradon(sinogram, angles; geometry=geometry_weight); +projection_weight = backproject(sinogram, angles; geometry=geometry_weight); simshow(projection_weight) ``` @@ -85,7 +85,7 @@ simshow(projection_weight) The ray gets some attenuation with `exp(-μ*x)` where `x` is the distance traveled to the entry point of the circle. `μ` is in units of pixel. ```julia -projected_exp = iradon(sinogram, angles; geometry=geometry_extreme, μ=0.04); +projected_exp = backproject(sinogram, angles; geometry=geometry_extreme, μ=0.04); simshow(projected_exp) ``` diff --git a/docs/src/index.md b/docs/src/index.md index a3be4de..c691d41 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 (iradon) transform implementation using [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl). +A simple yet sufficiently fast Radon and inverse Radon (backproject) transform implementation using [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl). ```@raw html @@ -11,13 +11,13 @@ A simple yet sufficiently fast Radon and inverse Radon (iradon) transform implem # Quick Overview * [x] For 2D and 3D arrays -* [x] parallel `radon` and `iradon` (`?RadonParallelCircle`) -* [x] attenuated `radon` and `iradon` (see the parameter `μ`) +* [x] parallel `radon` and `backproject` (`?RadonParallelCircle`) +* [x] attenuated `radon` and `backproject` (see the parameter `μ`) * [x] arbitrary 2D geometries where starting and endpoint of each ray can be specified (fan beam could be a special case if this) (`?RadonFlexibleCircle`) * [x] It is restricted to the incircle of radius `N ÷ 2 - 1` if the array has size `(N, N, N_z)` * [x] based on [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl) * [x] tested on `CPU()` and `CUDABackend` -* [x] registered adjoint rules for both `radon` and `iradon` +* [x] registered adjoint rules for both `radon` and `backproject` * [x] high performance however not ultra high performance * [x] simple API @@ -39,7 +39,7 @@ angles = range(0f0, 2f0π, 500)[begin:end-1] @time sinogram = radon(img, angles); # 0.268649 seconds (147 allocations: 1.015 MiB) -@time backproject = RadonKA.iradon(sinogram, angles); +@time backproject = RadonKA.backproject(sinogram, angles); simshow(sinogram) simshow(backproject) @@ -47,12 +47,12 @@ simshow(backproject) ```@raw html - + ``` # Examples See either the [documentation](https://roflmaostc.github.io/RadonKA.jl/dev/tutorial). -Otherwise, this [example](https://github.com/roflmaostc/RadonKA.jl/blob/main/examples/example_radon_iradon.jl) shows the main features, including CUDA support. +Otherwise, this [example](https://github.com/roflmaostc/RadonKA.jl/blob/main/examples/example_radon_backproject.jl) shows the main features, including CUDA support. There is one tutorial about [Gradient Descent optimization](https://github.com/roflmaostc/RadonKA.jl/blob/main/examples/CT_with_optimizer.jl). Another one covers how the Radon transform is used in [Volumetric Additive Manufacturing](https://github.com/roflmaostc/RadonKA.jl/blob/main/examples/volumetric_printing.jl). diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md index c7c892d..86aea1a 100644 --- a/docs/src/tutorial.md +++ b/docs/src/tutorial.md @@ -24,10 +24,10 @@ angles = range(0f0, 2f0π, 500)[begin:end-1] ``` ![](../assets/sinogram.png) -# inverse Radon (iradon) transform +# inverse Radon (backproject) transform ```julia # 0.268649 seconds (147 allocations: 1.015 MiB) -@time backproject = RadonKA.iradon(sinogram, angles); +@time backproject = RadonKA.backproject(sinogram, angles); # in Pluto or Jupyter simshow(sinogram) @@ -35,7 +35,7 @@ simshow(sinogram) [simshow(img) simshow(backproject)] ``` Left is the original sample and right the simple backprojection. -![](../assets/sample_iradon.png) +![](../assets/sample_backproject.png) # Filtered Backprojection diff --git a/examples/0_example_radon_iradon.jl b/examples/0_example_radon_iradon.jl index b7a6bb0..648e6e6 100644 --- a/examples/0_example_radon_iradon.jl +++ b/examples/0_example_radon_iradon.jl @@ -103,16 +103,16 @@ simshow(Array(sinogram_c[:, :, i_z])) md"# IRadon Transform" # ╔═╡ 7e27da5a-1b04-4d4c-8c62-eaffa7f4f9ce -@time backproject = RadonKA.iradon(sinogram, angles); +@time backproject = RadonKA.backproject(sinogram, angles); # ╔═╡ 4e367035-eb2f-4dfa-9646-7c182a111c49 -md"Use this slider to add more and more angles to the iradon transform" +md"Use this slider to add more and more angles to the backproject transform" # ╔═╡ b9bf49a0-7320-4269-9a6a-ac2533ab5fde @bind angle_limit Slider(1:size(sinogram, 2), default=size(sinogram, 2), show_value=true) # ╔═╡ 93a7ab4a-b2dc-4fc2-bf69-66e6d615103f -CUDA.@time CUDA.@sync backproject_cu = RadonKA.iradon(sinogram_c[:, begin:angle_limit, :], togoc(angles[begin:angle_limit])); +CUDA.@time CUDA.@sync backproject_cu = RadonKA.backproject(sinogram_c[:, begin:angle_limit, :], togoc(angles[begin:angle_limit])); # ╔═╡ 52a86ed8-4504-4d9e-9ea6-6aeaf8540406 @bind i_z3 Slider(1:size(sinogram, 3), show_value=true) diff --git a/examples/1_documentation_different_geometries.jl b/examples/1_documentation_different_geometries.jl index 81f12c5..a036532 100644 --- a/examples/1_documentation_different_geometries.jl +++ b/examples/1_documentation_different_geometries.jl @@ -38,7 +38,7 @@ end geometry_parallel = RadonParallelCircle(N, -(N-1)÷2:(N-1)÷2) # ╔═╡ 07bec511-7802-45e4-b331-0d35e8f850c1 -projection_parallel = iradon(sinogram, angles; geometry=geometry_parallel); +projection_parallel = backproject(sinogram, angles; geometry=geometry_parallel); # ╔═╡ 9e5a58be-0de5-4820-bc3b-6dff931271a9 simshow(projection_parallel) @@ -56,7 +56,7 @@ end geometry_small = RadonParallelCircle(200, -49:49) # ╔═╡ 27d5aac8-2528-46f2-a0dd-91c3fe5d275c -projection_small = iradon(sinogram_small, angles; geometry=geometry_small); +projection_small = backproject(sinogram_small, angles; geometry=geometry_small); # ╔═╡ 29c9d1ff-694c-4fff-b551-836cc9eaf347 simshow(projection_small) @@ -68,7 +68,7 @@ md"# Similar to fan Beam Tomography" geometry_fan = RadonFlexibleCircle(N, -(N-1)÷2:(N-1)÷2, range(-(N-1)÷4, (N-1)÷4, N-1)) # ╔═╡ 37d760fa-74e6-47d1-b8e6-3315c1747b4c -projected_fan = iradon(sinogram, angles; geometry=geometry_fan); +projected_fan = backproject(sinogram, angles; geometry=geometry_fan); # ╔═╡ 878121c5-4ad5-477b-88c7-f53df7510052 simshow(projected_fan, γ=0.01) @@ -80,7 +80,7 @@ md"# Extreme fan Beam Tomography" geometry_extreme = RadonFlexibleCircle(N, -(N-1)÷2:(N-1)÷2, zeros((199,))) # ╔═╡ 0756e11b-35ca-4eef-be48-52b8b5c098cd -projected_extreme = iradon(sinogram, angles; geometry=geometry_extreme); +projected_extreme = backproject(sinogram, angles; geometry=geometry_extreme); # ╔═╡ 1e22032d-d57d-42d3-a9c5-9f2e94531346 simshow(projected_extreme, γ=0.01) @@ -94,7 +94,7 @@ For example, if in your application some rays are stronger than others you can i geometry_weight = RadonParallelCircle(N, -(N-1)÷2:(N-1)÷2, abs.(-(N-1)÷2:(N-1)÷2)) # ╔═╡ 45b82171-c581-4fcf-9695-bc51464a2172 -projection_weight = iradon(sinogram, angles; geometry=geometry_weight); +projection_weight = backproject(sinogram, angles; geometry=geometry_weight); # ╔═╡ bbb51355-fa72-4064-b10e-46e29f4b2809 simshow(projection_weight) @@ -105,7 +105,7 @@ The ray gets some attenuation with `exp(-μ*x)` where `x` is the distance travel " # ╔═╡ b11cb860-3fab-4fc9-8921-bef32a8b7c12 -projected_exp = iradon(sinogram, angles; geometry=geometry_extreme, μ=0.04); +projected_exp = backproject(sinogram, angles; geometry=geometry_extreme, μ=0.04); # ╔═╡ 51561641-67a3-46d4-ae87-b896dc351a60 simshow(projected_exp) @@ -135,10 +135,10 @@ geometry_extreme2 = RadonFlexibleCircle(N2, -(N2-1)÷2:(N2-1)÷2, zeros((N2-1,)) simshow(sg_img) # ╔═╡ 29c25bd9-6359-489d-b0a3-ea767f832a02 -img_iradon = iradon(sg_img, angles2, geometry=geometry_extreme2, μ=0.008); +img_backproject = backproject(sg_img, angles2, geometry=geometry_extreme2, μ=0.008); # ╔═╡ 940d3b8b-1797-45f1-b493-02592a39eb19 -simshow(img_iradon) +simshow(img_backproject) # ╔═╡ c7dcc341-0afe-4913-b0a8-1869d8f819b4 md"# How to reconstruct? diff --git a/examples/2_CT_with_optimizer.jl b/examples/2_CT_with_optimizer.jl index 2416fb2..56c8fd0 100644 --- a/examples/2_CT_with_optimizer.jl +++ b/examples/2_CT_with_optimizer.jl @@ -73,10 +73,10 @@ md"""# Simple Backprojection""" img_bp = filtered_backprojection(measurement, angles); # ╔═╡ 49e59001-0e8d-4872-ae50-47c38486b3fd -img_iradon = iradon(measurement, angles); +img_backproject = backproject(measurement, angles); # ╔═╡ 6fac5606-2350-4f22-9f35-120936114d5a -[simshow(img_bp) simshow(img_iradon)] +[simshow(img_bp) simshow(img_backproject)] # ╔═╡ 1feccdec-cc35-4cc8-9a76-0e0e99bc7be3 md"# Optimization with gradient descent diff --git a/examples/4_Tomographic_Volumetric_Additive_Manufacturing_with_Refraction.jl b/examples/4_Tomographic_Volumetric_Additive_Manufacturing_with_Refraction.jl index ce9e78c..7bafe06 100644 --- a/examples/4_Tomographic_Volumetric_Additive_Manufacturing_with_Refraction.jl +++ b/examples/4_Tomographic_Volumetric_Additive_Manufacturing_with_Refraction.jl @@ -156,7 +156,7 @@ Revise.errors() geometry = RadonKA.RadonFlexibleCircle(size(target, 1), ray_startpoints, ray_endpoints) # ╔═╡ 85b6e54c-3256-4ba8-9889-56be82de0203 -simshow(Array(iradon(togoc(sinogram), togoc(range(0, 2π, kk + 1)[1:end-1]); geometry, μ=0.01)), γ=0.6) +simshow(Array(backproject(togoc(sinogram), togoc(range(0, 2π, kk + 1)[1:end-1]); geometry, μ=0.01)), γ=0.6) # ╔═╡ 5d55464f-c7d9-4e45-a1e7-907c2137be95 md"# Define Optimization algorithm @@ -198,7 +198,7 @@ function iter!(buffer, img, θs, μ; clip_sinogram=true, geometry) sinogram .= max.(sinogram, 0) end - img_recon = iradon(sinogram, θs; μ, geometry) + img_recon = backproject(sinogram, θs; μ, geometry) img_recon ./= maximum(img_recon) buffer .= max.(img_recon, 0) return buffer, sinogram @@ -224,13 +224,13 @@ function OSMO(img::AbstractArray{T}, θs, μ=nothing; guess[isobject] .+= max.(0, thresholds[2] .- tmp[isobject]) end - printed = iradon(s, θs; μ, geometry) + printed = backproject(s, θs; μ, geometry) printed ./= maximum(printed) return guess, s, printed end # ╔═╡ 12e7c6b7-fcf3-44e4-9d80-785704d93cb5 -fwd = x -> iradon(max.(0,x), angles; geometry) +fwd = x -> backproject(max.(0,x), angles; geometry) # ╔═╡ 942e1cf7-a47b-43af-a460-1ea896a2c9b8 f, g! = make_fg!(fwd, target, (0.65, 0.75)) @@ -266,19 +266,19 @@ simshow(Array(printed), cmap=:turbo) simshow(Array(fwd(res.minimizer)), cmap=:turbo) # ╔═╡ 7c998816-b4de-4dbf-977b-bdab9355be79 -simshow(Array(iradon(patterns[:, 96:96], togoc(angles[96:96]); μ=1/350f0, geometry)), cmap=:turbo) +simshow(Array(backproject(patterns[:, 96:96], togoc(angles[96:96]); μ=1/350f0, geometry)), cmap=:turbo) # ╔═╡ 28d8c378-b818-4d61-b9d4-a06ad7cc21e3 -simshow(Array(iradon(patterns[:, 306:306], togoc(angles[306:306]); μ=1/350f0, geometry)), cmap=:turbo) +simshow(Array(backproject(patterns[:, 306:306], togoc(angles[306:306]); μ=1/350f0, geometry)), cmap=:turbo) # ╔═╡ beb61bb1-3f7c-4301-92ca-5418bc0445cf -simshow(Array(iradon(patterns[:, 1:1], togoc(angles[1:1]); μ=1/350f0, geometry)), cmap=:turbo) +simshow(Array(backproject(patterns[:, 1:1], togoc(angles[1:1]); μ=1/350f0, geometry)), cmap=:turbo) # ╔═╡ 88895d8d-432e-4337-a51c-12ac4e4a4325 @bind i Slider(1:size(angles,1), show_value=true) # ╔═╡ 61048e10-8e23-41a6-8c49-0b1a8c0adc24 -simshow(Array(iradon(patterns[:, i:i], togoc(angles[i:i]); μ=1/350f0, geometry)), cmap=:turbo) +simshow(Array(backproject(patterns[:, i:i], togoc(angles[i:i]); μ=1/350f0, geometry)), cmap=:turbo) # ╔═╡ 21c0b369-a4b2-43ea-87b8-702474d88584 simshow(Array(patterns)) diff --git a/examples/99_docs_example.jl b/examples/99_docs_example.jl index 54d3981..be73b1f 100644 --- a/examples/99_docs_example.jl +++ b/examples/99_docs_example.jl @@ -36,7 +36,7 @@ angles = range(0f0, 2f0π, 500)[begin:end-1] simshow(sinogram) # ╔═╡ 451af1b2-c840-4ae8-935f-eb7965e7104b -@time backproject = RadonKA.iradon(sinogram, angles); +@time backproject = RadonKA.backproject(sinogram, angles); # ╔═╡ 42c42599-96b4-4a69-bbbd-5602cca700f4 [simshow(img) simshow(backproject)] diff --git a/examples/benchmarks/benchmark_REPL.jl b/examples/benchmarks/benchmark_REPL.jl index 491f7fe..7886f10 100644 --- a/examples/benchmarks/benchmark_REPL.jl +++ b/examples/benchmarks/benchmark_REPL.jl @@ -12,10 +12,10 @@ for sz in [(512,512,100), (1920, 1920)] if d === Array @btime radon($arr, $angles) - @btime iradon($iarr, $angles) + @btime backproject($iarr, $angles) else @btime CUDA.@sync radon($arr, $angles) - @btime CUDA.@sync iradon($iarr, $angles) + @btime CUDA.@sync backproject($iarr, $angles) end end end diff --git a/examples/benchmarks/benchmark_radon_iradon.jl b/examples/benchmarks/benchmark_radon_iradon.jl index 72399c6..21b5f74 100644 --- a/examples/benchmarks/benchmark_radon_iradon.jl +++ b/examples/benchmarks/benchmark_radon_iradon.jl @@ -60,13 +60,13 @@ simshow(sinogram) save("/tmp/radonka_sinogram.png", sinogram ./ maximum(sinogram)); # ╔═╡ a11fe70a-3425-4335-ba09-00efbae6c9af - sample_iradon = iradon(sinogram,angles); + sample_backproject = backproject(sinogram,angles); # ╔═╡ 29f7d88c-3bab-43a9-a253-d3bcfc98b706 -save("/tmp/radonka_iradon.png", sample_iradon ./ maximum(sample_iradon)); +save("/tmp/radonka_backproject.png", sample_backproject ./ maximum(sample_backproject)); # ╔═╡ da297929-f4f9-4fa2-b860-cfc159ac39f1 -@benchmark sample_iradon = iradon($sinogram, $angles) +@benchmark sample_backproject = backproject($sinogram, $angles) # ╔═╡ a61ea8f2-6329-4b6b-b4a2-1f5ac17bc618 sinogram_c = radon(sample_2D_c, angles_c, backend=CUDABackend()) @@ -75,7 +75,7 @@ sinogram_c = radon(sample_2D_c, angles_c, backend=CUDABackend()) @benchmark CUDA.@sync sinogram_c = radon($sample_2D_c, $angles_c, backend=CUDABackend()) # ╔═╡ 1b8e4e3f-d9db-4646-8ccf-840c21ecc328 -@benchmark CUDA.@sync sample_iradon_c = iradon($sinogram_c, $angles_c, backend=CUDABackend()) +@benchmark CUDA.@sync sample_backproject_c = backproject($sinogram_c, $angles_c, backend=CUDABackend()) # ╔═╡ c602a8df-d646-4f31-9d6f-902456f07ae5 md"# 3D" @@ -93,7 +93,7 @@ sample_3D_c = CuArray(sample_3D) @benchmark radon($sample_3D, $angles) # ╔═╡ b5825676-6600-4cbe-865b-788ef1567b85 -@benchmark iradon($sinogram_3D, $angles) +@benchmark backproject($sinogram_3D, $angles) # ╔═╡ 494190e4-e28a-42ad-b1c1-926a740ea2cb sinogram_3D_c = radon(sample_3D_c, angles_c, backend=CUDABackend()) @@ -102,7 +102,7 @@ sinogram_3D_c = radon(sample_3D_c, angles_c, backend=CUDABackend()) @benchmark CUDA.@sync radon($sample_3D_c, $angles_c, backend=CUDABackend()) # ╔═╡ cd166ef5-863f-48b9-bb89-a2a32daa166b -@benchmark CUDA.@sync iradon($sinogram_3D_c, $angles_c, backend=CUDABackend()) +@benchmark CUDA.@sync backproject($sinogram_3D_c, $angles_c, backend=CUDABackend()) # ╔═╡ Cell order: # ╠═e5a21d4e-928b-11ee-3909-53f97530eefa diff --git a/src/RadonKA.jl b/src/RadonKA.jl index f06a4b3..9f5bf57 100644 --- a/src/RadonKA.jl +++ b/src/RadonKA.jl @@ -10,7 +10,7 @@ using PrecompileTools include("geometry.jl") include("utils.jl") include("radon.jl") -include("iradon.jl") +include("backproject.jl") include("filtered_backprojection.jl") @@ -21,16 +21,16 @@ include("filtered_backprojection.jl") @compile_workload begin r = radon(Float32.(img), Float32.(angles)) - iradon(r, Float32.(angles)) + backproject(r, Float32.(angles)) RadonKA.filtered_backprojection(r, angles) r = radon(img, angles) - iradon(r, angles) + backproject(r, angles) RadonKA.filtered_backprojection(r, angles) r = radon(Float32.(img), Float32.(angles), μ=0.1f0) - iradon(r, Float32.(angles), μ=0.1f0) + backproject(r, Float32.(angles), μ=0.1f0) r = radon(img, angles, μ=0.1) - iradon(r, angles, μ=0.1) + backproject(r, angles, μ=0.1) end end diff --git a/src/iradon.jl b/src/backproject.jl similarity index 83% rename from src/iradon.jl rename to src/backproject.jl index fcf7417..edcbd97 100644 --- a/src/iradon.jl +++ b/src/backproject.jl @@ -1,18 +1,18 @@ -export iradon - +export backproject +@deprecate iradon backproject # handle 2D -function iradon(sinogram::AbstractArray{T, 2}, angles::AbstractArray{T2, 1}; +function backproject(sinogram::AbstractArray{T, 2}, angles::AbstractArray{T2, 1}; geometry=RadonParallelCircle(size(sinogram,1) + 1,-(size(sinogram,1))÷2:(size(sinogram,1))÷2), μ=nothing) where {T, T2} - view(iradon(reshape(sinogram, (size(sinogram)..., 1)), angles; geometry, μ), :, :, 1) + view(backproject(reshape(sinogram, (size(sinogram)..., 1)), angles; geometry, μ), :, :, 1) end """ - iradon(sinogram, θs; ) + backproject(sinogram, θs; ) -Conceptually the adjoint operation of [`radon`](@ref). Intuitively, the `iradon` smears rays back into the space. +Conceptually the adjoint operation of [`radon`](@ref). Intuitively, the `backproject` smears rays back into the space. See also [`radon`](@ref). # Example @@ -22,7 +22,7 @@ julia> arr = zeros((5,2)); arr[2,:] .= 1; arr[4, :] .= 1 1.0 1.0 -julia> iradon(arr, [0, π/2]) +julia> backproject(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.0 0.0 0.0 0.0 @@ -33,7 +33,7 @@ julia> iradon(arr, [0, π/2]) julia> arr = ones((2,1)); -julia> iradon(arr, [0], geometry=RadonFlexibleCircle(10, [-3, 3], [0,0])) +julia> backproject(arr, [0], geometry=RadonFlexibleCircle(10, [-3, 3], [0,0])) 10×10 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.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 @@ -47,13 +47,13 @@ julia> iradon(arr, [0], geometry=RadonFlexibleCircle(10, [-3, 3], [0,0])) 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ``` """ -function iradon(sinogram::AbstractArray{T, 3}, angles_T::AbstractVector; +function backproject(sinogram::AbstractArray{T, 3}, angles_T::AbstractVector; geometry=RadonParallelCircle(size(sinogram,1) + 1,-(size(sinogram,1))÷2:(size(sinogram,1))÷2), μ=nothing) where {T} - return _iradon(sinogram::AbstractArray{T, 3}, angles_T::AbstractVector, geometry, μ) + return _backproject(sinogram::AbstractArray{T, 3}, angles_T::AbstractVector, geometry, μ) end -function _iradon(sinogram::AbstractArray{T, 3}, angles_T::AbstractVector, geometry::Union{RadonParallelCircle, RadonFlexibleCircle}, μ) where T +function _backproject(sinogram::AbstractArray{T, 3}, angles_T::AbstractVector, geometry::Union{RadonParallelCircle, RadonFlexibleCircle}, μ) where T @assert size(sinogram, 2) == length(angles_T) "size of angles does not match sinogram size" @assert size(sinogram, 1) == size(geometry.in_height, 1) if geometry isa RadonFlexibleCircle @@ -96,14 +96,14 @@ function _iradon(sinogram::AbstractArray{T, 3}, angles_T::AbstractVector, geomet absorb_f = make_absorption_f(μ, T) # of the kernel goes - kernel! = iradon_kernel2!(backend) + kernel! = backproject_kernel2!(backend) kernel!(img, sinogram, weights, in_height, out_height, angles, mid, radius, absorb_f, ndrange=(size(sinogram, 1), N_angles, size(img, 3))) KernelAbstractions.synchronize(backend) return img end -@kernel function iradon_kernel2!(img::AbstractArray{T}, sinogram::AbstractArray{T}, weights, in_height, +@kernel function backproject_kernel2!(img::AbstractArray{T}, sinogram::AbstractArray{T}, weights, in_height, out_height, angles, mid, radius, absorb_f) where {T} i, iangle, i_z = @index(Global, NTuple) @@ -159,12 +159,12 @@ end # define adjoint rules -function ChainRulesCore.rrule(::typeof(_iradon), array::AbstractArray, angles, +function ChainRulesCore.rrule(::typeof(_backproject), array::AbstractArray, angles, geometry, μ) - res = _iradon(array, angles, geometry, μ) - function pb_iradon(ȳ) + res = _backproject(array, angles, geometry, μ) + function pb_backproject(ȳ) ad = _radon(unthunk(ȳ), angles, geometry, μ) return NoTangent(), ad, NoTangent(), NoTangent(), NoTangent() end - return res, pb_iradon + return res, pb_backproject end diff --git a/src/filtered_backprojection.jl b/src/filtered_backprojection.jl index 48e3b8d..9ecabb9 100644 --- a/src/filtered_backprojection.jl +++ b/src/filtered_backprojection.jl @@ -17,5 +17,5 @@ function filtered_backprojection(sinogram::AbstractArray{T}, θs::AbstractVector p = plan_fft(sinogram, (1,)) sinogram = real(inv(p) * (p * sinogram .* ifftshift(filter))) - return iradon(sinogram, θs; geometry, μ) + return backproject(sinogram, θs; geometry, μ) end diff --git a/src/geometry.jl b/src/geometry.jl index 113a374..8c6382f 100644 --- a/src/geometry.jl +++ b/src/geometry.jl @@ -3,13 +3,13 @@ export RadonGeometry, RadonParallelCircle, RadonFlexibleCircle """ abstract type RadonGeometry end -Abstract supertype for all geometries which are supported by `radon` and `iradon`. +Abstract supertype for all geometries which are supported by `radon` and `backproject`. List of geometries: * [`RadonParallelCircle`](@ref) * [`RadonFlexibleCircle`](@ref) -See [`radon`](@ref) and [`iradon`](@ref) how to apply. +See [`radon`](@ref) and [`backproject`](@ref) how to apply. """ abstract type RadonGeometry end @@ -25,7 +25,7 @@ So the resulting sinogram has the shape: `(9, length(angles), size(array, 3))`. - `weights` can weight each of the rays with different strength. Per default `weights = 0 .* in_height .+ 1` -See [`radon`](@ref) and [`iradon`](@ref) how to apply. +See [`radon`](@ref) and [`backproject`](@ref) how to apply. """ struct RadonParallelCircle{T, T2} <: RadonGeometry N::Int @@ -56,7 +56,7 @@ This is an extreme form of fan beam tomography. - `weights` can weight each of the rays with different strength. Per default `weights = 0 .* in_height .+ 1` -See [`radon`](@ref) and [`iradon`](@ref) how to apply. +See [`radon`](@ref) and [`backproject`](@ref) how to apply. """ struct RadonFlexibleCircle{T, T2, T3} <: RadonGeometry N::Int diff --git a/src/radon.jl b/src/radon.jl index a02d54a..a37c1f0 100644 --- a/src/radon.jl +++ b/src/radon.jl @@ -32,7 +32,7 @@ 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. -Both `radon` and [`iradon`](@ref) are differentiable with respect to `I`. +Both `radon` and [`backproject`](@ref) are differentiable with respect to `I`. # Keywords ## `μ=nothing` - Attenuated Radon Transform @@ -47,7 +47,7 @@ See `?RadonGeometries` for a full list of geometries. There is also the very fle -See also [`iradon`](@ref). +See also [`backproject`](@ref). # Example The reason the sinogram has the value `1.41421` for the diagonal ray `π/4` is, @@ -236,7 +236,7 @@ function ChainRulesCore.rrule(::typeof(_radon), array::AbstractArray, angles, geometry, μ) res = _radon(array, angles, geometry, μ) function pb_radon(ȳ) - ad = _iradon(unthunk(ȳ), angles, geometry, μ) + ad = _backproject(unthunk(ȳ), angles, geometry, μ) return NoTangent(), ad, NoTangent(), NoTangent(), NoTangent() end return res, pb_radon diff --git a/test/notebook.jl b/test/notebook.jl index 188c36a..f7810dd 100644 --- a/test/notebook.jl +++ b/test/notebook.jl @@ -27,19 +27,19 @@ end simshow(sinogram[:, :, 1]) # ╔═╡ 7fcb40e4-4530-4ce1-a267-84b46309ab2a -@show iradon(sinogram, [0f0, π*0.5]) +@show backproject(sinogram, [0f0, π*0.5]) # ╔═╡ c4cd76e9-1685-40b4-bde6-2d1b6c03fa09 -simshow(iradon(sinogram, [0f0, π*0.5]))[:, :, 1] +simshow(backproject(sinogram, [0f0, π*0.5]))[:, :, 1] # ╔═╡ 49629e15-c5f4-42ff-b196-24e13124f8f3 -@show I_r = iradon(sinogram, Float32[π/4 + π]) +@show I_r = backproject(sinogram, Float32[π/4 + π]) # ╔═╡ be118fec-3051-47a4-bf15-ac9d8d8c8e2e -sum(iradon(sinogram, Float32[2pi - pi/4])) +sum(backproject(sinogram, Float32[2pi - pi/4])) # ╔═╡ 80970ba6-a09f-46db-8b15-fa7226cdc400 -sum(iradon(sinogram, Float32[0])) +sum(backproject(sinogram, Float32[0])) # ╔═╡ c969f6d0-cfcc-4846-b970-9f6b734ac7cd simshow(I_r[:, :, 1]) @@ -64,16 +64,16 @@ begin end # ╔═╡ 8d9b6ea3-bbc0-48c1-a824-9749d6b8df0f -sum(iradon(sinogram_big, [0f0])[:, :, 1]) +sum(backproject(sinogram_big, [0f0])[:, :, 1]) # ╔═╡ 73a756d9-9984-4452-b90a-9cced7959813 -sum(iradon(sinogram_big, [pi/4f0])[:, :, 1]) +sum(backproject(sinogram_big, [pi/4f0])[:, :, 1]) # ╔═╡ b1607e05-d9ac-46cc-83b0-0f2b2c433059 -simshow(iradon(sinogram_big, [3.4f0])[:, :, 1]) +simshow(backproject(sinogram_big, [3.4f0])[:, :, 1]) # ╔═╡ 773e3961-e398-4d4b-a5e6-76b457561976 -simshow(iradon(sinogram_big, [0])[:, :, 1]) +simshow(backproject(sinogram_big, [0])[:, :, 1]) # ╔═╡ 5103d58e-34c8-4b9c-a0ab-c4f88f436186 md"# Radon" @@ -228,7 +228,7 @@ simshow(array3) simshow(sinogram2[:,:,1]) # ╔═╡ 2a3fe939-77b0-4388-9b9e-f9bee9e35e1b -I_2 = iradon(sinogram2, angles2); +I_2 = backproject(sinogram2, angles2); # ╔═╡ a33c5c99-00c3-437d-bfa7-2c0d12fb2339 simshow(I_2[:,:,1]) @@ -240,10 +240,10 @@ simshow(sinogram2) array_filtered = filtered_backprojection(sinogram2, angles2) # ╔═╡ ecb578ff-4dca-4260-9124-0bc17499d071 -array_iradon = iradon(sinogram2, angles2) +array_backproject = backproject(sinogram2, angles2) # ╔═╡ 2adbe3a7-c0bf-4e2d-bdec-4f62e29055af -sum(array_iradon) +sum(array_backproject) # ╔═╡ 65781501-0237-4caf-aa4e-4a83338e4658 simshow(array_filtered / sum(array_filtered)) diff --git a/test/runtests.jl b/test/runtests.jl index 12e04c6..99d7a76 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,13 +6,13 @@ using Zygote @testset "RadonKA.jl" begin - @testset "Simple iradon test" begin + @testset "Simple backproject test" begin sinogram = zeros(Float32, (9, 2, 1)) sinogram[5, :, 1] .= 1 - @test iradon(sinogram, [0.0f0, π * 0.5]) ≈ Float32[0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 1.0 1.0 1.0 1.0 2.0 1.0 1.0 1.0 1.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0;;;] + @test backproject(sinogram, [0.0f0, π * 0.5]) ≈ Float32[0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 1.0 1.0 1.0 1.0 2.0 1.0 1.0 1.0 1.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0;;;] - @test iradon(sinogram[:, 1, :], Float32[π / 4 + π]) == Float32[0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.96446574 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 1.4142138 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 1.4142135 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.4142132 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 1.4142137 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.4142133 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.3486991f-6 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0] + @test backproject(sinogram[:, 1, :], Float32[π / 4 + π]) == Float32[0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.96446574 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 1.4142138 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 1.4142135 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.4142132 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 1.4142137 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.4142133 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.3486991f-6 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0] end @testset "Simple API test" begin @@ -29,24 +29,24 @@ using Zygote sinogram = ones((3, 1)) - @test iradon(sinogram, [0], geometry = RadonFlexibleCircle(20, [0, 0, 0], [0, 0, 0])) ≈ [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0] + @test backproject(sinogram, [0], geometry = RadonFlexibleCircle(20, [0, 0, 0], [0, 0, 0])) ≈ [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0] - @test iradon(sinogram, [0], geometry = RadonFlexibleCircle(20, [-5, 0, 5], [0, 0, 0])) ≈ [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0848743343786964 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.8356890637711426 1.4134962068364105 0.8356890637711426 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0424371671893484 1.0 1.0424371671893484 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0424371671893484 1.0 1.0424371671893484 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.42219285693473696 0.620244310254611 1.0 0.620244310254611 0.42219285693473696 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.042437167189348 0.0 1.0 0.0 1.042437167189348 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0424371671893484 0.0 1.0 0.0 1.0424371671893484 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.008696650098325685 1.0337405170910223 0.0 1.0 0.0 1.0337405170910223 0.008696650098325685 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0424371671893484 0.0 0.0 1.0 0.0 0.0 1.0424371671893484 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.042437167189348 0.0 0.0 1.0 0.0 0.0 1.042437167189348 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0424371671893484 0.0 0.0 1.0 0.0 0.0 1.0424371671893484 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.6376376104512703 0.40479955673807827 0.0 0.0 1.0 0.0 0.0 0.40479955673807827 0.6376376104512703 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 1.0424371671893482 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0424371671893484 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 1.0424371671893482 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0424371671893484 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.22414140361486262 0.8182957635744856 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.8182957635744816 0.22414140361486629 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0424371671893482 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.042437167189348 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.5038252833980154 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.5038252833980154 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0] + @test backproject(sinogram, [0], geometry = RadonFlexibleCircle(20, [-5, 0, 5], [0, 0, 0])) ≈ [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0848743343786964 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.8356890637711426 1.4134962068364105 0.8356890637711426 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0424371671893484 1.0 1.0424371671893484 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0424371671893484 1.0 1.0424371671893484 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.42219285693473696 0.620244310254611 1.0 0.620244310254611 0.42219285693473696 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.042437167189348 0.0 1.0 0.0 1.042437167189348 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0424371671893484 0.0 1.0 0.0 1.0424371671893484 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.008696650098325685 1.0337405170910223 0.0 1.0 0.0 1.0337405170910223 0.008696650098325685 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0424371671893484 0.0 0.0 1.0 0.0 0.0 1.0424371671893484 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.042437167189348 0.0 0.0 1.0 0.0 0.0 1.042437167189348 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0424371671893484 0.0 0.0 1.0 0.0 0.0 1.0424371671893484 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.6376376104512703 0.40479955673807827 0.0 0.0 1.0 0.0 0.0 0.40479955673807827 0.6376376104512703 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 1.0424371671893482 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0424371671893484 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 1.0424371671893482 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0424371671893484 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.22414140361486262 0.8182957635744856 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.8182957635744816 0.22414140361486629 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0424371671893482 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.042437167189348 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.5038252833980154 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.5038252833980154 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0] - @test ≈(sum(iradon(sinogram, [0], geometry=RadonFlexibleCircle(20, 0.01 .+[-1,0,1], 0.01 .+[-1,0,1]))) + 1, sum(iradon(sinogram, [0], geometry=RadonFlexibleCircle(20, 0.00 .+[-1,0,1], 0.00 .+[-1,0,1]))), rtol=0.001) + @test ≈(sum(backproject(sinogram, [0], geometry=RadonFlexibleCircle(20, 0.01 .+[-1,0,1], 0.01 .+[-1,0,1]))) + 1, sum(backproject(sinogram, [0], geometry=RadonFlexibleCircle(20, 0.00 .+[-1,0,1], 0.00 .+[-1,0,1]))), rtol=0.001) end - @testset "Exponential iradon" begin + @testset "Exponential backproject" begin sinogram = zeros(Float64, (9, 1)) sinogram[5, :] .= 1 angles = [0] - arr = iradon(sinogram, angles, μ=0.1) + arr = backproject(sinogram, angles, μ=0.1) @test exp.(-(9:-1:1) .* 0.1) ≈ arr[begin+1:end, 6][:] angles = [π] - arr = iradon(sinogram, angles, μ=0.1) + arr = backproject(sinogram, angles, μ=0.1) @test exp.(-(1:1:8).* 0.1) ≈ arr[begin+1:end-1, 6][:] end @@ -112,11 +112,11 @@ using Zygote for geometry in [geometry, geometry1] test_rrule(RadonKA._radon, x, [0, π/4, π/2, 2π, 0.1, 0.00001] ⊢ ChainRulesTestUtils.NoTangent(), geometry ⊢ ChainRulesTestUtils.NoTangent(), nothing ⊢ ChainRulesTestUtils.NoTangent()) y = radon(x, [0, π/4, π/2, 2π, 0.1, 0.00001]) - test_rrule(RadonKA._iradon, y, [0, π/4, π/2, 2π, 0.1, 0.00001] ⊢ ChainRulesTestUtils.NoTangent(), geometry ⊢ ChainRulesTestUtils.NoTangent(), nothing ⊢ ChainRulesTestUtils.NoTangent()) + test_rrule(RadonKA._backproject, y, [0, π/4, π/2, 2π, 0.1, 0.00001] ⊢ ChainRulesTestUtils.NoTangent(), geometry ⊢ ChainRulesTestUtils.NoTangent(), nothing ⊢ ChainRulesTestUtils.NoTangent()) test_rrule(RadonKA._radon, x, [0, π/4, π/2, 2π, 0.1, 0.00001] ⊢ ChainRulesTestUtils.NoTangent(), geometry ⊢ ChainRulesTestUtils.NoTangent(), 0.1⊢ ChainRulesTestUtils.NoTangent()) y = radon(x, [0, π/4, π/2, 2π, 0.1, 0.00001]) - test_rrule(RadonKA._iradon, y, [0, π/4, π/2, 2π, 0.1, 0.00001] ⊢ ChainRulesTestUtils.NoTangent(), geometry ⊢ ChainRulesTestUtils.NoTangent(), 0.1 ⊢ ChainRulesTestUtils.NoTangent()) + test_rrule(RadonKA._backproject, y, [0, π/4, π/2, 2π, 0.1, 0.00001] ⊢ ChainRulesTestUtils.NoTangent(), geometry ⊢ ChainRulesTestUtils.NoTangent(), 0.1 ⊢ ChainRulesTestUtils.NoTangent()) end end end