Improve implementation of has_ghost_vertices and is_boundary_node

.gitignore

@@ -11,3 +11,4 @@ LocalPreferences.toml
 docs/src/applications/cell_simulation.mp4
 varying_weight_power.mp4
 varying_weight.mp4
+*settings.local.json

NEWS.md

@@ -1,5 +1,9 @@
 # Changelog
 
+## 1.6.6
+
+- The `Triangulation` struct now stores a `has_ghosts::Bool` field and a `boundary_vertex_to_ghost` map for efficient lookup. See [#240](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/240).
+
 ## 1.6.5
 
 - Clarified the counter-clockwise requirement for the `clip_polygon` argument in `voronoi`. See [#230](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/230).

Project.toml

@@ -1,7 +1,7 @@
 name = "DelaunayTriangulation"
 uuid = "927a84f5-c5f4-47a5-9785-b46e178433df"
 authors = ["Daniel VandenHeuvel <danj.vandenheuvel@gmail.com>"]
-version = "1.6.6"
+version = "1.6.7"
 
 [deps]
 AdaptivePredicates = "35492f91-a3bd-45ad-95db-fcad7dcfedb7"

src/algorithms/triangulation/basic_operations/add_boundary_information.jl

@@ -7,6 +7,8 @@ function add_boundary_information!(tri::Triangulation)
     I = integer_type(tri)
     ghost_vertex = I(𝒢)
     bn = get_boundary_nodes(tri)
+    # Clear any existing boundary_vertex_to_ghost entries before repopulating
+    empty!(get_boundary_vertex_to_ghost(tri))
     if has_multiple_curves(tri)
         add_boundary_curve_information!(tri, bn, ghost_vertex)
     elseif has_multiple_sections(tri)
@@ -18,14 +20,18 @@ function add_boundary_information!(tri::Triangulation)
 end
 function add_boundary_node_information!(tri::Triangulation, bn, ghost_vertex)
     n_edge = num_boundary_edges(bn)
+    n_edge == 0 && return ghost_vertex
     u = get_boundary_nodes(bn, n_edge + 1)
     for j in n_edge:-1:1
         v = get_boundary_nodes(bn, j)
         add_adjacent!(tri, u, v, ghost_vertex)
         add_adjacent2vertex!(tri, ghost_vertex, u, v)
         add_neighbour!(tri, ghost_vertex, u, v)
+        add_boundary_vertex_to_ghost!(tri, u, ghost_vertex)
         u = v
     end
+    # Add the last vertex
+    add_boundary_vertex_to_ghost!(tri, u, ghost_vertex)
     return ghost_vertex
 end
 function add_boundary_segment_information!(tri::Triangulation, bn, ghost_vertex)

src/algorithms/triangulation/basic_operations/add_ghost_triangles.jl

@@ -5,15 +5,31 @@ Adds all the ghost triangles to `tri`.
 """
 function add_ghost_triangles!(tri::Triangulation)
     T = get_triangles(tri)
-    for g in each_ghost_vertex(tri)
+    # If we have boundary nodes but the boundary_vertex_to_ghost map is empty,
+    # we need to populate it. This can happen when ghost triangles were deleted
+    # and are now being re-added.
+    if has_boundary_nodes(tri) && isempty(get_boundary_vertex_to_ghost(tri))
+        add_boundary_information!(tri)
+    end
+    # Use all_ghost_vertices instead of each_ghost_vertex because each_ghost_vertex
+    # depends on has_ghost_vertices(graph) which may be false if ghost vertices
+    # were removed from the graph by clear_empty_vertices!.
+    adj2v_dict = get_adjacent2vertex(get_adjacent2vertex(tri))
+    for g in all_ghost_vertices(tri)
+        # Skip ghost vertices that don't have any edges in adjacent2vertex.
+        # This can happen when ghost_vertex_ranges has a default entry but no
+        # actual boundary data was set up.
+        haskey(adj2v_dict, g) || continue
         for e in each_edge(get_adjacent2vertex(tri, g))
             u, v = edge_vertices(e)
             add_adjacent!(tri, v, g, u)
             add_adjacent!(tri, g, u, v)
             add_adjacent2vertex!(tri, u, v, g)
             add_adjacent2vertex!(tri, v, g, u)
+            add_neighbour!(tri, g, u, v)
             add_triangle!(T, u, v, g)
         end
     end
+    set_has_ghosts!(tri, true)
     return tri
 end

src/algorithms/triangulation/basic_operations/delete_ghost_triangles.jl

@@ -5,13 +5,13 @@ Deletes all the ghost triangles from `tri`.
 
 !!! warning "Ghost vertices"
 
-    Ghost vertices are still used in the `keys` of the [`Adjacent2Vertex`](@ref) 
-    of `tri`, and are still present in the [`Graph`](@ref). If you want to delete the 
-    ghost vertex `keys` from the [`Adjacent2Vertex`](@ref), you need to use 
-    [`delete_adjacent2vertex!`](@ref). For deleting the ghost vertices from the 
-    [`Graph`](@ref), you need [`delete_ghost_vertices_from_graph!`](@ref). Additionally, 
-    edges in [`Adjacent`](@ref) can still map to ghost vertices. If you also want to delete 
-    those, you need to filter through the `values` of the [`Adjacent`](@ref) map 
+    Ghost vertices are still used in the `keys` of the [`Adjacent2Vertex`](@ref)
+    of `tri`, and are still present in the [`Graph`](@ref). If you want to delete the
+    ghost vertex `keys` from the [`Adjacent2Vertex`](@ref), you need to use
+    [`delete_adjacent2vertex!`](@ref). For deleting the ghost vertices from the
+    [`Graph`](@ref), you need [`delete_ghost_vertices_from_graph!`](@ref). Additionally,
+    edges in [`Adjacent`](@ref) can still map to ghost vertices. If you also want to delete
+    those, you need to filter through the `values` of the [`Adjacent`](@ref) map
     that are ghost vertices, and use [`delete_adjacent!`](@ref).
 """
 function delete_ghost_triangles!(tri::Triangulation)
@@ -26,5 +26,7 @@ function delete_ghost_triangles!(tri::Triangulation)
             delete_triangle!(T, u, v, g)
         end
     end
+    empty!(get_boundary_vertex_to_ghost(tri))
+    set_has_ghosts!(tri, false)
     return tri
 end

src/algorithms/triangulation/basic_operations/lock_convex_hull.jl

@@ -30,16 +30,24 @@ function lock_convex_hull!(tri::Triangulation; rng::Random.AbstractRNG = Random.
     interior_segments = get_interior_segments(tri)
     interior_segments_on_hull = get_interior_segments_on_hull(get_cache(tri))
     empty!(interior_segments_on_hull)
+    ghost_vertex = I(𝒢)
     for i in 1:ne
         u = get_boundary_nodes(bn, i)
         v = get_boundary_nodes(bn, i + 1)
         e = construct_edge(E, u, v)
         bnn_map[e] = (bn, i)
+        # Populate the boundary_vertex_to_ghost map to maintain consistency
+        # with the boundary_nodes. This prevents add_ghost_triangles! from
+        # needing to repopulate it later.
+        add_boundary_vertex_to_ghost!(tri, u, ghost_vertex)
         if contains_unoriented_edge(e, interior_segments)
             delete_unoriented_edge!(interior_segments, e)
             add_edge!(interior_segments_on_hull, e)
         end
     end
+    # Add the last vertex (closed boundary)
+    u = get_boundary_nodes(bn, ne + 1)
+    add_boundary_vertex_to_ghost!(tri, u, ghost_vertex)
     for e in keys(bnn_map)
         add_segment!(tri, e; rng, predicates)
     end

src/algorithms/triangulation/basic_operations/unlock_convex_hull.jl

@@ -25,6 +25,9 @@ function unlock_convex_hull!(tri::Triangulation; reconstruct = false)
     bnn_map = get_boundary_edge_map(tri)
     empty!(bn)
     bn_map[I(𝒢)] = bn
+    # Clear the boundary_vertex_to_ghost map since we're unlocking the convex hull
+    # and these are no longer constrained boundary nodes
+    empty!(get_boundary_vertex_to_ghost(tri))
     segments = get_interior_segments(tri)
     all_segments = get_all_segments(tri)
     for e in keys(bnn_map)

src/algorithms/triangulation/constrained_triangulation.jl

@@ -332,6 +332,8 @@ function remake_triangulation_with_constraints(tri::Triangulation, segments, bou
     polygon_hierarchy = get_polygon_hierarchy(tri)
     boundary_enricher = get_boundary_enricher(tri)
     cache = get_cache(tri)
+    has_ghosts_flag = has_ghosts(tri)
+    boundary_vertex_to_ghost_map = get_boundary_vertex_to_ghost(tri)
     return new_ghost_vertex_map, new_ghost_vertex_ranges, Triangulation(
             points,
             triangles,
@@ -351,6 +353,8 @@ function remake_triangulation_with_constraints(tri::Triangulation, segments, bou
             polygon_hierarchy,
             boundary_enricher,
             cache,
+            has_ghosts_flag,
+            boundary_vertex_to_ghost_map,
         )
 end
 
@@ -385,6 +389,8 @@ function replace_ghost_vertex_information(tri::Triangulation, ghost_vertex_map,
     representative_point_list = get_representative_point_list(tri)
     boundary_enricher = get_boundary_enricher(tri)
     cache = get_cache(tri)
+    has_ghosts_flag = has_ghosts(tri)
+    boundary_vertex_to_ghost_map = get_boundary_vertex_to_ghost(tri)
     return Triangulation(
         points,
         triangles,
@@ -404,6 +410,8 @@ function replace_ghost_vertex_information(tri::Triangulation, ghost_vertex_map,
         polygon_hierarchy,
         boundary_enricher,
         cache,
+        has_ghosts_flag,
+        boundary_vertex_to_ghost_map,
     )
 end
 

src/algorithms/triangulation/triangulate_convex.jl

@@ -159,27 +159,36 @@ function postprocess_triangulate_convex!(tri::Triangulation, S; delete_ghosts, d
     append!(hull, S)
     push!(hull, S[begin])
     I = integer_type(tri)
+    ghost_vertex = I(𝒢)
+    # Only populate boundary_vertex_to_ghost map for triangulations with actual constrained boundary nodes.
+    # Convex hull vertices (from unconstrained triangulations) should not be in this map.
+    # This matches the behavior of add_ghost_triangles!
+    populate_map = has_boundary_nodes(tri) && !delete_ghosts
     for i in firstindex(S):(lastindex(S) - 1)
         u = S[i]
         v = S[i + 1]
         if !delete_ghosts
-            add_triangle!(tri, v, u, I(𝒢); protect_boundary = true, update_ghost_edges = false)
+            add_triangle!(tri, v, u, ghost_vertex; protect_boundary = true, update_ghost_edges = false)
         else
             # Still want the ghost vertex in the graph, adjacent2vertex map, and the adjacent map.
-            add_neighbour!(tri, I(𝒢), u)
-            add_adjacent2vertex!(tri, I(𝒢), v, u)
-            add_adjacent!(tri, v, u, I(𝒢))
+            add_neighbour!(tri, ghost_vertex, u)
+            add_adjacent2vertex!(tri, ghost_vertex, v, u)
+            add_adjacent!(tri, v, u, ghost_vertex)
         end
+        # Add boundary vertex to ghost mapping only for constrained boundaries
+        populate_map && add_boundary_vertex_to_ghost!(tri, u, ghost_vertex)
     end
     u = S[end]
     v = S[begin]
     if !delete_ghosts
-        add_triangle!(tri, v, u, I(𝒢); protect_boundary = true, update_ghost_edges = false)
+        add_triangle!(tri, v, u, ghost_vertex; protect_boundary = true, update_ghost_edges = false)
     else
-        add_neighbour!(tri, I(𝒢), u)
-        add_adjacent2vertex!(tri, I(𝒢), v, u)
-        add_adjacent!(tri, v, u, I(𝒢))
+        add_neighbour!(tri, ghost_vertex, u)
+        add_adjacent2vertex!(tri, ghost_vertex, v, u)
+        add_adjacent!(tri, v, u, ghost_vertex)
     end
+    populate_map && add_boundary_vertex_to_ghost!(tri, u, ghost_vertex)
+    set_has_ghosts!(tri, !delete_ghosts)
     delete_empty_features && clear_empty_features!(tri)
     empty_representative_points!(tri)
     cx, cy = mean_points(get_points(tri), S)

src/algorithms/triangulation/unconstrained_triangulation.jl

@@ -83,6 +83,7 @@ function initialise_bowyer_watson!(tri::Triangulation, insertion_order, predicat
     add_triangle!(tri, v, u, g; protect_boundary = true, update_ghost_edges = false)
     add_triangle!(tri, w, v, g; protect_boundary = true, update_ghost_edges = false)
     add_triangle!(tri, u, w, g; protect_boundary = true, update_ghost_edges = false)
+    set_has_ghosts!(tri, true)
     new_representative_point!(tri, I(1))
     for i in triangle_vertices(initial_triangle)
         p = get_point(tri, i)

src/data_structures/triangulation/methods/boundary_nodes.jl

@@ -126,11 +126,22 @@ function split_boundary_edge!(tri::Triangulation, i, j, node)
     delete_unoriented_edge!(segments, construct_edge(E, i, j))
     !contains_edge(construct_edge(E, node, i), segments) && add_edge!(segments, construct_edge(E, i, node))
     !contains_edge(construct_edge(E, j, node), segments) && add_edge!(segments, construct_edge(E, node, j))
-    # Sometimes, the user might accidentally also put an interior segment in that forms part of the boundary. Let's remove it. 
+    # Sometimes, the user might accidentally also put an interior segment in that forms part of the boundary. Let's remove it.
     interior_segments = get_interior_segments(tri)
     if contains_unoriented_edge(construct_edge(E, i, j), interior_segments)
         delete_unoriented_edge!(interior_segments, construct_edge(E, i, j))
     end
+    # Update boundary_vertex_to_ghost map for the new boundary node
+    # We need to find the ghost vertex for the section being split, not just use
+    # vertex i's ghost (which may be from a different section in multiply-connected domains)
+    section_path = pos[1]
+    ghost_vertex_map = get_ghost_vertex_map(tri)
+    for (g, path) in ghost_vertex_map
+        if path == section_path
+            add_boundary_vertex_to_ghost!(tri, node, g)
+            break
+        end
+    end
     return tri
 end
 
@@ -149,6 +160,7 @@ function merge_boundary_edge!(tri::Triangulation, i, j, node)
     node_pos = (pos[1], pos[2] + 1)
     bnn = get_boundary_edge_map(tri)
     delete_boundary_node!(tri, node_pos)
+    delete_boundary_vertex_from_ghost_map!(tri, node)
     E = edge_type(tri)
     delete!(bnn, construct_edge(E, i, node))
     delete!(bnn, construct_edge(E, node, j))

src/data_structures/triangulation/methods/checks.jl

@@ -23,15 +23,10 @@ unoriented_edge_exists(tri::Triangulation, i, j) = unoriented_edge_exists(tri, c
     has_ghost_triangles(tri::Triangulation) -> Bool
 
 Returns `true` if `tri` has ghost triangles, and `false` otherwise.
+
+This is an alias for [`has_ghosts`](@ref).
 """
-function has_ghost_triangles(tri::Triangulation)
-    num_ghost_vertices(tri) == 0 && return false
-    I = integer_type(tri)
-    some_outer_boundary_edges = get_adjacent2vertex(tri, I(𝒢))
-    isempty(some_outer_boundary_edges) && return false
-    e = (first ∘ each_edge)(some_outer_boundary_edges)
-    return edge_exists(tri, terminal(e), I(𝒢))
-end
+has_ghost_triangles(tri::Triangulation) = has_ghosts(tri)
 
 """
     is_constrained(tri::Triangulation) -> Bool