From 59688900180ead256cfea06075593de9e443164f Mon Sep 17 00:00:00 2001 From: roflmaostc Date: Mon, 27 Mar 2023 17:54:15 +0200 Subject: [PATCH] Add threading for radon --- src/ImageReconstruction.jl | 31 +++++++++++++++++-------------- 1 file changed, 17 insertions(+), 14 deletions(-) diff --git a/src/ImageReconstruction.jl b/src/ImageReconstruction.jl index 755e7bb..f65d58b 100644 --- a/src/ImageReconstruction.jl +++ b/src/ImageReconstruction.jl @@ -13,30 +13,33 @@ Radon transform of a image `I` producing a sinogram from view angles `θ` in radians and detector sampling `t`. """ function radon(I::AbstractMatrix, θ::AbstractRange, t::AbstractRange) - P = zeros(eltype(I), length(t), length(θ)) + Ps = [zeros(eltype(I), length(t), length(θ)) for _ = 1:nthreads()] ax1, ax2 = axes(I) nax1, nax2 = length(ax1), length(ax2) - for j in ax2, i in ax1 - x = j - nax2 / 2 + 0.5 - y = i - nax1 / 2 + 0.5 - @inbounds for (k, θₖ) in enumerate(θ) - t′ = x * cos(θₖ) + y * sin(θₖ) + Threads.@threads for j in ax2 + for i in ax1 + tid = threadid() + x = j - nax2 / 2 + 0.5 + y = i - nax1 / 2 + 0.5 + @inbounds for (k, θₖ) in enumerate(θ) + t′ = x * cos(θₖ) + y * sin(θₖ) - a = convert(Int, round((t′ - minimum(t)) / step(t) + 1)) + a = convert(Int, round((t′ - minimum(t)) / step(t) + 1)) - (a < 1 || a > length(t)) && continue - α = abs(t′ - t[a]) + (a < 1 || a > length(t)) && continue + α = abs(t′ - t[a]) - I′ = I[i, j] - P[a, k] += (1 - α) * I′ + I′ = I[i, j] + Ps[tid][a, k] += (1 - α) * I′ - (a > length(t) + 1) && continue - P[a+1, k] += α * I′ + (a > length(t) + 1) && continue + Ps[tid][a+1, k] += α * I′ + end end end - P + return sum(Ps) end function _ramp_spatial(N::Int, τ, Npad::Int = N)