N-D ray reconstruction

2 methods:
- integration of the backtracking equations (generalization of the former 2D
  method)
- jump to minimum-time neighbor (new; seems more stable in complex cases)

src/Eikonal.jl

@@ -59,7 +59,7 @@ end
 function Base.show(io::IO, fs::FastSweeping)
     (m,n) = size(fs.v)
     grid_size = join(string.(size(fs.v)), "×")
-    println(io, "FastSweeping solver on a $grid_size grid")
+    print(io, "FastSweeping solver on a $grid_size grid")
 end
 
 
@@ -73,6 +73,8 @@ end
 
 function init!(fs :: FastSweeping, source)
     fs.t[source...] = 0
+    fs.change .= Inf
+    fs
 end
 
 """
@@ -101,19 +103,18 @@ end
     c = sum(square, tᵢ) - v^2
     Δ = b^2 - 4*a*c
 
+    t = typemax(T)
     if Δ >= 0
         t = (-b + √(Δ)) / 2a
-        t <= maximum(tᵢ) && return typemax(T)
-        return t
-    else
-        # ∇t actually belongs to the boundary of the orthant
-        # => perform an N-1 dimensional update
-        t = typemax(T)
-        for tᵢ′ in subtuples(tᵢ)
-            t = min(t, update(tᵢ′, v))
-        end
-        return t
+        t <= maximum(tᵢ) && (t = typemax(T))
+    end
+
+    # Hyp: ∇t lies on the boundary of the orthant
+    # => perform N-1 dimensional updates
+    for tᵢ′ in subtuples(tᵢ)
+        t = min(t, update(tᵢ′, v))
     end
+    return t
 end
 
 @inline function update(tᵢ::NTuple{1, T}, v) where {T}
@@ -185,10 +186,6 @@ function sweep!(fs :: FastSweeping{T};
     error("Max number of iterations reached!")
 end
 
-
-
-
-
 function update(t::AbstractMatrix{T}, v, (i,j), (δᵢ, δⱼ)) where {T}
     sᵢ = δᵢ<0 ? 0 : 1   # Offset between vertex-based indices (in `t`)
     sⱼ = δⱼ<0 ? 0 : 1   # and cell-based indices (in `v`)
@@ -299,35 +296,70 @@ end
 
 
 # Trajectory reconstruction using a gradient descent algorithm
-function ray(t::AbstractMatrix{T}, pos; ρ=T(0.1)) where {T}
+ray(t::AbstractArray{T,D}, pos) where {T,D} = ray(t, pos, Integration(0.1))
+
+struct Integration
+    rho :: Float64
+end
+
+function ray(t::AbstractArray{T,D}, pos, method::Integration) where {T,D}
     # ρ is the step of the gradient descent (in terms of cell size)
     # should be O(1) but can be decreased in cases where the gradient can be very large
+    ρ = method.rho
 
     res = [pos]
 
-    (i,j) = pos               # We need to keep track of integer cell indices
-    (x,y) = eltype(t).(pos)   # as well as floating-point position
+    # We need to keep track of floating point coordinates
+    x = eltype(t).(pos)
 
     # Gradient descent
     while true
-        tn = get(t, CartesianIndex(i+1, j), 2*t[i,j])
-        ts = get(t, CartesianIndex(i-1, j), 2*t[i,j])
-        tw = get(t, CartesianIndex(i, j+1), 2*t[i,j])
-        te = get(t, CartesianIndex(i, j-1), 2*t[i,j])
-        ∇ᵢ = (tn-ts) / 2   # Centered differences
-        ∇ⱼ = (tw-te) / 2   # (inconsistent with the FSM: is it a problem?)
-        ∇  = (∇ᵢ, ∇ⱼ)
-
-        # t is not differentiable near the origin
-        # => stop a bit before arriving
-        t[i,j] <= ρ * norm(∇) && break
-
-        # Fixed step length ρ
-        x -= ρ * ∇ᵢ / norm(∇); i = round(Int, x)
-        y -= ρ * ∇ⱼ / norm(∇); j = round(Int, y)
+        I = CartesianIndex(round.(Int, x)) # Integer coordinates
+        ∇ = ntuple(Val(D)) do k
+            δ = CartesianIndex(ntuple(l -> k==l ? 1 : 0, Val(D)))
+            t1 = min(t[I + δ], 2*t[I])
+            t2 = min(t[I - δ], 2*t[I])
+            (t1-t2) / 2
+        end
+
+        t_min = ntuple(Val(D)) do k
+            δ = CartesianIndex(ntuple(l -> k==l ? 1 : 0, Val(D)))
+            min(t[I + δ], t[I - δ])
+        end
+
+        # Stop when we reached a local minimum
+        minimum(t_min) >= t[I] && break
+
+        # Update position with a fixed step length ρ
+        x = x .- ρ .* ∇ ./ norm(∇)
 
         # Only record point if it is new
-        (i,j) != last(res) && push!(res, (i,j))
+        pos = Tuple(I)
+        pos != last(res) && push!(res, pos)
+    end
+
+    res
+end
+
+struct NearestMin end
+function ray(t::AbstractArray{T,D}, pos, ::NearestMin) where {T, D}
+    I = CartesianIndex(pos)
+    res = [pos]
+
+    while true
+        # Find minimal time in a box surrounding the current position
+        box = I-oneunit(I):I+oneunit(I)
+        t_min, ibox = findmin(box) do I′
+            I′ ∈ CartesianIndices(t) ? t[I′] : Inf
+        end
+        I′ = box[ibox]
+
+        # Stop when we reached a local minimum
+        I == I′ && break
+
+        # Update current position
+        I = I′
+        push!(res, Tuple(I))
     end
 
     res

test/Manifest.toml

@@ -1,7 +1,8 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.7.1"
+julia_version = "1.9.2"
 manifest_format = "2.0"
+project_hash = "71d91126b5a1fb1020e1098d9d492de2a4438fd2"
 
 [[deps.Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
@@ -23,6 +24,7 @@ uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [[deps.SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+version = "0.7.0"
 
 [[deps.Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

test/runtests.jl

@@ -1,29 +1,36 @@
 using Eikonal
 using Test
 
+using Eikonal: NearestMin
+
+dist(x, y) = maximum(x.-y)
+
 @testset "Eikonal.jl" begin
     @testset "FastSweeping" begin
-        siz = (1000, 1200)
+        siz = (1100, 1200)
         fsm = FastSweeping(Float64, siz)
-        @test sprint(show, fsm) == "FastSweeping solver on a 1000×1200 grid\n"
+        @test sprint(show, fsm) == "FastSweeping solver on a 1100×1200 grid"
 
         fsm.v .= 1
 
-        x0, y0 = 500, 550
-        init!(fsm, (x0, y0))
+        source = 500, 550
+        init!(fsm, source)
         sweep!(fsm, verbose=true)
 
-        x1, y1 = 50, 50
-        @test fsm.t[x1, y1] ≈ sqrt((x1-x0)^2 + (y1-y0)^2)  rtol=1e-2
+        pos = 50, 50
+        @test fsm.t[pos...] ≈ sqrt(sum((pos.-source).^2))  rtol=1e-2
+
+        r = ray(fsm.t, pos)
+        @test dist(last(r), source) <= 1
 
-        r = ray(fsm.t, (x1, y1))
-        @test last(r) == (x0, y0)
+        r = ray(fsm.t, pos, NearestMin())
+        @test dist(last(r), source) <= 1
     end
 
     @testset "FastSweeping 3D" begin
         siz = (100, 110, 120)
         fsm = FastSweeping(Float64, siz)
-        @test sprint(show, fsm) == "FastSweeping solver on a 100×110×120 grid\n"
+        @test sprint(show, fsm) == "FastSweeping solver on a 100×110×120 grid"
 
         fsm.v .= 1
 
@@ -33,22 +40,28 @@ using Test
 
         pos = 5, 5, 6
         @test fsm.t[pos...] ≈ sqrt(sum((pos.-source).^2))  rtol=5e-2
+
+        r = ray(fsm.t, pos)
+        @test dist(last(r), source) <= 1
+
+        r = ray(fsm.t, pos, NearestMin())
+        @test dist(last(r), source) <= 1
     end
 
     @testset "FastMarching" begin
         m, n = 1000, 1200
         fmm = FastMarching(m, n)
         fmm.v .= 1
 
-        x0, y0 = 500, 550
-        init!(fmm, (x0, y0))
+        source = 500, 550
+        init!(fmm, source)
         march!(fmm, verbose=true)
 
-        x1, y1 = 50, 50
-        @test fmm.t[x1, y1] ≈ sqrt((x1-x0)^2 + (y1-y0)^2)  rtol=1e-2
+        pos = 50, 50
+        @test fmm.t[pos...] ≈ sqrt(sum((pos.-source).^2))  rtol=1e-2
 
-        r = ray(fmm.t, (x1, y1))
-        @test last(r) == (x0, y0)
+        r = ray(fmm.t, pos)
+        @test dist(last(r), source) <= 1
     end
 
     @testset "Utilities" begin
@@ -64,8 +77,10 @@ using Test
         @test fsm.t[goal...] ≈ 3937        rtol=1e-2
         @test fsm.t[goal...] ≈ t[goal...]  rtol=1e-2
 
-        fsm.t .= min.(fsm.t, 4100)
         r = ray(fsm.t, goal)
-        @test last(r) == entrance
+        @test dist(last(r), entrance) <= 1
+
+        r = ray(fsm.t, goal, NearestMin())
+        @test dist(last(r), entrance) <= 1
     end
 end