Skip to content

Commit

Permalink
Merge pull request #9 from roflmaostc/rename_iradon
Browse files Browse the repository at this point in the history
Rename iradon to backproject
  • Loading branch information
roflmaostc committed Feb 9, 2024
2 parents c4f9521 + 84ed86c commit 65165e9
Show file tree
Hide file tree
Showing 25 changed files with 121 additions and 121 deletions.
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "RadonKA"
uuid = "86de8297-835b-47df-b249-c04e8db91db5"
authors = ["Felix Wechsler (roflmaostc) <fxw+git@mailbox.org>"]
version = "0.3.2"
version = "0.3.3"

[deps]
Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458"
Expand Down
20 changes: 10 additions & 10 deletions README.md
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# RadonKA.jl
A simple yet sufficiently fast Radon and inverse Radon (iradon) transform implementation using [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl).
A simple yet sufficiently fast Radon and inverse Radon (backproject) transform implementation using [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl).
It offers multithreading and CUDA support and outperforms any existing Julia Radon transforms (at least the ones we are aware of).
On CUDA it is faster much than Matlab and it offers the same or faster speed than ASTRA.

Expand All @@ -14,13 +14,13 @@ On CUDA it is faster much than Matlab and it offers the same or faster speed tha

# Quick Overview
* [x] For 2D and 3D arrays
* [x] parallel `radon` and `iradon` (`?RadonParallelCircle`)
* [x] attenuated `radon` and `iradon` (see the parameter `μ`) and see this [paper](https://iopscience.iop.org/article/10.1088/0266-5611/17/1/309/meta) as reference)
* [x] parallel `radon` and `backproject` (`?RadonParallelCircle`)
* [x] attenuated `radon` and `backproject` (see the parameter `μ`) and see this [paper](https://iopscience.iop.org/article/10.1088/0266-5611/17/1/309/meta) as reference)
* [x] arbitrary 2D geometries where starting and endpoint of each ray can be specified (fan beam could be a special case if this) (`?RadonFlexibleCircle`)
* [x] different strength weighting of rays
* [x] based on [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl)
* [x] tested on `CPU()` and `CUDABackend()`
* [x] registered adjoint rules for both `radon` and `iradon`
* [x] registered adjoint rules for both `radon` and `backproject`
* [x] high performance however not ultra high performance. On par with ASTRA, on CUDA faster than Matlab.
* [x] simple and extensible API

Expand All @@ -40,7 +40,7 @@ angles = range(0f0, 2f0π, 500)[begin:end-1]
# 0.085398 seconds (260 allocations: 1.006 MiB)
@time sinogram = radon(img, angles);
# 0.127043 seconds (251 allocations: 1.036 MiB)
@time backproject = RadonKA.iradon(sinogram, angles);
@time backproject = RadonKA.backproject(sinogram, angles);

simshow(sinogram)
simshow(backproject)
Expand All @@ -50,10 +50,10 @@ img_c = CuArray(img)
# 0.003363 seconds (244 CPU allocations: 18.047 KiB) (7 GPU allocations: 1007.934 KiB, 0.96% memmgmt time)
CUDA.@time sinogram = radon(img_c, angles);
# 0.005928 seconds (218 CPU allocations: 16.109 KiB) (7 GPU allocations: 1.012 MiB, 0.49% memmgmt time)
CUDA.@time backproject = RadonKA.iradon(sinogram, angles);
CUDA.@time backproject = RadonKA.backproject(sinogram, angles);
```
<a href="docs/src/assets/sinogram.png"><img src="docs/src/assets/sinogram.png" width="300"></a>
<a href="docs/src/assets/radonka_iradon.png"><img src="docs/src/assets/radonka_iradon.png" width="308"></a>
<a href="docs/src/assets/radonka_backproject.png"><img src="docs/src/assets/radonka_backproject.png" width="308"></a>

# Examples
See the [documentation](https://roflmaostc.github.io/RadonKA.jl/dev/tutorial).
Expand All @@ -69,8 +69,8 @@ julia> using Pluto; Pluto.run()
```
A browser should open.
The following examples show case the ability of this package:
* Simple `radon` and `iradon`: [Pluto notebook](examples/0_example_radon_iradon.jl)
* Different geometries: [Pluto notebook](examples/0_example_radon_iradon.jl)
* Simple `radon` and `backproject`: [Pluto notebook](examples/0_example_radon_backproject.jl)
* Different geometries: [Pluto notebook](examples/0_example_radon_backproject.jl)
* Reconstruction of a CT dataset with an optimizer: [Pluto notebook](examples/2_CT_with_optimizer.jl)
* How this package is used in Tomographic Volumetric Additive Manufacturing (3D printing): [Pluto notebook](examples/4_Tomographic_Volumetric_Additive_Manufacturing_with_Refraction.jl)

Expand All @@ -93,4 +93,4 @@ There exists [Sinograms.jl](https://github.com/JuliaImageRecon/Sinograms.jl) and
Again, no arbitrary geometries can be specified. And also no attenuated Radon transform is possible.

## Matlab
Matlab has built-in a `radon` and `iradon` transform which is similar to our lightweight API. However, no CUDA acceleration, no 3D arrays and no attenuated Radon transform.
Matlab has built-in a `radon` and `backproject` transform which is similar to our lightweight API. However, no CUDA acceleration, no 3D arrays and no attenuated Radon transform.
File renamed without changes
File renamed without changes
File renamed without changes
File renamed without changes
20 changes: 10 additions & 10 deletions docs/src/benchmark.md
Original file line number Diff line number Diff line change
Expand Up @@ -30,26 +30,26 @@ Tested on a AMD Ryzen 9 5900X 12-Core Processor with 24 Threads and a NVIDIA GeF
sinogram = radon(sample_2D, angles);
@btime sinogram = radon($sample_2D, $angles);
simshow(sinogram)
@btime sample_iradon = iradon($sinogram, $angles);
@btime sample_backproject = backproject($sinogram, $angles);

@btime CUDA.@sync sinogram_c = radon($sample_2D_c, $angles_c);
sinogram_c = radon(sample_2D_c, angles_c);
@btime CUDA.@sync sample_iradon_c = iradon($sinogram_c, $angles_c);
@btime CUDA.@sync sample_backproject_c = backproject($sinogram_c, $angles_c);


sample_3D = randn(Float32, (512, 512, 100));
sample_3D_c = CuArray(sample_3D)

sinogram_3D = radon(sample_3D, angles);
@btime radon($sample_3D, $angles)
@btime iradon($sinogram_3D, $angles)
@btime backproject($sinogram_3D, $angles)

sinogram_3D_c = radon(sample_3D_c, angles_c)
@btime CUDA.@sync radon($sample_3D_c, $angles_c)
@btime CUDA.@sync iradon($sinogram_3D_c, $angles_c)
@btime CUDA.@sync backproject($sinogram_3D_c, $angles_c)
```
![](../assets/radonka_sinogram.png)
![](../assets/radonka_iradon.png)
![](../assets/radonka_backproject.png)


## Matlab (R2023a)
Expand All @@ -65,10 +65,10 @@ R = R / max(R(:));
imwrite(R, "/tmp/matlab_sinogram.png");
tic;
iR = iradon(R, theta, "linear", "none");
iR = backproject(R, theta, "linear", "none");
toc;
iR = iR / max(iR(:));
imwrite(iR, "/tmp/matlab_iradon.png");
imwrite(iR, "/tmp/matlab_backproject.png");
Expand All @@ -87,13 +87,13 @@ toc;
tic;
for i = 1:size(y, 1)
ix(i, :, :) = iradon(squeeze(y(i, :, :)), theta);
ix(i, :, :) = backproject(squeeze(y(i, :, :)), theta);
end
toc;
```

![](../assets/matlab_sinogram.png)
![](../assets/matlab_iradon.png)
![](../assets/matlab_backproject.png)

## Astra
Some of the benchmarks did not run properly or were providing non-sense.
Expand Down Expand Up @@ -174,5 +174,5 @@ np.copy(rec)
```

![](../assets/astra_sinogram.png)
![](../assets/astra_iradon.png)
![](../assets/astra_backproject.png)

2 changes: 1 addition & 1 deletion docs/src/functions.md
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ radon

## IRadon
```@docs
iradon
backproject
filtered_backprojection
```

Expand Down
12 changes: 6 additions & 6 deletions docs/src/geometries.md
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ sinogram[1:5:end] .= 1

geometry_parallel = RadonParallelCircle(N, -(N-1)÷2:(N-1)÷2)

projection_parallel = iradon(sinogram, angles; geometry=geometry_parallel);
projection_parallel = backproject(sinogram, angles; geometry=geometry_parallel);

simshow(projection_parallel)
```
Expand All @@ -37,7 +37,7 @@ sinogram_small[1:3:end] .= 1

geometry_small = RadonParallelCircle(200, -49:49)

projection_small = iradon(sinogram_small, angles; geometry=geometry_small);
projection_small = backproject(sinogram_small, angles; geometry=geometry_small);

simshow(projection_small)
```
Expand All @@ -54,7 +54,7 @@ The second range indicates the position upon exit of the circle.
```julia
geometry_fan = RadonFlexibleCircle(N, -(N-1)÷2:(N-1)÷2, range(-(N-1)÷4, (N-1)÷4, N-1))

projected_fan = iradon(sinogram, angles; geometry=geometry_fan);
projected_fan = backproject(sinogram, angles; geometry=geometry_fan);

simshow(projected_fan, γ=0.01)
```
Expand All @@ -63,7 +63,7 @@ simshow(projected_fan, γ=0.01)
```julia
geometry_extreme = RadonFlexibleCircle(N, -(N-1)÷2:(N-1)÷2, zeros((199,)))

projected_extreme = iradon(sinogram, angles; geometry=geometry_extreme);
projected_extreme = backproject(sinogram, angles; geometry=geometry_extreme);

simshow(projected_extreme, γ=0.01)
```
Expand All @@ -74,7 +74,7 @@ For example, if in your application some rays are stronger than others you can i

```julia
geometry_weight = RadonParallelCircle(N, -(N-1)÷2:(N-1)÷2, abs.(-(N-1)÷2:(N-1)÷2))
projection_weight = iradon(sinogram, angles; geometry=geometry_weight);
projection_weight = backproject(sinogram, angles; geometry=geometry_weight);

simshow(projection_weight)
```
Expand All @@ -85,7 +85,7 @@ simshow(projection_weight)
The ray gets some attenuation with `exp(-μ*x)` where `x` is the distance traveled to the entry point of the circle. `μ` is in units of pixel.

```julia
projected_exp = iradon(sinogram, angles; geometry=geometry_extreme, μ=0.04);
projected_exp = backproject(sinogram, angles; geometry=geometry_extreme, μ=0.04);

simshow(projected_exp)
```
Expand Down
14 changes: 7 additions & 7 deletions docs/src/index.md
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# RadonKA.jl
A simple yet sufficiently fast Radon and inverse Radon (iradon) transform implementation using [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl).
A simple yet sufficiently fast Radon and inverse Radon (backproject) transform implementation using [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl).

```@raw html
<a href="assets/RadonKA_logo.png"><img src="assets/RadonKA_logo.png" width="200"></a>
Expand All @@ -11,13 +11,13 @@ A simple yet sufficiently fast Radon and inverse Radon (iradon) transform implem

# Quick Overview
* [x] For 2D and 3D arrays
* [x] parallel `radon` and `iradon` (`?RadonParallelCircle`)
* [x] attenuated `radon` and `iradon` (see the parameter `μ`)
* [x] parallel `radon` and `backproject` (`?RadonParallelCircle`)
* [x] attenuated `radon` and `backproject` (see the parameter `μ`)
* [x] arbitrary 2D geometries where starting and endpoint of each ray can be specified (fan beam could be a special case if this) (`?RadonFlexibleCircle`)
* [x] It is restricted to the incircle of radius `N ÷ 2 - 1` if the array has size `(N, N, N_z)`
* [x] based on [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl)
* [x] tested on `CPU()` and `CUDABackend`
* [x] registered adjoint rules for both `radon` and `iradon`
* [x] registered adjoint rules for both `radon` and `backproject`
* [x] high performance however not ultra high performance
* [x] simple API

Expand All @@ -39,20 +39,20 @@ angles = range(0f0, 2f0π, 500)[begin:end-1]
@time sinogram = radon(img, angles);

# 0.268649 seconds (147 allocations: 1.015 MiB)
@time backproject = RadonKA.iradon(sinogram, angles);
@time backproject = RadonKA.backproject(sinogram, angles);

simshow(sinogram)
simshow(backproject)
```

```@raw html
<a href="assets/sinogram.png"><img src="assets/sinogram.png" width="300"></a>
<a href="assets/radonka_iradon.png"><img src="assets/radonka_iradon.png" width="308"></a>
<a href="assets/radonka_backproject.png"><img src="assets/radonka_backproject.png" width="308"></a>
```

# Examples
See either the [documentation](https://roflmaostc.github.io/RadonKA.jl/dev/tutorial).
Otherwise, this [example](https://github.com/roflmaostc/RadonKA.jl/blob/main/examples/example_radon_iradon.jl) shows the main features, including CUDA support.
Otherwise, this [example](https://github.com/roflmaostc/RadonKA.jl/blob/main/examples/example_radon_backproject.jl) shows the main features, including CUDA support.
There is one tutorial about [Gradient Descent optimization](https://github.com/roflmaostc/RadonKA.jl/blob/main/examples/CT_with_optimizer.jl).
Another one covers how the Radon transform is used in [Volumetric Additive Manufacturing](https://github.com/roflmaostc/RadonKA.jl/blob/main/examples/volumetric_printing.jl).

Expand Down
6 changes: 3 additions & 3 deletions docs/src/tutorial.md
Original file line number Diff line number Diff line change
Expand Up @@ -24,18 +24,18 @@ angles = range(0f0, 2f0π, 500)[begin:end-1]
```
![](../assets/sinogram.png)

# inverse Radon (iradon) transform
# inverse Radon (backproject) transform
```julia
# 0.268649 seconds (147 allocations: 1.015 MiB)
@time backproject = RadonKA.iradon(sinogram, angles);
@time backproject = RadonKA.backproject(sinogram, angles);

# in Pluto or Jupyter
simshow(sinogram)

[simshow(img) simshow(backproject)]
```
Left is the original sample and right the simple backprojection.
![](../assets/sample_iradon.png)
![](../assets/sample_backproject.png)


# Filtered Backprojection
Expand Down
6 changes: 3 additions & 3 deletions examples/0_example_radon_iradon.jl
Original file line number Diff line number Diff line change
Expand Up @@ -103,16 +103,16 @@ simshow(Array(sinogram_c[:, :, i_z]))
md"# IRadon Transform"

# ╔═╡ 7e27da5a-1b04-4d4c-8c62-eaffa7f4f9ce
@time backproject = RadonKA.iradon(sinogram, angles);
@time backproject = RadonKA.backproject(sinogram, angles);

# ╔═╡ 4e367035-eb2f-4dfa-9646-7c182a111c49
md"Use this slider to add more and more angles to the iradon transform"
md"Use this slider to add more and more angles to the backproject transform"

# ╔═╡ b9bf49a0-7320-4269-9a6a-ac2533ab5fde
@bind angle_limit Slider(1:size(sinogram, 2), default=size(sinogram, 2), show_value=true)

# ╔═╡ 93a7ab4a-b2dc-4fc2-bf69-66e6d615103f
CUDA.@time CUDA.@sync backproject_cu = RadonKA.iradon(sinogram_c[:, begin:angle_limit, :], togoc(angles[begin:angle_limit]));
CUDA.@time CUDA.@sync backproject_cu = RadonKA.backproject(sinogram_c[:, begin:angle_limit, :], togoc(angles[begin:angle_limit]));

# ╔═╡ 52a86ed8-4504-4d9e-9ea6-6aeaf8540406
@bind i_z3 Slider(1:size(sinogram, 3), show_value=true)
Expand Down
16 changes: 8 additions & 8 deletions examples/1_documentation_different_geometries.jl
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ end
geometry_parallel = RadonParallelCircle(N, -(N-1)÷2:(N-1)÷2)

# ╔═╡ 07bec511-7802-45e4-b331-0d35e8f850c1
projection_parallel = iradon(sinogram, angles; geometry=geometry_parallel);
projection_parallel = backproject(sinogram, angles; geometry=geometry_parallel);

# ╔═╡ 9e5a58be-0de5-4820-bc3b-6dff931271a9
simshow(projection_parallel)
Expand All @@ -56,7 +56,7 @@ end
geometry_small = RadonParallelCircle(200, -49:49)

# ╔═╡ 27d5aac8-2528-46f2-a0dd-91c3fe5d275c
projection_small = iradon(sinogram_small, angles; geometry=geometry_small);
projection_small = backproject(sinogram_small, angles; geometry=geometry_small);

# ╔═╡ 29c9d1ff-694c-4fff-b551-836cc9eaf347
simshow(projection_small)
Expand All @@ -68,7 +68,7 @@ md"# Similar to fan Beam Tomography"
geometry_fan = RadonFlexibleCircle(N, -(N-1)÷2:(N-1)÷2, range(-(N-1)÷4, (N-1)÷4, N-1))

# ╔═╡ 37d760fa-74e6-47d1-b8e6-3315c1747b4c
projected_fan = iradon(sinogram, angles; geometry=geometry_fan);
projected_fan = backproject(sinogram, angles; geometry=geometry_fan);

# ╔═╡ 878121c5-4ad5-477b-88c7-f53df7510052
simshow(projected_fan, γ=0.01)
Expand All @@ -80,7 +80,7 @@ md"# Extreme fan Beam Tomography"
geometry_extreme = RadonFlexibleCircle(N, -(N-1)÷2:(N-1)÷2, zeros((199,)))

# ╔═╡ 0756e11b-35ca-4eef-be48-52b8b5c098cd
projected_extreme = iradon(sinogram, angles; geometry=geometry_extreme);
projected_extreme = backproject(sinogram, angles; geometry=geometry_extreme);

# ╔═╡ 1e22032d-d57d-42d3-a9c5-9f2e94531346
simshow(projected_extreme, γ=0.01)
Expand All @@ -94,7 +94,7 @@ For example, if in your application some rays are stronger than others you can i
geometry_weight = RadonParallelCircle(N, -(N-1)÷2:(N-1)÷2, abs.(-(N-1)÷2:(N-1)÷2))

# ╔═╡ 45b82171-c581-4fcf-9695-bc51464a2172
projection_weight = iradon(sinogram, angles; geometry=geometry_weight);
projection_weight = backproject(sinogram, angles; geometry=geometry_weight);

# ╔═╡ bbb51355-fa72-4064-b10e-46e29f4b2809
simshow(projection_weight)
Expand All @@ -105,7 +105,7 @@ The ray gets some attenuation with `exp(-μ*x)` where `x` is the distance travel
"

# ╔═╡ b11cb860-3fab-4fc9-8921-bef32a8b7c12
projected_exp = iradon(sinogram, angles; geometry=geometry_extreme, μ=0.04);
projected_exp = backproject(sinogram, angles; geometry=geometry_extreme, μ=0.04);

# ╔═╡ 51561641-67a3-46d4-ae87-b896dc351a60
simshow(projected_exp)
Expand Down Expand Up @@ -135,10 +135,10 @@ geometry_extreme2 = RadonFlexibleCircle(N2, -(N2-1)÷2:(N2-1)÷2, zeros((N2-1,))
simshow(sg_img)

# ╔═╡ 29c25bd9-6359-489d-b0a3-ea767f832a02
img_iradon = iradon(sg_img, angles2, geometry=geometry_extreme2, μ=0.008);
img_backproject = backproject(sg_img, angles2, geometry=geometry_extreme2, μ=0.008);

# ╔═╡ 940d3b8b-1797-45f1-b493-02592a39eb19
simshow(img_iradon)
simshow(img_backproject)

# ╔═╡ c7dcc341-0afe-4913-b0a8-1869d8f819b4
md"# How to reconstruct?
Expand Down
4 changes: 2 additions & 2 deletions examples/2_CT_with_optimizer.jl
Original file line number Diff line number Diff line change
Expand Up @@ -73,10 +73,10 @@ md"""# Simple Backprojection"""
img_bp = filtered_backprojection(measurement, angles);

# ╔═╡ 49e59001-0e8d-4872-ae50-47c38486b3fd
img_iradon = iradon(measurement, angles);
img_backproject = backproject(measurement, angles);

# ╔═╡ 6fac5606-2350-4f22-9f35-120936114d5a
[simshow(img_bp) simshow(img_iradon)]
[simshow(img_bp) simshow(img_backproject)]

# ╔═╡ 1feccdec-cc35-4cc8-9a76-0e0e99bc7be3
md"# Optimization with gradient descent
Expand Down

2 comments on commit 65165e9

@roflmaostc
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

Release notes:

  • rename iradon to backproject. Added a deprecation notice.

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/100571

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.3.3 -m "<description of version>" 65165e968ff76ac339beac841b4d39ca3fd7aed2
git push origin v0.3.3

Please sign in to comment.