Add benchmark section to docs

docs/src/benchmark.md

@@ -0,0 +1,164 @@
+# 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.
+
+# 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        |
+
+
+
+## RadonKA.jl
+See this [benchmark example](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());
+```
+![](../assets/radonka_sinogram.png)
+![](../assets/radonka_iradon.png)
+
+
+## Matlab (R2023a)
+```matlab
+arr = single(imread("/tmp/sample.jpg")); 
+%arr = rand(1920, 1920, "single");
+theta = linspace(0, 360, 500);
+
+tic;
+R = radon(arr, theta);
+toc;
+R = R / max(R(:));
+imwrite(R, "/tmp/matlab_sinogram.png");
+
+tic; 
+iR = iradon(R, theta, "linear", "none");
+toc;
+iR = iR / max(iR(:));
+imwrite(iR, "/tmp/matlab_iradon.png");
+
+
+
+x = (rand(100, 512, 512, 'single'));
+
+theta = linspace(0, 360, 500);
+y = (zeros(size(x, 1), size(radon(squeeze(x(1, :, :)), theta), 1), size(radon(squeeze(x(1, :, :)), theta), 2), 'single'));
+ix = zeros(100, 514, 514, "single");
+
+tic;
+for i = 1:size(x, 1)
+    y(i, :, :) = radon(squeeze(x(i, :, :)), theta);
+end
+toc;
+
+
+tic;
+for i = 1:size(y, 1)
+    ix(i, :, :) = iradon(squeeze(y(i, :, :)), theta);
+end
+toc;
+```
+
+![](../assets/matlab_sinogram.png)
+![](../assets/matlab_iradon.png)
+
+## Astra
+Some of the benchmarks did not run properly or were providing non-sense.
+Astra's docs are unfortunately slightly confusing...
+```python
+import numpy as np
+import matplotlib.pyplot as plt
+import imageio.v3 as iio
+import astra
+
+im = np.array(iio.imread('/tmp/sample.png'), dtype=np.float32)
+
+plt.imshow(im)
+
+
+vol_geom = astra.create_vol_geom(1920, 1920)
+angles =  np.linspace(0,2 * np.pi,500)
+
+proj_geom = astra.create_proj_geom('parallel', 1.0, 1920, angles)
+proj_id = astra.create_projector('line', proj_geom,vol_geom)
+
+
+%%time
+sinogram_id, sinogram = astra.create_sino(im, proj_id)
+np.copy(sinogram);
+
+plt.imshow(sinogram)
+
+proj_geom = astra.create_proj_geom('parallel', 1.0, 1920, angles)
+proj_id = astra.create_projector('cuda', proj_geom,vol_geom)
+
+%%time
+sinogram_id, sinogram = astra.create_sino(im, proj_id)
+np.copy(sinogram);
+
+rec_id = astra.data2d.create("-vol", vol_geom)
+cfg = astra.astra_dict('BP')
+cfg["ReconstructionDataId"] = rec_id
+cfg["ProjectionDataId"] = sinogram_id
+cfg["ProjectorId"] = proj_id
+# Create the algorithm object from the configuration structure
+alg_id = astra.algorithm.create(cfg)
+# Run back-projection and get the reconstruction
+
+%%time
+astra.algorithm.run(alg_id)
+f_rec = astra.data2d.get(rec_id)
+#np.copy(f_rec)
+
+
+im_3D = np.random.rand(100, 512, 512)
+
+vol_geom = astra.create_vol_geom(512, 512, 100)
+
+proj_geom = astra.create_proj_geom('parallel3d', 1.0, 1.0, 512, 512, angles)
+#proj_id = astra.create_projector('line', proj_geom,vol_geom)
+
+%%time
+proj_id, proj_data = astra.create_sino3d_gpu(im_3D, proj_geom, vol_geom)
+np.copy(proj_data)
+
+
+rec_id = astra.data3d.create('-vol', vol_geom)
+
+# Set up the parameters for a reconstruction algorithm using the GPU
+cfg = astra.astra_dict('BP3D_CUDA')
+cfg['ReconstructionDataId'] = rec_id
+cfg['ProjectionDataId'] = proj_id
+
+
+# Create the algorithm object from the configuration structure
+alg_id = astra.algorithm.create(cfg)
+
+%time
+astra.algorithm.run(alg_id, 1)
+rec = astra.data3d.get(rec_id)
+np.copy(rec)
+```
+
+![](../assets/astra_sinogram.png)
+![](../assets/astra_iradon.png)
+

examples/benchmark_radon_iradon.jl

@@ -0,0 +1,139 @@
+### A Pluto.jl notebook ###
+# v0.19.32
+
+using Markdown
+using InteractiveUtils
+
+# ╔═╡ e5a21d4e-928b-11ee-3909-53f97530eefa
+begin
+	using Pkg
+	Pkg.activate(".")
+	using Revise
+end
+
+# ╔═╡ 3757ba3d-8ff2-4bc8-93bc-850592ac0ec9
+Pkg.add("FileIO")
+
+# ╔═╡ 114e830d-9bd0-4546-ad6d-2e62b4a1194b
+using IndexFunArrays, ImageShow, Plots, ImageIO, PlutoUI, PlutoTest, TestImages
+
+# ╔═╡ dc48ab00-c044-4292-a908-a155af022fde
+using FileIO
+
+# ╔═╡ a4dc63b8-4ba3-42fe-9743-3d7169604415
+using RadonKA
+
+# ╔═╡ c9b6dc38-dcfd-4b1f-91a8-2d105cda9c40
+using BenchmarkTools
+
+# ╔═╡ 1fddfd52-e964-4e36-a2f6-53302aa7e217
+using CUDA, CUDA.CUDAKernels
+
+# ╔═╡ a96e7332-dab7-41f0-b979-9dcf39c9e0bf
+sample_2D = Float32.(testimage("resolution_test_1920"));
+
+# ╔═╡ 50565c52-fa1a-4ed2-96d8-15ca67b44081
+save("/tmp/sample.png", sample_2D)
+
+# ╔═╡ 8a7a0f9d-fa4f-4a7f-b385-a02d441aa974
+sample_2D_c = CuArray(sample_2D);
+
+# ╔═╡ 3c11a680-a6f2-4f81-9ae0-091ae95002e5
+simshow(sample_2D)
+
+# ╔═╡ df8f25ee-dc19-442d-97fe-3100fb31c5d0
+angles = range(0, 2π, 500);
+
+# ╔═╡ 4adb8e5d-e6ca-4978-bb3f-98788b38defe
+angles_c = CuArray(angles);
+
+# ╔═╡ 6672cff8-e1af-4425-8c19-fbd2448f40ee
+@benchmark sinogram = radon($sample_2D, $angles)
+
+# ╔═╡ 3aaf53a3-1087-4d8f-9533-1963319ece01
+sinogram = radon(sample_2D, angles)
+
+# ╔═╡ 1fb7d9bd-daa8-4898-80e2-282d03bb53a8
+simshow(sinogram)
+
+# ╔═╡ 76df1a12-2e52-4440-8c39-9f126715e77c
+save("/tmp/radonka_sinogram.png", sinogram ./ maximum(sinogram));
+
+# ╔═╡ a11fe70a-3425-4335-ba09-00efbae6c9af
+ sample_iradon = iradon(sinogram,angles);
+
+# ╔═╡ 29f7d88c-3bab-43a9-a253-d3bcfc98b706
+save("/tmp/radonka_iradon.png", sample_iradon ./ maximum(sample_iradon));
+
+# ╔═╡ da297929-f4f9-4fa2-b860-cfc159ac39f1
+@benchmark sample_iradon = iradon($sinogram, $angles)
+
+# ╔═╡ a61ea8f2-6329-4b6b-b4a2-1f5ac17bc618
+sinogram_c = radon(sample_2D_c, angles_c, backend=CUDABackend())
+
+# ╔═╡ c7bb742c-5f2e-44e1-9304-b14db275043d
+@benchmark CUDA.@sync sinogram_c = radon($sample_2D_c, $angles_c, backend=CUDABackend())
+
+# ╔═╡ 1b8e4e3f-d9db-4646-8ccf-840c21ecc328
+@benchmark CUDA.@sync sample_iradon_c = iradon($sinogram_c, $angles_c, backend=CUDABackend())
+
+# ╔═╡ c602a8df-d646-4f31-9d6f-902456f07ae5
+md"# 3D"
+
+# ╔═╡ 3570e170-77e5-449d-8e97-5ed6090ebc55
+sample_3D = randn(Float32, (512, 512, 100));
+
+# ╔═╡ d5e9a965-1e93-4998-b953-57ca0fdd1a50
+sample_3D_c = CuArray(sample_3D)
+
+# ╔═╡ 137f6c8b-c774-49e4-b163-eede8b60abfd
+@time sinogram_3D = radon(sample_3D, angles);
+
+# ╔═╡ a5bf8b01-aefb-42e6-8ec5-588190bb3fa2
+@benchmark radon($sample_3D, $angles)
+
+# ╔═╡ b5825676-6600-4cbe-865b-788ef1567b85
+@benchmark iradon($sinogram_3D, $angles)
+
+# ╔═╡ 494190e4-e28a-42ad-b1c1-926a740ea2cb
+sinogram_3D_c = radon(sample_3D_c, angles_c, backend=CUDABackend())
+
+# ╔═╡ 52c781ab-c08a-4b10-9568-011abec63bda
+@benchmark CUDA.@sync radon($sample_3D_c, $angles_c, backend=CUDABackend())
+
+# ╔═╡ cd166ef5-863f-48b9-bb89-a2a32daa166b
+@benchmark CUDA.@sync iradon($sinogram_3D_c, $angles_c, backend=CUDABackend())
+
+# ╔═╡ Cell order:
+# ╠═e5a21d4e-928b-11ee-3909-53f97530eefa
+# ╠═114e830d-9bd0-4546-ad6d-2e62b4a1194b
+# ╠═dc48ab00-c044-4292-a908-a155af022fde
+# ╠═3757ba3d-8ff2-4bc8-93bc-850592ac0ec9
+# ╠═a4dc63b8-4ba3-42fe-9743-3d7169604415
+# ╠═c9b6dc38-dcfd-4b1f-91a8-2d105cda9c40
+# ╠═1fddfd52-e964-4e36-a2f6-53302aa7e217
+# ╠═a96e7332-dab7-41f0-b979-9dcf39c9e0bf
+# ╠═50565c52-fa1a-4ed2-96d8-15ca67b44081
+# ╠═8a7a0f9d-fa4f-4a7f-b385-a02d441aa974
+# ╠═3c11a680-a6f2-4f81-9ae0-091ae95002e5
+# ╠═df8f25ee-dc19-442d-97fe-3100fb31c5d0
+# ╠═4adb8e5d-e6ca-4978-bb3f-98788b38defe
+# ╠═6672cff8-e1af-4425-8c19-fbd2448f40ee
+# ╠═3aaf53a3-1087-4d8f-9533-1963319ece01
+# ╠═1fb7d9bd-daa8-4898-80e2-282d03bb53a8
+# ╠═76df1a12-2e52-4440-8c39-9f126715e77c
+# ╠═29f7d88c-3bab-43a9-a253-d3bcfc98b706
+# ╠═a11fe70a-3425-4335-ba09-00efbae6c9af
+# ╠═da297929-f4f9-4fa2-b860-cfc159ac39f1
+# ╠═a61ea8f2-6329-4b6b-b4a2-1f5ac17bc618
+# ╠═c7bb742c-5f2e-44e1-9304-b14db275043d
+# ╠═1b8e4e3f-d9db-4646-8ccf-840c21ecc328
+# ╠═c602a8df-d646-4f31-9d6f-902456f07ae5
+# ╠═3570e170-77e5-449d-8e97-5ed6090ebc55
+# ╠═d5e9a965-1e93-4998-b953-57ca0fdd1a50
+# ╠═137f6c8b-c774-49e4-b163-eede8b60abfd
+# ╠═a5bf8b01-aefb-42e6-8ec5-588190bb3fa2
+# ╠═b5825676-6600-4cbe-865b-788ef1567b85
+# ╠═494190e4-e28a-42ad-b1c1-926a740ea2cb
+# ╠═52c781ab-c08a-4b10-9568-011abec63bda
+# ╠═cd166ef5-863f-48b9-bb89-a2a32daa166b