Fix jump_and_march for domains with concave boundaries and for triangulations of disjoint domains

src/point_location/brute_force.jl

@@ -1,10 +1,10 @@
 """
-    brute_force_search(tri::Triangulation, q)
+    brute_force_search(tri::Triangulation, q; itr = each_triangle(tri)
 
 Returns the triangle in `tri` containing `q` using brute force search.
 """
-function brute_force_search(tri::Triangulation, q)
-    for V in each_triangle(tri)
+function brute_force_search(tri::Triangulation, q; itr = each_triangle(tri))
+    for V in itr
         cert = point_position_relative_to_triangle(tri, V, q)
         !is_outside(cert) && return V 
     end

src/point_location/jump_and_march.jl

@@ -71,7 +71,8 @@ end
         store_history::F=Val(false),
         history=nothing,
         rng::AbstractRNG=Random.default_rng(),
-        exterior_curve_index=1) where {C,F}
+        exterior_curve_index=1,
+        maxiters=2 + length(exterior_curve_index) - num_solid_vertices(tri) + num_solid_edges(tri)) where {C,F}
 
 Returns the triangle containing `q` using the jump-and-march algorithm.
 
@@ -89,6 +90,7 @@ Returns the triangle containing `q` using the jump-and-march algorithm.
 - `rng::AbstractRNG=Random.default_rng()`: The random number generator to use.
 - `exterior_curve_index=1`: The curve (or curves) corresponding to the outermost boundary.
 - `maxiters = num_triangles(tri)`: Maximum number of iterations to perform before restarting the algorithm at a new initial point. 
+- `concavity_protection=false`: When your triangulation has concave boundaries, it is possible that a ghost triangle is incorrectly classified as containing the point. By setting this to `true`, this will be protected against.
 
 !!! note 
 
