Skip to content

Commit

Permalink
Add exponential Radon transform
Browse files Browse the repository at this point in the history
  • Loading branch information
@roflmaostc
roflmaostc committed Dec 3, 2023
1 parent 591e3a6 commit 6b15882
Show file tree
Hide file tree
Showing 3 changed files with 130 additions and 29 deletions.
70 changes: 62 additions & 8 deletions examples/volumetric_printing.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,8 +55,8 @@ See this paper: Rackson, Charles M., et al. "Object-space optimization of tomogr
"""

# ╔═╡ 0c1abf70-8004-460a-b008-750050f6e718
function iter!(buffer, img, θs; filtered=true, clip_sinogram=true, backend=CPU())
sinogram = radon(img, θs; backend)
function iter!(buffer, img, θs, μ; filtered=true, clip_sinogram=true, backend=CPU())
sinogram = radon(img, θs, μ; backend)
if filtered
sinogram = ramp_filter(sinogram)
end
Expand All @@ -65,14 +65,14 @@ function iter!(buffer, img, θs; filtered=true, clip_sinogram=true, backend=CPU(
sinogram .= max.(sinogram, 0)
end

img_recon = iradon(sinogram, θs; backend)
img_recon = iradon(sinogram, θs, μ; backend)

buffer .= max.(img_recon, 0)
return buffer, sinogram
end

# ╔═╡ 5ac76842-f5a3-42c0-9764-9b20f98ab1d5
function optimize(img, θs; thresholds=(0.65, 0.75), N_iter = 2, backend=CPU())
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)
Expand All @@ -87,11 +87,11 @@ function optimize(img, θs; thresholds=(0.65, 0.75), N_iter = 2, backend=CPU())
isobject = isone.(img)

buffer = copy(img)
tmp, s = iter!(buffer, guess, θs; filtered=false, clip_sinogram=true, backend)
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)
tmp, s = iter!(buffer, guess, θs, μ; backend, filtered=false, clip_sinogram=true)

guess[isobject] .+= max.(0, thresholds[2] .- tmp[isobject])
end
Expand All @@ -109,7 +109,7 @@ angles = range(0, 2π, 350)
thresholds = (0.7, 0.8)

# ╔═╡ 96f6eba3-1700-4c29-aff6-5f774d7d4b3b
@time virtual_object, patterns_optimized = optimize(img, angles; N_iter=50, thresholds)
@time virtual_object, patterns_optimized = optimize(img, angles; N_iter=100, thresholds)

# ╔═╡ 5568ca5b-2aad-4024-8cf0-d68260a4ac36
simshow(patterns_optimized)
Expand All @@ -128,7 +128,8 @@ simshow([object_printed object_printed .> thresh2 abs.((object_printed .> thresh

# ╔═╡ b5ae35db-bf7d-4c00-a069-8572af9bc1de
begin
histogram(object_printed[:], bins=(0.0:0.01:1), xlim=(0.0, 1.0), label="dose distribution", ylabel="voxel count", xlabel="normalized intensity", ylim=(0, 1000))
histogram(object_printed[img .== 0], bins=(0.0:0.01:1), xlim=(0.0, 1.0), label="dose distribution", 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", 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")
Expand All @@ -152,6 +153,45 @@ 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"

Expand All @@ -170,7 +210,10 @@ simshow(volume[:, :, 1])
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));
Expand Down Expand Up @@ -205,6 +248,17 @@ simshow(object_printed_c2 .> 0.74)[:, :, 1]
# ╠═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
Expand Down
47 changes: 39 additions & 8 deletions src/iradon.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,24 +24,54 @@ iradon(sinogram::AbstractArray{T, 3}, angles::AbstractArray{T2, 1}, μ=nothing;
iradon(sinogram, T.(angles), μ; backend)

"""
iradon(sinogram, θs; backend=CPU())
iradon(sinogram, θs, μ=nothing; backend=CPU())
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 a odd number. And `size(sinogram, 2)` has to be equal to
`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!
Works either with a `AbstractArray{T, 3}` or `AbstractArray{T, 2}`.
`θs` is a vector or range storing the angles in radians.
`backend` can be either `CPU()` for multithreaded CPU execution or
`CUDABackend()` for CUDA execution. In principle, all backends of KernelAbstractions.jl should work
but are not tested.
# 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.
# Keywords
* `backend` can be either `CPU()` for multithreaded CPU execution or `CUDABackend()` for CUDA execution. In principle, all backends of KernelAbstractions.jl should work but are not tested.
See also [`radon`](@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(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(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;
backend=CPU()) where T
Expand All @@ -51,7 +81,7 @@ function iradon(sinogram::AbstractArray{T, 3}, angles::AbstractArray{T, 1}, μ=n
if isnothing(μ)
(x, y, x_start, y_start) -> one(T)
else
(x, y, x_start, y_start) -> exp(-μ * sqrt((x - x_start)^2 + (y - y_start)^2))
(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
Expand Down Expand Up @@ -127,7 +157,8 @@ end
distance = sqrt((a_old - a) ^2 +
(b_old - b) ^2)
# cell value times distance travelled through
Atomix.@atomic img[cell_i, cell_j, i_z] += distance * sinogram[k, r, i_z] * absorption_f(a, b, a0, b0)
Atomix.@atomic img[cell_i, cell_j, i_z] +=
distance * sinogram[k, r, i_z] * absorption_f(a, b, a0, b0)
end
end

42 changes: 29 additions & 13 deletions src/radon.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,12 +2,12 @@ export radon



radon(img::AbstractArray{T, 3}, angles::AbstractArray{T2, 1}; backend=CPU()) where {T, T2} =
radon(img, T.(angles); backend)
radon(img::AbstractArray{T, 3}, angles::AbstractArray{T2, 1}, μ=nothing; backend=CPU()) where {T, T2} =
radon(img, T.(angles), μ; backend)

# handle 2D
radon(img::AbstractArray{T, 2}, angles::AbstractArray{T2, 1}; backend=CPU()) where {T, T2} =
view(radon(reshape(img, (size(img)..., 1)), angles; backend), :, :, 1)
radon(img::AbstractArray{T, 2}, angles::AbstractArray{T2, 1}, μ=nothing; backend=CPU()) where {T, T2} =
view(radon(reshape(img, (size(img)..., 1)), angles, μ; backend), :, :, 1)

"""
radon(I, θs; backend=CPU())
Expand All @@ -21,9 +21,13 @@ Works either with a `AbstractArray{T, 3}` or `AbstractArray{T, 2}`.
`θs` is a vector or range storing the angles in radians.
`backend` can be either `CPU()` for multithreaded CPU execution or
`CUDABackend()` for CUDA execution. In principle, all backends of KernelAbstractions.jl should work
but are not tested
# 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.
# Keywords
* `backend` can be either `CPU()` for multithreaded CPU execution or `CUDABackend()` for CUDA execution. In principle, all backends of KernelAbstractions.jl should work but are not tested.
Please note: the implementation is not quite optimized for cache efficiency and
Expand All @@ -44,9 +48,19 @@ julia> radon(arr, [0, π/4, π/2])
0.0 0.0 0.0
```
"""
function radon(img::AbstractArray{T, 3}, angles::AbstractArray{T, 1};
function radon(img::AbstractArray{T, 3}, angles::AbstractArray{T, 1}, μ=nothing;
backend=CPU()) 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"


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)
Expand All @@ -67,18 +81,20 @@ function radon(img::AbstractArray{T, 3}, angles::AbstractArray{T, 1};
#@show typeof(sinogram), typeof(img), typeof(y_dists), typeof(angles)
kernel! = radon_kernel!(backend)
kernel!(sinogram::AbstractArray{T}, img, y_dists, angles, mid, radius,
ndrange=(N, N_angles, size(img, 3)))
absorption_f,
ndrange=(N, N_angles, size(img, 3)))

return sinogram
end

"""
radon_kernel!(sinogram, img,
y_dists, angles, mid, radius)
y_dists, angles, mid, radius,
absorption_f)
"""
@kernel function radon_kernel!(sinogram::AbstractArray{T},
img, y_dists, angles, mid, radius) where 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
Expand All @@ -89,7 +105,7 @@ end

# 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
Expand Down Expand Up @@ -117,7 +133,7 @@ end
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]
@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

0 comments on commit 6b15882

Please sign in to comment.