Mercurial > repos > PredictPDPS / file diff

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Radon.jl	Fri Apr 19 17:00:37 2024 +0300
@@ -0,0 +1,85 @@
+# From https://github.com/LMescheder/Radon.jl/blob/master/src/Radon.jl
+# Updated to modern Julia. Fixed radon! to zero out `radonim`.
+
+__precompile__()
+
+
+module Radon
+export radon!, backproject!
+# package code goes here
+
+function radon(im::AbstractArray{T, 2}) where T<:AbstractFloat
+    Nx, Ny = size(im)
+    radonim = similar(im, min(Nx, Ny), max(Nx, Ny))
+    radon!(radonim, im)
+    radonim
+end
+
+function radon!(radonim::AbstractArray{T, 2}, im::AbstractArray{T, 2}) where T<:AbstractFloat
+    Nφ, Ns = size(radonim)
+    Nx, Ny = size(im)
+    Ldiag = hypot(Nx, Ny)
+    K = round(Int, Ldiag)
+
+    dLz = Ldiag/(K-1)
+    dLs = Ldiag/(Ns-1)
+    dφ = π/Nφ
+
+    @inbounds for is = 1:Ns, iφ = 1:Nφ
+        φ = (iφ-1)*dφ
+        s = -Ldiag/2 + (is-1)*dLs
+        tx = cos(φ)
+        ty = sin(φ)
+
+        radonim[iφ, is] = zero(T)
+        for k=1:K
+            z = -Ldiag/2 + (k-1)*dLz
+            x = round(Int, Nx/2 + z*ty + s*tx)
+            y = round(Int, Ny/2 - z*tx + s*ty)
+
+            if 1 <= x <= Nx && 1 <= y <= Ny
+                radonim[iφ, is] += im[x, y]*dLz
+            end
+        end
+    end
+    radonim
+end
+
+function backproject(radonim::AbstractArray{T, 2}) where T<:AbstractFloat
+    Nφ, Ns = size(radonim)
+    im = similar(radonim, Ns, Ns)
+    backproject!(im, radonim)
+    im
+end
+
+function backproject!(im::AbstractArray{T, 2}, radonim::AbstractArray{T, 2}) where T<:AbstractFloat
+    Nφ, Ns = size(radonim)
+    Nx, Ny = size(im)
+    Ldiag = hypot(Nx, Ny)
+    K = round(Int, Ldiag)
+
+    dLz = Ldiag/(K-1)
+    dLs = Ldiag/(Ns-1)
+    dφ = π/Nφ
+
+    im .= 0.
+    @inbounds for is = 1:Ns, iφ = 1:Nφ
+        φ = (iφ-1)*dφ
+        s = -Ldiag/2 + (is-1)*dLs
+        tx = cos(φ)
+        ty = sin(φ)
+
+        for k=1:K
+            z = -Ldiag/2 + (k-1)*dLz
+            x = round(Int, Nx/2 + z*ty + s*tx)
+            y = round(Int, Ny/2 - z*tx + s*ty)
+
+            if 1 <= x <= Nx && 1 <= y <= Ny
+                im[x, y] += radonim[iφ, is]*dLz
+            end
+        end
+    end
+    im
+end
+
+end # module