diff --git a/Project.toml b/Project.toml index d87a37d..8382d4d 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.2.1" +version = "0.3" [deps] Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458" diff --git a/README.md b/README.md index acfe949..5affa0b 100644 --- a/README.md +++ b/README.md @@ -9,13 +9,15 @@ 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` -* [x] arbitrary (2D) geometries can be specified with `ray_startpoints` and `ray_endpoints` -* [x] attenuated `radon` and `iradon` (see this [paper](https://iopscience.iop.org/article/10.1088/0266-5611/17/1/309/meta) as reference) +* [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] arbitrary 2D geometries where starting and endpoint of each ray can be specified (cone 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] high performance however not ultra high performance +* [x] simple API # Installation Requires Julia at least 1.9 @@ -23,7 +25,6 @@ Requires Julia at least 1.9 julia> ]add RadonKA ``` - # Simple use ```julia using RadonKA, ImageShow, ImageIO, TestImages diff --git a/docs/Manifest.toml b/docs/Manifest.toml index ab1987e..e29e7b1 100644 --- a/docs/Manifest.toml +++ b/docs/Manifest.toml @@ -1,6 +1,6 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.9.4" +julia_version = "1.10.0" manifest_format = "2.0" project_hash = "eddb71d2ac352012abcf8111e85f5032b7826df9" @@ -27,9 +27,9 @@ version = "0.4.4" [[deps.Adapt]] deps = ["LinearAlgebra", "Requires"] -git-tree-sha1 = "02f731463748db57cc2ebfbd9fbc9ce8280d3433" +git-tree-sha1 = "0fb305e0253fd4e833d486914367a2ee2c2e78d0" uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -version = "3.7.1" +version = "4.0.1" weakdeps = ["StaticArrays"] [deps.Adapt.extensions] @@ -52,25 +52,25 @@ version = "0.1.0" uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" [[deps.CEnum]] -git-tree-sha1 = "eb4cb44a499229b3b8426dcfb5dd85333951ff90" +git-tree-sha1 = "389ad5c84de1ae7cf0e28e381131c98ea87d54fc" uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82" -version = "0.4.2" +version = "0.5.0" [[deps.ChainRulesCore]] deps = ["Compat", "LinearAlgebra"] -git-tree-sha1 = "e0af648f0692ec1691b5d094b8724ba1346281cf" +git-tree-sha1 = "1287e3872d646eed95198457873249bd9f0caed2" uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -version = "1.18.0" +version = "1.20.1" weakdeps = ["SparseArrays"] [deps.ChainRulesCore.extensions] ChainRulesCoreSparseArraysExt = "SparseArrays" [[deps.Compat]] -deps = ["UUIDs"] -git-tree-sha1 = "886826d76ea9e72b35fcd000e535588f7b60f21d" +deps = ["TOML", "UUIDs"] +git-tree-sha1 = "75bd5b6fc5089df449b5d35fa501c846c9b6549b" uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "4.10.1" +version = "4.12.0" weakdeps = ["Dates", "LinearAlgebra"] [deps.Compat.extensions] @@ -79,7 +79,7 @@ weakdeps = ["Dates", "LinearAlgebra"] [[deps.CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" -version = "1.0.5+0" +version = "1.0.5+1" [[deps.Dates]] deps = ["Printf"] @@ -93,9 +93,9 @@ version = "0.9.3" [[deps.Documenter]] deps = ["ANSIColoredPrinters", "AbstractTrees", "Base64", "Dates", "DocStringExtensions", "Downloads", "Git", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "MarkdownAST", "Pkg", "PrecompileTools", "REPL", "RegistryInstances", "SHA", "Test", "Unicode"] -git-tree-sha1 = "43aa88b72dffff46b1b19f66483ea3e2f907c4fa" +git-tree-sha1 = "2613dbec8f4748273bbe30ba71fd5cb369966bac" uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -version = "1.2.0" +version = "1.2.1" [[deps.Downloads]] deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] @@ -110,9 +110,9 @@ version = "2.5.0+0" [[deps.FFTW]] deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"] -git-tree-sha1 = "b4fbdd20c889804969571cc589900803edda16b7" +git-tree-sha1 = "4820348781ae578893311153d69049a93d05f39d" uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" -version = "1.7.1" +version = "1.8.0" [[deps.FFTW_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -131,27 +131,27 @@ version = "1.3.0" [[deps.Git_jll]] deps = ["Artifacts", "Expat_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "Libiconv_jll", "OpenSSL_jll", "PCRE2_jll", "Zlib_jll"] -git-tree-sha1 = "bb8f7cc77ec1152414b2af6db533d9471cfbb2d1" +git-tree-sha1 = "b30c473c97fcc1e1e44fab8f3e88fd1b89c9e9d1" uuid = "f8c6e375-362e-5223-8a59-34ff63f689eb" -version = "2.42.0+0" +version = "2.43.0+0" [[deps.IOCapture]] deps = ["Logging", "Random"] -git-tree-sha1 = "d75853a0bdbfb1ac815478bacd89cd27b550ace6" +git-tree-sha1 = "8b72179abc660bfab5e28472e019392b97d0985c" uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89" -version = "0.2.3" +version = "0.2.4" [[deps.IndexFunArrays]] deps = ["ChainRulesCore", "LinearAlgebra"] -git-tree-sha1 = "e66a2aeb6d5814015004080e5203dfff44d2856f" +git-tree-sha1 = "6f78703c7a4ba06299cddd8694799c91de0157ac" uuid = "613c443e-d742-454e-bfc6-1d7f8dd76566" -version = "0.2.6" +version = "0.2.7" [[deps.IntelOpenMP_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "ad37c091f7d7daf900963171600d7c1c5c3ede32" +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "5fdf2fe6724d8caabf43b557b84ce53f3b7e2f6b" uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0" -version = "2023.2.0+0" +version = "2024.0.2+0" [[deps.InteractiveUtils]] deps = ["Markdown"] @@ -171,9 +171,9 @@ version = "0.21.4" [[deps.KernelAbstractions]] deps = ["Adapt", "Atomix", "InteractiveUtils", "LinearAlgebra", "MacroTools", "PrecompileTools", "Requires", "SparseArrays", "StaticArrays", "UUIDs", "UnsafeAtomics", "UnsafeAtomicsLLVM"] -git-tree-sha1 = "b0737cbbe1c8da6f1139d1c23e35e7cea129c0af" +git-tree-sha1 = "4e0cb2f5aad44dcfdc91088e85dee4ecb22c791c" uuid = "63c18a36-062a-441e-b654-da1e3ab1ce7c" -version = "0.9.13" +version = "0.9.16" [deps.KernelAbstractions.extensions] EnzymeExt = "EnzymeCore" @@ -183,9 +183,9 @@ version = "0.9.13" [[deps.LLVM]] deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Preferences", "Printf", "Requires", "Unicode"] -git-tree-sha1 = "c879e47398a7ab671c782e02b51a4456794a7fa3" +git-tree-sha1 = "cb4619f7353fc62a1a22ffa3d7ed9791cfb47ad8" uuid = "929cbde3-209d-540e-8aea-75f648917ca0" -version = "6.4.0" +version = "6.4.2" [deps.LLVM.extensions] BFloat16sExt = "BFloat16s" @@ -219,9 +219,14 @@ uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" version = "8.4.0+0" [[deps.LibGit2]] -deps = ["Base64", "NetworkOptions", "Printf", "SHA"] +deps = ["Base64", "LibGit2_jll", "NetworkOptions", "Printf", "SHA"] uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" +[[deps.LibGit2_jll]] +deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll"] +uuid = "e37daf67-58a4-590a-8e99-b0245dd2ffc5" +version = "1.6.4+0" + [[deps.LibSSH2_jll]] deps = ["Artifacts", "Libdl", "MbedTLS_jll"] uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" @@ -244,16 +249,16 @@ uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" [[deps.MKL_jll]] -deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"] -git-tree-sha1 = "eb006abbd7041c28e0d16260e50a24f8f9104913" +deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl"] +git-tree-sha1 = "72dc3cf284559eb8f53aa593fe62cb33f83ed0c0" uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7" -version = "2023.2.0+0" +version = "2024.0.0+0" [[deps.MacroTools]] deps = ["Markdown", "Random"] -git-tree-sha1 = "9ee1618cbf5240e6d4e0371d6f24065083f60c48" +git-tree-sha1 = "2fa9ee3e63fd3a4f7a9a4f4744a52f4856de82df" uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" -version = "0.5.11" +version = "0.5.13" [[deps.Markdown]] deps = ["Base64"] @@ -268,14 +273,14 @@ version = "0.1.2" [[deps.MbedTLS_jll]] deps = ["Artifacts", "Libdl"] uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" -version = "2.28.2+0" +version = "2.28.2+1" [[deps.Mmap]] uuid = "a63ad114-7e13-5084-954f-fe012c677804" [[deps.MozillaCACerts_jll]] uuid = "14a3606d-f60d-562e-9121-12d972cd8159" -version = "2022.10.11" +version = "2023.1.10" [[deps.NetworkOptions]] uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" @@ -284,29 +289,29 @@ version = "1.2.0" [[deps.OpenBLAS_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" -version = "0.3.21+4" +version = "0.3.23+2" [[deps.OpenSSL_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "cc6e1927ac521b659af340e0ca45828a3ffc748f" +git-tree-sha1 = "60e3045590bd104a16fefb12836c00c0ef8c7f8c" uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" -version = "3.0.12+0" +version = "3.0.13+0" [[deps.PCRE2_jll]] deps = ["Artifacts", "Libdl"] uuid = "efcefdf7-47ab-520b-bdef-62a2eaa19f15" -version = "10.42.0+0" +version = "10.42.0+1" [[deps.Parsers]] deps = ["Dates", "PrecompileTools", "UUIDs"] -git-tree-sha1 = "a935806434c9d4c506ba941871b327b96d41f2bf" +git-tree-sha1 = "8489905bcdbcfac64d1daa51ca07c0d8f0283821" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "2.8.0" +version = "2.8.1" [[deps.Pkg]] deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -version = "1.9.2" +version = "1.10.0" [[deps.PrecompileTools]] deps = ["Preferences"] @@ -329,13 +334,13 @@ deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" [[deps.RadonKA]] -deps = ["Adapt", "Atomix", "FFTW", "IndexFunArrays", "KernelAbstractions"] +deps = ["Atomix", "ChainRulesCore", "FFTW", "IndexFunArrays", "KernelAbstractions", "PrecompileTools"] path = ".." uuid = "86de8297-835b-47df-b249-c04e8db91db5" -version = "0.1.0" +version = "0.3.0" [[deps.Random]] -deps = ["SHA", "Serialization"] +deps = ["SHA"] uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [[deps.Reexport]] @@ -368,17 +373,20 @@ uuid = "6462fe0b-24de-5631-8697-dd941f90decc" [[deps.SparseArrays]] deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"] uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +version = "1.10.0" [[deps.StaticArrays]] deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"] -git-tree-sha1 = "5ef59aea6f18c25168842bded46b16662141ab87" +git-tree-sha1 = "7b0e9c14c624e435076d19aea1e5cbdec2b9ca37" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.7.0" +version = "1.9.2" [deps.StaticArrays.extensions] + StaticArraysChainRulesCoreExt = "ChainRulesCore" StaticArraysStatisticsExt = "Statistics" [deps.StaticArrays.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" [[deps.StaticArraysCore]] @@ -387,9 +395,9 @@ uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" version = "1.4.2" [[deps.SuiteSparse_jll]] -deps = ["Artifacts", "Libdl", "Pkg", "libblastrampoline_jll"] +deps = ["Artifacts", "Libdl", "libblastrampoline_jll"] uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c" -version = "5.10.1+6" +version = "7.2.1+1" [[deps.TOML]] deps = ["Dates"] @@ -426,12 +434,12 @@ version = "0.1.3" [[deps.Zlib_jll]] deps = ["Libdl"] uuid = "83775a58-1f1d-513f-b197-d71354ab007a" -version = "1.2.13+0" +version = "1.2.13+1" [[deps.libblastrampoline_jll]] deps = ["Artifacts", "Libdl"] uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" -version = "5.8.0+0" +version = "5.8.0+1" [[deps.nghttp2_jll]] deps = ["Artifacts", "Libdl"] @@ -441,4 +449,4 @@ version = "1.52.0+1" [[deps.p7zip_jll]] deps = ["Artifacts", "Libdl"] uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" -version = "17.4.0+0" +version = "17.4.0+2" diff --git a/docs/src/benchmark.md b/docs/src/benchmark.md index 0466e3d..441728d 100644 --- a/docs/src/benchmark.md +++ b/docs/src/benchmark.md @@ -1,51 +1,52 @@ # Benchmark and Comparison with Matlab and Astra -Tested on a AMD Ryzen 9 5900X 12-Core Processor with 24 Threads and a NVIDIA GeForce RTX 3060 with Julia 1.9.4 on Ubuntu 22.04. +Tested on a AMD Ryzen 9 5900X 12-Core Processor with 24 Threads and a NVIDIA GeForce RTX 3060 with Julia 1.10.0 on Ubuntu 22.04. # Results | |RadonKA.jl CPU | RadonKA.jl GPU | Matlab CPU | Astra CPU | Astra GPU | |-------------------|---------------|-------------------|------------|-----------|-----------| -|2D sample - Radon |1.2s |0.091s |0.39s |7.0s |0.025s | -|2D sample - IRadon |1.8s |0.52s |0.37s |6.4s |N/A | -|3D sample - Radon |8.4s |0.47s |9.01s |N/A |1.12s | -|3D sample - IRadon |10.5s |0.59s |3.24s |N/A |N/A | +|2D sample - Radon |1.1s |0.07s |0.39s |7.0s |0.025s | +|2D sample - IRadon |1.4s |0.50s |0.37s |6.4s |N/A | +|3D sample - Radon |7.5s |0.29s |9.01s |N/A |1.12s | +|3D sample - IRadon |7.9s |0.53s |3.24s |N/A |N/A | ## RadonKA.jl -See this [benchmark example](https://github.com/roflmaostc/RadonKA.jl/blob/main/examples/benchmark_radon_iradon.jl) ```julia -using IndexFunArrays, ImageShow, Plots, ImageIO, PlutoUI, PlutoTest, TestImages -using RadonKA -using CUDA, CUDA.CUDAKernels - -sample_2D = Float32.(testimage("resolution_test_1920")); -sample_2D_c = CuArray(sample_2D); -simshow(sample_2D) - -angles = range(0, 2π, 500); -angles_c = CuArray(angles); - -# run those cells multiple times -@time sinogram = radon(sample_2D, angles); -simshow(sinogram) -@time sample_iradon = iradon(sinogram, angles); - -CUDA.@time CUDA.@sync sinogram_c = radon(sample_2D_c, angles_c, backend=CUDABackend()); -CUDA.@time CUDA.@sync sample_iradon_c = iradon(sinogram_c, angles_c, backend=CUDABackend()); - - -sample_3D = randn(Float32, (512, 512, 100)); -sample_3D_c = CuArray(sample_3D) - -@time sinogram_3D = radon(sample_3D, angles); -@benchmark radon($sample_3D, $angles) -@benchmark iradon($sinogram_3D, $angles) - -sinogram_3D_c = radon(sample_3D_c, angles_c, backend=CUDABackend()) -@benchmark CUDA.@sync radon($sample_3D_c, $angles_c, backend=CUDABackend()) -@benchmark CUDA.@sync iradon($sinogram_3D_c, $angles_c, backend=CUDABackend()) - + using IndexFunArrays, ImageShow, Plots, ImageIO, PlutoUI, PlutoTest, TestImages + using RadonKA + using CUDA, CUDA.CUDAKernels + using BenchmarkTools + + sample_2D = Float32.(testimage("resolution_test_1920")); + sample_2D_c = CuArray(sample_2D); + simshow(sample_2D) + + angles = range(0, 2π, 500); + angles_c = CuArray(angles); + + # run those cells multiple times + sinogram = radon(sample_2D, angles); + @btime sinogram = radon($sample_2D, $angles); + simshow(sinogram) + @btime sample_iradon = iradon($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); + + + sample_3D = randn(Float32, (512, 512, 100)); + sample_3D_c = CuArray(sample_3D) + + sinogram_3D = radon(sample_3D, angles); + @btime radon($sample_3D, $angles) + @btime iradon($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) ``` ![](../assets/radonka_sinogram.png) ![](../assets/radonka_iradon.png) diff --git a/docs/src/index.md b/docs/src/index.md index 3ff6866..527ae9d 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,5 +1,9 @@ -# RadonKA -A simple but still decently fast Radon and IRadon transform based on [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl). +# RadonKA.jl +A simple yet sufficiently fast Radon and inverse Radon (iradon) transform implementation using [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl). + +```@raw html + +``` [![Build Status](https://github.com/roflmaostc/RadonKA.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/roflmaostc/RadonKA.jl/actions/workflows/CI.yml?query=branch%3Amain) [![Coverage](https://codecov.io/gh/roflmaostc/RadonKA.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/roflmaostc/RadonKA.jl) [![Documentation for stable version](https://img.shields.io/badge/docs-stable-blue.svg)](https://roflmaostc.github.io/RadonKA.jl/stable) [![Documentation for development version](https://img.shields.io/badge/docs-main-blue.svg)](https://roflmaostc.github.io/RadonKA.jl/dev) @@ -7,12 +11,15 @@ A simple but still decently fast Radon and IRadon transform based on [KernelAbst # Quick Overview * [x] For 2D and 3D arrays -* [x] parallel `radon` and `iradon` -* [x] parallel exponential `radon` and `iradon` +* [x] parallel `radon` and `iradon` (`?RadonParallelCircle`) +* [x] attenuated `radon` and `iradon` (see the parameter `μ`) +* [x] arbitrary 2D geometries where starting and endpoint of each ray can be specified (cone 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] high performance however not ultra high performance +* [x] simple API # Installation Requires Julia 1.9 @@ -20,5 +27,36 @@ Requires Julia 1.9 julia> ]add RadonKA ``` +# Simple use +```julia +using RadonKA, ImageShow, ImageIO, TestImages + +img = Float32.(testimage("resolution_test_512")) + +angles = range(0f0, 2f0π, 500)[begin:end-1] + +# 0.196049 seconds (145 allocations: 1009.938 KiB) +@time sinogram = radon(img, angles); + +# 0.268649 seconds (147 allocations: 1.015 MiB) +@time backproject = RadonKA.iradon(sinogram, angles); + +simshow(sinogram) +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. +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). + # Development File an [issue](https://github.com/roflmaostc/RadonKA.jl/issues) on [GitHub](https://github.com/roflmaostc/RadonKA.jl) if you encounter any problems. + + diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md index 751068d..c7c892d 100644 --- a/docs/src/tutorial.md +++ b/docs/src/tutorial.md @@ -52,7 +52,7 @@ RadonKA.jl supports CUDA and usually provides a 10-20x speedup on typical RTX 3x Pay attention that the element type of the array/img should be `Float32` for good speedup. In my case we used a AMD Ryzen 5 5600X 6-Core Processor and a CUDA RTX 3060 Super. ```julia -using CUDA, CUDA.CUDAKernels +using CUDA img_c = CuArray(img); angles_c = CuArray(angles); diff --git a/examples/CT_reconstruction.jl b/examples/CT_reconstruction.jl deleted file mode 100644 index c771f20..0000000 --- a/examples/CT_reconstruction.jl +++ /dev/null @@ -1,176 +0,0 @@ -### A Pluto.jl notebook ### -# v0.19.36 - -using Markdown -using InteractiveUtils - -# ╔═╡ 5cc65690-91d8-11ee-07f9-557c9a76e478 -begin - using Pkg - Pkg.activate(".") - using Revise -end - -# ╔═╡ c4b5d8f3-7ccd-4888-8b83-4642c1e81916 -using TestImages, RadonKA, ImageShow, ImageIO, Noise, PlutoUI - -# ╔═╡ 4458d362-1aae-4a24-9bf7-2bed6286579a -using RegularizedLeastSquares, IterativeSolvers, LinearOperators - -# ╔═╡ b881ba80-7b02-4f2a-9fad-b61c4ee218c1 -using CUDA, CUDA.CUDAKernels - -# ╔═╡ f95dfc57-060f-4173-aedc-25e36b67fcb4 -md"# Load Testimage" - -# ╔═╡ 162a9e47-985c-406d-b8db-914a8ee17047 -N = 256 - -# ╔═╡ 2e58525e-3c9b-4bd4-a903-149a27939af0 -sample = Float32.(TestImages.shepp_logan(N)); - -# ╔═╡ 8503c029-c37f-44eb-8e66-3acc01dd678a -simshow(sample) - -# ╔═╡ b31c0d18-4338-450f-bf58-c825ec813519 -md"# Make Sinogram" - -# ╔═╡ e90ac7f4-504b-47ea-b087-9a25afe8ec22 -K = 300 - -# ╔═╡ 29ca822c-fe80-463c-8d6e-efcb178d8d67 -angles = Float32.(range(0, π, K)) - -# ╔═╡ 0378cea6-8a61-4e16-a99a-e2160c83518b -@time sinogram = poisson(radon(sample, angles), 400); - -# ╔═╡ c93ecf01-1a06-4c3e-9c4d-ae1516f38749 -size(sinogram) - -# ╔═╡ f2573a91-c0a6-421d-8cb9-8745b5c5d276 -simshow(sinogram) - -# ╔═╡ c40638f2-9bcc-4954-b80a-126c6afbd344 -md"# Backproject with iradon" - -# ╔═╡ fa6f64e8-557f-4c4b-81ba-c1bae5b98f36 -@time sample_iradon = iradon(sinogram, angles); - -# ╔═╡ c9ffc233-0e0b-452d-8f53-09b8654c8bca -simshow(sample_iradon) - -# ╔═╡ 7b7579a8-8799-4ba9-8931-ee4642a36124 -md"# Filtered Backprojection" - -# ╔═╡ ca4fdba3-2504-4f26-860e-e95c6aae8f12 -sample_filtered = filtered_backprojection(sinogram, angles); - -# ╔═╡ eff77ed3-a9ea-4459-b64e-c4459b31619b -simshow(sample_filtered) - -# ╔═╡ a303b94f-8807-4fe5-974a-7914db10f8eb -md"# Iterative Reconstruction -For that we define a linear `radon_map` which defines the efficient forward operation and the adjoint. We use [LinearOperators](https://github.com/JuliaSmoothOptimizers/LinearOperators.jl) for that. We adapted [this example](https://jso.dev/tutorials/introduction-to-linear-operators/#using_functions). - -For the Regularized Least Squares we use [RegularizedLeastSquares.jl](https://github.com/JuliaImageRecon/RegularizedLeastSquares.jl). - -" - -# ╔═╡ 104b78de-022b-4419-867a-1d99a2663f8c -function radon_f!(res, v, α, β) - if β == 0 - res .= α .* vec(radon(reshape(v, (N, N)), angles)) - else - res .= α .*vec(radon(reshape(v, (N, N)), angles)) .+ β .* res - end -end - -# ╔═╡ 3873e3c9-b9cf-48d7-9ef5-b1372778c138 -function iradon_f!(res, w, α, β) - if β == 0 - res .= α .* vec(iradon(reshape(w, (N-1, K)), angles)) - else - res .= α .* vec(iradon(reshape(w, (N-1, K)), angles)) .+ β .* res - end -end - -# ╔═╡ d3b27172-3373-4176-a8bd-475acd825e8d -radon_map = LinearOperator(Float32, (N - 1) * K, N*N, false, false, - radon_f!, - nothing, # will be inferred - iradon_f!) - -# ╔═╡ 90ef6756-60ac-4b5d-8620-73566a465d22 -# seem to work -simshow(reshape(radon_map * vec(sample), (N - 1, K))) - -# ╔═╡ e2207cd2-9471-460f-9eb4-fcfed59cbe7e -adjoint(radon_map) * vec(sinogram) ≈ sample_iradon[:] - -# ╔═╡ ad5b1b23-2187-411c-89c2-ea4a9c8bd2b0 -guess = vec(copy(sample_filtered)); - -# ╔═╡ 48bbb97c-d12d-4282-92f9-02525d21f3ee -reg = TVRegularization(0.001f0; shape=(N,N), startVector=guess) - -# ╔═╡ ea6d5aa7-09de-4c1e-9a46-00e7d5267433 -solver = createLinearSolver(ADMM, radon_map; reg=reg, ρ=0.1, iterations=20) - -# ╔═╡ 8917be57-4a1b-4f1a-bcf2-c29b32e0c726 -@time Ireco = reshape(solve(solver, vec(sinogram)), (N, N)) - -# ╔═╡ 924f005c-3b18-4ac7-9723-fda33068a90d -[simshow(sample) simshow(Ireco) simshow(sample_filtered)] - -# ╔═╡ 38eca17b-deca-47e2-a7fb-6b5fb93e1b3e -cg!(guess, radon_map, vec(sinogram)) - -# ╔═╡ 01b6465b-eabf-48b2-95f9-a689784a2516 -md"# Accelerate with CUDA" - -# ╔═╡ 846099aa-c8fe-4559-ac53-204146f7854d -sinogram_c = CuArray(sinogram); - -# ╔═╡ 5bdf05d9-fbc7-417b-b26e-2d645151fc60 -guess_c = CuArray(guess); - -# ╔═╡ 0166e1a8-03b5-480f-951f-4c1b73ebb954 -Ireco_c = reshape(solve(solver, vec(sinogram_c), startVector=guess_c), (N, N)) - -# ╔═╡ Cell order: -# ╠═5cc65690-91d8-11ee-07f9-557c9a76e478 -# ╠═c4b5d8f3-7ccd-4888-8b83-4642c1e81916 -# ╟─f95dfc57-060f-4173-aedc-25e36b67fcb4 -# ╠═162a9e47-985c-406d-b8db-914a8ee17047 -# ╠═2e58525e-3c9b-4bd4-a903-149a27939af0 -# ╠═8503c029-c37f-44eb-8e66-3acc01dd678a -# ╟─b31c0d18-4338-450f-bf58-c825ec813519 -# ╠═e90ac7f4-504b-47ea-b087-9a25afe8ec22 -# ╠═29ca822c-fe80-463c-8d6e-efcb178d8d67 -# ╠═0378cea6-8a61-4e16-a99a-e2160c83518b -# ╠═c93ecf01-1a06-4c3e-9c4d-ae1516f38749 -# ╠═f2573a91-c0a6-421d-8cb9-8745b5c5d276 -# ╟─c40638f2-9bcc-4954-b80a-126c6afbd344 -# ╠═fa6f64e8-557f-4c4b-81ba-c1bae5b98f36 -# ╠═c9ffc233-0e0b-452d-8f53-09b8654c8bca -# ╟─7b7579a8-8799-4ba9-8931-ee4642a36124 -# ╠═ca4fdba3-2504-4f26-860e-e95c6aae8f12 -# ╠═eff77ed3-a9ea-4459-b64e-c4459b31619b -# ╟─a303b94f-8807-4fe5-974a-7914db10f8eb -# ╠═4458d362-1aae-4a24-9bf7-2bed6286579a -# ╠═104b78de-022b-4419-867a-1d99a2663f8c -# ╠═3873e3c9-b9cf-48d7-9ef5-b1372778c138 -# ╠═d3b27172-3373-4176-a8bd-475acd825e8d -# ╠═90ef6756-60ac-4b5d-8620-73566a465d22 -# ╠═e2207cd2-9471-460f-9eb4-fcfed59cbe7e -# ╠═ad5b1b23-2187-411c-89c2-ea4a9c8bd2b0 -# ╠═ea6d5aa7-09de-4c1e-9a46-00e7d5267433 -# ╠═48bbb97c-d12d-4282-92f9-02525d21f3ee -# ╠═8917be57-4a1b-4f1a-bcf2-c29b32e0c726 -# ╠═924f005c-3b18-4ac7-9723-fda33068a90d -# ╠═38eca17b-deca-47e2-a7fb-6b5fb93e1b3e -# ╟─01b6465b-eabf-48b2-95f9-a689784a2516 -# ╠═b881ba80-7b02-4f2a-9fad-b61c4ee218c1 -# ╠═846099aa-c8fe-4559-ac53-204146f7854d -# ╠═5bdf05d9-fbc7-417b-b26e-2d645151fc60 -# ╠═0166e1a8-03b5-480f-951f-4c1b73ebb954 diff --git a/examples/CT_with_optimizer.jl b/examples/CT_with_optimizer.jl index 3f280dd..bf8e77d 100644 --- a/examples/CT_with_optimizer.jl +++ b/examples/CT_with_optimizer.jl @@ -1,5 +1,5 @@ ### A Pluto.jl notebook ### -# v0.19.36 +# v0.19.38 using Markdown using InteractiveUtils diff --git a/examples/Manifest.toml b/examples/Manifest.toml index 53b38e8..4bcf654 100644 --- a/examples/Manifest.toml +++ b/examples/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.10.0" manifest_format = "2.0" -project_hash = "c7b2d79018ada0dcd95eadb1b7f6353fbad9fd04" +project_hash = "ba21296ff49d3ce1d48865df6c57170e38033bac" [[deps.ADTypes]] git-tree-sha1 = "41c37aa88889c171f1300ceac1313c06e891d245" @@ -39,9 +39,9 @@ version = "0.4.4" [[deps.Adapt]] deps = ["LinearAlgebra", "Requires"] -git-tree-sha1 = "cde29ddf7e5726c9fb511f340244ea3481267608" +git-tree-sha1 = "0fb305e0253fd4e833d486914367a2ee2c2e78d0" uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -version = "3.7.2" +version = "4.0.1" weakdeps = ["StaticArrays"] [deps.Adapt.extensions] @@ -107,9 +107,9 @@ version = "0.4.2" [[deps.BangBang]] deps = ["Compat", "ConstructionBase", "InitialValues", "LinearAlgebra", "Requires", "Setfield", "Tables"] -git-tree-sha1 = "e28912ce94077686443433c2800104b061a827ed" +git-tree-sha1 = "7aa7ad1682f3d5754e3491bb59b8103cae28e3a3" uuid = "198e06fe-97b7-11e9-32a5-e1d131e6ad66" -version = "0.3.39" +version = "0.3.40" [deps.BangBang.extensions] BangBangChainRulesCoreExt = "ChainRulesCore" @@ -152,9 +152,9 @@ version = "0.1.5" [[deps.Bzip2_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2" +git-tree-sha1 = "9e2a6b69137e6969bab0152632dcb3bc108c8bdd" uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0" -version = "1.0.8+0" +version = "1.0.8+1" [[deps.CEnum]] git-tree-sha1 = "389ad5c84de1ae7cf0e28e381131c98ea87d54fc" @@ -168,10 +168,10 @@ uuid = "2a0fbf3d-bb9c-48f3-b0a9-814d99fd7ab9" version = "0.2.4" [[deps.CUDA]] -deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CUDA_Driver_jll", "CUDA_Runtime_Discovery", "CUDA_Runtime_jll", "Crayons", "DataFrames", "ExprTools", "GPUArrays", "GPUCompiler", "KernelAbstractions", "LLVM", "LLVMLoopInfo", "LazyArtifacts", "Libdl", "LinearAlgebra", "Logging", "NVTX", "Preferences", "PrettyTables", "Printf", "Random", "Random123", "RandomNumbers", "Reexport", "Requires", "SparseArrays", "Statistics", "UnsafeAtomicsLLVM"] -git-tree-sha1 = "95ac638373ac40e29c1b6d086a3698f5026ff6a6" +deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CUDA_Driver_jll", "CUDA_Runtime_Discovery", "CUDA_Runtime_jll", "Crayons", "DataFrames", "ExprTools", "GPUArrays", "GPUCompiler", "KernelAbstractions", "LLVM", "LLVMLoopInfo", "LazyArtifacts", "Libdl", "LinearAlgebra", "Logging", "NVTX", "Preferences", "PrettyTables", "Printf", "Random", "Random123", "RandomNumbers", "Reexport", "Requires", "SparseArrays", "StaticArrays", "Statistics"] +git-tree-sha1 = "baa8ea7a1ea63316fa3feb454635215773c9c845" uuid = "052768ef-5323-5732-b1bb-66c8b64840ba" -version = "5.1.2" +version = "5.2.0" weakdeps = ["ChainRulesCore", "SpecialFunctions"] [deps.CUDA.extensions] @@ -186,15 +186,15 @@ version = "0.7.0+1" [[deps.CUDA_Runtime_Discovery]] deps = ["Libdl"] -git-tree-sha1 = "bcc4a23cbbd99c8535a5318455dcf0f2546ec536" +git-tree-sha1 = "2cb12f6b2209f40a4b8967697689a47c50485490" uuid = "1af6417a-86b4-443c-805f-a4643ffb695f" -version = "0.2.2" +version = "0.2.3" [[deps.CUDA_Runtime_jll]] deps = ["Artifacts", "CUDA_Driver_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"] -git-tree-sha1 = "9704e50c9158cf8896c2776b8dbc5edd136caf80" +git-tree-sha1 = "8e25c009d2bf16c2c31a70a6e9e8939f7325cc84" uuid = "76a88914-d11a-5bdc-97e0-2f5a05c973a2" -version = "0.10.1+0" +version = "0.11.1+0" [[deps.Cairo_jll]] deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"] @@ -204,15 +204,15 @@ version = "1.16.1+1" [[deps.ChainRules]] deps = ["Adapt", "ChainRulesCore", "Compat", "Distributed", "GPUArraysCore", "IrrationalConstants", "LinearAlgebra", "Random", "RealDot", "SparseArrays", "SparseInverseSubset", "Statistics", "StructArrays", "SuiteSparse"] -git-tree-sha1 = "0aa0a3dd7b9bacbbadf1932ccbdfa938985c5561" +git-tree-sha1 = "213f001d1233fd3b8ef007f50c8cab29061917d8" uuid = "082447d4-558c-5d27-93f4-14fc19e9eca2" -version = "1.58.1" +version = "1.61.0" [[deps.ChainRulesCore]] deps = ["Compat", "LinearAlgebra"] -git-tree-sha1 = "2118cb2765f8197b08e5958cdd17c165427425ee" +git-tree-sha1 = "1287e3872d646eed95198457873249bd9f0caed2" uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -version = "1.19.0" +version = "1.20.1" weakdeps = ["SparseArrays"] [deps.ChainRulesCore.extensions] @@ -224,11 +224,17 @@ git-tree-sha1 = "c0216e792f518b39b22212127d4a84dc31e4e386" uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" version = "1.3.5" +[[deps.CodecBzip2]] +deps = ["Bzip2_jll", "Libdl", "TranscodingStreams"] +git-tree-sha1 = "9b1ca1aa6ce3f71b3d1840c538a8210a043625eb" +uuid = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd" +version = "0.8.2" + [[deps.CodecZlib]] deps = ["TranscodingStreams", "Zlib_jll"] -git-tree-sha1 = "cd67fc487743b2f0fd4380d4cbd3a24660d0eec8" +git-tree-sha1 = "59939d8a997469ee05c4b4944560a820f9ba0d73" uuid = "944b1d66-785c-5afd-91f1-9de20f533193" -version = "0.7.3" +version = "0.7.4" [[deps.ColorSchemes]] deps = ["ColorTypes", "ColorVectorSpace", "Colors", "FixedPointNumbers", "PrecompileTools", "Random"] @@ -270,10 +276,10 @@ uuid = "bbf7d656-a473-5ed7-a52c-81e309532950" version = "0.3.0" [[deps.Compat]] -deps = ["UUIDs"] -git-tree-sha1 = "886826d76ea9e72b35fcd000e535588f7b60f21d" +deps = ["TOML", "UUIDs"] +git-tree-sha1 = "75bd5b6fc5089df449b5d35fa501c846c9b6549b" uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "4.10.1" +version = "4.12.0" weakdeps = ["Dates", "LinearAlgebra"] [deps.Compat.extensions] @@ -347,9 +353,9 @@ uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f" version = "4.1.1" [[deps.DataAPI]] -git-tree-sha1 = "8da84edb865b0b5b0100c0666a9bc9a0b71c553c" +git-tree-sha1 = "abe83f3a2f1b857aac70ef8b269080af17764bbe" uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" -version = "1.15.0" +version = "1.16.0" [[deps.DataFrames]] deps = ["Compat", "DataAPI", "DataStructures", "Future", "InlineStrings", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrecompileTools", "PrettyTables", "Printf", "REPL", "Random", "Reexport", "SentinelArrays", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"] @@ -473,9 +479,9 @@ version = "4.4.2+2" [[deps.FFTW]] deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"] -git-tree-sha1 = "ec22cbbcd01cba8f41eecd7d44aac1f23ee985e3" +git-tree-sha1 = "4820348781ae578893311153d69049a93d05f39d" uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" -version = "1.7.2" +version = "1.8.0" [[deps.FFTW_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -616,15 +622,15 @@ version = "3.3.9+0" [[deps.GPUArrays]] deps = ["Adapt", "GPUArraysCore", "LLVM", "LinearAlgebra", "Printf", "Random", "Reexport", "Serialization", "Statistics"] -git-tree-sha1 = "85d7fb51afb3def5dcb85ad31c3707795c8bccc1" +git-tree-sha1 = "47e4686ec18a9620850bad110b79966132f14283" uuid = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7" -version = "9.1.0" +version = "10.0.2" [[deps.GPUArraysCore]] deps = ["Adapt"] -git-tree-sha1 = "2d6ca471a6c7b536127afccfa7564b5b39227fe0" +git-tree-sha1 = "ec632f177c0d990e64d955ccc1b8c04c485a0950" uuid = "46192b85-c4d5-4398-a991-12ede77f4527" -version = "0.1.5" +version = "0.1.6" [[deps.GPUCompiler]] deps = ["ExprTools", "InteractiveUtils", "LLVM", "Libdl", "Logging", "Scratch", "TimerOutputs", "UUIDs"] @@ -687,9 +693,9 @@ version = "0.1.16" [[deps.Hyperscript]] deps = ["Test"] -git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9" +git-tree-sha1 = "179267cfa5e712760cd43dcae385d7ea90cc25a4" uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91" -version = "0.0.4" +version = "0.0.5" [[deps.HypertextLiteral]] deps = ["Tricks"] @@ -699,15 +705,15 @@ version = "0.9.5" [[deps.IOCapture]] deps = ["Logging", "Random"] -git-tree-sha1 = "d75853a0bdbfb1ac815478bacd89cd27b550ace6" +git-tree-sha1 = "8b72179abc660bfab5e28472e019392b97d0985c" uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89" -version = "0.2.3" +version = "0.2.4" [[deps.IRTools]] deps = ["InteractiveUtils", "MacroTools", "Test"] -git-tree-sha1 = "8aa91235360659ca7560db43a7d57541120aa31d" +git-tree-sha1 = "5d8c5713f38f7bc029e26627b687710ba406d0dd" uuid = "7869d1d1-7146-5819-86e3-90919afe41df" -version = "0.4.11" +version = "0.4.12" [[deps.IfElse]] git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1" @@ -727,10 +733,10 @@ uuid = "c817782e-172a-44cc-b673-b171935fbb9e" version = "0.1.7" [[deps.ImageCore]] -deps = ["AbstractFFTs", "ColorVectorSpace", "Colors", "FixedPointNumbers", "MappedArrays", "MosaicViews", "OffsetArrays", "PaddedViews", "PrecompileTools", "Reexport"] -git-tree-sha1 = "fc5d1d3443a124fde6e92d0260cd9e064eba69f8" +deps = ["ColorVectorSpace", "Colors", "FixedPointNumbers", "MappedArrays", "MosaicViews", "OffsetArrays", "PaddedViews", "PrecompileTools", "Reexport"] +git-tree-sha1 = "b2a7eaa169c13f5bcae8131a83bc30eff8f71be0" uuid = "a09fc81d-aa75-5fe9-8630-4744c3626534" -version = "0.10.1" +version = "0.10.2" [[deps.ImageIO]] deps = ["FileIO", "IndirectArrays", "JpegTurbo", "LazyModules", "Netpbm", "OpenEXR", "PNGFiles", "QOI", "Sixel", "TiffImages", "UUIDs"] @@ -770,9 +776,9 @@ version = "3.1.7+0" [[deps.IndexFunArrays]] deps = ["ChainRulesCore", "LinearAlgebra"] -git-tree-sha1 = "e66a2aeb6d5814015004080e5203dfff44d2856f" +git-tree-sha1 = "6f78703c7a4ba06299cddd8694799c91de0157ac" uuid = "613c443e-d742-454e-bfc6-1d7f8dd76566" -version = "0.2.6" +version = "0.2.7" [[deps.IndirectArrays]] git-tree-sha1 = "012e604e1c7458645cb8b436f8fba789a51b257f" @@ -883,9 +889,9 @@ version = "3.0.1+0" [[deps.JuliaInterpreter]] deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"] -git-tree-sha1 = "e49bce680c109bc86e3e75ebcb15040d6ad9e1d3" +git-tree-sha1 = "04663b9e1eb0d0eabf76a6d0752e0dac83d53b36" uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a" -version = "0.9.27" +version = "0.9.28" [[deps.JuliaNVTXCallbacks_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -901,9 +907,9 @@ version = "0.2.4" [[deps.KernelAbstractions]] deps = ["Adapt", "Atomix", "InteractiveUtils", "LinearAlgebra", "MacroTools", "PrecompileTools", "Requires", "SparseArrays", "StaticArrays", "UUIDs", "UnsafeAtomics", "UnsafeAtomicsLLVM"] -git-tree-sha1 = "653e0824fc9ab55b3beec67a6dbbe514a65fb954" +git-tree-sha1 = "4e0cb2f5aad44dcfdc91088e85dee4ecb22c791c" uuid = "63c18a36-062a-441e-b654-da1e3ab1ce7c" -version = "0.9.15" +version = "0.9.16" [deps.KernelAbstractions.extensions] EnzymeExt = "EnzymeCore" @@ -1100,9 +1106,9 @@ uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" [[deps.LinearMaps]] deps = ["LinearAlgebra"] -git-tree-sha1 = "9df2ab050ffefe870a09c7b6afdb0cde381703f2" +git-tree-sha1 = "9948d6f8208acfebc3e8cf4681362b2124339e7e" uuid = "7a12625a-238d-50fd-b39a-03d52299707e" -version = "3.11.1" +version = "3.11.2" weakdeps = ["ChainRulesCore", "SparseArrays", "Statistics"] [deps.LinearMaps.extensions] @@ -1169,9 +1175,9 @@ version = "1.0.3" [[deps.LoweredCodeUtils]] deps = ["JuliaInterpreter"] -git-tree-sha1 = "38756922d32476c8f41f73560b910fc805a5a103" +git-tree-sha1 = "20ce1091ba18bcdae71ad9b71ee2367796ba6c48" uuid = "6f1432cf-f94c-5a45-995e-cdbf5db27b0b" -version = "2.4.3" +version = "2.4.4" [[deps.MIMEs]] git-tree-sha1 = "65f28ad4b594aebe22157d6fac869786a255b7eb" @@ -1191,9 +1197,9 @@ version = "0.4.17" [[deps.MacroTools]] deps = ["Markdown", "Random"] -git-tree-sha1 = "b211c553c199c111d998ecdaf7623d1b89b69f93" +git-tree-sha1 = "2fa9ee3e63fd3a4f7a9a4f4744a52f4856de82df" uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" -version = "0.5.12" +version = "0.5.13" [[deps.Malt]] deps = ["Distributed", "Logging", "RelocatableFolders", "Serialization", "Sockets"] @@ -1215,6 +1221,12 @@ version = "0.4.2" deps = ["Base64"] uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" +[[deps.MathOptInterface]] +deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "DataStructures", "ForwardDiff", "JSON", "LinearAlgebra", "MutableArithmetics", "NaNMath", "OrderedCollections", "PrecompileTools", "Printf", "SparseArrays", "SpecialFunctions", "Test", "Unicode"] +git-tree-sha1 = "8b40681684df46785a0012d352982e22ac3be59e" +uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" +version = "1.25.2" + [[deps.MbedTLS]] deps = ["Dates", "MbedTLS_jll", "MozillaCACerts_jll", "NetworkOptions", "Random", "Sockets"] git-tree-sha1 = "c067a280ddc25f196b5e7df3877c6b226d390aaf" @@ -1262,6 +1274,12 @@ git-tree-sha1 = "f5db02ae992c260e4826fe78c942954b48e1d9c2" uuid = "99f44e22-a591-53d1-9472-aa23ef4bd671" version = "1.2.1" +[[deps.MutableArithmetics]] +deps = ["LinearAlgebra", "SparseArrays", "Test"] +git-tree-sha1 = "806eea990fb41f9b36f1253e5697aa645bf6a9f8" +uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" +version = "1.4.0" + [[deps.NDTools]] deps = ["LinearAlgebra", "OffsetArrays", "PaddedViews", "Random", "Statistics"] git-tree-sha1 = "4d5fc006e0a006875f57be883c81d9c4a5d56bc6" @@ -1276,9 +1294,9 @@ version = "7.8.3" [[deps.NVTX]] deps = ["Colors", "JuliaNVTXCallbacks_jll", "Libdl", "NVTX_jll"] -git-tree-sha1 = "8bc9ce4233be3c63f8dcd78ccaf1b63a9c0baa34" +git-tree-sha1 = "53046f0483375e3ed78e49190f1154fa0a4083a1" uuid = "5da4648a-3479-48b8-97b9-01cb529c0a1f" -version = "0.3.3" +version = "0.3.4" [[deps.NVTX_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -1370,10 +1388,10 @@ uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" version = "0.5.5+0" [[deps.Optim]] -deps = ["Compat", "FillArrays", "ForwardDiff", "LineSearches", "LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "PositiveFactorizations", "Printf", "SparseArrays", "StatsBase"] -git-tree-sha1 = "01f85d9269b13fedc61e63cc72ee2213565f7a72" +deps = ["Compat", "FillArrays", "ForwardDiff", "LineSearches", "LinearAlgebra", "MathOptInterface", "NLSolversBase", "NaNMath", "Parameters", "PositiveFactorizations", "Printf", "SparseArrays", "StatsBase"] +git-tree-sha1 = "47fea72de134f75b105a5d4a1abe5c6aec89d390" uuid = "429524aa-4258-5aef-a3af-852621145aeb" -version = "1.7.8" +version = "1.9.1" [[deps.Optimisers]] deps = ["ChainRulesCore", "Functors", "LinearAlgebra", "Random", "Statistics"] @@ -1382,10 +1400,10 @@ uuid = "3bd65402-5787-11e9-1adc-39752487f4e2" version = "0.3.1" [[deps.Optimization]] -deps = ["ADTypes", "ArrayInterface", "ConsoleProgressMonitor", "DocStringExtensions", "LinearAlgebra", "Logging", "LoggingExtras", "Pkg", "Printf", "ProgressLogging", "Reexport", "SciMLBase", "SparseArrays", "SymbolicIndexingInterface", "TerminalLoggers"] -git-tree-sha1 = "a38741a9f0ef4847b892151957b6d01a5da49214" +deps = ["ADTypes", "ArrayInterface", "ConsoleProgressMonitor", "DocStringExtensions", "LinearAlgebra", "Logging", "LoggingExtras", "Pkg", "Printf", "ProgressLogging", "Reexport", "SciMLBase", "SparseArrays", "TerminalLoggers"] +git-tree-sha1 = "145baedf71770d84bf590d307b4686a9ea619c4b" uuid = "7f7a1694-90dd-40f0-9382-eb1efda571ba" -version = "3.21.0" +version = "3.22.0" [deps.Optimization.extensions] OptimizationEnzymeExt = "Enzyme" @@ -1410,15 +1428,15 @@ version = "3.21.0" [[deps.OptimizationOptimJL]] deps = ["Optim", "Optimization", "Reexport", "SparseArrays"] -git-tree-sha1 = "31083728a4a9066edd81b97835df130de0e47d79" +git-tree-sha1 = "6cdc08135141c153b49372de47416ac15fed42de" uuid = "36348300-93cb-4f02-beb5-3c3902f8871e" -version = "0.2.0" +version = "0.2.2" [[deps.OptimizationOptimisers]] deps = ["Optimisers", "Optimization", "Printf", "ProgressLogging", "Reexport"] -git-tree-sha1 = "2b63d33faf7a850564c9a4343daa37744e0f2f82" +git-tree-sha1 = "6c12ab297301d706b542d8876e55bfadfe773b00" uuid = "42dfb2eb-d2b4-4451-abcd-913932933ac1" -version = "0.2.0" +version = "0.2.1" [[deps.Opus_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -1495,10 +1513,10 @@ uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043" version = "1.4.0" [[deps.Plots]] -deps = ["Base64", "Contour", "Dates", "Downloads", "FFMPEG", "FixedPointNumbers", "GR", "JLFzf", "JSON", "LaTeXStrings", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "Pkg", "PlotThemes", "PlotUtils", "PrecompileTools", "Preferences", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "RelocatableFolders", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs", "UnicodeFun", "UnitfulLatexify", "Unzip"] -git-tree-sha1 = "ccee59c6e48e6f2edf8a5b64dc817b6729f99eb5" +deps = ["Base64", "Contour", "Dates", "Downloads", "FFMPEG", "FixedPointNumbers", "GR", "JLFzf", "JSON", "LaTeXStrings", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "Pkg", "PlotThemes", "PlotUtils", "PrecompileTools", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "RelocatableFolders", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs", "UnicodeFun", "UnitfulLatexify", "Unzip"] +git-tree-sha1 = "38a748946dca52a622e79eea6ed35c6737499109" uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" -version = "1.39.0" +version = "1.40.0" [deps.Plots.extensions] FileIOExt = "FileIO" @@ -1516,9 +1534,9 @@ version = "1.39.0" [[deps.Pluto]] deps = ["Base64", "Configurations", "Dates", "Downloads", "ExpressionExplorer", "FileWatching", "FuzzyCompletions", "HTTP", "HypertextLiteral", "InteractiveUtils", "Logging", "LoggingExtras", "MIMEs", "Malt", "Markdown", "MsgPack", "Pkg", "PrecompileSignatures", "PrecompileTools", "REPL", "RegistryInstances", "RelocatableFolders", "Scratch", "Sockets", "TOML", "Tables", "URIs", "UUIDs"] -git-tree-sha1 = "e6a92bf27d9e8eda41b672772ad05f6652513e02" +git-tree-sha1 = "449f468cbb80c3eec6e6d8443a0913d8bbad4d0d" uuid = "c3e4b0f8-55cb-11ea-2926-15256bba5781" -version = "0.19.36" +version = "0.19.38" [[deps.PlutoTest]] deps = ["HypertextLiteral", "InteractiveUtils", "Markdown", "Test"] @@ -1528,9 +1546,9 @@ version = "0.2.2" [[deps.PlutoUI]] deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "FixedPointNumbers", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "MIMEs", "Markdown", "Random", "Reexport", "URIs", "UUIDs"] -git-tree-sha1 = "bd7c69c7f7173097e7b5e1be07cee2b8b7447f51" +git-tree-sha1 = "68723afdb616445c6caaef6255067a8339f91325" uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8" -version = "0.7.54" +version = "0.7.55" [[deps.PoissonRandom]] deps = ["Random"] @@ -1671,21 +1689,25 @@ version = "0.6.12" [[deps.RecursiveArrayTools]] deps = ["Adapt", "ArrayInterface", "DocStringExtensions", "GPUArraysCore", "IteratorInterfaceExtensions", "LinearAlgebra", "RecipesBase", "SparseArrays", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables"] -git-tree-sha1 = "4873672a6d7990c683b4badf8671c50043ad321b" +git-tree-sha1 = "2bd309f5171a628efdf5309361cd8a779b9e63a9" uuid = "731186ca-8d62-57ce-b412-fbd966d074cd" -version = "3.4.2" +version = "3.8.0" [deps.RecursiveArrayTools.extensions] RecursiveArrayToolsFastBroadcastExt = "FastBroadcast" + RecursiveArrayToolsForwardDiffExt = "ForwardDiff" RecursiveArrayToolsMeasurementsExt = "Measurements" RecursiveArrayToolsMonteCarloMeasurementsExt = "MonteCarloMeasurements" + RecursiveArrayToolsReverseDiffExt = ["ReverseDiff", "Zygote"] RecursiveArrayToolsTrackerExt = "Tracker" RecursiveArrayToolsZygoteExt = "Zygote" [deps.RecursiveArrayTools.weakdeps] FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" + ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca" + ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" @@ -1702,9 +1724,9 @@ version = "0.1.0" [[deps.RegularizedLeastSquares]] deps = ["FLoops", "InteractiveUtils", "IterativeSolvers", "LinearAlgebra", "LinearOperatorCollection", "LinearOperators", "Random", "SparseArrays", "Statistics", "StatsBase", "VectorizationBase"] -git-tree-sha1 = "a910c2559719b06d64e8e53b16c0bc7350279972" +git-tree-sha1 = "5cfd4be08e8fe095b9d333d5287698b7224b32fc" uuid = "1e9c538a-f78c-5de5-8ffb-0b6dbe892d23" -version = "0.12.0" +version = "0.13.0" [[deps.RelocatableFolders]] deps = ["SHA", "Scratch"] @@ -1741,9 +1763,9 @@ version = "0.1.0" [[deps.SciMLBase]] deps = ["ADTypes", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "FillArrays", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "PrecompileTools", "Preferences", "Printf", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLOperators", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables", "TruncatedStacktraces"] -git-tree-sha1 = "c3a0fe8eb1cfcda4c94b073f4d456c2641d99d3e" +git-tree-sha1 = "75bae786dc8b07ec3c2159d578886691823bcb42" uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462" -version = "2.16.5" +version = "2.23.1" [deps.SciMLBase.extensions] SciMLBaseChainRulesCoreExt = "ChainRulesCore" @@ -1877,9 +1899,9 @@ weakdeps = ["OffsetArrays", "StaticArrays"] [[deps.StaticArrays]] deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"] -git-tree-sha1 = "4e17a790909b17f7bf1496e3aec138cf01b60b3b" +git-tree-sha1 = "7b0e9c14c624e435076d19aea1e5cbdec2b9ca37" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.9.0" +version = "1.9.2" weakdeps = ["ChainRulesCore", "Statistics"] [deps.StaticArrays.extensions] @@ -1922,9 +1944,9 @@ version = "0.3.4" [[deps.StructArrays]] deps = ["Adapt", "ConstructionBase", "DataAPI", "GPUArraysCore", "StaticArraysCore", "Tables"] -git-tree-sha1 = "0a3db38e4cce3c54fe7a71f831cd7b6194a54213" +git-tree-sha1 = "1b0b1205a56dc288b71b1961d48e351520702e24" uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" -version = "0.6.16" +version = "0.6.17" [[deps.SuiteSparse]] deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"] @@ -1936,9 +1958,9 @@ uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c" version = "7.2.1+1" [[deps.SymbolicIndexingInterface]] -git-tree-sha1 = "74502f408d99fc217a9d7cd901d9ffe45af892b1" +git-tree-sha1 = "b3103f4f50a3843e66297a2456921377c78f5e31" uuid = "2efcf032-c050-4f8e-a9bb-153293bab1f5" -version = "0.3.3" +version = "0.3.5" [[deps.TOML]] deps = ["Dates"] @@ -1997,9 +2019,9 @@ uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" version = "0.5.23" [[deps.TranscodingStreams]] -git-tree-sha1 = "1fbeaaca45801b4ba17c251dd8603ef24801dd84" +git-tree-sha1 = "54194d92959d8ebaa8e26227dbe3cdefcdcd594f" uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" -version = "0.10.2" +version = "0.10.3" weakdeps = ["Random", "Test"] [deps.TranscodingStreams.extensions] @@ -2007,9 +2029,9 @@ weakdeps = ["Random", "Test"] [[deps.Transducers]] deps = ["Adapt", "ArgCheck", "BangBang", "Baselet", "CompositionsBase", "ConstructionBase", "DefineSingletons", "Distributed", "InitialValues", "Logging", "Markdown", "MicroCollections", "Requires", "Setfield", "SplittablesBase", "Tables"] -git-tree-sha1 = "e579d3c991938fecbb225699e8f611fa3fbf2141" +git-tree-sha1 = "3064e780dbb8a9296ebb3af8f440f787bb5332af" uuid = "28d57a85-8fef-5791-bfe6-a80928e7c999" -version = "0.4.79" +version = "0.4.80" [deps.Transducers.extensions] TransducersBlockArraysExt = "BlockArrays" @@ -2288,9 +2310,9 @@ version = "1.5.5+0" [[deps.Zygote]] deps = ["AbstractFFTs", "ChainRules", "ChainRulesCore", "DiffRules", "Distributed", "FillArrays", "ForwardDiff", "GPUArrays", "GPUArraysCore", "IRTools", "InteractiveUtils", "LinearAlgebra", "LogExpFunctions", "MacroTools", "NaNMath", "PrecompileTools", "Random", "Requires", "SparseArrays", "SpecialFunctions", "Statistics", "ZygoteRules"] -git-tree-sha1 = "30c1b8bfc2b3c7c5d8bba7cd32e8b6d5f968e7c3" +git-tree-sha1 = "4ddb4470e47b0094c93055a3bcae799165cc68f1" uuid = "e88e6eb3-aa80-5325-afca-941959d7151f" -version = "0.6.68" +version = "0.6.69" [deps.Zygote.extensions] ZygoteColorsExt = "Colors" @@ -2304,9 +2326,9 @@ version = "0.6.68" [[deps.ZygoteRules]] deps = ["ChainRulesCore", "MacroTools"] -git-tree-sha1 = "9d749cd449fb448aeca4feee9a2f4186dbb5d184" +git-tree-sha1 = "27798139afc0a2afa7b1824c206d5e87ea587a00" uuid = "700de1a5-db45-46bc-99cf-38207098b444" -version = "0.2.4" +version = "0.2.5" [[deps.fzf_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] diff --git a/examples/Project.toml b/examples/Project.toml index b44d446..87417ba 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -23,6 +23,8 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Pluto = "c3e4b0f8-55cb-11ea-2926-15256bba5781" PlutoTest = "cb4044da-4d16-4ffa-a6a3-8cad7f73ebdc" PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8" +ProgressLogging = "33c8b6b6-d38a-422a-b730-caa89a2f386c" +ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" RadonKA = "86de8297-835b-47df-b249-c04e8db91db5" RegularizedLeastSquares = "1e9c538a-f78c-5de5-8ffb-0b6dbe892d23" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" diff --git a/examples/Tomographic_Volumetric_Additive_Manufacturing_with_Refraction.jl b/examples/Tomographic_Volumetric_Additive_Manufacturing_with_Refraction.jl new file mode 100644 index 0000000..ce9e78c --- /dev/null +++ b/examples/Tomographic_Volumetric_Additive_Manufacturing_with_Refraction.jl @@ -0,0 +1,358 @@ +### A Pluto.jl notebook ### +# v0.19.38 + +using Markdown +using InteractiveUtils + +# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error). +macro bind(def, element) + quote + local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end + local el = $(esc(element)) + global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el) + el + end +end + +# ╔═╡ 179d107e-b3be-11ee-0a6c-49f1bf2a10fd +begin + using Pkg + Pkg.activate(".") + Pkg.instantiate() + using Revise +end + +# ╔═╡ e3a3dc83-6536-445f-b56b-d8ad62cd9a85 +using TestImages, RadonKA, ImageShow, ImageIO, Noise, PlutoUI, BenchmarkTools, CUDA, Zygote, IndexFunArrays, FileIO, NDTools, Plots, ProgressLogging, Optim + +# ╔═╡ eb4b7177-4f93-4728-8c91-4ff6a86f40cf +md"# Load packages and check CUDA" + +# ╔═╡ f3277f00-4162-43b2-98e0-53bcb508df1e +TableOfContents() + +# ╔═╡ 42f31ceb-0f5d-46cc-b868-0314b41a5407 +begin + use_CUDA = Ref(CUDA.functional()) + var"@mytime" = use_CUDA[] ? CUDA.var"@time" : var"@time" + togoc(x) = use_CUDA[] ? CuArray(x) : x + togoc(x::AbstractArray{Float64}) = use_CUDA[] ? CuArray(Float32.(x)) : x +end + +# ╔═╡ 78746412-00f2-40ad-9d24-e884d8dd7e77 +md"""# CUDA is functional: $(CUDA.functional() ? "yes" : "no")""" + +# ╔═╡ 980e9656-202d-4a47-837f-b2197fe5d3b4 +md"# Load JOSS logo as target" + +# ╔═╡ f0d86933-5eab-4d86-8909-c51ab9ca99d9 +background = select_region(ImageShow.Gray.(load(download("https://upload.wikimedia.org/wikipedia/commons/thumb/9/94/JOSS.png/459px-JOSS.png"))) .< 0.001, new_size=(660,660)); + +# ╔═╡ 36eab63f-5643-4b43-8e12-0f7d971c1040 +img = select_region(.-ImageShow.Gray.(load(download("https://upload.wikimedia.org/wikipedia/commons/thumb/9/94/JOSS.png/459px-JOSS.png"))), new_size=(660,660)); + +# ╔═╡ fc049fbc-4ae7-4256-b39f-27fdff98d45a +target = togoc(Float32.((img .- background) .> 0.1)); + +# ╔═╡ a6bdc90f-3c10-4524-9892-6e976d556a95 +simshow(Array(target)) + +# ╔═╡ 066664f7-9ae5-4a4d-b541-663067cea361 +md"# Specify Angles and ray endpoints + +We consider that we have a vial with outer radius 7mm and a inner radius of 6.5mm. The vial has a refractive index of 1.5 whereas the resin has one of 1.48. + +The full derivation of `distort_rays_vial` is linked in this [PDF](https://github.com/roflmaostc/RadonKA.jl/blob/main/docs/src/angle_derivation/refraction_inside_vial_resin.tex). +" + +# ╔═╡ 5ef75bb5-5ab0-4e98-9026-a8e2c5b7833d +angles = togoc(range(0, 2f0 * π, 400)); + +# ╔═╡ 5aa03590-c552-410e-93aa-1e47dabd2f9e +N_s = size(target,1) - 1 + +# ╔═╡ 9d582f0d-7da4-41db-9369-fef7a6dbe590 +range(-N_s/2 - 1, N_s / 2 - 1, N_s) + +# ╔═╡ 6b91d2bd-1170-4ce5-bb20-10d69f8cd6d7 +md" +The `ray_startpoints` and `ray_endpoints` indicate each the beginning and end position of the ray. +" + +# ╔═╡ 8f212882-5a4f-4915-a01a-9c7f342ea1ef +Rvial = 7f-3 + +# ╔═╡ 07ff2197-836d-4e9e-b716-4e0eacd89e6e +Rresin = 6.5f-3 + +# ╔═╡ 6deeaf05-1008-4f1e-bb33-2b74196c9d80 +nvial = 1.5f0 + +# ╔═╡ b94f82e8-1c24-46a0-b3e9-f98a0e3e1c12 +nresin = 1.48f0 + +# ╔═╡ 66a09fd3-545b-4456-be74-6593d2d15a28 +radius_pixel = N_s/2f0 + +# ╔═╡ af83cb81-49b8-41f4-a24e-c1e33e76d925 +function distort_rays_vial(y::T, Rₒ::T, Rᵢ::T, nvial::T, nresin::T) where T + y, Rₒ, Rᵢ, nvial, nresin = Float64.((y, Rₒ, Rᵢ, nvial, nresin)) + + if iszero(y) + return zero(T), zero(T) + end + + α = asin(y / Rₒ) + β = asin(sin(α) / nvial) + + x = Rₒ * cos(β) - sqrt(Rₒ^2*(cos(β)^2 - 1) + Rᵢ^2) + + ϵ = acos((-Rₒ^2 + x^2 + Rᵢ^2) / 2 / Rᵢ / x) - Float64(π) / 2 + + β_ = sign(y) * (Float64(π) / 2 - ϵ) + γ = asin(nvial * sin(β_) / nresin) + + δ₁ = α - β + δ₂ = β_ - γ + δ_ges = δ₁ + δ₂ + + y_ = abs(Rᵢ * sin(γ)) + p = sqrt(Rₒ^2 - y_^2) + + η = - (asin(p / Rₒ) - sign(y) * (Float64(π)/2 - δ_ges)) + yf = Rₒ * sin(η) + + yi = 2 * p * sin(δ_ges) + yf + + return T(yi), T(yf) +end + +# ╔═╡ 6721b65f-48ed-4dc7-99dd-c460bb124c1e +ray_points = distort_rays_vial.(range(-N_s/2f0 + 1f0, N_s / 2f0 - 1f0, N_s), + radius_pixel, + radius_pixel / Rvial * Rresin, + nvial, + nresin) + +# ╔═╡ 76b4c0de-af7a-456e-8d49-04f10a0a17fb +ray_startpoints, ray_endpoints = map(first, ray_points), map(x->x[2], ray_points) + +# ╔═╡ 4cdc5592-3b92-4e98-8718-d77dd939ebf1 +md"## Ray distribution inside vial" + +# ╔═╡ 494c5b20-0af2-4903-959c-812459877c31 +kk = 1 + +# ╔═╡ dc22dead-1b83-4787-a147-9811136af085 +begin + sinogram = zeros((length(ray_startpoints), kk)) + sinogram[1:1:end,:] .= 1 +end; + +# ╔═╡ e10549a7-0f31-48fa-8d38-820907d8554e +Revise.errors() + +# ╔═╡ e83be634-dde4-4ed4-9226-2ba895a84d15 +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) + +# ╔═╡ 5d55464f-c7d9-4e45-a1e7-907c2137be95 +md"# Define Optimization algorithm +See this paper: +[*Charles M. Rackson, Kyle M. Champley, Joseph T. Toombs, Erika J. Fong, Vishal Bansal, Hayden K. Taylor, Maxim Shusteff, Robert R. McLeod, Object-space optimization of tomographic reconstructions for additive manufacturing*](https://www.sciencedirect.com/science/article/abs/pii/S2214860421005212) +" + +# ╔═╡ 834cd0f6-9be0-4589-82c0-ca0021911dd5 +function make_fg!(fwd_operator, target, threshold) + + isobject = target .== 1 + notobject = target .== 0 + f = let fwd_operator=fwd_operator + # apply sqrt for anscombe + f(x) = begin + fwd = fwd_operator(x) + return (sum(abs2, max.(0, fwd[notobject] .- threshold[1])) + + sum(abs2, max.(0, threshold[2] .- fwd[isobject])) + + sum(abs2, max.(0, fwd[isobject] .- 1))) + end + end + + g! = let f=f + g!(G, x) = begin + if !isnothing(G) + G .= Zygote.gradient(f, x)[1] + end + end + end + + return f, g! +end + +# ╔═╡ 697d4c7d-697b-452c-8a6a-25f2a1089529 +function iter!(buffer, img, θs, μ; clip_sinogram=true, geometry) + sinogram = radon(img, θs; μ, geometry) + + if clip_sinogram + sinogram .= max.(sinogram, 0) + end + + img_recon = iradon(sinogram, θs; μ, geometry) + img_recon ./= maximum(img_recon) + buffer .= max.(img_recon, 0) + return buffer, sinogram +end + +# ╔═╡ 624858de-e7da-4e3c-91b5-6ab7b5e645d5 +# see https://www.sciencedirect.com/science/article/abs/pii/S2214860421005212 +function OSMO(img::AbstractArray{T}, θs, μ=nothing; + thresholds=(0.65, 0.75), N_iter = 2, + geometry) where T + N = size(img, 1) + guess = copy(img) + + notobject = togoc(iszero.(img)) + isobject = togoc(isone.(img)) + + losses = T[] + buffer = copy(img) + tmp, s = iter!(buffer, guess, θs, μ; clip_sinogram=true, geometry) + @progress for i in 1:N_iter + guess[notobject] .-= max.(0, tmp[notobject] .- thresholds[1]) + tmp, s = iter!(buffer, guess, θs, μ; clip_sinogram=true, geometry) + guess[isobject] .+= max.(0, thresholds[2] .- tmp[isobject]) + end + + printed = iradon(s, θs; μ, geometry) + printed ./= maximum(printed) + return guess, s, printed +end + +# ╔═╡ 12e7c6b7-fcf3-44e4-9d80-785704d93cb5 +fwd = x -> iradon(max.(0,x), angles; geometry) + +# ╔═╡ 942e1cf7-a47b-43af-a460-1ea896a2c9b8 + f, g! = make_fg!(fwd, target, (0.65, 0.75)) + +# ╔═╡ d80f6c47-1c55-4a7e-9374-104806873a9a +md"# Run optimization +With a decent GPU (such as RTX 3060) it should take around ~20-30s. +With a 12 core multithreaded CPU around ~15min. +" + +# ╔═╡ be45bfe1-e91e-45eb-8902-62e8df0d8ac7 +begin + init0 = radon(target, angles); + init0 .= 0 +end; + +# ╔═╡ e14a65c3-ece0-4d1e-96a7-9a9ab00001ad +@time res = Optim.optimize(f, g!, init0, LBFGS(), + Optim.Options(iterations = 15, + store_trace=true)) + +# ╔═╡ 2ca8698e-3ff8-445c-a5e3-3a20d2e2afc3 +@mytime a_object, patterns, printed = OSMO(target, angles, 1/350f0; N_iter=200, + geometry) + +# ╔═╡ ec528b63-95eb-430b-b683-04e207dd61f0 +simshow((Array(printed) .> 0.7) .- 0 .* Array(target), cmap=:turbo) + +# ╔═╡ 43e4b941-b838-483d-9a8f-8897720f8b90 +simshow(Array(printed), cmap=:turbo) + +# ╔═╡ 45884ebb-a02d-446a-a7ca-05b8b4abc244 +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) + +# ╔═╡ 28d8c378-b818-4d61-b9d4-a06ad7cc21e3 +simshow(Array(iradon(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) + +# ╔═╡ 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) + +# ╔═╡ 21c0b369-a4b2-43ea-87b8-702474d88584 +simshow(Array(patterns)) + +# ╔═╡ b53c8016-177f-4b1b-b07d-48cbd2b2fc65 +function plot_histogram(target, object_printed, thresholds; yscale=:log10) + # :stephist vs :barhist + + plot_font = "Computer Modern" + default(fontfamily=plot_font, + linewidth=2, framestyle=:box, label=nothing, grid=false) + plot(object_printed[target .== 0], seriestype=:stephist, bins=(0.0:0.01:1), xlim=(0.0, 1.0), label="dose distribution void", ylabel="voxel count", xlabel="normalized intensity", ylim=(10, 500000), linewidth=1, legend=:topleft, yscale=yscale, size=(500, 300)) + plot!(object_printed[target .== 1], seriestype=:stephist, bins=(0.0:0.01:1), xlim=(0.0, 1.0), label="dose distribution object", ylabel="voxel count", xlabel="normalized intensity", ylim=(10, 500000), linewidth=1, legend=:topleft, yscale=yscale, size=(500, 300)) + plot!([thresholds[1], thresholds[1]], [1, 10000_000], label="lower threshold", linewidth=3) + plot!([thresholds[2], thresholds[2]], [1, 10000_000], label="upper threshold", linewidth=3) + #plot!([chosen_threshold, chosen_threshold], [1, 30000000], label="chosen threshold", linewidth=3) +end + +# ╔═╡ 4e3274ed-ad71-4568-a930-c50abe7c6bd5 +plot_histogram(Array(target), Array(fwd(res.minimizer)), (0.65, 0.75)) + +# ╔═╡ 03223eb4-3c15-4899-93bf-5440fb6b5c7a +plot_histogram(Array(target), Array(printed), (0.65, 0.75)) + +# ╔═╡ Cell order: +# ╠═179d107e-b3be-11ee-0a6c-49f1bf2a10fd +# ╟─eb4b7177-4f93-4728-8c91-4ff6a86f40cf +# ╠═e3a3dc83-6536-445f-b56b-d8ad62cd9a85 +# ╠═f3277f00-4162-43b2-98e0-53bcb508df1e +# ╠═42f31ceb-0f5d-46cc-b868-0314b41a5407 +# ╟─78746412-00f2-40ad-9d24-e884d8dd7e77 +# ╟─980e9656-202d-4a47-837f-b2197fe5d3b4 +# ╠═f0d86933-5eab-4d86-8909-c51ab9ca99d9 +# ╠═36eab63f-5643-4b43-8e12-0f7d971c1040 +# ╠═fc049fbc-4ae7-4256-b39f-27fdff98d45a +# ╠═a6bdc90f-3c10-4524-9892-6e976d556a95 +# ╟─066664f7-9ae5-4a4d-b541-663067cea361 +# ╠═5ef75bb5-5ab0-4e98-9026-a8e2c5b7833d +# ╠═5aa03590-c552-410e-93aa-1e47dabd2f9e +# ╠═9d582f0d-7da4-41db-9369-fef7a6dbe590 +# ╟─6b91d2bd-1170-4ce5-bb20-10d69f8cd6d7 +# ╠═8f212882-5a4f-4915-a01a-9c7f342ea1ef +# ╠═07ff2197-836d-4e9e-b716-4e0eacd89e6e +# ╠═6deeaf05-1008-4f1e-bb33-2b74196c9d80 +# ╠═b94f82e8-1c24-46a0-b3e9-f98a0e3e1c12 +# ╠═66a09fd3-545b-4456-be74-6593d2d15a28 +# ╠═6721b65f-48ed-4dc7-99dd-c460bb124c1e +# ╠═76b4c0de-af7a-456e-8d49-04f10a0a17fb +# ╠═af83cb81-49b8-41f4-a24e-c1e33e76d925 +# ╟─4cdc5592-3b92-4e98-8718-d77dd939ebf1 +# ╠═dc22dead-1b83-4787-a147-9811136af085 +# ╠═494c5b20-0af2-4903-959c-812459877c31 +# ╠═e10549a7-0f31-48fa-8d38-820907d8554e +# ╠═e83be634-dde4-4ed4-9226-2ba895a84d15 +# ╠═85b6e54c-3256-4ba8-9889-56be82de0203 +# ╟─5d55464f-c7d9-4e45-a1e7-907c2137be95 +# ╠═834cd0f6-9be0-4589-82c0-ca0021911dd5 +# ╠═697d4c7d-697b-452c-8a6a-25f2a1089529 +# ╠═624858de-e7da-4e3c-91b5-6ab7b5e645d5 +# ╠═12e7c6b7-fcf3-44e4-9d80-785704d93cb5 +# ╠═942e1cf7-a47b-43af-a460-1ea896a2c9b8 +# ╟─d80f6c47-1c55-4a7e-9374-104806873a9a +# ╠═be45bfe1-e91e-45eb-8902-62e8df0d8ac7 +# ╠═4e3274ed-ad71-4568-a930-c50abe7c6bd5 +# ╠═e14a65c3-ece0-4d1e-96a7-9a9ab00001ad +# ╠═2ca8698e-3ff8-445c-a5e3-3a20d2e2afc3 +# ╠═03223eb4-3c15-4899-93bf-5440fb6b5c7a +# ╠═ec528b63-95eb-430b-b683-04e207dd61f0 +# ╠═43e4b941-b838-483d-9a8f-8897720f8b90 +# ╠═45884ebb-a02d-446a-a7ca-05b8b4abc244 +# ╠═7c998816-b4de-4dbf-977b-bdab9355be79 +# ╠═28d8c378-b818-4d61-b9d4-a06ad7cc21e3 +# ╠═beb61bb1-3f7c-4301-92ca-5418bc0445cf +# ╠═61048e10-8e23-41a6-8c49-0b1a8c0adc24 +# ╠═88895d8d-432e-4337-a51c-12ac4e4a4325 +# ╠═21c0b369-a4b2-43ea-87b8-702474d88584 +# ╠═b53c8016-177f-4b1b-b07d-48cbd2b2fc65 diff --git a/examples/benchmark_REPL.jl b/examples/benchmark_REPL.jl new file mode 100644 index 0000000..6bf8e71 --- /dev/null +++ b/examples/benchmark_REPL.jl @@ -0,0 +1,21 @@ +using RadonKA, CUDA, BenchmarkTools + + + + +for sz in [(512,512,100), (256, 256)] + for d in [Array, CuArray] + @show d, sz + angles = d(range(0f0, 360f0, 360)) + arr = d(randn(Float32, sz)) + iarr = radon(arr, angles) + + if d === Array + @btime radon($arr, $angles) + @btime iradon($iarr, $angles) + else + @btime CUDA.@sync radon($arr, $angles) + @btime CUDA.@sync iradon($iarr, $angles) + end + end +end diff --git a/examples/example_radon_iradon.jl b/examples/example_radon_iradon.jl index 190b008..b7a6bb0 100644 --- a/examples/example_radon_iradon.jl +++ b/examples/example_radon_iradon.jl @@ -1,5 +1,5 @@ ### A Pluto.jl notebook ### -# v0.19.36 +# v0.19.38 using Markdown using InteractiveUtils diff --git a/examples/playground.jl b/examples/playground.jl new file mode 100644 index 0000000..becc264 --- /dev/null +++ b/examples/playground.jl @@ -0,0 +1,386 @@ +### A Pluto.jl notebook ### +# v0.19.38 + +using Markdown +using InteractiveUtils + +# ╔═╡ 4eb3148e-8f8b-11ee-3cfe-854d3bd5cc80 +begin + using Pkg + Pkg.activate(".") + Pkg.instantiate() + using Revise +end + +# ╔═╡ b336e55e-0be4-422f-b48a-0a2242cb0915 +using RadonKA + +# ╔═╡ 1311e853-c4cd-42bb-8bf3-5e0d564bf9c5 +using IndexFunArrays, ImageShow, Plots, ImageIO, PlutoUI, PlutoTest, TestImages + +# ╔═╡ e90a05dd-5781-44ef-9b39-5b1ee85ca477 +using BenchmarkTools + +# ╔═╡ 03bccb92-b47f-477a-9bdb-74cc404da690 +using KernelAbstractions, CUDA, CUDA.CUDAKernels + +# ╔═╡ 6f6c8d28-5a54-440e-9b7a-52e1fce959ab +md"# Activate example environment" + +# ╔═╡ f5e2024b-deaf-4344-b610-a4b956abbfaa +md"# Load CUDA +In case you have no CUDA card available, it will run on CPU :) +" + +# ╔═╡ ef92457a-87c0-43bf-a046-9fe82afbbe13 +begin + use_CUDA = Ref(true && CUDA.functional()) + var"@mytime" = use_CUDA[] ? CUDA.var"@time" : var"@time" + togoc(x) = use_CUDA[] ? CuArray(x) : x +end + +# ╔═╡ 810aebe4-2c6e-4ba6-916b-9e4306df33c9 +TableOfContents() + +# ╔═╡ 23288c14-cd27-4bca-b28a-8fb02bd58bd8 + + +# ╔═╡ f80cbda7-1efd-4323-913b-fb7cdcf7fcff +# ╠═╡ disabled = true +#=╠═╡ +begin + img2 = zeros((10,10,1)) + img2[6,5,1] = 3 + img2[6,6,1] = 1 +end + ╠═╡ =# + +# ╔═╡ 59c69bee-e15b-4e43-b824-7c588c3aabf7 +begin + img2 = zeros((8,8,1)) + img2[5,4,1] = 3 + img2[5,5,1] = 1 +end + +# ╔═╡ b8c43ce0-9a08-489c-a8bc-0f9d18473369 +img2 + +# ╔═╡ 51d38dd3-fc3b-4c14-9b80-e6b39051ac51 +simshow(img2[:, :, 1]) + +# ╔═╡ c508389f-a605-4457-9360-410a8021b46c +Revise.errors() + +# ╔═╡ c4a249b8-ac6d-4389-91c0-7df12954871f +@time RadonKA.radon2(randn(Float32, (100, 100)), range(0, 2*π, 100)) + +# ╔═╡ 29a1afa0-7225-4851-98b1-17b06606ba3a +Revise.errors() + +# ╔═╡ aab7ddad-8c54-4c6c-9096-0f043dc2716a +x = abs.(ones((100,100))) + +# ╔═╡ e575b9c3-5410-4e9f-9ccf-7809edae4721 +a = [0] + +# ╔═╡ 80532842-b487-4cc1-9798-f0aa3762034a +simshow(x) + +# ╔═╡ bbbc01c6-c70d-42bf-969b-7b4fc4bbdfd6 +g = RadonParallelCircle(size(x,1), -(size(x,1)-1)÷2:(size(x,1)-1)÷2) + +# ╔═╡ 46956618-d12f-4096-a91d-e4fdfa32a7c7 +simshow(iradon2(radon2(x, a; geometry=g), a; geometry=g, μ=1/30)) + +# ╔═╡ 57fca1d7-cf47-4a59-afc8-bf27497d27e4 +Revise.errors() + +# ╔═╡ 29ea08cf-f148-4ad6-aad1-f11d51cdfe82 + + +# ╔═╡ ae0912af-af3a-4049-ba7e-6a4d6d61af35 +Revise.errors() + +# ╔═╡ c5a80f13-736e-4650-8d70-0c0d28914a37 + + +# ╔═╡ badba484-fedc-4251-8095-01860c826ee2 +Revise.errors() + +# ╔═╡ ea813d04-865e-4f3e-af1c-b50c3f40da54 +Revise.errors() + +# ╔═╡ 28ee6c01-a5dd-47f0-a7b0-177523183c3a +RadonKA.next_cell_intersection + +# ╔═╡ 6ca44a82-2104-4959-b2be-d1d222685353 +# ╠═╡ disabled = true +#=╠═╡ +simshow(radon(img,range(0f0, 2f0*π, 100))[:, :, 1]) + ╠═╡ =# + +# ╔═╡ 6f639b7e-4d50-4696-b566-37b28072a168 +# ╠═╡ disabled = true +#=╠═╡ +@btime sinogram2 = RadonKA.radon2($img2, range(0, 2π, 20)) + ╠═╡ =# + +# ╔═╡ 71789bcf-9c42-4b78-897b-fa254cdc98e4 +# ╠═╡ disabled = true +#=╠═╡ +@btime sinogram2 = RadonKA.radon($img2, range(0, 2π, 20)) + ╠═╡ =# + +# ╔═╡ df04736c-17fe-4fd4-b53e-46ac5afbc6b4 +# ╠═╡ disabled = true +#=╠═╡ +@time radon(img2, range(0, 2π, 20)); + ╠═╡ =# + +# ╔═╡ ceda298a-43b3-4012-a9cc-0d5276678532 +# ╠═╡ disabled = true +#=╠═╡ +simshow(radon(img2, range(0, 2π, 50))[:, :, 1]) + ╠═╡ =# + +# ╔═╡ 190129d1-984c-4582-9d44-880f40b53fcf +Revise.errors() + +# ╔═╡ d0e59091-7a0e-4828-8c17-7c2e643ef4d5 +1 ≤ 3 ≤ 5 + +# ╔═╡ 7d45d80d-d7f4-4db8-9326-1d036ee04881 +RadonKA.RadonCircle(nothing,[1,2]) + +# ╔═╡ d25c1381-baf1-429b-8150-622b8f731d83 +md"# Example Image" + +# ╔═╡ ecd4662c-639a-49fb-b947-b759a733ba22 + + +# ╔═╡ d8d230e5-f870-46cf-a254-47248e3953cd + + +# ╔═╡ bd921606-320b-43d8-8355-e7014037d085 + + +# ╔═╡ 6932f9ea-65fc-408a-b604-65213a17a776 + + +# ╔═╡ 80535445-0bc8-4afe-a497-c9f764428749 + + +# ╔═╡ 04ccf689-b0f6-4fc3-bcb0-9dfc9697d0ec + + +# ╔═╡ 2ccd4d32-269a-43e3-9967-5bf173b18e45 + + +# ╔═╡ 7e40f887-ae05-40f1-b417-5536c8f97b09 + + +# ╔═╡ bfbba881-ee79-4ba1-8aa8-9c94955e7d87 + + +# ╔═╡ c0b3a91c-10f0-45cc-a227-107a7b5d05f2 + + +# ╔═╡ ed576ff5-1b51-4ff2-b3bf-16cb62553f86 + + +# ╔═╡ 54208d78-cf55-41d7-b4bf-6d1ab4927bbb +begin + N = 512 + N_z = 5 + img = box(Float32, (N, N, N_z), (N ÷4, N ÷ 4, 20), offset=(N ÷ 2 + 60, N ÷ 2 -50, N_z ÷ 2)) |> collect + + #img = box(Float32, (N, N, N_z), (1, 1, 1)) |> collect + #img = box(Float32, (N, N, N_z), (N ÷2, N ÷ 2, 50), offset=(N ÷ 2 - 50, N ÷ 2 + 50, N_z ÷ 2)) |> collect + + img .+= 0.0f0 .+ (rr2(Float32, (N, N, N_z), offset=(180, 210, N_z÷2)) .< 30 .^2) + img[:, :, 1] .= testimage("resolution_test_512") + #img = box(Float32, (100, 100), (3, 3), offset=(51, 51)) |> collect +end; + +# ╔═╡ dccc67e0-411f-40d8-b28b-33a1701230bc +@benchmark RadonKA.radon($img, Float32.(range(0, 2*π, 100))) + +# ╔═╡ 72074d52-9ce0-4a1c-9b02-cbc593fadf66 +sinogram = @time RadonKA.radon2(img, Float32.(range(0, 2*π, 100))); + +# ╔═╡ b5f14855-5f26-490e-92cc-e5345d024a70 +simshow(iradon(sinogram, range(0, 2*π, 100)))[:, :, 1] + +# ╔═╡ 7185a544-8a84-426b-bb3f-7b227ca865ff +simshow(iradon2(sinogram, range(0, 2*π, 100)))[:, :, 1] + +# ╔═╡ d12b4d0f-dea0-448e-9982-6d1d8c410bf8 +@benchmark RadonKA.radon2($img, Float32.(range(0, 2*π, 100))) + +# ╔═╡ ce5aa941-b7bb-41f9-96b4-b6e1fe052bc2 +@benchmark RadonKA.radon($img, Float32.(range(0, 2*π, 100))) + +# ╔═╡ 7ad93d26-6041-46ed-9fb9-d116ea2e5608 +@benchmark RadonKA.radon2($img, Float32.(range(0, 2*π, 100))) + +# ╔═╡ 14074b7b-c937-4dfc-896c-fbe9d4f226de +@benchmark RadonKA.radon2($img, Float32.(range(0, 2*π, 100)), μ=0.1f0) + +# ╔═╡ b2c6a7fb-0d8c-45c3-8de4-0b6ae09450cb +@benchmark RadonKA.radon2($img, Float32.(range(0, 2*π, 100))) + +# ╔═╡ 117f1b48-c96b-48d3-8940-d6639e29ef76 +@time sinogram3 = RadonKA.radon(img, range(0, 2*π, 1000)) + +# ╔═╡ f7dbad0d-ba37-473c-82e5-0fa441fc66a6 +simshow(sinogram3[:, :, 1]) + +# ╔═╡ 9b3ac362-3194-4b56-ad41-8848eb17b2e2 +@time sinogram2 = RadonKA.radon2(img, Float32.(range(0, 2*π, 1000)), + geometry=RadonParallelCircle(-250:1:250)); + +# ╔═╡ 9315be14-1822-4091-a1ed-f570bc379414 +simshow(sinogram2[:, :, 1]) + +# ╔═╡ 969fa003-3050-46d3-92b6-6df44e7f903f +simshow(sinogram2[:, :, 1]) + +# ╔═╡ c946a431-a3dc-4dd5-92e1-a5d24bd0c55e +sinogram2 + +# ╔═╡ 1393d029-66be-40aa-a2f9-f31317222575 +img_c = togoc(img); + +# ╔═╡ aecb76ad-d6a2-4d73-b8d5-8327e2448d5b +img_c + +# ╔═╡ b87ee00b-d071-41be-a2c1-f7c32b6105e6 +@benchmark CUDA.@sync RadonKA.radon($img_c, CuArray(range(0, 2*π, 100))) + +# ╔═╡ 4e6924ba-c877-4e31-9e92-425dfbd0c3c6 +sinogramc = RadonKA.radon2(img_c, CuArray(range(0f0, 2f0*π, 100))) + +# ╔═╡ c0cd7d00-ae6e-4ba8-aad6-5b1c9c63736e +@benchmark CUDA.@sync RadonKA.iradon2($sinogramc, CuArray(range(0f0, 2f0*π, 100))) + +# ╔═╡ 0cd70b79-b526-46c5-9f34-8056d1a3478a +@benchmark CUDA.@sync RadonKA.radon2($img_c, CuArray(range(0f0, 2f0*π, 100))) + +# ╔═╡ 08e670a7-4d64-4a01-acc8-35a2f4fc7401 +@benchmark CUDA.@sync RadonKA.radon2($img_c, range(0f0, 2f0*π, 100)) + +# ╔═╡ 1414d566-24d5-4d79-8806-9f3d42f9bbff +CUDA.@time RadonKA.radon2(img_c, collect(range(0f0, 2f0*π, 100))); + +# ╔═╡ f9ad870a-204f-4e0d-8a60-2ddfa42d1506 +CUDA.@time RadonKA.radon2(img_c, CuArray(range(0f0, 2f0*π, 100))); + +# ╔═╡ 43da57cb-f777-4d0b-9754-1e863743f72f +CUDA.@time RadonKA.radon2(img_c, CuArray(range(0f0, 2f0*π, 100)), μ=0.10f0); + +# ╔═╡ f994ed14-67a7-4ee1-878f-61927882f136 +CUDA.@time RadonKA.radon2(img_c, CuArray(range(0f0, 2f0*π, 100))); + +# ╔═╡ 1611bd14-6b59-4e37-8bfb-09c737cd6035 +CUDA.@time RadonKA.radon(img_c, CuArray(range(0f0, 2f0*π, 100))); + +# ╔═╡ 1e027ca9-5050-4b74-af5a-b92f28ee7205 +@benchmark RadonKA.radon2($img_c, CuArray(range(0f0, 2f0*π, 100)), +geometry=RadonParallelCircle(-250.0:1.0:250.0)) + +# ╔═╡ 18f1fd94-bc04-4602-a7f4-b88e0dab2904 +@benchmark RadonKA.radon2($img_c, CuArray(range(0f0, 2f0*π, 100)), +geometry=RadonParallelCircle(-250.0f0:1.0f0:250f0)) + +# ╔═╡ 419f16ba-f4fe-4421-91fe-412d92d51df6 +@code_warntype RadonKA.radon2(img_c, CuArray(range(0f0, 2f0*π, 100))); + +# ╔═╡ 3ce983fc-fd0a-48bf-b2ed-a39457f877d0 +size(img_c) + +# ╔═╡ 01b4b8f8-37d5-425f-975e-ebb3890d8624 +simshow(img[:, :, 1]) + +# ╔═╡ Cell order: +# ╠═6f6c8d28-5a54-440e-9b7a-52e1fce959ab +# ╠═4eb3148e-8f8b-11ee-3cfe-854d3bd5cc80 +# ╠═b336e55e-0be4-422f-b48a-0a2242cb0915 +# ╠═1311e853-c4cd-42bb-8bf3-5e0d564bf9c5 +# ╠═e90a05dd-5781-44ef-9b39-5b1ee85ca477 +# ╟─f5e2024b-deaf-4344-b610-a4b956abbfaa +# ╠═03bccb92-b47f-477a-9bdb-74cc404da690 +# ╠═ef92457a-87c0-43bf-a046-9fe82afbbe13 +# ╟─810aebe4-2c6e-4ba6-916b-9e4306df33c9 +# ╠═23288c14-cd27-4bca-b28a-8fb02bd58bd8 +# ╠═f80cbda7-1efd-4323-913b-fb7cdcf7fcff +# ╠═59c69bee-e15b-4e43-b824-7c588c3aabf7 +# ╠═b8c43ce0-9a08-489c-a8bc-0f9d18473369 +# ╠═51d38dd3-fc3b-4c14-9b80-e6b39051ac51 +# ╠═c508389f-a605-4457-9360-410a8021b46c +# ╠═aecb76ad-d6a2-4d73-b8d5-8327e2448d5b +# ╠═dccc67e0-411f-40d8-b28b-33a1701230bc +# ╠═c4a249b8-ac6d-4389-91c0-7df12954871f +# ╠═29a1afa0-7225-4851-98b1-17b06606ba3a +# ╠═aab7ddad-8c54-4c6c-9096-0f043dc2716a +# ╠═e575b9c3-5410-4e9f-9ccf-7809edae4721 +# ╠═80532842-b487-4cc1-9798-f0aa3762034a +# ╠═bbbc01c6-c70d-42bf-969b-7b4fc4bbdfd6 +# ╠═46956618-d12f-4096-a91d-e4fdfa32a7c7 +# ╠═b5f14855-5f26-490e-92cc-e5345d024a70 +# ╠═7185a544-8a84-426b-bb3f-7b227ca865ff +# ╠═72074d52-9ce0-4a1c-9b02-cbc593fadf66 +# ╠═57fca1d7-cf47-4a59-afc8-bf27497d27e4 +# ╠═d12b4d0f-dea0-448e-9982-6d1d8c410bf8 +# ╠═ce5aa941-b7bb-41f9-96b4-b6e1fe052bc2 +# ╠═b87ee00b-d071-41be-a2c1-f7c32b6105e6 +# ╠═29ea08cf-f148-4ad6-aad1-f11d51cdfe82 +# ╠═4e6924ba-c877-4e31-9e92-425dfbd0c3c6 +# ╠═c0cd7d00-ae6e-4ba8-aad6-5b1c9c63736e +# ╠═0cd70b79-b526-46c5-9f34-8056d1a3478a +# ╠═08e670a7-4d64-4a01-acc8-35a2f4fc7401 +# ╠═1414d566-24d5-4d79-8806-9f3d42f9bbff +# ╠═f9ad870a-204f-4e0d-8a60-2ddfa42d1506 +# ╠═43da57cb-f777-4d0b-9754-1e863743f72f +# ╠═f994ed14-67a7-4ee1-878f-61927882f136 +# ╠═1611bd14-6b59-4e37-8bfb-09c737cd6035 +# ╠═1e027ca9-5050-4b74-af5a-b92f28ee7205 +# ╠═18f1fd94-bc04-4602-a7f4-b88e0dab2904 +# ╠═419f16ba-f4fe-4421-91fe-412d92d51df6 +# ╠═ae0912af-af3a-4049-ba7e-6a4d6d61af35 +# ╠═3ce983fc-fd0a-48bf-b2ed-a39457f877d0 +# ╠═7ad93d26-6041-46ed-9fb9-d116ea2e5608 +# ╠═14074b7b-c937-4dfc-896c-fbe9d4f226de +# ╠═c5a80f13-736e-4650-8d70-0c0d28914a37 +# ╠═badba484-fedc-4251-8095-01860c826ee2 +# ╠═ea813d04-865e-4f3e-af1c-b50c3f40da54 +# ╠═b2c6a7fb-0d8c-45c3-8de4-0b6ae09450cb +# ╠═9315be14-1822-4091-a1ed-f570bc379414 +# ╠═117f1b48-c96b-48d3-8940-d6639e29ef76 +# ╠═9b3ac362-3194-4b56-ad41-8848eb17b2e2 +# ╠═28ee6c01-a5dd-47f0-a7b0-177523183c3a +# ╠═f7dbad0d-ba37-473c-82e5-0fa441fc66a6 +# ╠═969fa003-3050-46d3-92b6-6df44e7f903f +# ╠═6ca44a82-2104-4959-b2be-d1d222685353 +# ╠═6f639b7e-4d50-4696-b566-37b28072a168 +# ╠═71789bcf-9c42-4b78-897b-fa254cdc98e4 +# ╠═df04736c-17fe-4fd4-b53e-46ac5afbc6b4 +# ╠═ceda298a-43b3-4012-a9cc-0d5276678532 +# ╠═c946a431-a3dc-4dd5-92e1-a5d24bd0c55e +# ╠═190129d1-984c-4582-9d44-880f40b53fcf +# ╠═d0e59091-7a0e-4828-8c17-7c2e643ef4d5 +# ╠═7d45d80d-d7f4-4db8-9326-1d036ee04881 +# ╠═d25c1381-baf1-429b-8150-622b8f731d83 +# ╠═ecd4662c-639a-49fb-b947-b759a733ba22 +# ╠═d8d230e5-f870-46cf-a254-47248e3953cd +# ╠═bd921606-320b-43d8-8355-e7014037d085 +# ╠═6932f9ea-65fc-408a-b604-65213a17a776 +# ╠═80535445-0bc8-4afe-a497-c9f764428749 +# ╠═04ccf689-b0f6-4fc3-bcb0-9dfc9697d0ec +# ╠═2ccd4d32-269a-43e3-9967-5bf173b18e45 +# ╠═7e40f887-ae05-40f1-b417-5536c8f97b09 +# ╠═bfbba881-ee79-4ba1-8aa8-9c94955e7d87 +# ╠═c0b3a91c-10f0-45cc-a227-107a7b5d05f2 +# ╠═ed576ff5-1b51-4ff2-b3bf-16cb62553f86 +# ╠═54208d78-cf55-41d7-b4bf-6d1ab4927bbb +# ╠═1393d029-66be-40aa-a2f9-f31317222575 +# ╠═01b4b8f8-37d5-425f-975e-ebb3890d8624 diff --git a/examples/volumetric_additive_manufacturing_optimization.jl b/examples/volumetric_additive_manufacturing_optimization.jl deleted file mode 100644 index 68baeee..0000000 --- a/examples/volumetric_additive_manufacturing_optimization.jl +++ /dev/null @@ -1,531 +0,0 @@ -### A Pluto.jl notebook ### -# v0.19.36 - -using Markdown -using InteractiveUtils - -# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error). -macro bind(def, element) - quote - local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end - local el = $(esc(element)) - global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el) - el - end -end - -# ╔═╡ 453f6cca-91ea-11ee-0c48-b923a785e914 -begin - using Pkg - Pkg.activate(".") - using Revise -end - -# ╔═╡ d956ede0-5a1f-4ed6-a819-1730f52a3536 -using TestImages, RadonKA, ImageShow, ImageIO, Noise, PlutoUI, IndexFunArrays, FFTW - -# ╔═╡ 9fe3b45a-3901-40e4-bad9-32ee3735d971 -using Plots, FileIO - -# ╔═╡ 57d629e6-4de8-4131-bd97-a505c31ab475 -using NDTools, Statistics - -# ╔═╡ 20f68941-4141-4d81-bae8-5ee9661dd80a -using CUDA, CUDA.CUDAKernels, KernelAbstractions - -# ╔═╡ 310da375-4964-4165-aa6b-c50c151ea65e -begin - plot_font = "Computer Modern" - default(fontfamily=plot_font, - linewidth=2, framestyle=:box, label=nothing, grid=false) - scalefontsizes(1.3) -end - -# ╔═╡ a264aa76-2057-42d4-8dd6-254c3d0cdf44 -# use CUDA if functional -use_CUDA = Ref(true && CUDA.functional()) - -# ╔═╡ 4678b348-f61e-4720-af21-7313c79365f7 -var"@mytime" = use_CUDA[] ? CUDA.var"@time" : var"@time" - -# ╔═╡ fbd565e1-ea6c-441a-8ed4-39e65071b798 -togoc(x) = use_CUDA[] ? CuArray(x) : x - -# ╔═╡ 21f890e4-2bb7-4b06-bd79-5e8bf5cce6f5 -md"# A sample we want to print" - -# ╔═╡ 358dda4a-d589-4993-9d6f-d42bcddda1a3 -# ╠═╡ disabled = true -#=╠═╡ -begin - img = zeros(Float32, (200,200)) - - #img .+= rr2(img, offset=(70, 100)) .< 30^2 - #img .-= rr2(img, offset=(70, 100)) .< 20^2 - img .+= box(img, (10, 70), offset=(130, 80)) - #img .+= box(img, (1, 1), offset=(130, 100)) - #img .+= box(img, (1, 1), offset=(110, 100)) - img .+= rr2(img, offset=(70, 100)) .< 30^2 -end; - ╠═╡ =# - -# ╔═╡ b6ed1974-d562-47d4-a342-66eb474bc556 -# ╠═╡ disabled = true -#=╠═╡ -function make_grid_hole() - sz = (512, 512) - - img = zeros(Bool, sz) - img = img .|| Bool.(box(Int, sz, (600, 70), offset=(251, 100))) - img = img .|| Bool.(box(Int, sz, (70, 600), offset=(100, 251))) - - img = img .|| Bool.(box(Int, sz, (30, 600), offset=(170, 251))) - img = img .|| Bool.(box(Int, sz, (600, 30), offset=(251, 170))) - - - img = img .|| Bool.(box(Int, sz, (30, 600), offset=(420, 251))) - img = img .|| Bool.(box(Int, sz, (600, 30), offset=(251, 420))) - - img = img .|| Bool.(box(Int, sz, (15, 600), offset=(200, 251))) - img = img .|| Bool.(box(Int, sz, (600, 15), offset=(251, 200))) - - - - img = img .|| Bool.(box(Int, sz, (5, 600), offset=(240, 251))) - img = img .|| Bool.(box(Int, sz, (600, 5), offset=(251, 240))) - - - img = img .|| Bool.(rr2(Int, sz, offset=(300, 300)) .<= 80 .^2) - img = img .- Bool.(rr2(Int, sz, offset=(300, 300)) .<= 30 .^2) - - # img = img .|| Bool.(box(Int, sz, (30, 50), offset=(50, 50))) - - return Float32.(img) .* (rr2(Int,sz) .<= 250 ^2) -end - ╠═╡ =# - -# ╔═╡ 79c1b042-2b9d-43db-9525-d66418db5da2 -img = Float32.((Float32.(Gray.(select_region(load("/home/felix/Documents/code/Printing_Siemens_Star_Target/siemens_star_512px.png"), new_size=(512, 512)))) .* (rr2(Float32, (512, 512)) .<= 255 .^2)) .> 0.5); - -# ╔═╡ aa1cccec-f846-4b1e-9aae-81e939cd3b8c -md"""# Algorithm - -See this paper: Rackson, Charles M., et al. "Object-space optimization of tomographic reconstructions for additive manufacturing." Additive Manufacturing 48 (2021): 102367. -""" - -# ╔═╡ 0c1abf70-8004-460a-b008-750050f6e718 -function iter!(buffer, img, θs, μ; filtered=true, clip_sinogram=true) - sinogram = radon(img, θs, μ) - if filtered - sinogram = ramp_filter(sinogram) - end - - if clip_sinogram - sinogram .= max.(sinogram, 0) - end - - img_recon = iradon(sinogram, θs, μ) - img_recon ./= maximum(img_recon) - buffer .= max.(img_recon, 0) - return buffer, sinogram -end - -# ╔═╡ 5ac76842-f5a3-42c0-9764-9b20f98ab1d5 -function optimize(img::AbstractArray{T}, θs, μ=nothing; thresholds=(0.65, 0.75), N_iter = 2) where T - N = size(img, 1) - fx = (-N / 2):1:(N /2 -1) - R2D = similar(img) - R2D .= sqrt.(fx'.^2 .+ fx.^2) - - p = plan_fft(img, (1,2)) - guess = max.(0, real.(inv(p) * ((p * img) .* ifftshift(R2D, (1,2))))) - guess ./= maximum(guess) - - loss(x) = (sum(max.(0,thresholds[2] .- x[isobject])) + sum(max.(0, x[notobject] .- thresholds[1]))) / length(x) - #guess = copy(img) - notobject = iszero.(img) - isobject = isone.(img) - - losses = T[] - buffer = copy(img) - tmp, s = iter!(buffer, guess, θs, μ; filtered=false, clip_sinogram=true) - for i in 1:N_iter - guess[notobject] .-= max.(0, tmp[notobject] .- thresholds[1]) - - tmp, s = iter!(buffer, guess, θs, μ; filtered=false, clip_sinogram=true) - - guess[isobject] .+= max.(0, thresholds[2] .- tmp[isobject]) - - push!(losses, loss(tmp)) - end - - printed = iradon(s, θs, μ) - printed ./= maximum(printed) - return guess, s, printed, losses -end - -# ╔═╡ e19a1fb6-d237-44c8-a1b1-45e38aebf948 -function errors(printed, isobject, notobject, thresholds) - mid_thresh = (thresholds[2] + thresholds[1]) / 2 - W_not = sum(printed[notobject] .> thresholds[1]) - W_not_is = sum(printed[notobject] .> mid_thresh) - W_is = sum(printed[isobject] .< thresholds[2]) - - N_not = sum(notobject) - N_is = sum(isobject) - - voxels_object_wrong_printed = sum(abs.((printed .> mid_thresh)[isobject] .- img[isobject])) - voxels_void_wrong_printed = sum(abs.((printed .> mid_thresh)[notobject] .- img[notobject])) - - #voxels_object_wrong_printed / N_is, W_not_is / N_not, W_not / N_not, W_is / N_is - - @info "Object pixels not printed $(round(voxels_object_wrong_printed / N_is * 100, digits=4))%" - @info "Void pixels falsely printed $(round(voxels_void_wrong_printed / N_not * 100, digits=4))%" -end - -# ╔═╡ dceec7d4-706c-4a3b-8cf4-fa73f9800e3d - - -# ╔═╡ 5290976e-f60e-4a02-b58d-f91a870a5f89 -function plot_histogram(img, object_printed, thresholds, chosen_threshold) - plot(object_printed[img .== 0], bins=(0.0:0.01:1), seriestype=:barhist, xlim=(0.0, 1.0), label="dose distribution void", ylabel="voxel count", xlabel="normalized intensity", ylim=(10, 1000000), yscale=:log10, linewidth=1, legend=:topleft) - plot!(object_printed[img .== 1], seriestype=:barhist, bins=(0.0:0.01:1), xlim=(0.0, 1.0), label="dose distribution object", ylabel="voxel count", xlabel="normalized intensity", ylim=(10, 1000000), linewidth=1, yscale=:log10,) - plot!([thresholds[1], thresholds[1]], [1, 10000_000], label="lower threshold", linewidth=3) - plot!([thresholds[2], thresholds[2]], [1, 10000_000], label="upper threshold", linewidth=3) - #plot!([chosen_threshold, chosen_threshold], [1, 30000000], label="chosen threshold", linewidth=3) -end - -# ╔═╡ 7d694ff4-d0f6-4f10-a26b-d3c43930dd4f -md"# Optimize" - -# ╔═╡ 578f360d-66b2-434e-bdf2-469673116711 -angles = range(0, 2π, 805) - -# ╔═╡ 448f7462-6a87-4ca3-a2a3-0f76934e97a1 -N_iter = 500 - -# ╔═╡ aff2fc67-5569-40e6-9ed3-936f16de1aba -thresholds = (0.65, 0.75) - -# ╔═╡ 96f6eba3-1700-4c29-aff6-5f774d7d4b3b -@time _, patterns_065_075, printed_065_075, losses_065_075 = Array.(optimize(togoc(img), togoc(angles); N_iter=20, thresholds=(0.65, 0.75))) - -# ╔═╡ 198375e7-6857-4971-92c8-6aec7cb25090 -size(printed_065_075) - -# ╔═╡ 4b04002a-449c-46a6-921e-b3fe72395018 -size(patterns_065_075) - -# ╔═╡ e55b5e16-1ea1-4cfa-a0df-c32e03848c9f -errors(printed_065_075, isone.(img), iszero.(img), thresholds) - -# ╔═╡ 1aa687fd-9647-4bce-a800-5cfe6730b28b -@bind thresh1 Slider(0.0:0.005:1, show_value=true) - -# ╔═╡ 0ac01881-af91-4067-a216-6cef66d43bc0 -md"Threshold =$thresh1" - -# ╔═╡ f6b67921-cf3f-4144-bbdb-3c423d6ee81d -sum(printed_065_075 .> 0.5) - -# ╔═╡ 0c7375cb-2cd1-4ebd-9227-5b561783e37e -simshow([printed_065_075 printed_065_075 .> thresh1 abs.((printed_065_075 .> thresh1) .- img)]) - -# ╔═╡ c81c9248-29b5-494a-95e1-276a1f8ecfa2 -begin - save("/tmp/printed_065_075.png", simshow(printed_065_075)) - save("/tmp/printed_065_075_thresh.png", simshow(printed_065_075 .> 0.7)) - save("/tmp/printed_065_075_thresh_diff.png", simshow(abs.((printed_065_075 .> 0.7) .- img))) -end - -# ╔═╡ bf5a5074-ad69-4d17-b8de-892688a76891 -begin - p = plot_histogram(img[1:1:end], printed_065_075[1:1:end], thresholds, 0.7) - savefig(p, "/tmp/histogram_printed_065_075.pdf") - p -end - -# ╔═╡ ec51fa20-725c-46a4-b5ab-eaefdbdda501 - - -# ╔═╡ ef196630-d7fd-4481-b4e9-ce410dec8e73 -@time _, patterns_05_08, printed_05_08, losses_05_08 = Array.(optimize(togoc(img), togoc(angles); N_iter=N_iter, thresholds=(0.5, 0.8))); - -# ╔═╡ 5e1f0420-45d6-48b5-a249-8e5b7048bd65 -errors(printed_05_08, isone.(img), iszero.(img), thresholds) - -# ╔═╡ 118cbae4-04af-461a-9301-dafc9b1efbf7 -begin - save("/tmp/printed_05_08.png", simshow(printed_05_08)) - save("/tmp/printed_05_08_thresh_2.png", simshow(printed_05_08 .> 0.65)) - save("/tmp/printed_05_08_thresh_diff_2.png", simshow(abs.((printed_05_08 .> 0.65) .- img))) -end - -# ╔═╡ 6dd162d4-f099-49f7-860b-dc1e0088090c -begin - p2 = plot_histogram(img, printed_05_08, (0.5, 0.8), 0.65) - savefig(p2, "/tmp/histogram_printed_05_08.pdf") - p2 -end - -# ╔═╡ 0ca23230-da00-43fc-8ad7-c93a0f56e123 - - -# ╔═╡ f2bfcd8f-2e65-4fd4-903b-667171d685ff -@time _, patterns_065_075_3mu, printed_065_075_3mu, losses_065_075_3mu = Array.(optimize(togoc(img), togoc(angles), 3/ size(img, 1); N_iter=N_iter, thresholds=(0.65, 0.75))) - -# ╔═╡ d77c1ac9-9e08-47ca-89db-3680e550eb4e -simshow(printed_065_075_3mu) - -# ╔═╡ d53fcc24-5fee-4e4f-a144-f9fa2b0ee5cb -errors(printed_065_075_3mu, isone.(img), iszero.(img), thresholds) - -# ╔═╡ 9351709f-7831-433b-839f-84c5b8b315c1 -begin - save("/tmp/printed_065_075_3mu.png", simshow(printed_065_075_3mu)) - save("/tmp/printed_065_075_3mu_thresh.png", simshow(printed_065_075_3mu .> 0.7)) - save("/tmp/printed_065_075_3mu_thresh_diff.png", simshow(abs.((printed_065_075_3mu .> 0.7) .- img))) -end - -# ╔═╡ 5670277d-65b8-40f1-95a2-78358f082bed -begin - p3 = plot_histogram(img, printed_065_075_3mu, thresholds, 0.7) - savefig(p3, "/tmp/histogram_printed_065_075_3mu.pdf") - p3 -end - -# ╔═╡ ac932771-c8e1-4ecf-ab09-008350a39b28 - - -# ╔═╡ 19def696-ef29-466b-9ace-a7ed8f9d8c28 -@time _, patterns_065_075_1mu, printed_065_075_1mu, losses_065_075_1mu = Array.(optimize(togoc(img), togoc(angles), 1/ size(img, 1); N_iter=N_iter, thresholds=(0.65, 0.75))) - -# ╔═╡ 51950c3e-bfb0-443d-89d3-fe52e0b720ac -errors(printed_065_075_1mu, isone.(img), iszero.(img), thresholds) - -# ╔═╡ f28558d3-bab7-4c0d-9a48-e470df8d6b2a -begin - save("/tmp/printed_065_075_1mu.png", simshow(printed_065_075_1mu)) - save("/tmp/printed_065_075_1mu_thresh.png", simshow(printed_065_075_1mu .> 0.7)) - save("/tmp/printed_065_075_1mu_thresh_diff.png", simshow(abs.((printed_065_075_1mu .> 0.7) .- img))) -end - -# ╔═╡ fae06fab-a707-487c-aeea-6f0c8e5ca3f6 -begin - p4 = plot_histogram(img, printed_065_075_1mu, thresholds, 0.7) - savefig(p4, "/tmp/histogram_printed_065_075_1mu.pdf") - p4 -end - -# ╔═╡ 48aafe12-1e21-4960-b0af-37468ad613f8 - - -# ╔═╡ b1309600-9d64-47cb-aaed-86a1fa3a591e -# ╠═╡ disabled = true -#=╠═╡ -@time _, patterns_065_075_1mu_long, printed_065_075_1mu_long, losses_065_075_1mu_long = Array.(optimize(togoc(img), togoc(angles), 1/ size(img, 1); N_iter=10_000, thresholds=(0.65, 0.75))) - ╠═╡ =# - -# ╔═╡ e6329988-a919-4569-9288-b8c4620a2d67 -#=╠═╡ -errors(printed_065_075_1mu_long, isone.(img), iszero.(img), thresholds) - ╠═╡ =# - -# ╔═╡ a7cc2f68-5592-498c-8891-328f54e35019 -#=╠═╡ -begin - p5 = plot_histogram(img, printed_065_075_1mu_long, thresholds, 0.7) - savefig(p5, "/tmp/histogram_printed_065_075_1mu_long.pdf") - p5 -end - ╠═╡ =# - -# ╔═╡ 3ef9a22c-ce55-4a90-9a5d-d5a5b835253e -#=╠═╡ -save("/tmp/patterns_065_075_1mu_long.png",simshow(patterns_065_075_1mu_long, cmap=:turbo)) - ╠═╡ =# - -# ╔═╡ 2a9f6404-50ca-473a-89f7-a55f163e5f61 -save("/tmp/patterns_065_075_1mu.png",simshow(patterns_065_075_1mu, cmap=:turbo)) - -# ╔═╡ ba576d5e-3f1f-4f1b-aa1c-ef0f425ca14a - - -# ╔═╡ 5fc73726-c51d-41ef-b258-511b435d7ee8 - - -# ╔═╡ adb5d5bb-a925-433b-ba66-c5b7ec4f4d73 -#=╠═╡ -sum(patterns_065_075_1mu_long ./ maximum(patterns_065_075_1mu_long)/ length(patterns_065_075_1mu_long)) - ╠═╡ =# - -# ╔═╡ ef4f6fa4-ee46-47d1-9587-65125609acb2 -#=╠═╡ -sum(patterns_065_075_1mu ./ maximum(patterns_065_075_1mu) / length(patterns_065_075_1mu_long)) - ╠═╡ =# - -# ╔═╡ 708e26fa-89d1-4821-a9b9-6efce3b27b2b -#=╠═╡ -sum(patterns_065_075_1mu / length(patterns_065_075_1mu_long)) - ╠═╡ =# - -# ╔═╡ fa62003a-2748-45fe-9b37-7c32d215c49b -#=╠═╡ -extrema(patterns_065_075_1mu_long) - ╠═╡ =# - -# ╔═╡ 04bc9aed-dd0d-46bb-a1f8-87a385806127 -maximum(patterns_05_08) / mean(patterns_05_08) - -# ╔═╡ 8dbbfc81-7df2-481b-b758-64ed23134f2d -#=╠═╡ -maximum(patterns_065_075_1mu_long) / mean(patterns_065_075_1mu_long) - ╠═╡ =# - -# ╔═╡ c41d4159-9c31-4ece-98d2-a38007f1049b -#=╠═╡ -simshow(patterns_065_075_1mu_long) - ╠═╡ =# - -# ╔═╡ e99e8aaa-f538-4647-82af-8e85397d43c8 - - -# ╔═╡ 16384727-05a8-48bb-8755-8c47bd28b97c -# ╠═╡ disabled = true -#=╠═╡ -begin - p4 = plot(losses, yscale=:log10, label="\$T_L =0.65\$, \$T_U=0.75\$", ylabel="normalized loss", xlabel="iterations", yticks=[10^2, 10^3, 10^4, 10^5], grid=true) - plot!(losses_2, label="\$T_L =0.5\$, \$T_U=0.8\$") - #plot!(losses_3, label="\$T_L =0.75\$, \$T_U=0.8\$") - savefig(p4, "/tmp/loss.pdf") - p4 -end - ╠═╡ =# - -# ╔═╡ b974c854-fcf7-4cbc-b3c8-1a8ab34f9190 - - -# ╔═╡ fb10aac8-2068-4452-83e1-1624a8375d2c - - -# ╔═╡ 40cd392b-8020-4cbb-ac7c-ea32cddd4410 - - -# ╔═╡ 621f86b9-0856-4aae-8ad7-4f1d557ebf81 - - -# ╔═╡ ec67bddb-5fd1-4d26-8e76-41930d22f09c - - -# ╔═╡ bcab5404-1816-4e5a-b0d2-f46735826a65 - - -# ╔═╡ f5d516a4-dd22-4971-977b-1526db3571c8 -md"# Large optimization" - -# ╔═╡ 51dbb38f-ad62-4664-ba95-2f943a33ce60 -begin - volume = zeros(Float32, (512, 512, 100)) - volume .+= box(size(volume)[1:2], (30, 100), offset=(230, 380)) - volume .-= rr2(size(volume)[1:2], offset=(170, 200)) .< 30^2 - volume .+= rr2(size(volume)[1:2], offset=(170, 200)) .< 100^2 -end; - -# ╔═╡ bd6bca52-a797-4a12-a24d-7989e9d941c8 -simshow(volume[:, :, 1]) - -# ╔═╡ 600dcf9d-bcf3-42ea-bc9b-bd399483f028 -angles2 = range(0, π, 200) - -# ╔═╡ 2fd28802-e7dd-4eec-903a-3b94b87cbc16 -# ╠═╡ disabled = true -#=╠═╡ -CUDA.@sync CUDA.@time virtual_object_c2, patterns_optimized_c2 = optimize(CuArray(volume), CuArray(angles2); N_iter=50, thresholds, backend=CUDABackend()) - ╠═╡ =# - -# ╔═╡ da12770d-4e4e-457b-bbaa-affad72eac33 -object_printed_c2 = Array(iradon(Array(patterns_optimized_c2), angles2)); - -# ╔═╡ 9bd41858-b191-4935-b9cc-24fdf6be8cbe -simshow(object_printed_c2 .> 0.75)[:, :, 1] - -# ╔═╡ Cell order: -# ╠═453f6cca-91ea-11ee-0c48-b923a785e914 -# ╠═d956ede0-5a1f-4ed6-a819-1730f52a3536 -# ╠═9fe3b45a-3901-40e4-bad9-32ee3735d971 -# ╠═57d629e6-4de8-4131-bd97-a505c31ab475 -# ╠═310da375-4964-4165-aa6b-c50c151ea65e -# ╠═20f68941-4141-4d81-bae8-5ee9661dd80a -# ╠═a264aa76-2057-42d4-8dd6-254c3d0cdf44 -# ╠═4678b348-f61e-4720-af21-7313c79365f7 -# ╠═fbd565e1-ea6c-441a-8ed4-39e65071b798 -# ╟─21f890e4-2bb7-4b06-bd79-5e8bf5cce6f5 -# ╠═358dda4a-d589-4993-9d6f-d42bcddda1a3 -# ╠═b6ed1974-d562-47d4-a342-66eb474bc556 -# ╠═79c1b042-2b9d-43db-9525-d66418db5da2 -# ╟─aa1cccec-f846-4b1e-9aae-81e939cd3b8c -# ╠═0c1abf70-8004-460a-b008-750050f6e718 -# ╠═5ac76842-f5a3-42c0-9764-9b20f98ab1d5 -# ╠═e19a1fb6-d237-44c8-a1b1-45e38aebf948 -# ╠═dceec7d4-706c-4a3b-8cf4-fa73f9800e3d -# ╠═5290976e-f60e-4a02-b58d-f91a870a5f89 -# ╟─7d694ff4-d0f6-4f10-a26b-d3c43930dd4f -# ╠═578f360d-66b2-434e-bdf2-469673116711 -# ╠═448f7462-6a87-4ca3-a2a3-0f76934e97a1 -# ╠═aff2fc67-5569-40e6-9ed3-936f16de1aba -# ╠═96f6eba3-1700-4c29-aff6-5f774d7d4b3b -# ╠═198375e7-6857-4971-92c8-6aec7cb25090 -# ╠═4b04002a-449c-46a6-921e-b3fe72395018 -# ╠═e55b5e16-1ea1-4cfa-a0df-c32e03848c9f -# ╟─0ac01881-af91-4067-a216-6cef66d43bc0 -# ╠═1aa687fd-9647-4bce-a800-5cfe6730b28b -# ╠═f6b67921-cf3f-4144-bbdb-3c423d6ee81d -# ╠═0c7375cb-2cd1-4ebd-9227-5b561783e37e -# ╠═c81c9248-29b5-494a-95e1-276a1f8ecfa2 -# ╠═bf5a5074-ad69-4d17-b8de-892688a76891 -# ╠═ec51fa20-725c-46a4-b5ab-eaefdbdda501 -# ╠═ef196630-d7fd-4481-b4e9-ce410dec8e73 -# ╠═5e1f0420-45d6-48b5-a249-8e5b7048bd65 -# ╠═118cbae4-04af-461a-9301-dafc9b1efbf7 -# ╠═6dd162d4-f099-49f7-860b-dc1e0088090c -# ╠═0ca23230-da00-43fc-8ad7-c93a0f56e123 -# ╠═f2bfcd8f-2e65-4fd4-903b-667171d685ff -# ╠═d77c1ac9-9e08-47ca-89db-3680e550eb4e -# ╠═d53fcc24-5fee-4e4f-a144-f9fa2b0ee5cb -# ╠═9351709f-7831-433b-839f-84c5b8b315c1 -# ╠═5670277d-65b8-40f1-95a2-78358f082bed -# ╠═ac932771-c8e1-4ecf-ab09-008350a39b28 -# ╠═19def696-ef29-466b-9ace-a7ed8f9d8c28 -# ╠═51950c3e-bfb0-443d-89d3-fe52e0b720ac -# ╠═f28558d3-bab7-4c0d-9a48-e470df8d6b2a -# ╠═fae06fab-a707-487c-aeea-6f0c8e5ca3f6 -# ╠═48aafe12-1e21-4960-b0af-37468ad613f8 -# ╠═b1309600-9d64-47cb-aaed-86a1fa3a591e -# ╠═e6329988-a919-4569-9288-b8c4620a2d67 -# ╠═a7cc2f68-5592-498c-8891-328f54e35019 -# ╠═3ef9a22c-ce55-4a90-9a5d-d5a5b835253e -# ╠═2a9f6404-50ca-473a-89f7-a55f163e5f61 -# ╠═ba576d5e-3f1f-4f1b-aa1c-ef0f425ca14a -# ╠═5fc73726-c51d-41ef-b258-511b435d7ee8 -# ╠═adb5d5bb-a925-433b-ba66-c5b7ec4f4d73 -# ╠═ef4f6fa4-ee46-47d1-9587-65125609acb2 -# ╠═708e26fa-89d1-4821-a9b9-6efce3b27b2b -# ╠═fa62003a-2748-45fe-9b37-7c32d215c49b -# ╠═04bc9aed-dd0d-46bb-a1f8-87a385806127 -# ╠═8dbbfc81-7df2-481b-b758-64ed23134f2d -# ╠═c41d4159-9c31-4ece-98d2-a38007f1049b -# ╠═e99e8aaa-f538-4647-82af-8e85397d43c8 -# ╠═16384727-05a8-48bb-8755-8c47bd28b97c -# ╠═b974c854-fcf7-4cbc-b3c8-1a8ab34f9190 -# ╠═fb10aac8-2068-4452-83e1-1624a8375d2c -# ╠═40cd392b-8020-4cbb-ac7c-ea32cddd4410 -# ╠═621f86b9-0856-4aae-8ad7-4f1d557ebf81 -# ╠═ec67bddb-5fd1-4d26-8e76-41930d22f09c -# ╠═bcab5404-1816-4e5a-b0d2-f46735826a65 -# ╟─f5d516a4-dd22-4971-977b-1526db3571c8 -# ╠═51dbb38f-ad62-4664-ba95-2f943a33ce60 -# ╠═bd6bca52-a797-4a12-a24d-7989e9d941c8 -# ╠═600dcf9d-bcf3-42ea-bc9b-bd399483f028 -# ╠═2fd28802-e7dd-4eec-903a-3b94b87cbc16 -# ╠═da12770d-4e4e-457b-bbaa-affad72eac33 -# ╠═9bd41858-b191-4935-b9cc-24fdf6be8cbe diff --git a/examples/volumetric_printing.jl b/examples/volumetric_printing.jl deleted file mode 100644 index 4ce0567..0000000 --- a/examples/volumetric_printing.jl +++ /dev/null @@ -1,268 +0,0 @@ -### A Pluto.jl notebook ### -# v0.19.32 - -using Markdown -using InteractiveUtils - -# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error). -macro bind(def, element) - quote - local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end - local el = $(esc(element)) - global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el) - el - end -end - -# ╔═╡ 453f6cca-91ea-11ee-0c48-b923a785e914 -begin - using Pkg - Pkg.activate(".") - using Revise -end - -# ╔═╡ d956ede0-5a1f-4ed6-a819-1730f52a3536 -using TestImages, RadonKA, ImageShow, ImageIO, Noise, PlutoUI, IndexFunArrays, FFTW - -# ╔═╡ 9fe3b45a-3901-40e4-bad9-32ee3735d971 -using Plots - -# ╔═╡ 29703cc9-d5ba-4ff3-bc0c-2ad325c7aa7d -using CUDA, CUDA.CUDAKernels, KernelAbstractions - -# ╔═╡ 21f890e4-2bb7-4b06-bd79-5e8bf5cce6f5 -md"# A sample we want to print" - -# ╔═╡ 358dda4a-d589-4993-9d6f-d42bcddda1a3 -begin - img = zeros(Float32, (200,200)) - - #img .+= rr2(img, offset=(70, 100)) .< 30^2 - #img .-= rr2(img, offset=(70, 100)) .< 20^2 - img .+= box(img, (10, 70), offset=(130, 80)) - #img .+= box(img, (1, 1), offset=(130, 100)) - #img .+= box(img, (1, 1), offset=(110, 100)) - img .+= rr2(img, offset=(70, 100)) .< 30^2 -end; - -# ╔═╡ 57da750b-ec5c-4636-bac3-549e7fadc3c8 -simshow(img) - -# ╔═╡ aa1cccec-f846-4b1e-9aae-81e939cd3b8c -md"""# Algorithm - -See this paper: Rackson, Charles M., et al. "Object-space optimization of tomographic reconstructions for additive manufacturing." Additive Manufacturing 48 (2021): 102367. -""" - -# ╔═╡ 0c1abf70-8004-460a-b008-750050f6e718 -function iter!(buffer, img, θs, μ; filtered=true, clip_sinogram=true, backend=CPU()) - sinogram = radon(img, θs, μ; backend) - if filtered - sinogram = ramp_filter(sinogram) - end - - if clip_sinogram - sinogram .= max.(sinogram, 0) - end - - img_recon = iradon(sinogram, θs, μ; backend) - - buffer .= max.(img_recon, 0) - return buffer, sinogram -end - -# ╔═╡ 5ac76842-f5a3-42c0-9764-9b20f98ab1d5 -function optimize(img, θs, μ=nothing; thresholds=(0.65, 0.75), N_iter = 2, backend=CPU()) - N = size(img, 1) - fx = (-N / 2):1:(N /2 -1) - R2D = similar(img) - R2D .= sqrt.(fx'.^2 .+ fx.^2) - - p = plan_fft(img, (1,2)) - guess = max.(0, real.(inv(p) * ((p * img) .* ifftshift(R2D, (1,2))))) - guess ./= maximum(guess) - - guess = copy(img) - notobject = iszero.(img) - isobject = isone.(img) - - buffer = copy(img) - tmp, s = iter!(buffer, guess, θs, μ; filtered=false, clip_sinogram=true, backend) - for i in 1:N_iter - guess[notobject] .-= max.(0, tmp[notobject] .- thresholds[1]) - - tmp, s = iter!(buffer, guess, θs, μ; backend, filtered=false, clip_sinogram=true) - - guess[isobject] .+= max.(0, thresholds[2] .- tmp[isobject]) - end - - return guess, s#radon(guess, θs, det) -end - -# ╔═╡ 7d694ff4-d0f6-4f10-a26b-d3c43930dd4f -md"# Optimize" - -# ╔═╡ 578f360d-66b2-434e-bdf2-469673116711 -angles = range(0, 2π, 350) - -# ╔═╡ aff2fc67-5569-40e6-9ed3-936f16de1aba -thresholds = (0.7, 0.8) - -# ╔═╡ 96f6eba3-1700-4c29-aff6-5f774d7d4b3b -@time virtual_object, patterns_optimized = optimize(img, angles; N_iter=100, thresholds) - -# ╔═╡ 5568ca5b-2aad-4024-8cf0-d68260a4ac36 -simshow(patterns_optimized) - -# ╔═╡ 26613767-5e05-4349-8bc2-2cfe357749c8 -object_printed = iradon(patterns_optimized, angles); - -# ╔═╡ 2d602a28-c28b-4267-966b-81a7de2e2fbc -simshow(virtual_object) - -# ╔═╡ 0ac01881-af91-4067-a216-6cef66d43bc0 -md"Threshold =$(@bind thresh2 Slider(0.0:0.01:1, show_value=true))" - -# ╔═╡ 0c7375cb-2cd1-4ebd-9227-5b561783e37e -simshow([object_printed object_printed .> thresh2 abs.((object_printed .> thresh2) .- img)]) - -# ╔═╡ b5ae35db-bf7d-4c00-a069-8572af9bc1de -begin - histogram(object_printed[img .== 0], bins=(0.0:0.01:1), xlim=(0.0, 1.0), label="dose distribution void", ylabel="voxel count", xlabel="normalized intensity", ylim=(0, 1000)) - histogram!(object_printed[img .== 1], bins=(0.0:0.01:1), xlim=(0.0, 1.0), label="dose distribution object", ylabel="voxel count", xlabel="normalized intensity", ylim=(0, 1000)) - plot!([thresholds[1], thresholds[1]], [0, 100_000], label="Lower threshold") - plot!([thresholds[2], thresholds[2]], [0, 100_000], label="Upper threshold") - plot!([thresh2, thresh2], [0, 3000], label="Real threshold") -end - -# ╔═╡ b277afc7-bd21-4031-866c-3fe0f0747adb -md"# Speed Up with CUDA" - -# ╔═╡ 4dd89f19-68b2-4e01-9ad9-de8fcbcd00ba -radon(CuArray(img), CuArray(angles), backend=CUDABackend()); - -# ╔═╡ e08da75f-d3a9-4aab-acb0-e4a1d5371fae -CUDA.@sync CUDA.@time virtual_object_c, patterns_optimized_c = optimize(CuArray(img), CuArray(angles); N_iter=50, thresholds, backend=CUDABackend()) - -# ╔═╡ 2515881d-be85-4e5a-87fb-881aaa860464 -object_printed_c = Array(iradon(Array(patterns_optimized_c), angles)); - -# ╔═╡ 178846ee-5678-4332-9cb9-f900090f06fa -Array(object_printed_c) ≈ object_printed - -# ╔═╡ 8be0b6ee-5a54-4d7f-a7bf-e71b0a637439 -simshow([object_printed object_printed_c .> thresh2 abs.((object_printed .> thresh2) .- img)]) - -# ╔═╡ ebc5cce3-180c-412d-8b09-0004115a6e78 - - -# ╔═╡ 3767200c-6a46-4afa-8390-d566a36d4ebd -md"# Exponential Absorption" - -# ╔═╡ 8ca1db87-9a1f-45af-840b-e10458047cfa -thresholds2 = (0.6, 0.7) - -# ╔═╡ 96c079e9-0676-4115-bf5e-f6084293b884 -μ=1/70f0 - -# ╔═╡ ed4495a1-89aa-4415-8604-53a4bd05046c -CUDA.@sync CUDA.@time virtual_object_c_abs, patterns_optimized_c_abs = optimize(CuArray(img), CuArray(angles), μ; N_iter=1000, thresholds=thresholds2, backend=CUDABackend()) - -# ╔═╡ 56400f42-8f18-4bb2-a224-3b2513a6c35b -object_printed_abs = Array(iradon(Array(patterns_optimized_c_abs), angles, μ)); - -# ╔═╡ ed116dd3-5d1c-4a5e-869c-101c44e0c42b -md"Threshold =$(@bind thresh3 Slider(0.0:0.005:1, show_value=true))" - -# ╔═╡ 9b39771c-92f0-43db-9880-dcf26f99e455 -begin - histogram(object_printed_abs[img .== 0], bins=(0.0:0.01:1), xlim=(0.0, 1.0), label="dose distribution", ylabel="voxel count", xlabel="normalized intensity", ylim=(0, 500)) - histogram!(object_printed_abs[img .== 1], bins=(0.0:0.01:1), label="dose distribution", ylim=(0, 500)) - plot!([thresholds2[1], thresholds2[1]], [0, 100_000], label="Lower threshold") - plot!([thresholds2[2], thresholds2[2]], [0, 100_000], label="Upper threshold") - plot!([thresh3, thresh3], [0, 3000], label="Real threshold") -end - -# ╔═╡ bf450d7c-536b-452f-96c6-d07aa78e9f31 -simshow([object_printed_abs object_printed_abs .> thresh3 abs.((object_printed_abs .> thresh3) .- img)]) - -# ╔═╡ 8e843705-f4b5-46b1-80f5-d60753f9a05a - - -# ╔═╡ fb10aac8-2068-4452-83e1-1624a8375d2c - - -# ╔═╡ f5d516a4-dd22-4971-977b-1526db3571c8 -md"# Large optimization" - -# ╔═╡ 51dbb38f-ad62-4664-ba95-2f943a33ce60 -begin - volume = zeros(Float32, (512, 512, 100)) - volume .+= box(size(volume)[1:2], (30, 100), offset=(230, 380)) - volume .-= rr2(size(volume)[1:2], offset=(170, 200)) .< 30^2 - volume .+= rr2(size(volume)[1:2], offset=(170, 200)) .< 100^2 -end; - -# ╔═╡ bd6bca52-a797-4a12-a24d-7989e9d941c8 -simshow(volume[:, :, 1]) - -# ╔═╡ 600dcf9d-bcf3-42ea-bc9b-bd399483f028 -angles2 = range(0, π, 200) - -# ╔═╡ 2fd28802-e7dd-4eec-903a-3b94b87cbc16 -# ╠═╡ disabled = true -#=╠═╡ -CUDA.@sync CUDA.@time virtual_object_c2, patterns_optimized_c2 = optimize(CuArray(volume), CuArray(angles2); N_iter=50, thresholds, backend=CUDABackend()) - ╠═╡ =# - -# ╔═╡ da12770d-4e4e-457b-bbaa-affad72eac33 -object_printed_c2 = Array(iradon(Array(patterns_optimized_c2), angles2)); - -# ╔═╡ 9bd41858-b191-4935-b9cc-24fdf6be8cbe -simshow(object_printed_c2 .> 0.74)[:, :, 1] - -# ╔═╡ Cell order: -# ╠═453f6cca-91ea-11ee-0c48-b923a785e914 -# ╠═d956ede0-5a1f-4ed6-a819-1730f52a3536 -# ╠═9fe3b45a-3901-40e4-bad9-32ee3735d971 -# ╟─21f890e4-2bb7-4b06-bd79-5e8bf5cce6f5 -# ╠═358dda4a-d589-4993-9d6f-d42bcddda1a3 -# ╠═57da750b-ec5c-4636-bac3-549e7fadc3c8 -# ╟─aa1cccec-f846-4b1e-9aae-81e939cd3b8c -# ╠═0c1abf70-8004-460a-b008-750050f6e718 -# ╠═5ac76842-f5a3-42c0-9764-9b20f98ab1d5 -# ╟─7d694ff4-d0f6-4f10-a26b-d3c43930dd4f -# ╠═578f360d-66b2-434e-bdf2-469673116711 -# ╠═aff2fc67-5569-40e6-9ed3-936f16de1aba -# ╠═96f6eba3-1700-4c29-aff6-5f774d7d4b3b -# ╠═5568ca5b-2aad-4024-8cf0-d68260a4ac36 -# ╠═26613767-5e05-4349-8bc2-2cfe357749c8 -# ╠═2d602a28-c28b-4267-966b-81a7de2e2fbc -# ╟─0ac01881-af91-4067-a216-6cef66d43bc0 -# ╠═0c7375cb-2cd1-4ebd-9227-5b561783e37e -# ╠═b5ae35db-bf7d-4c00-a069-8572af9bc1de -# ╟─b277afc7-bd21-4031-866c-3fe0f0747adb -# ╠═29703cc9-d5ba-4ff3-bc0c-2ad325c7aa7d -# ╠═4dd89f19-68b2-4e01-9ad9-de8fcbcd00ba -# ╠═e08da75f-d3a9-4aab-acb0-e4a1d5371fae -# ╠═2515881d-be85-4e5a-87fb-881aaa860464 -# ╠═178846ee-5678-4332-9cb9-f900090f06fa -# ╠═8be0b6ee-5a54-4d7f-a7bf-e71b0a637439 -# ╠═ebc5cce3-180c-412d-8b09-0004115a6e78 -# ╟─3767200c-6a46-4afa-8390-d566a36d4ebd -# ╠═8ca1db87-9a1f-45af-840b-e10458047cfa -# ╠═96c079e9-0676-4115-bf5e-f6084293b884 -# ╠═ed4495a1-89aa-4415-8604-53a4bd05046c -# ╠═56400f42-8f18-4bb2-a224-3b2513a6c35b -# ╠═9b39771c-92f0-43db-9880-dcf26f99e455 -# ╟─ed116dd3-5d1c-4a5e-869c-101c44e0c42b -# ╟─bf450d7c-536b-452f-96c6-d07aa78e9f31 -# ╠═8e843705-f4b5-46b1-80f5-d60753f9a05a -# ╠═fb10aac8-2068-4452-83e1-1624a8375d2c -# ╟─f5d516a4-dd22-4971-977b-1526db3571c8 -# ╠═51dbb38f-ad62-4664-ba95-2f943a33ce60 -# ╠═bd6bca52-a797-4a12-a24d-7989e9d941c8 -# ╠═600dcf9d-bcf3-42ea-bc9b-bd399483f028 -# ╠═2fd28802-e7dd-4eec-903a-3b94b87cbc16 -# ╠═da12770d-4e4e-457b-bbaa-affad72eac33 -# ╠═9bd41858-b191-4935-b9cc-24fdf6be8cbe diff --git a/src/RadonKA.jl b/src/RadonKA.jl index 83bab24..642a664 100644 --- a/src/RadonKA.jl +++ b/src/RadonKA.jl @@ -8,13 +8,15 @@ using ChainRulesCore using PrecompileTools +include("geometry.jl") include("utils.jl") include("radon.jl") include("iradon.jl") +include("filtered_backprojection.jl") +include("radon_legacy.jl") - - +# PrecompileTools @setup_workload begin img = randn((2,2)) angles = range(0, π, 2) @@ -27,10 +29,10 @@ include("iradon.jl") iradon(r, angles) RadonKA.filtered_backprojection(r, angles) - r = radon(Float32.(img), Float32.(angles), 0.1f0) - iradon(r, Float32.(angles), 0.1f0) - r = radon(img, angles, 0.1) - iradon(r, angles, 0.1) + r = radon(Float32.(img), Float32.(angles), μ=0.1f0) + iradon(r, Float32.(angles), μ=0.1f0) + r = radon(img, angles, μ=0.1) + iradon(r, angles, μ=0.1) end end diff --git a/src/filtered_backprojection.jl b/src/filtered_backprojection.jl new file mode 100644 index 0000000..48e3b8d --- /dev/null +++ b/src/filtered_backprojection.jl @@ -0,0 +1,21 @@ +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 iradon(sinogram, θs; geometry, μ) +end diff --git a/src/geometry.jl b/src/geometry.jl new file mode 100644 index 0000000..87f5d81 --- /dev/null +++ b/src/geometry.jl @@ -0,0 +1,38 @@ +export RadonGeometry, RadonParallelCircle, RadonFlexibleCircle + +abstract type RadonGeometry end + +""" + RadonParallelCircle(N, in_height) + +`N` is the size of the first and second dimension of the array. +`in_height` is a vector or range indicating the positions of the detector +with respect to the midpoint which is located at `N ÷ 2 + 1`. + +So an array of size `N=10` the default definition is: `RadonParallelCircle(10, -4:4)` +So the resulting sinogram has the shape: `(9, length(angles), size(array, 3))` +""" +struct RadonParallelCircle{T} <: RadonGeometry + N::Int + in_height::AbstractVector{T} +end + + +""" + RadonFlexibleCircle(N, in_height, out_height) + +`N` is the size of the first and second dimension of the array. +`in_height` is a vector or range indicating the vertical positions of the rays entering the circle +with respect to the midpoint which is located at `N ÷ 2 + 1`. +`out_height` is a vector or range indicating the vertical positions of the rays exiting the circle +with respect to the midpoint which is located at `N ÷ 2 + 1`. + +One definition could be: `RadonFlexibleCircle(10, -4:4, zeros((9,)))` +It would describe rays which enter the circle at positions `-4:4` but all of them would focus at the position 0 when leaving the circle. +This is an extreme form of cone beam tomography. +""" +struct RadonFlexibleCircle{T} <: RadonGeometry + N::Int + in_height::AbstractVector{T} + out_height::AbstractVector{T} +end diff --git a/src/iradon.jl b/src/iradon.jl index c7b41f5..eb42741 100644 --- a/src/iradon.jl +++ b/src/iradon.jl @@ -1,51 +1,22 @@ -export iradon, filtered_backprojection +export iradon -""" - filtered_backprojection(sinogram, θs) - -Calculates the simple Filtered Backprojection in CT with applying a ramp filter -in Fourier space. -""" -function filtered_backprojection(sinogram::AbstractArray{T}, θs::AbstractVector, μ=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 iradon(sinogram, θs, μ) -end # handle 2D -iradon(sinogram::AbstractArray{T, 2}, angles::AbstractArray{T2, 1}, μ=nothing) where {T, T2} = - view(iradon(reshape(sinogram, (size(sinogram)..., 1)), T.(angles), μ), :, :, 1) +function iradon(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) +end -iradon(sinogram::AbstractArray{T, 3}, angles::AbstractArray{T2, 1}, μ=nothing) where {T, T2} = - iradon(sinogram, T.(angles), μ) """ - iradon(sinogram, θs, μ=nothing) - -Calculates the parallel inverse Radon transform of the `sinogram`. -The first two dimensions are y and x. The third dimension is z, the rotational axis. -Works either with a `AbstractArray{T, 3}` or `AbstractArray{T, 2}`. -`size(sinogram, 1)` has to be an odd number. And `size(sinogram, 2)` has to be equal to -`length(angles)`. -The inverse Radon transform is rotated around the pixel `size(sinogram, 1) ÷ 2`, so there -is always a real center pixel! - -`θs` is a vector or range storing the angles in radians. - -# Exponential IRadon Transform -If `μ != nothing`, then the rays are attenuated with `exp(-μ * dist)` where `dist` -is the distance to the circular boundary of the field of view. -`μ` is in units of pixel length. So `μ=1` corresponds to an attenuation of `exp(-1)` if propagated through one pixel. + iradon(sinogram, θs; ) -In principle, all backends of KernelAbstractions.jl should work but are not tested. CUDA and CPU arrays are actively tested. See also [`radon`](@ref). -# Examples +# Example +# # Examples ```jldoctest julia> arr = zeros((5,2)); arr[2,:] .= 1; arr[4, :] .= 1 2-element view(::Matrix{Float64}, 4, :) with eltype Float64: @@ -60,111 +31,113 @@ julia> iradon(arr, [0, π/2]) 0.0 0.0 0.1 0.0 0.1 0.0 0.0 0.1 0.2 0.1 0.2 0.0232051 0.0 0.0 0.0232051 0.0 0.0232051 0.0 - -julia> iradon(arr, [0, π/2], 1) # exponential -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.00145226 0.0 0.00145226 0.0 - 0.0 0.00145226 0.00789529 0.0107308 0.033117 0.0183994 - 0.0 0.0 0.0107308 0.0 0.0107308 0.0 - 0.0 0.00145226 0.033117 0.0107308 0.0583388 0.0183994 - 0.0 0.0 0.0183994 0.0 0.0183994 0.0 ``` """ -function iradon(sinogram::AbstractArray{T, 3}, angles::AbstractArray{T, 1}, μ=nothing) where T - @assert isodd(size(sinogram, 1)) && size(sinogram, 2) == length(angles) - backend=KernelAbstractions.get_backend(sinogram) - absorption_f = let μ=μ - if isnothing(μ) - (x, y, x_start, y_start) -> one(T) - else - (x, y, x_start, y_start) -> exp(-T(μ) * sqrt((x - x_start)^2 + (y - y_start)^2)) - end - end - # this is the actual size we are using for calculation - N = size(sinogram, 1) - N_angles = size(angles, 1) +function iradon(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, μ) +end + + +function _iradon(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) + backend = KernelAbstractions.get_backend(sinogram) + + # angles_T might be a normal vector instead of Cuvector. fix it. + angles = similar(sinogram, (size(angles_T, 1),)) + angles .= typeof(angles)(angles_T) - # the only significant allocation - img = similar(sinogram, (N + 1, N + 1, size(sinogram, 3))) + in_height = similar(sinogram, (size(geometry.in_height, 1),)) + in_height .= typeof(in_height)(geometry.in_height) + + if geometry isa RadonFlexibleCircle + out_height = similar(sinogram, (size(geometry.out_height, 1),)) + out_height .= typeof(out_height)(geometry.out_height) + else + out_height = in_height + end + + # geometry can be very densely sampled, hence the sinogram depends on geometry size + N = geometry.N + # we only propagate inside this in circle + # convert radius to correct float type, very important for performance! + radius = T((N - 1) ÷ 2) + # the midpoint of the array + # convert to good type + mid = T(N ÷ 2 + 1 + 1 // 2) + N_angles = length(angles) + + img = similar(sinogram, (N, N, size(sinogram, 3))) fill!(img, 0) - - # radius of the cylinder we are projecting through - radius = size(img, 1) ÷ 2 - 1 - # mid point, it is actually N ÷ 2 + 1 - # but because of how adress the indices, we need 1.5 instead of +1 - mid = size(img, 1) ÷ 2 + T(1.5) - # the y_dists samplings, in principle we can add this as function parameter - y_dists = similar(img, (size(img, 1) - 1, )) - y_dists .= -radius:radius - - #@show typeof(sinogram), typeof(img), typeof(y_dists), typeof(angles) - kernel! = iradon_kernel!(backend) - kernel!(sinogram::AbstractArray{T}, img, y_dists, angles, mid, radius, - absorption_f, - ndrange=(N, N_angles, size(img, 3))) + + # create an absorption function, maps just to 1 in case isnothing(μ) + absorb_f = make_absorption_f(μ, T) + + # of the kernel goes + kernel! = iradon_kernel2!(backend) + kernel!(img, sinogram, in_height, out_height, angles, mid, radius, absorb_f, + ndrange=(size(sinogram, 1), N_angles, size(img, 3))) KernelAbstractions.synchronize(backend) - return img + return img end -""" - iradon_kernel!(sinogram, img, - y_dists, angles, mid, radius, - absorption_f) -""" -@kernel function iradon_kernel!(sinogram::AbstractArray{T}, - img, y_dists, angles, mid, radius, - absorption_f) where T - # r is the index of the angles - # k is the index of the detector spatial coordinate - # i_z is the index for the z dimension - k, r, i_z = @index(Global, NTuple) - - angle = angles[r] - sinα, cosα = sincos(angle) +@kernel function iradon_kernel2!(img::AbstractArray{T}, sinogram::AbstractArray{T}, in_height, out_height, angles, mid, radius, absorb_f) where {T} + i, iangle, i_z = @index(Global, NTuple) - # x0, y0, x1, y1 beginning and end point of each ray - a, b, c, d = next_ray_on_circle(img, angle, y_dists[k], mid, radius, sinα, cosα) - a0, b0, c0, d0 = a, b, c, d - # different comparisons depending which direction the ray is propagating - cac = a <= c ? (a,c) -> a<=c : (a,c) -> a>=c - cbd = b <= d ? (b,d) -> b<=d : (b,d) -> b>=d - - l = 1 - # acculumators of the intensity + @inbounds sinα, cosα = sincos(angles[iangle]) + + @inbounds ybegin = T(in_height[i]) + @inbounds yend = T(out_height[i]) + + # map the detector positions on a circle + # also rotate according to the angle + x_dist_rot, y_dist_rot, x_dist_end_rot, y_dist_end_rot = + next_ray_on_circle(ybegin, yend, mid, radius, sinα, cosα) + + # new will be always the current coordinates + # end the final destination + # and old is the previous one + xnew = x_dist_rot + ynew = y_dist_rot + xend = x_dist_end_rot + yend = y_dist_end_rot + xold, yold = xnew, ynew + + + tmp = zero(T) - while cac(a, c) && cbd(b, d) - a_old, b_old = a, b - # would be good to move this branch outside of the while loop - # but maybe branch prediction is doing a good job here - if a ≈ c && b ≈ d - break - end - - # find the next intersection for the ray - a, b = find_next_intersection(a,b,c,d) - l += 1 + # we store old_direction + # if the sign of old_direction changes, we have to stop tracing + # because then we hit the end point + _, _, sxold, syold = next_cell_intersection(xnew, ynew, xend, yend) + while true + # find next intersection point with integer grid + xnew, ynew, sx, sy = next_cell_intersection(xnew, ynew, xend, yend) + + # if we leave the circle or the direction of marching changes, this is the end + inside_circ(xnew, ynew, mid, radius + T(0.5)) && + (sx == sxold && sy == syold) || break - # find the cell it is cutting through - cell_i, cell_j = find_cell((a_old, b_old), - (a, b)) - - # distance travelled through that cell - distance = sqrt((a_old - a) ^2 + - (b_old - b) ^2) - # cell value times distance travelled through - @inbounds Atomix.@atomic img[cell_i, cell_j, i_z] += - distance * sinogram[k, r, i_z] * absorption_f(a, b, a0, b0) - end + # switch of i and j intentional to keep it consistent with existing code + icell, jcell = find_cell(xnew, ynew, xold, yold) + + # calculate intersection distance + distance = sqrt((xnew - xold)^2 + (ynew - yold) ^2) + # add value to ray, potentially attenuate by attenuated exp factor + @inbounds Atomix.@atomic img[icell, jcell, i_z] += distance * sinogram[i, iangle, i_z] * absorb_f(xnew, ynew, x_dist_rot, y_dist_rot) + xold, yold = xnew, ynew + end end # define adjoint rules -function ChainRulesCore.rrule(::typeof(iradon), sinogram::AbstractArray, angles, μ=nothing) - res = iradon(sinogram, angles, μ) +function ChainRulesCore.rrule(::typeof(_iradon), array::AbstractArray, angles, + geometry, μ) + res = _iradon(array, angles, geometry, μ) function pb_iradon(ȳ) - ad = radon(unthunk(ȳ), angles, μ) - return NoTangent(), ad, NoTangent(), NoTangent() + ad = _radon(unthunk(ȳ), angles, geometry, μ) + return NoTangent(), ad, NoTangent(), NoTangent(), NoTangent() end return res, pb_iradon end diff --git a/src/radon.jl b/src/radon.jl index 330db7f..be83b2f 100644 --- a/src/radon.jl +++ b/src/radon.jl @@ -1,32 +1,47 @@ export radon +function make_absorption_f(μ, ::Type{T}) where T + if isnothing(μ) + @inline (x, y, x_start, y_start) -> one(T) + else + @inline (x, y, x_start, y_start) -> exp(-T(μ) * sqrt((x - x_start)^2 + (y - y_start)^2)) + end +end -radon(img::AbstractArray{T, 3}, angles::AbstractArray{T2, 1}, μ=nothing) where {T, T2} = - radon(img, T.(angles), μ) # handle 2D -radon(img::AbstractArray{T, 2}, angles::AbstractArray{T2, 1}, μ=nothing) where {T, T2} = - view(radon(reshape(img, (size(img)..., 1)), T.(angles), μ), :, :, 1) +function radon(img::AbstractArray{T, 2}, angles::AbstractArray{T2, 1}; + geometry=RadonParallelCircle(size(img, 1), -(size(img,1)-1)÷2:(size(img,1)-1)÷2), μ=nothing) where {T, T2} + view(radon(reshape(img, (size(img)..., 1)), angles; geometry, μ), :, :, 1) +end + """ - radon(I, θs, μ=nothing) + radon(img, θs; ) Calculates the parallel Radon transform of the AbstractArray `I`. The first two dimensions are y and x. The third dimension is z, the rotational axis. -`size(I, 1)` and `size(I, 2)` have to be equal and a even number. +`size(I, 1)` and `size(I, 2)` have to be equally sized. The Radon transform is rotated around the pixel `size(I, 1) ÷ 2 + 1`, so there -is always a real center pixel! +is always an integer center pixel! Works either with a `AbstractArray{T, 3}` or `AbstractArray{T, 2}`. `θs` is a vector or range storing the angles in radians. -# Exponential IRadon Transform + +In principle, all backends of KernelAbstractions.jl should work but are not tested. CUDA and CPU arrays are actively tested. + + +# Keywords +## `μ=nothing` - Exponential IRadon Transform If `μ != nothing`, then the rays are attenuated with `exp(-μ * dist)` where `dist` is the distance to the circular boundary of the field of view. `μ` is in units of pixel length. So `μ=1` corresponds to an attenuation of `exp(-1)` if propagated through one pixel. -In principle, all backends of KernelAbstractions.jl should work but are not tested. CUDA and CPU arrays are actively tested. +## `geometry=geometry=RadonParallelCircle(-(size(img,1)-1)÷2:(size(img,1)-1)÷2)` +This corresponds to a parallel Radon transform. +See `?RadonGeometries` for a full list of geometries. There is also the very flexible `RadonFlexibleCircle`. Please note: the implementation is not quite optimized for cache efficiency and @@ -47,103 +62,121 @@ julia> radon(arr, [0, π/4, π/2]) 0.0 0.0 0.0 ``` """ -function radon(img::AbstractArray{T, 3}, angles::AbstractArray{T, 1}, μ=nothing) where T - @assert iseven(size(img, 1)) && iseven(size(img, 2)) && size(img, 1) == size(img, 2) "Arrays has to be quadratic and even sized shape" - - backend = KernelAbstractions.get_backend(img) +function radon(img::AbstractArray{T, 3}, angles_T::AbstractVector; + geometry=RadonParallelCircle(size(img, 1), -(size(img,1)-1)÷2:(size(img,1)-1)÷2), + μ=nothing) where {T} + return _radon(img::AbstractArray{T, 3}, angles_T::AbstractVector, geometry, μ) +end - absorption_f = let μ=μ - if isnothing(μ) - (x, y, x_start, y_start) -> one(T) - else - (x, y, x_start, y_start) -> exp(-T(μ) * sqrt((x - x_start)^2 + (y - y_start)^2)) - end - end +""" - # this is the actual size we are using for calculation - N = size(img, 1) - 1 - N_angles = size(angles, 1) - - # the only significant allocation - sinogram = similar(img, (N, N_angles, size(img, 3))) - fill!(sinogram, 0) +internal method which handles the different dispatches. +""" +function _radon(img::AbstractArray{T, 3}, angles_T::AbstractVector, + geometry::Union{RadonParallelCircle, RadonFlexibleCircle}, μ) where {T} + + @assert size(img, 1) == size(img, 2) "Arrays has to be quadratically shaped" + @assert size(img, 1) == geometry.N + backend = KernelAbstractions.get_backend(img) + + # angles_T might be a normal vector instead of Cuvector. fix it. + angles = similar(img, (size(angles_T, 1),)) + angles .= typeof(angles)(angles_T) - # radius of the cylinder we are projecting through - radius = size(img, 1) ÷ 2 - 1 - # mid point, it is actually N ÷ 2 + 1 - # but because of how adress the indices, we need 1.5 instead of +1 - mid = size(img, 1) ÷ 2 + T(1.5) - # the y_dists samplings, in principle we can add this as function parameter - y_dists = similar(img, (size(img, 1) - 1, )) - y_dists .= -radius:radius + in_height = similar(img, (size(geometry.in_height, 1),)) + in_height .= typeof(in_height)(geometry.in_height) - #@show typeof(sinogram), typeof(img), typeof(y_dists), typeof(angles) + if geometry isa RadonFlexibleCircle + out_height = similar(img, (size(geometry.out_height, 1),)) + out_height .= typeof(out_height)(geometry.out_height) + else + out_height = in_height + end + + # geometry can be very densely sampled, hence the sinogram depends on geometry size + N_sinogram = length(geometry.in_height) + # we only propagate inside this in circle + # convert radius to correct float type, very important for performance! + radius = T((size(img, 1) - 1) ÷ 2) + # the midpoint of the array + # convert to good type + mid = T(size(img, 1) ÷ 2 + 1 + 1 // 2) + N_angles = length(angles) + + # final sinogram + sinogram = similar(img, (N_sinogram, N_angles, size(img, 3))) + fill!(sinogram, 0) + + # create an absorption function, maps just to 1 in case isnothing(μ) + absorb_f = make_absorption_f(μ, T) + + # of the kernel goes kernel! = radon_kernel!(backend) - kernel!(sinogram::AbstractArray{T}, img, y_dists, angles, mid, radius, - absorption_f, - ndrange=(N, N_angles, size(img, 3))) + kernel!(sinogram::AbstractArray{T}, img, in_height, out_height, angles, mid, radius, absorb_f, + ndrange=(N_sinogram, N_angles, size(img, 3))) KernelAbstractions.synchronize(backend) return sinogram end -""" - radon_kernel!(sinogram, img, - y_dists, angles, mid, radius, - absorption_f) +@inline inside_circ(ii, jj, mid, radius) = (ii - mid)^2 + (jj - mid)^2 ≤ radius ^2 -""" -@kernel function radon_kernel!(sinogram::AbstractArray{T}, - img, y_dists, angles, mid, radius, absorption_f) where T - # r is the index of the angles - # k is the index of the detector spatial coordinate - # i_z is the index for the z dimension - k, r, i_z = @index(Global, NTuple) - - angle = angles[r] - sinα, cosα = sincos(angle) - - # x0, y0, x1, y1 beginning and end point of each ray - a, b, c, d = next_ray_on_circle(img, angle, y_dists[k], mid, radius, sinα, cosα) - a0, b0, c0, d0 = a, b, c, d - # different comparisons depending which direction the ray is propagating - cac = a <= c ? (a,c) -> a<=c : (a,c) -> a>=c - cbd = b <= d ? (b,d) -> b<=d : (b,d) -> b>=d +@kernel function radon_kernel!(sinogram::AbstractArray{T}, img::AbstractArray{T}, in_height, out_height, angles, mid, radius, absorb_f) where {T} + i, iangle, i_z = @index(Global, NTuple) - l = 1 - # acculumators of the intensity + @inbounds sinα, cosα = sincos(angles[iangle]) + + @inbounds ybegin = T(in_height[i]) + @inbounds yend = T(out_height[i]) + + # map the detector positions on a circle + # also rotate according to the angle + x_dist_rot, y_dist_rot, x_dist_end_rot, y_dist_end_rot = + next_ray_on_circle(ybegin, yend, mid, radius, sinα, cosα) + + # new will be always the current coordinates + # end the final destination + # and old is the previous one + xnew = x_dist_rot + ynew = y_dist_rot + xend = x_dist_end_rot + yend = y_dist_end_rot + xold, yold = xnew, ynew + + tmp = zero(T) - while cac(a, c) && cbd(b, d) - a_old, b_old = a, b - # would be good to move this branch outside of the while loop - # but maybe branch prediction is doing a good job here - if a ≈ c && b ≈ d - break - end + # we store old_direction + # if the sign of old_direction changes, we have to stop tracing + # because then we hit the end point + _, _, sxold, syold = next_cell_intersection(xnew, ynew, xend, yend) + while true + # find next intersection point with integer grid + xnew, ynew, sx, sy = next_cell_intersection(xnew, ynew, xend, yend) + + # if we leave the circle or the direction of marching changes, this is the end + inside_circ(xnew, ynew, mid, radius + T(0.5)) && + (sx == sxold && sy == syold) || break - # find the next intersection for the ray - a, b = find_next_intersection(a,b,c,d) - l += 1 - - # find the cell it is cutting through - cell_i, cell_j = find_cell((a_old, b_old), - (a, b)) - - # distance travelled through that cell - distance = sqrt((a_old - a) ^2 + - (b_old - b) ^2) - # cell value times distance travelled through - @inbounds tmp += distance * img[cell_i, cell_j, i_z] * absorption_f(a, b, a0, b0) - end - @inbounds sinogram[k, r, i_z] = tmp + # switch of i and j intentional to keep it consistent with existing code + icell, jcell = find_cell(xnew, ynew, xold, yold) + + # calculate intersection distance + distance = sqrt((xnew - xold)^2 + (ynew - yold) ^2) + # add value to ray, potentially attenuate by attenuated exp factor + @inbounds tmp += distance * img[icell, jcell, i_z] * absorb_f(xnew, ynew, x_dist_rot, y_dist_rot) + xold, yold = xnew, ynew + end + @inbounds sinogram[i, iangle, i_z] = tmp end + # define adjoint rules -function ChainRulesCore.rrule(::typeof(radon), array::AbstractArray, angles, μ=nothing) - res = radon(array, angles, μ) +function ChainRulesCore.rrule(::typeof(_radon), array::AbstractArray, angles, + geometry, μ) + res = _radon(array, angles, geometry, μ) function pb_radon(ȳ) - ad = iradon(unthunk(ȳ), angles, μ) - return NoTangent(), ad, NoTangent(), NoTangent() + ad = _iradon(unthunk(ȳ), angles, geometry, μ) + return NoTangent(), ad, NoTangent(), NoTangent(), NoTangent() end return res, pb_radon end diff --git a/src/radon_legacy.jl b/src/radon_legacy.jl new file mode 100644 index 0000000..491e47c --- /dev/null +++ b/src/radon_legacy.jl @@ -0,0 +1,374 @@ +export radon_old +export iradon_old + + +# handle 2D +iradon_old(sinogram::AbstractArray{T, 2}, angles::AbstractArray{T2, 1}, μ=nothing; + ray_startpoints=nothing, ray_endpoints=nothing) where {T, T2} = begin + angles_T = similar(sinogram, (size(angles, 1),)) + angles_T .= angles + view(iradon_old(reshape(sinogram, (size(sinogram)..., 1)), angles_T, μ; ray_startpoints, ray_endpoints), :, :, 1) +end + + +iradon_old(sinogram::AbstractArray{T, 3}, angles::AbstractArray{T2, 1}, μ=nothing; + ray_startpoints=nothing, ray_endpoints=nothing) where {T, T2} = begin + angles_T = similar(sinogram, (size(angles, 1),)) + angles_T .= angles + iradon_old(sinogram, angles_T, μ; ray_startpoints, ray_endpoints) +end +""" + iradon_old(sinogram, θs, μ=nothing; ) + +Calculates the parallel inverse radon_old transform of the `sinogram`. +The first two dimensions are y and x. The third dimension is z, the rotational axis. +Works either with a `AbstractArray{T, 3}` or `AbstractArray{T, 2}`. +`size(sinogram, 1)` has to be an odd number. And `size(sinogram, 2)` has to be equal to +`length(angles)`. +The inverse radon_old transform is rotated around the pixel `size(sinogram, 1) ÷ 2`, so there +is always a real center pixel! + +`θs` is a vector or range storing the angles in radians. + +# Exponential Iradon_old Transform +If `μ != nothing`, then the rays are attenuated with `exp(-μ * dist)` where `dist` +is the distance to the circular boundary of the field of view. +`μ` is in units of pixel length. So `μ=1` corresponds to an attenuation of `exp(-1)` if propagated through one pixel. + + +# Keywords +* `ray_startpoints = -size(img, 1) ÷ 2 - 1: size(img, 1) ÷ 2 - 1`: this indicates the starting position of a ray at the entrance on the array +* `ray_endpoints = -size(img, 1) ÷ 2 - 1: size(img, 1) ÷ 2 - 1`: this indicates the end position of a ray at the exit of the array + +Even though we trace everything inside the circle only, start and endpoint of the ray are all the first and last index position in the array. +So it's not the y position on the circle but rather the y position at the smallest and greatest x value. + + +In principle, all backends of KernelAbstractions.jl should work but are not tested. CUDA and CPU arrays are actively tested. + +See also [`radon_old`](@ref). + +# Examples +```jldoctest +julia> arr = zeros((5,2)); arr[2,:] .= 1; arr[4, :] .= 1 +2-element view(::Matrix{Float64}, 4, :) with eltype Float64: + 1.0 + 1.0 + +julia> iradon_old(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.1 0.0 0.1 0.0 + 0.0 0.1 0.2 0.1 0.2 0.0232051 + 0.0 0.0 0.1 0.0 0.1 0.0 + 0.0 0.1 0.2 0.1 0.2 0.0232051 + 0.0 0.0 0.0232051 0.0 0.0232051 0.0 + +julia> iradon_old(arr, [0, π/2], 1) # exponential +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.00145226 0.0 0.00145226 0.0 + 0.0 0.00145226 0.00789529 0.0107308 0.033117 0.0183994 + 0.0 0.0 0.0107308 0.0 0.0107308 0.0 + 0.0 0.00145226 0.033117 0.0107308 0.0583388 0.0183994 + 0.0 0.0 0.0183994 0.0 0.0183994 0.0 +``` +""" +function iradon_old(sinogram::AbstractArray{T, 3}, angles::AbstractArray{T, 1}, μ=nothing; + ray_startpoints=nothing, ray_endpoints=nothing) where T + @assert isodd(size(sinogram, 1)) && size(sinogram, 2) == length(angles) + backend=KernelAbstractions.get_backend(sinogram) + absorption_f = let μ=μ + if isnothing(μ) + (x, y, x_start, y_start) -> one(T) + else + (x, y, x_start, y_start) -> exp(-T(μ) * sqrt((x - x_start)^2 + (y - y_start)^2)) + end + end + # this is the actual size we are using for calculation + N = size(sinogram, 1) + N_angles = size(angles, 1) + + # the only significant allocation + img = similar(sinogram, (N + 1, N + 1, size(sinogram, 3))) + fill!(img, 0) + + # radius of the cylinder we are projecting through + radius = size(img, 1) ÷ 2 - 1 + # mid point, it is actually N ÷ 2 + 1 + # but because of how adress the indices, we need 1.5 instead of +1 + mid = size(img, 1) ÷ 2 + T(1.5) + + if isnothing(ray_endpoints) + # the y_dists samplings, in principle we can add this as function parameter + y_dists_end = similar(img, (size(img, 1) - 1, )) + y_dists_end .= -radius:radius + else + y_dists_end = similar(img, (size(img, 1) - 1, )) + y_dists_end .= ray_endpoints + end + + if isnothing(ray_startpoints) + y_dists = similar(img, (size(img, 1) - 1, )) + y_dists .= -radius:radius + else + y_dists = similar(img, (size(img, 1) - 1, )) + y_dists .= ray_startpoints + end + + #@show typeof(sinogram), typeof(img), typeof(y_dists), typeof(angles) + kernel! = iradon_old_kernel!(backend) + kernel!(sinogram::AbstractArray{T}, img, y_dists, y_dists_end, angles, mid, radius, + absorption_f, + ndrange=(N, N_angles, size(img, 3))) + KernelAbstractions.synchronize(backend) + return img +end + +""" + iradon_old_kernel!(sinogram, img, + y_dists, angles, mid, radius, + absorption_f) +""" +@kernel function iradon_old_kernel!(sinogram::AbstractArray{T}, + img, y_dists, y_dists_end, angles, mid, radius, + absorption_f) where T + # r is the index of the angles + # k is the index of the detector spatial coordinate + # i_z is the index for the z dimension + k, r, i_z = @index(Global, NTuple) + + angle = angles[r] + sinα, cosα = sincos(angle) + + # x0, y0, x1, y1 beginning and end point of each ray + a, b, c, d = next_ray_on_circle(y_dists[k], y_dists_end[k], mid, radius, sinα, cosα) + a0, b0, c0, d0 = a, b, c, d + # different comparisons depending which direction the ray is propagating + cac = a <= c ? (a,c) -> a<=c : (a,c) -> a>=c + cbd = b <= d ? (b,d) -> b<=d : (b,d) -> b>=d + + l = 1 + # acculumators of the intensity + tmp = zero(T) + while cac(a, c) && cbd(b, d) + a_old, b_old = a, b + # would be good to move this branch outside of the while loop + # but maybe branch prediction is doing a good job here + if a ≈ c && b ≈ d + break + end + + # find the next intersection for the ray + a, b = find_next_intersection(a,b,c,d) + l += 1 + + # find the cell it is cutting through + cell_i, cell_j = find_cell((a_old, b_old), + (a, b)) + + # distance travelled through that cell + distance = sqrt((a_old - a) ^2 + + (b_old - b) ^2) + # cell value times distance travelled through + @inbounds Atomix.@atomic img[cell_i, cell_j, i_z] += + distance * sinogram[k, r, i_z] * absorption_f(a, b, a0, b0) + end +end + + + # define adjoint rules +function ChainRulesCore.rrule(::typeof(iradon_old), sinogram::AbstractArray, angles, μ=nothing; ray_endpoints=nothing) + res = iradon_old(sinogram, angles, μ; ray_endpoints) + function pb_iradon_old(ȳ) + ad = radon_old(unthunk(ȳ), angles, μ; ray_endpoints) + return NoTangent(), ad, NoTangent(), NoTangent() + end + return res, pb_iradon_old +end + + +radon_old(img::AbstractArray{T, 3}, angles::AbstractArray{T2, 1}, μ=nothing; + ray_startpoints=nothing, ray_endpoints=nothing) where {T, T2} = begin + angles_T = similar(img, (size(angles, 1),)) + angles_T .= angles + radon_old(img, angles_T, μ; ray_startpoints, ray_endpoints) +end + +# handle 2D +radon_old(img::AbstractArray{T, 2}, angles::AbstractArray{T2, 1}, μ=nothing; + ray_startpoints=nothing, ray_endpoints=nothing) where {T, T2} = begin + angles_T = similar(img, (size(angles, 1),)) + angles_T .= angles + view(radon_old(reshape(img, (size(img)..., 1)), angles_T, μ; + ray_startpoints, ray_endpoints), :, :, 1) +end +""" + radon_old(I, θs, μ=nothing; ) + +Calculates the parallel radon_old transform of the AbstractArray `I`. +The first two dimensions are y and x. The third dimension is z, the rotational axis. +`size(I, 1)` and `size(I, 2)` have to be equal and a even number. +The radon_old transform is rotated around the pixel `size(I, 1) ÷ 2 + 1`, so there +is always a real center pixel! +Works either with a `AbstractArray{T, 3}` or `AbstractArray{T, 2}`. + +`θs` is a vector or range storing the angles in radians. + +# Exponential Iradon_old Transform +If `μ != nothing`, then the rays are attenuated with `exp(-μ * dist)` where `dist` +is the distance to the circular boundary of the field of view. +`μ` is in units of pixel length. So `μ=1` corresponds to an attenuation of `exp(-1)` if propagated through one pixel. + +In principle, all backends of KernelAbstractions.jl should work but are not tested. CUDA and CPU arrays are actively tested. + + +# Keywords +* `ray_startpoints = -size(img, 1) ÷ 2 - 1: size(img, 1) ÷ 2 - 1`: this indicates the starting position of a ray at the entrance on the array +* `ray_endpoints = -size(img, 1) ÷ 2 - 1: size(img, 1) ÷ 2 - 1`: this indicates the end position of a ray at the exit of the array + +Even though we trace everything inside the circle only, start and endpoint of the ray are all the first and last index position in the array. +So it's not the y position on the circle but rather the y position at the smallest and greatest x value. + + + +Please note: the implementation is not quite optimized for cache efficiency and +it is a very naive algorithm. But still, it is fast! + +See also [`iradon_old`](@ref). + +# Example +The reason the sinogram has the value `1.41421` for the diagonal ray `π/4` is, +that such a diagonal travels a longer distance through the pixel. +```jldoctest +julia> arr = zeros((4,4)); arr[3,3] = 1; + +julia> radon_old(arr, [0, π/4, π/2]) +3×3 view(::Array{Float64, 3}, :, :, 1) with eltype Float64: + 0.0 0.0 0.0 + 1.0 1.41421 1.0 + 0.0 0.0 0.0 +``` +""" +function radon_old(img::AbstractArray{T, 3}, angles::AbstractArray{T, 1}, + μ=nothing; ray_startpoints=nothing, ray_endpoints=nothing) where T + @assert iseven(size(img, 1)) && iseven(size(img, 2)) && size(img, 1) == size(img, 2) "Arrays has to be quadratic and even sized shape" + + backend = KernelAbstractions.get_backend(img) + + absorption_f = let μ=μ + if isnothing(μ) + (x, y, x_start, y_start) -> one(T) + else + (x, y, x_start, y_start) -> exp(-T(μ) * sqrt((x - x_start)^2 + (y - y_start)^2)) + end + end + + # this is the actual size we are using for calculation + N = size(img, 1) - 1 + N_angles = size(angles, 1) + + # the only significant allocation + sinogram = similar(img, (N, N_angles, size(img, 3))) + fill!(sinogram, 0) + + # radius of the cylinder we are projecting through + radius = size(img, 1) ÷ 2 - 1 + # mid point, it is actually N ÷ 2 + 1 + # but because of how adress the indices, we need 1.5 instead of +1 + mid = size(img, 1) ÷ 2 + T(1.5) + + + if isnothing(ray_endpoints) + # the y_dists samplings, in principle we can add this as function parameter + y_dists_end = similar(img, (size(img, 1) - 1, )) + y_dists_end .= -radius:radius + else + y_dists_end = similar(img, (size(img, 1) - 1, )) + y_dists_end .= ray_endpoints + end + + if isnothing(ray_startpoints) + y_dists = similar(img, (size(img, 1) - 1, )) + y_dists .= -radius:radius + else + y_dists = similar(img, (size(img, 1) - 1, )) + y_dists .= ray_startpoints + end + + + # the y_dists samplings, in principle we can add this as function parameter + y_dists = similar(img, (size(img, 1) - 1, )) + y_dists .= -radius:radius + + #@show typeof(sinogram), typeof(img), typeof(y_dists), typeof(angles) + kernel! = radon_old_kernel!(backend) + kernel!(sinogram::AbstractArray{T}, img, y_dists, y_dists_end, angles, mid, radius, + absorption_f, + ndrange=(N, N_angles, size(img, 3))) + KernelAbstractions.synchronize(backend) + return sinogram +end + +""" + radon_old_kernel!(sinogram, img, + y_dists, angles, mid, radius, + absorption_f) + +""" +@kernel function radon_old_kernel!(sinogram::AbstractArray{T}, + img, y_dists, y_dists_end, angles, mid, radius, absorption_f) where T + # r is the index of the angles + # k is the index of the detector spatial coordinate + # i_z is the index for the z dimension + k, r, i_z = @index(Global, NTuple) + + angle = angles[r] + sinα, cosα = sincos(angle) + + # x0, y0, x1, y1 beginning and end point of each ray + a, b, c, d = next_ray_on_circle(y_dists[k], y_dists_end[k], mid, radius, sinα, cosα) + a0, b0, c0, d0 = a, b, c, d + # different comparisons depending which direction the ray is propagating + cac = a <= c ? (a,c) -> a<=c : (a,c) -> a>=c + cbd = b <= d ? (b,d) -> b<=d : (b,d) -> b>=d + + l = 1 + # acculumators of the intensity + tmp = zero(T) + while cac(a, c) && cbd(b, d) + a_old, b_old = a, b + # would be good to move this branch outside of the while loop + # but maybe branch prediction is doing a good job here + if a ≈ c && b ≈ d + break + end + + # find the next intersection for the ray + a, b = find_next_intersection(a,b,c,d) + l += 1 + + # find the cell it is cutting through + cell_i, cell_j = find_cell((a_old, b_old), + (a, b)) + + # distance travelled through that cell + distance = sqrt((a_old - a) ^2 + + (b_old - b) ^2) + # cell value times distance travelled through + @inbounds tmp += distance * img[cell_i, cell_j, i_z] * absorption_f(a, b, a0, b0) + end + @inbounds sinogram[k, r, i_z] = tmp +end + + + # define adjoint rules +function ChainRulesCore.rrule(::typeof(radon_old), array::AbstractArray, angles, μ=nothing; + ray_startpoints=nothing, ray_endpoints=nothing) + res = radon_old(array, angles, μ; ray_startpoints, ray_endpoints) + function pb_radon_old(ȳ) + ad = iradon_old(unthunk(ȳ), angles, μ; ray_startpoints, ray_endpoints) + return NoTangent(), ad, NoTangent(), NoTangent() + end + return res, pb_radon_old +end diff --git a/src/utils.jl b/src/utils.jl index 1e49a63..97e4afd 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -138,11 +138,64 @@ corner of the cell. return x, y end +""" + find_cell(x, y, xnew, ynew) + +Find the correct intersected cell. +`a` and `b` are a Vector of length 2 + +Cells are identified by the coordinates in the lower left +corner of the cell. + + +|---------------------------| +| | | | +|(3,1) |(3,2) |(3,3) | +|--------|---------|--------| +| | | | +|(2,1) |(2,2) |(2,3) | +|--------|---------|--------| +| | | | +|(1,1) |(1,2) |(1,3) | +|--------|---------|--------| +""" +@inline function find_cell(x, y, xnew, ynew) + x = floor(Int, (x + xnew) / 2) + y = floor(Int, (y + ynew) / 2) + return x, y +end + + +@inline function next_cell_intersection(x₀, y₀, x₁, y₁) + # inspired by + floor_or_ceilx = x₁ ≥ x₀ ? ceil : floor + floor_or_ceily = y₁ ≥ y₀ ? ceil : floor + # https://cs.stackexchange.com/questions/132887/determining-the-intersections-of-a-line-segment-and-grid + txx = (floor_or_ceilx(x₀)-x₀) + xdiff = x₁ - x₀ + ydiff = y₁ - y₀ + txx = ifelse(txx != 0, txx, txx + sign(xdiff)) + tyy = (floor_or_ceily(y₀)-y₀) + tyy = ifelse(tyy != 0, tyy, tyy + sign(ydiff)) + + tx = txx / (xdiff) + ty = tyy / (ydiff) + # decide which t is smaller, and hence the next step + t = ifelse(tx > ty, ty, tx) + + # calculate new coordinates + x = x₀ + t * (xdiff) + y = y₀ + t * (ydiff) + return x, y, sign(xdiff), sign(ydiff) +end + + + -function next_ray_on_circle(img, angle, y_dist, mid, radius, sinα, cosα) +@inline function next_ray_on_circle(y_dist, y_dist_end, mid, radius, sinα, cosα) x_dist = sqrt(radius^2 - y_dist^2) - x_dist_end = -x_dist - y_dist_end = y_dist + x_dist_end = -sqrt(radius^2 - y_dist_end^2) + y_dist_end = y_dist_end x_dist_rot = mid + cosα * x_dist - sinα * y_dist y_dist_rot = mid + sinα * x_dist + cosα * y_dist diff --git a/test/runtests.jl b/test/runtests.jl index abd3d02..58ec225 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,8 +11,8 @@ using Zygote sinogram[5, :, 1] .= 1 @test iradon(sinogram, [0.0f0, π * 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 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 0.5; 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 0.5 0.0 0.0 0.0 0.0;;;] - @test iradon(sinogram[:, 1, :], Float32[π / 4 + π]) ≈ [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.46446568 6.7434956f-7 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.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.4142135 0.0; 0.0 0.0 0.0 0.0 0.0 0.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.46446568 6.7434956f-7 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.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 "Exponential iradon" begin @@ -21,9 +21,9 @@ using Zygote angles = [0] - arr = iradon(sinogram, angles, 0.1) + arr = iradon(sinogram, angles, μ=0.1) - exp.(-(8.5:-1:1.5) .* 0.1) ≈ arr[6, begin+1:end-1] + @test exp.(-(8.5:-1:1.5) .* 0.1) ≈ arr[begin+1:end-1, 6][:] end @@ -71,19 +71,25 @@ using Zygote sinogram2 = radon(array3, angles2) array_filtered = filtered_backprojection(sinogram2, angles2) - @test ≈(array_filtered / sum(array_filtered) .+ 0.1, array3 / sum(array3) .+ 0.1, rtol=0.05) + @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) + @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 @testset "Test gradients" begin - x = randn((10, 10)) - test_rrule(radon, x, [0, π/4, π/2, 2π, 0.1, 0.00001] ⊢ ChainRulesTestUtils.NoTangent(), nothing ⊢ ChainRulesTestUtils.NoTangent()) + x = randn((10, 10, 1)) + geometry=RadonParallelCircle(size(x, 1), -(size(x,1)-1)÷2:(size(x,1)-1)÷2) + 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(iradon, y, [0, π/4, π/2, 2π, 0.1, 0.00001] ⊢ ChainRulesTestUtils.NoTangent(), nothing ⊢ ChainRulesTestUtils.NoTangent()) + test_rrule(RadonKA._iradon, y, [0, π/4, π/2, 2π, 0.1, 0.00001] ⊢ ChainRulesTestUtils.NoTangent(), geometry ⊢ ChainRulesTestUtils.NoTangent(), nothing ⊢ ChainRulesTestUtils.NoTangent()) - x = randn((10, 10)) - test_rrule(radon, x, [0, π/4, π/2, 2π, 0.1, 0.00001] ⊢ ChainRulesTestUtils.NoTangent(), 0.1⊢ ChainRulesTestUtils.NoTangent()) + x = randn((10, 10, 1)) + 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(iradon, y, [0, π/4, π/2, 2π, 0.1, 0.00001] ⊢ ChainRulesTestUtils.NoTangent(), 0.1 ⊢ ChainRulesTestUtils.NoTangent()) + test_rrule(RadonKA._iradon, y, [0, π/4, π/2, 2π, 0.1, 0.00001] ⊢ ChainRulesTestUtils.NoTangent(), geometry ⊢ ChainRulesTestUtils.NoTangent(), 0.1 ⊢ ChainRulesTestUtils.NoTangent()) end