@@ -110,15 +112,16 @@ Returns `V`, the triangle in `tri` containing `q`.
     when it is outside of the triangulation unless ghost triangles are present. 
 """
 function jump_and_march(tri::Triangulation{P,Ts,I}, q;
-    point_indices =each_solid_vertex(tri),
-    m = default_num_samples(num_vertices(point_indices)),
-    try_points=(), 
+    point_indices=each_solid_vertex(tri),
+    m=default_num_samples(num_vertices(point_indices)),
+    try_points=(),
     rng::AbstractRNG=Random.default_rng(),
     k=select_initial_point(tri, q; point_indices, m, try_points, rng),
     store_history::F=Val(false),
     history=nothing,
     exterior_curve_index=1,
-    maxiters=2 + length(exterior_curve_index) - num_solid_vertices(tri) + num_solid_edges(tri)) where {P,Ts,I,F}
+    maxiters=2 + length(exterior_curve_index) - num_solid_vertices(tri) + num_solid_edges(tri),
+    concavity_protection=false) where {P,Ts,I,F}
     return _jump_and_march(
         tri,
         q,
@@ -129,6 +132,8 @@ function jump_and_march(tri::Triangulation{P,Ts,I}, q;
         exterior_curve_index,
         maxiters,
         zero(maxiters),
+        concavity_protection,
+        zero(maxiters)
     )
 end
 
@@ -141,18 +146,17 @@ function _jump_and_march(
     rng::AbstractRNG=Random.default_rng(),
     exterior_curve_index=1,
     maxiters=2 + length(exterior_curve_index) - num_vertices(graph) + num_edges(graph),
-    cur_iter=0) where {P,Ts,I,E,F}
+    cur_iter=0,
+    concavity_protection=false,
+    num_restarts=0) where {P,Ts,I,E,F}
     is_bnd, bnd_idx = is_boundary_node(tri, k)
     if !(is_bnd && get_curve_index(tri, bnd_idx) ∈ exterior_curve_index) || !has_ghost_triangles(tri)
         # If k is not a boundary node, then we can rotate around the point k to find an initial triangle 
         # to start the search. Additionally, if there are no ghost triangles, we can only hope that q 
         # is inside the interior, meaning we should only search for the initial triangle here anyway.
         p, i, j, pᵢ, pⱼ = select_initial_triangle_interior_node(tri, k, q, store_history, history, rng)
         if !edge_exists(i) && !edge_exists(j) # When we find a possible infinite loop, we use i==j==DefaultAdjacentValue. Let's reinitialise. 
-            m = 2num_solid_vertices(tri)
-            point_indices = each_solid_vertex(tri)
-            k = select_initial_point(tri, q; m=2m, point_indices, rng)
-            return _jump_and_march(tri, q, k, store_history, history, rng, exterior_curve_index, maxiters, cur_iter)
+            return restart_jump_and_march(tri, q, store_history, history, rng, exterior_curve_index, maxiters, cur_iter, concavity_protection, num_restarts + 1)
         end
         if is_true(store_history)
             add_triangle!(history, j, i, k)
@@ -168,7 +172,12 @@ function _jump_and_march(
                 return construct_triangle(triangle_type(Ts), u, v, w)
             else
                 u, v = exterior_jump_and_march(tri, u, q, bnd_idx)
-                return construct_triangle(triangle_type(Ts), u, v, get_adjacent(tri, u, v)) # Can't just use I(BoundaryIndex) here since there could be multiple - just use get_adjacent
+                V = construct_triangle(triangle_type(Ts), u, v, get_adjacent(tri, u, v)) # Can't just use I(BoundaryIndex) here since there could be multiple - just use get_adjacent
+                if _concavity_protection_check(concavity_protection, tri, V, q)
+                    return restart_jump_and_march(tri, q, store_history, history, rng, exterior_curve_index, maxiters, cur_iter, concavity_protection, num_restarts + 1)
+                else
+                    return V
+                end
             end
         end
         # If we did not find anything from the neighbouring boundary edges, we can search the neighbouring interior edges
@@ -178,17 +187,34 @@ function _jump_and_march(
             return construct_triangle(triangle_type(Ts), i, j, k)
         elseif is_none(edge_cert)
             u, v = exterior_jump_and_march(tri, k, q, bnd_idx)
-            return construct_triangle(triangle_type(Ts), u, v, get_adjacent(tri, u, v))
+            V = construct_triangle(triangle_type(Ts), u, v, get_adjacent(tri, u, v))
+            if _concavity_protection_check(concavity_protection, tri, V, q)
+                return restart_jump_and_march(tri, q, store_history, history, rng, exterior_curve_index, maxiters, cur_iter, concavity_protection, num_restarts + 1)
+            else
+                return V
+            end
         else
             p, pᵢ, pⱼ = get_point(tri, k, i, j)
         end
     end
     if q == p || q == pᵢ || q == pⱼ
-        orientation = triangle_orientation(pᵢ, pⱼ, p)
-        if is_positively_oriented(orientation)
-            return construct_triangle(triangle_type(Ts), i, j, k)
+        # Just return where we currently are. We do need to be careful, though: 
+        # If k was a boundary index, then one of pᵢ or pⱼ could come from the 
+        # representative point list, which could mean that q is equal to one of the 
+        # vertices, but without meaning that it is actually in that triangle. So, 
+        # we need to also check for the type of indices we have. 
+        safety_check = (q == p && !is_boundary_index(k)) ||
+                       (q == pᵢ && !is_boundary_index(i)) ||
+                       (q == pⱼ && !is_boundary_index(j))
+        if safety_check
+            orientation = triangle_orientation(pᵢ, pⱼ, p)
+            if is_positively_oriented(orientation)
+                return construct_triangle(triangle_type(Ts), i, j, k)
+            else
+                return construct_triangle(triangle_type(Ts), j, i, k)
+            end
         else
-            return construct_triangle(triangle_type(Ts), j, i, k)
+            return restart_jump_and_march(tri, q, store_history, history, rng, exterior_curve_index, maxiters, cur_iter, concavity_protection, num_restarts + 1)
         end
     end
     original_k = k
@@ -220,9 +246,14 @@ function _jump_and_march(
             # triangles there have the same orientation, so we can find them as normal.
             if has_ghost_triangles(tri)
                 i, j = exterior_jump_and_march(tri, last_changed == i ? j : i, q, k) # use last_changed to ensure we get the boundary point
-                return construct_triangle(triangle_type(Ts), i, j, k)
+                V = construct_triangle(triangle_type(Ts), i, j, get_adjacent(tri, i, j))
+                if _concavity_protection_check(concavity_protection, tri, V, q)
+                    return restart_jump_and_march(tri, q, store_history, history, rng, exterior_curve_index, maxiters, cur_iter, concavity_protection, num_restarts + 1)
+                else
+                    return V
+                end
             else
-                return _jump_and_march(tri, q, k, store_history, history, rng, exterior_curve_index, maxiters, cur_iter)
+                return _jump_and_march(tri, q, k, store_history, history, rng, exterior_curve_index, maxiters, cur_iter, concavity_protection, num_restarts + 1)
             end
         end
         # Now we can finally move forward. We use check_existence to protect against the issue mentioned above.
@@ -281,7 +312,12 @@ function _jump_and_march(
                 add_triangle!(history, i, j, k)
             end
             if !is_outside(in_cert)
-                return construct_triangle(triangle_type(Ts), i, j, k)
+                V = construct_triangle(triangle_type(Ts), i, j, k)
+                if _concavity_protection_check(concavity_protection, tri, V, q)
+                    return restart_jump_and_march(tri, q, store_history, history, rng, exterior_curve_index, maxiters, cur_iter, concavity_protection, num_restarts + 1)
+                else
+                    return V
+                end
             end
             # To decide which direction this collinear point is away from the line, we can just use the last changed variable:
             # If last_changed = i, this means that the left direction was what caused the collinear point, so make k go on the left. 
@@ -296,11 +332,7 @@ function _jump_and_march(
         end
         arrangement = triangle_orientation(pᵢ, pⱼ, q)
         if cur_iter ≥ maxiters
-            # Restart
-            m = 2num_solid_vertices(tri)
-            point_indices = each_solid_vertex(tri)
-            k = select_initial_point(tri, q; m=2m, point_indices, rng)
-            return _jump_and_march(tri, q, k, store_history, history, rng, exterior_curve_index, maxiters, zero(cur_iter))
+            return restart_jump_and_march(tri, q, store_history, history, rng, exterior_curve_index, maxiters, cur_iter, concavity_protection, num_restarts + 1)
         end
     end
     # We can finish the above loop even if q is not in the triangle, in which case pᵢpⱼq was a straight line. 
@@ -315,10 +347,49 @@ function _jump_and_march(
             add_right_vertex!(history, j)
         end
         if is_outside(in_cert)
-            return _jump_and_march(tri, q, last_changed == I(DefaultAdjacentValue) ? i : last_changed, store_history, history, rng, exterior_curve_index, maxiters, cur_iter)
+            return _jump_and_march(tri, q, last_changed == I(DefaultAdjacentValue) ? i : last_changed, store_history, history, rng, exterior_curve_index, maxiters, cur_iter, concavity_protection, num_restarts + 1)
         end
     end
     # Swap the orientation to get a positively oriented triangle, remembering that we kept pᵢ on the left of pq and pⱼ on the right 
     k = get_adjacent(tri, j, i)
-    return construct_triangle(triangle_type(Ts), j, i, k)
-end
+    V = construct_triangle(triangle_type(Ts), j, i, k)
+    if _concavity_protection_check(concavity_protection, tri, V, q)
+        return restart_jump_and_march(tri, q, store_history, history, rng, exterior_curve_index, maxiters, cur_iter, concavity_protection, num_restarts + 1)
+    end
+    return V
+end
+
+function _concavity_protection_check(concavity_protection, tri, V, q)
+    !concavity_protection && return false
+    points = get_points(tri)
+    bn = get_boundary_nodes(tri)
+    δ = distance_to_polygon(q, points, bn)
+    cert = point_position_relative_to_triangle(tri, V, q)
+    if is_outside(cert)
+        return true # need to restart
+    end
+    is_ghost = is_ghost_triangle(V)
+    need_to_restart = (is_ghost && δ > 0.0) || (!is_ghost && δ < 0.0)
+    return need_to_restart
+end
+
+const RESTART_LIMIT = 25
+function restart_jump_and_march(tri, q, store_history, history, rng, exterior_curve_index, maxiters, cur_iter, concavity_protection, num_restarts)
+    if num_restarts < RESTART_LIMIT
+        m = num_solid_vertices(tri)
+        point_indices = each_solid_vertex(tri)
+        k = select_initial_point(tri, q; m=(m ÷ 2) + 1, point_indices, rng) # don't want to try all points, still want to give the algorithm a chance
+        return _jump_and_march(tri, q, k, store_history, history, rng, exterior_curve_index, maxiters, zero(cur_iter), concavity_protection, num_restarts)
+    else
+        V = brute_force_search(tri, q)
+        V_is_bad = _concavity_protection_check(concavity_protection, tri, V, q)
+        if V_is_bad
+            if is_ghost_triangle(V)
+                V = brute_force_search(tri, q; itr=each_solid_triangle(tri))
+            else
+                V = brute_force_search(tri, q; itr=each_ghost_triangle(tri))
+            end
+        end
+        return V
+    end
+end

test/point_location/jump_and_march.jl

@@ -285,9 +285,129 @@ end
     i, j, k = 2, 25, 8
     T = (i, j, k)
     r = 5
-    for i in 1:100
-        @show i
-        V = jump_and_march(tri, get_point(tri, k); point_indices=nothing, m=nothing, try_points=nothing, k=r)
+    for i in 1:10000
+        V = jump_and_march(tri, get_point(tri, k); point_indices=nothing, m=nothing, try_points=nothing, k=r, concavity_protection=true)
         @test !DT.is_outside(DT.point_position_relative_to_triangle(tri, V, get_point(tri, k)))
     end
+end
+
+@testset "Finding points in a triangulation with concave boundaries" begin
+    a, b, c = (0.0, 8.0), (0.0, 6.0), (0.0, 4.0)
+    d, e, f = (0.0, 2.0), (0.0, 0.0), (2.0, 0.0)
+    g, h, i = (4.0, 0.0), (6.0, 0.0), (8.0, 0.0)
+    j, k, ℓ = (8.0, 1.0), (7.0, 2.0), (5.0, 2.0)
+    m, n, o = (3.0, 2.0), (2.0, 3.0), (2.0, 5.0)
+    p, q, r = (2.0, 7.0), (1.0, 8.0), (1.0, 2.2)
+    s, t, u = (0.4, 1.4), (1.2, 1.8), (2.8, 0.6)
+    v, w, z = (3.4, 1.2), (1.6, 1.4), (1.6, 2.2)
+    outer = [[a, b, c, d, e], [e, f, g, h, i, j, k, ℓ], [ℓ, m, n, o, p, q, a]]
+    inner = [[r, z, v, u, w, t, s, r]]
+    boundary_nodes, points = convert_boundary_points_to_indices([outer, inner])
+    rng = StableRNG(123)
+    tri = triangulate(points; rng, boundary_nodes, delete_ghosts=false)
+    refine!(tri; max_area=0.01get_total_area(tri), rng)
+    qs = [
+        (4.0, 5.0), (1.0, 5.6), (0.2, 5.0),
+        (0.0, -1.0), (0.5, 3.5), (2.5, 1.5),
+        (1.0, 2.0), (4.5, 1.0), (6.0, 1.5),
+        (0.5, 8.5), (1.0, 7.5), (1.2, 1.6)
+    ]
+    δs = [DelaunayTriangulation.distance_to_polygon(q, get_points(tri), get_boundary_nodes(tri)) for q in qs]
+    for i in 1:1000
+        @show i
+        Vs = [jump_and_march(tri, q, concavity_protection=true, rng=StableRNG(i + j)) for (j, q) in enumerate(qs)]
+        for (q, δ, V) in zip(qs, δs, Vs)
+            cert = DelaunayTriangulation.point_position_relative_to_triangle(tri, V, q)
+            if δ ≥ 0.0
+                @test !DelaunayTriangulation.is_outside(cert)
+                @test !DelaunayTriangulation.is_ghost_triangle(V)
+            else
+                @test !DelaunayTriangulation.is_outside(cert)
+                @test DelaunayTriangulation.is_ghost_triangle(V)
+            end
+        end
+    end
+    a, b, c, d = -1, 10, -1, 10
+    for i in 1:1000
+        rng = StableRNG(i)
+        qs = [((b - a) * rand(rng) + a, (d - c) * rand(rng) + c) for _ in 1:1000]
+        δs = [DelaunayTriangulation.distance_to_polygon(q, get_points(tri), get_boundary_nodes(tri)) for q in qs]
+        Vs = [jump_and_march(tri, q, concavity_protection=true, rng=StableRNG(i + j)) for (j, q) in enumerate(qs)]
+        for (q, δ, V) in zip(qs, δs, Vs)
+            cert = DelaunayTriangulation.point_position_relative_to_triangle(tri, V, q)
+            if δ ≥ 0.0
+                @test !DelaunayTriangulation.is_outside(cert)
+                @test !DelaunayTriangulation.is_ghost_triangle(V)
+            else
+                @test !DelaunayTriangulation.is_outside(cert)
+                @test DelaunayTriangulation.is_ghost_triangle(V)
+            end
+        end
+    end
+end
+
+@testset "Finding points in a triangulation with concave boundaries and disjoint domains" begin
+    a, b, c = (0.0, 8.0), (0.0, 6.0), (0.0, 4.0)
+    d, e, f = (0.0, 2.0), (0.0, 0.0), (2.0, 0.0)
+    g, h, i = (4.0, 0.0), (6.0, 0.0), (8.0, 0.0)
+    j, k, ℓ = (8.0, 1.0), (7.0, 2.0), (5.0, 2.0)
+    m, n, o = (3.0, 2.0), (2.0, 3.0), (2.0, 5.0)
+    p, q, r = (2.0, 7.0), (1.0, 8.0), (1.0, 2.2)
+    s, t, u = (0.4, 1.4), (1.2, 1.8), (2.8, 0.6)
+    v, w, z = (3.4, 1.2), (1.6, 1.4), (1.6, 2.2)
+    outer = [[a, b, c, d, e], [e, f, g, h, i, j, k, ℓ], [ℓ, m, n, o, p, q, a]]
+    inner = [[r, z, v, u, w, t, s, r]]
+    m₁, n₁, o₁ = (6.0, 8.0), (8.0, 8.0), (8.0, 4.0)
+    p₁, q₁, r₁ = (10.0, 4.0), (6.0, 6.0), (8.0, 6.0)
+    s₁, t₁, u₁ = (9.0, 7.0), (4.0, 4.0), (5.0, 4.0)
+    v₁, w₁ = (5.0, 3.0), (4.0, 3.0)
+    new_domain₁ = [[m₁, q₁, o₁, p₁, r₁, s₁, n₁, m₁]]
+    new_domain₂ = [[t₁, w₁, v₁, u₁, t₁]]
+    boundary_nodes, points = convert_boundary_points_to_indices([outer, inner, new_domain₁, new_domain₂])
+    rng = StableRNG(125123)
+    tri = triangulate(points; rng, boundary_nodes, check_arguments=false, delete_ghosts=false)
+    refine!(tri; max_area=0.001get_total_area(tri), rng)
+    qs = [
+        (0.6, 6.4), (1.4, 0.8), (3.1, 2.9),
+        (6.3, 4.9), (4.6, 3.5), (7.0, 7.0),
+        (8.9, 5.1), (5.8, 0.8), (1.0, 1.5),
+        (1.5, 2.0), (8.15, 6.0)
+    ]
+    q = (7.0, 7.0) # When you have a point that is exactly the same as a representative point, you need to be careful of (1) the point being exactly collinear with a ghost edge, and (2) the algorithm mistakenly regarding q as if it were equal to a triangle's vertices (since this is where the ghost vertices map). This is now fixed, but we need this isolated to avoid regressions.
+    V = jump_and_march(tri, q, rng=StableRNG(268), concavity_protection=true)
+    @test DelaunayTriangulation.compare_triangles(V, (377, 378, 68))
+    δs = [DelaunayTriangulation.distance_to_polygon(q, get_points(tri), get_boundary_nodes(tri)) for q in qs]
+    for i in 1:1000
+        @show i
+        Vs = [jump_and_march(tri, q; concavity_protection=true, rng=StableRNG(i)) for q in qs]
+        for (q, δ, V) in zip(qs, δs, Vs)
+            cert = DelaunayTriangulation.point_position_relative_to_triangle(tri, V, q)
+            if δ > 0.0
+                @test !DelaunayTriangulation.is_outside(cert)
+                @test !DelaunayTriangulation.is_ghost_triangle(V)
+            elseif δ < 0.0
+                @test !DelaunayTriangulation.is_outside(cert)
+                @test DelaunayTriangulation.is_ghost_triangle(V)
+            else
+                @test !DelaunayTriangulation.is_outside(cert)
+            end
+        end
+    end
+    a, b, c, d = -1, 10, -1, 10
+    for i in 1:100
+        rng = StableRNG(i)
+        qs = [((b - a) * rand(rng) + a, (d - c) * rand(rng) + c) for _ in 1:100]
+        δs = [DelaunayTriangulation.distance_to_polygon(q, get_points(tri), get_boundary_nodes(tri)) for q in qs]
+        Vs = [jump_and_march(tri, q, concavity_protection=true, rng=StableRNG(i + j)) for (j, q) in enumerate(qs)]
+        for (q, δ, V) in zip(qs, δs, Vs)
+            cert = DelaunayTriangulation.point_position_relative_to_triangle(tri, V, q)
+            if δ ≥ 0.0
+                @test !DelaunayTriangulation.is_outside(cert)
+                @test !DelaunayTriangulation.is_ghost_triangle(V)
+            else
+                @test !DelaunayTriangulation.is_outside(cert)
+                @test DelaunayTriangulation.is_ghost_triangle(V)
+            end
+        end
+    end
 end