From 8a55f7ac5e493daa7ae91efde56de8e3d6527551 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20F=C3=A9votte?= Date: Fri, 11 Aug 2023 17:14:58 +0200 Subject: [PATCH] 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 | 100 ++++++++++++++++++++++++++++++--------------- test/Manifest.toml | 4 +- test/runtests.jl | 49 ++++++++++++++-------- 3 files changed, 101 insertions(+), 52 deletions(-) diff --git a/src/Eikonal.jl b/src/Eikonal.jl index bee5c0b..c53a864 100644 --- a/src/Eikonal.jl +++ b/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 diff --git a/test/Manifest.toml b/test/Manifest.toml index e07c83f..974a54d 100644 --- a/test/Manifest.toml +++ b/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" diff --git a/test/runtests.jl b/test/runtests.jl index c29f139..b704398 100644 --- a/test/runtests.jl +++ b/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,6 +40,12 @@ 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 @@ -40,15 +53,15 @@ using Test 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