Change types to make them more stable

Project.toml

@@ -1,7 +1,7 @@
 name = "RadonKA"
 uuid = "86de8297-835b-47df-b249-c04e8db91db5"
 authors = ["Felix Wechsler (roflmaostc) <fxw+git@mailbox.org>"]
-version = "0.5.0"
+version = "0.6.0"
 
 [deps]
 Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458"

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
 ```

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

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 

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
-