diff --git a/Project.toml b/Project.toml index 9279ccc..2e7e38e 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.5.0" +version = "0.6.0" [deps] Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458" diff --git a/README.md b/README.md index f2a543b..18d0e71 100644 --- a/README.md +++ b/README.md @@ -23,7 +23,8 @@ On CUDA it is faster much than Matlab and it offers the same or faster speed tha * [x] simple and extensible API # Installation -Requires Julia at least 1.9 +This toolbox runs with CUDA support on Linux, Windows and MacOS! +Requires at least Julia 1.9 ```julia julia> ]add RadonKA ``` diff --git a/src/geometry.jl b/src/geometry.jl index 8c6382f..afd5d8b 100644 --- a/src/geometry.jl +++ b/src/geometry.jl @@ -29,15 +29,16 @@ See [`radon`](@ref) and [`backproject`](@ref) how to apply. """ struct RadonParallelCircle{T, T2} <: RadonGeometry N::Int - in_height::AbstractVector{T} - weights::AbstractVector{T2} + in_height::T + weights::T2 function RadonParallelCircle(N, in_height) - return new{eltype(in_height),eltype(in_height)}(N, in_height, in_height .* 0 .+ 1) + weights = in_height .* 0 .+ 1 + return new{typeof(in_height), typeof(weights)}(N, in_height, weights) end function RadonParallelCircle(N, in_height, weights) - return new{eltype(in_height),eltype(weights)}(N, in_height, weights) + return new{typeof(in_height),typeof(weights)}(N, in_height, weights) end end @@ -60,14 +61,15 @@ See [`radon`](@ref) and [`backproject`](@ref) how to apply. """ struct RadonFlexibleCircle{T, T2, T3} <: RadonGeometry N::Int - in_height::AbstractVector{T} - out_height::AbstractVector{T2} - weights::AbstractVector{T3} + in_height::T + out_height::T2 + weights::T3 function RadonFlexibleCircle(N, in_height, out_height) - return new{eltype(in_height), eltype(out_height), eltype(in_height)}(N, in_height, out_height, in_height .* 0 .+ 1) + weights = in_height .* 0 .+ 1 + return new{typeof(in_height), typeof(out_height), typeof(weights)}(N, in_height, out_height, weights) end function RadonFlexibleCircle(N, in_height, out_height, weights) - return new{eltype(in_height), eltype(out_height), eltype(weights)}(N, in_height, out_height, weights) + return new{typeof(in_height), typeof(out_height), typeof(weights)}(N, in_height, out_height, weights) end end diff --git a/src/radon.jl b/src/radon.jl index a37c1f0..856aad5 100644 --- a/src/radon.jl +++ b/src/radon.jl @@ -170,7 +170,7 @@ function _radon(img::AbstractArray{T, 3}, angles_T::AbstractVector, out_height, angles, mid, radius, absorb_f, ndrange=(N_sinogram, N_angles, size(img, 3))) KernelAbstractions.synchronize(backend) - return sinogram + return sinogram::typeof(img) end @inline inside_circ(ii, jj, mid, radius) = (ii - mid)^2 + (jj - mid)^2 ≤ radius ^2 diff --git a/src/utils.jl b/src/utils.jl index 2ffb6f4..33e7fc1 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -20,46 +20,45 @@ corner of the cell. |--------|---------|--------| """ @inline function find_cell(x, y, xnew, ynew) - x = floor(Int, (x + xnew) / 2) - y = floor(Int, (y + ynew) / 2) - return x, y + x = floor(Int32, (x + xnew) / Int32(2)) + y = floor(Int32, (y + ynew) / Int32(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₀) + # 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 = (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) + + 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 @inline function next_ray_on_circle(y_dist::T, y_dist_end, mid, radius, sinα, cosα) where T - x_dist = sqrt(radius^2 - y_dist^2) + T(0.5) + x_dist = sqrt(radius^2 - y_dist^2) + T(0.5) 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 - - x_dist_end_rot = mid + cosα * x_dist_end - sinα * y_dist_end - y_dist_end_rot = mid + sinα * x_dist_end + cosα * y_dist_end - return x_dist_rot, y_dist_rot, x_dist_end_rot, y_dist_end_rot + 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 + + x_dist_end_rot = mid + cosα * x_dist_end - sinα * y_dist_end + y_dist_end_rot = mid + sinα * x_dist_end + cosα * y_dist_end + return x_dist_rot, y_dist_rot, x_dist_end_rot, y_dist_end_rot end -