Fix references to duplicate points in boundary_nodes/segments within check_args

NEWS.md

@@ -1,5 +1,9 @@
 # Changelog
 
+## 1.6.5
+- (Internal) A function `set_boundary_node!` has been added for setting a specific boundary node to another inside the provided `boundary_nodes`. See [#224](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/224).
+- (Fix) Currently when checking for duplicate points, any extra points get skipped. However we did not correctly make sure those point's vertices inside `segments` or `boundary_nodes` were replaced with the first instance of the duplicate. This has been fixed. See [#224](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/224).
+
 ## 1.6.4 
 
 - An error is no longer thrown for inputs with duplicate points. Instead, a warning is thrown and any duplicates are merged into the `skip_points` keyword argument. With this, `DuplicatePointsError` has been removed. To silence the new warning, use `DelaunayTriangulation.toggle_warn_on_dupes!()`.

Project.toml

@@ -1,7 +1,7 @@
 name = "DelaunayTriangulation"
 uuid = "927a84f5-c5f4-47a5-9785-b46e178433df"
 authors = ["Daniel VandenHeuvel <[email protected]>"]
-version = "1.6.4"
+version = "1.6.5"
 
 [deps]
 AdaptivePredicates = "35492f91-a3bd-45ad-95db-fcad7dcfedb7"

docs/Project.toml

@@ -34,3 +34,6 @@ SimpleGraphs = "0.8"
 StableRNGs = "1.0"
 StatsBase = "0.34"
 Test = "1.11"
+
+[sources]
+DelaunayTriangulation = { path = ".." }

src/algorithms/triangulation/check_args.jl

@@ -1,5 +1,5 @@
 """
-    check_args(points, boundary_nodes, hierarchy::PolygonHierarchy, boundary_curves = (); skip_points = Set{Int}()) -> Bool 
+    check_args(points, boundary_nodes, segments, hierarchy::PolygonHierarchy, boundary_curves = (); skip_points = Set{Int}()) -> Bool 
 
 Check that the arguments `points` and `boundary_nodes` to [`triangulate`](@ref), and a constructed 
 [`PolygonHierarchy`](@ref) given by `hierarchy`, are valid. In particular, the function checks:
@@ -18,17 +18,24 @@
    so that e.g. the exterior boundary curves are all counter-clockwise (relative to just themselves), the next exterior-most curves inside those 
    exteriors are all clockwise (again, relative to just themselves), and so on.
 
+The arguments `boundary_nodes` and `segments` are also used when checking for duplicate points. Any duplicate points that are also referenced in `boundary_nodes`
+and `segments` are updated so the vertex refers to the first instance of the duplicate point. This is done in-place, so that the original `boundary_nodes` and `segments` are modified.
+
+!!! danger "Mutation"
+
+    If indeed duplicate points are found, the function modifies the `boundary_nodes` and `segments` in-place. This means that the original
+    `boundary_nodes` and `segments` are modified, and the original points are not modified. The indices of the duplicates are merged into `skip_points` in-place.
+
 !!! danger "Intersecting boundaries"
 
     Another requirement for [`triangulate`](@ref) is that none of the boundaries intersect in their interior, which also prohibits 
     interior self-intersections. This is NOT checked. Similarly, segments should not intersect in their interior, which is not checked.
 """
-function check_args(points, boundary_nodes, hierarchy, boundary_curves = (); skip_points = Set{Int}())
+function check_args(points, boundary_nodes, segments, hierarchy, boundary_curves=(); skip_points=Set{Int}())
     check_dimension(points)
-    has_unique_points!(skip_points, points)
+    has_unique_points!(skip_points, points, boundary_nodes, segments)
     has_enough_points(points)
-    has_bnd = has_boundary_nodes(boundary_nodes)
-    if has_bnd
+    if has_boundary_nodes(boundary_nodes)
         has_consistent_connections(boundary_nodes)
         has_consistent_orientations(hierarchy, boundary_nodes, is_curve_bounded(boundary_curves))
     end
@@ -38,7 +45,7 @@
 struct InsufficientPointsError{P} <: Exception
     points::P
 end
-struct InconsistentConnectionError{I, J} <: Exception
+struct InconsistentConnectionError{I,J} <: Exception
     curve_index::I
     segment_index₁::I
     segment_index₂::I
@@ -81,8 +88,8 @@
         # by a combination of multiple AbstractParametricCurves and possibly a PiecewiseLinear part. Thus, the above advice
         # might not be wrong.
         str2 = "\nIf this curve is defined by an AbstractParametricCurve, you may instead need to reverse the order of the control points defining" *
-            " the sections of the curve; the `positive` keyword may also be of interest for CircularArcs and EllipticalArcs. Alternatively, for individual" *
-            " AbstractParametricCurves, note that `reverse` can be used to reverse the orientation of the curve directly instead of the control points."
+               " the sections of the curve; the `positive` keyword may also be of interest for CircularArcs and EllipticalArcs. Alternatively, for individual" *
+               " AbstractParametricCurves, note that `reverse` can be used to reverse the orientation of the curve directly instead of the control points."
         str *= str2
     end
     sign = err.should_be_positive ? "positive" : "negative"
@@ -92,28 +99,82 @@
 end
 
 function check_dimension(points)
-    valid = is_planar(points) 
-    if !valid 
-        @warn "The provided points are not in the plane. All but the first two coordinates of each point will be ignored." maxlog=1
+    valid = is_planar(points)
+    if !valid
+        @warn "The provided points are not in the plane. All but the first two coordinates of each point will be ignored." maxlog = 1
     end
     return valid
 end
 
-function has_unique_points!(skip_points, points)
+function has_unique_points!(skip_points, points, boundary_nodes, segments)
     all_unique = points_are_unique(points)
-    if !all_unique 
+    if !all_unique
         dup_seen = find_duplicate_points(points)
         if WARN_ON_DUPES[]
             io = IOBuffer()
-            println(io, "There were duplicate points. Only one of each duplicate will be used, and all other duplicates will be skipped. The indices of the duplicates are:")
+            println(io, "There were duplicate points. Only one of each duplicate will be used (the first vertex encountered in the order that follows), and all other duplicates will be skipped. The indices of the duplicates are:")
         end
-        for (p, ivec) in dup_seen 
-            for j in 2:lastindex(ivec)
+        for (p, ivec) in dup_seen
+            # Skip all but the first duplicate point for each point.
+            for j in (firstindex(ivec)+1):lastindex(ivec)
                 push!(skip_points, ivec[j])
             end
             if WARN_ON_DUPES[]
                 println(io, "  ", p, " at indices ", ivec)
             end
+            # We need to be careful about the duplicate points in the boundary nodes and segments. 
+            # https://github.com/JuliaGeometry/DelaunayTriangulation.jl/issues/220
+            ref = first(ivec)
+            if !(isnothing(segments) || num_edges(segments) == 0)
+                # We can't delete the segments while iterating, so we need to build a list 
+                E = edge_type(segments)
+                segment_map = Dict{E,E}()
+                for e in each_edge(segments)
+                    u, v = edge_vertices(e)
+                    if u ∈ ivec
+                        segment_map[e] = construct_edge(E, ref, v)
+                    elseif v ∈ ivec
+                        segment_map[e] = construct_edge(E, u, ref)
+                    end
+                end
+                for (e, new_e) in segment_map
+                    delete_edge!(segments, e)
+                    add_edge!(segments, new_e)
+                end
+            end 
+            if has_boundary_nodes(boundary_nodes)
+                if has_multiple_curves(boundary_nodes)
+                    for k in 1:num_curves(boundary_nodes)
+                        curve = get_boundary_nodes(boundary_nodes, k)
+                        for j in 1:num_sections(curve)
+                            segment = get_boundary_nodes(curve, j)
+                            for i in 1:(num_boundary_edges(segment)+1)
+                                v = get_boundary_nodes(segment, i)
+                                if v ∈ ivec
+                                    set_boundary_node!(boundary_nodes, ((k, j), i), ref)
+                                end
+                            end
+                        end
+                    end
+                elseif has_multiple_sections(boundary_nodes)
+                    for j in 1:num_sections(boundary_nodes)
+                        segment = get_boundary_nodes(boundary_nodes, j)
+                        for i in 1:(num_boundary_edges(segment)+1)
+                            v = get_boundary_nodes(segment, i)
+                            if v ∈ ivec
+                                set_boundary_node!(boundary_nodes, (j, i), ref)
+                            end
+                        end
+                    end
+                else
+                    for i in 1:(num_boundary_edges(boundary_nodes)+1)
+                        v = get_boundary_nodes(boundary_nodes, i)
+                        if v ∈ ivec
+                            set_boundary_node!(boundary_nodes, (boundary_nodes, i), ref)
+                        end
+                    end
+                end
+            end
         end
         if WARN_ON_DUPES[]
             println(io, "To suppress this warning, call `DelaunayTriangulation.toggle_warn_on_dupes!()`.")
@@ -125,7 +186,7 @@
     return true
 end
 
-function has_enough_points(points)
+function has_enough_points(points) 
     has_enough = num_points(points) ≥ 3
     !has_enough && throw(InsufficientPointsError(points))
     return true
@@ -171,7 +232,7 @@
     end
     return true
 end
-function has_consistent_connections_multiple_sections(boundary_nodes, curve_index = 0)
+function has_consistent_connections_multiple_sections(boundary_nodes, curve_index=0)
     ns = num_sections(boundary_nodes)
     segmentⱼ₋₁ = get_boundary_nodes(boundary_nodes, 1)
     nn = num_boundary_edges(segmentⱼ₋₁) + 1

src/algorithms/triangulation/main.jl

@@ -86,9 +86,10 @@
 
 !!! warning "Mutation"
 
-    The `segments` may get mutated in two ways: (1) Segments may get rotated so that `(i, j)` becomes `(j, i)`. (2) If there are
+    The `segments` may get mutated in three ways: (1) Segments may get rotated so that `(i, j)` becomes `(j, i)`. (2) If there are
     segments that are collinear with other segments, then they may get split into chain of non-overlapping connecting segments (also see below). For 
-    curve-bounded domains, segments are also split so that no subsegment's diametral circle contains any other point.
+    curve-bounded domains, segments are also split so that no subsegment's diametral circle contains any other point. (3) If there are 
+    duplicate points in `points` that are also referenced in `segments`, then the duplicated vertex will be changed to the first occurrence of the vertex in `points`.
 
 !!! warning "Intersecting segments"
 
@@ -99,6 +100,10 @@
 
 - `boundary_nodes`
 
+!!! warning "Mutation"
+
+    If duplicate points are found, and they are also referenced in `boundary_nodes`, then the duplicated vertex will be changed to the first occurrence of the vertex in `points`.
+
 !!! warning "Points outside of boundary curves"
 
     While for standard domains with piecewise linear boundaries (or no boundaries) it is fine for points to be 
@@ -151,9 +156,9 @@
     if isnothing(full_polygon_hierarchy)
         full_polygon_hierarchy = construct_polygon_hierarchy(points, boundary_nodes; IntegerType)
     end
-    skip_points_set = Set{IntegerType}(skip_points)
+    skip_points_set = Set{IntegerType}(skip_points) 
     n = length(skip_points_set)
-    check_arguments && check_args(points, boundary_nodes, full_polygon_hierarchy, boundary_curves; skip_points = skip_points_set)
+    check_arguments && check_args(points, boundary_nodes, segments, full_polygon_hierarchy, boundary_curves; skip_points = skip_points_set)
     if length(skip_points_set) > n 
         setdiff!(insertion_order, skip_points_set) 
     end

src/data_structures/mesh_refinement/boundary_enricher.jl

@@ -388,13 +388,13 @@
 Base.copy(enricher::BoundaryEnricher) = enrcopy(enricher)
 enrcopy(::Nothing; kwargs...) = nothing
 
-function enrcopy(enricher::BoundaryEnricher; 
-        points = copy(get_points(enricher)), 
-        boundary_nodes  = copy(get_boundary_nodes(enricher)),
-        segments = copy(get_segments(enricher)),
-        boundary_curves = _plcopy.(get_boundary_curves(enricher); points),
-        polygon_hierarchy = copy(get_polygon_hierarchy(enricher)),
-        boundary_edge_map = copy(get_boundary_edge_map(enricher))) 
+function enrcopy(enricher::BoundaryEnricher;
+    points=copy(get_points(enricher)),
+    boundary_nodes=copy(get_boundary_nodes(enricher)),
+    segments=copy(get_segments(enricher)),
+    boundary_curves=_plcopy.(get_boundary_curves(enricher); points),
+    polygon_hierarchy=copy(get_polygon_hierarchy(enricher)),
+    boundary_edge_map=copy(get_boundary_edge_map(enricher)))
     parent_map = copy(get_parent_map(enricher))
     curve_index_map = copy(get_curve_index_map(enricher))
     spatial_tree = copy(get_spatial_tree(enricher))
@@ -454,7 +454,8 @@
     boundary_nodes = get_boundary_nodes(enricher)
     hierarchy = get_polygon_hierarchy(enricher)
     boundary_curves = get_boundary_curves(enricher)
-    return check_args(points, boundary_nodes, hierarchy, boundary_curves)
+    segments = get_segments(enricher)
+    return check_args(points, boundary_nodes, segments, hierarchy, boundary_curves)
 end
 
 """

src/geometric_primitives/boundary_nodes.jl

@@ -759,3 +759,25 @@
 @inline function _get_skeleton_contiguous(boundary_nodes, ::Type{I}) where {I}
     return I[]
 end
+
+"""
+    set_boundary_node!(boundary_nodes, pos, node)
+
+Given a set of `boundary_nodes`, sets the boundary node at position `pos` to `node`.
+Here, `pos[1]` is such that `get_boundary_nodes(boundary_nodes, pos[1])`
+is the section that the node will be set onto, and `pos[2]` gives the position
+of the array to set `node` into. In particular, 
+
+    set_boundary_node!(boundary_nodes, pos, node)
+
+is the same as
+
+    get_boundary_nodes(boundary_nodes, pos[1])[pos[2]] = node
+
+assuming `setindex!` is defined for the type of `boundary_nodes`.
+"""
+function set_boundary_node!(boundary_nodes, pos, node)
+    nodes = get_boundary_nodes(boundary_nodes, pos[1])
+    nodes[pos[2]] = node
+    return boundary_nodes
+end

test/Project.toml

@@ -37,3 +37,6 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 Aqua = "0.8.7"
+
+[sources]
+DelaunayTriangulation = {path = ".."}

test/helper_functions.jl

@@ -10,6 +10,7 @@ using OrderedCollections
 using Distributions
 using InteractiveUtils
 using Test
+
 using DelaunayTriangulation
 import SpatialIndexing as SI
 getxy((0.0, 0.0)) # avoid shadow

test/interfaces/boundary_nodes.jl

@@ -66,7 +66,7 @@ end
     map3 = DT.construct_ghost_vertex_map(bn3)
     idx = DT.𝒢
     @test map1 ==
-        Dict(
+          Dict(
         idx => (1, 1), idx - 1 => (1, 2), idx - 2 => (1, 3), idx - 3 => (1, 4),
         idx - 4 => (2, 1), idx - 5 => (2, 2),
         idx - 6 => (3, 1), idx - 7 => (3, 2),
@@ -157,7 +157,7 @@ end
     end
     cbn = copy(bn)
     _bn_map = DT._bemcopy(bn_map; boundary_nodes=cbn)
-    @test bn_map == _bn_map 
+    @test bn_map == _bn_map
     _bn = first.(values(_bn_map))
     @test all(x -> x === cbn, _bn)
     @test all(x -> !(x === bn), _bn)
@@ -185,7 +185,7 @@ end
     @test _bn_map == bn_map && !(bn_map === _bn_map)
     bn = Int[]
     bn_map = DT.construct_boundary_edge_map(bn)
-    @test bn_map == Dict{Tuple{Int32, Int32}, Tuple{Vector{Int}, Int}}()
+    @test bn_map == Dict{Tuple{Int32,Int32},Tuple{Vector{Int},Int}}()
 end
 
 @testset "insert_boundary_node!" begin
@@ -243,3 +243,46 @@ end
         [[13, 14, 16, 17], [17, 18, 20], [20]],
     ]
 end
+
+@testset "set_boundary_node!" begin
+    # Single curve
+    bn = [1, 2, 3, 4, 5, 6, 7, 1]
+    DT.set_boundary_node!(bn, (bn, 5), 17)
+    DT.set_boundary_node!(bn, (bn, 1), 13)
+    @test bn == [13, 2, 3, 4, 17, 6, 7, 1]
+
+    # Multiple sections
+    bn = [[1, 2, 3, 4], [4, 5, 6, 7, 8], [8, 9, 10, 1]]
+    DT.set_boundary_node!(bn, (1, 2), 22)
+    DT.set_boundary_node!(bn, (1, 4), 33)
+    DT.set_boundary_node!(bn, (2, 4), 44)
+    DT.set_boundary_node!(bn, (3, 1), 55)
+    @test bn == [[1, 22, 3, 33], [4, 5, 6, 44, 8], [55, 9, 10, 1]]
+
+    # Multiple curves
+    bn = [
+        [[1, 2, 3, 4, 5], [5, 6, 7], [7, 8], [8, 9, 10, 1]],
+        [[13, 14, 15, 16, 17], [17, 18, 19, 20], [20, 13]],
+    ]
+    DT.set_boundary_node!(bn, ((1, 1), 1), 99)
+    DT.set_boundary_node!(bn, ((1, 2), 3), 88)
+    DT.set_boundary_node!(bn, ((1, 3), 2), 77)
+    DT.set_boundary_node!(bn, ((1, 4), 3), 66)
+    DT.set_boundary_node!(bn, ((2, 1), 3), 55)
+    DT.set_boundary_node!(bn, ((2, 2), 3), 44)
+    DT.set_boundary_node!(bn, ((2, 3), 2), 33)
+    @test bn == [
+        [[99, 2, 3, 4, 5], [5, 6, 88], [7, 77], [8, 9, 66, 1]],
+        [[13, 14, 55, 16, 17], [17, 18, 44, 20], [20, 33]],
+    ]
+
+    # Invalid index should throw
+    @test_throws MethodError DT.set_boundary_node!([1, 2, 3], (1, 10), 99)
+    @test_throws BoundsError DT.set_boundary_node!([1, 2, 3], ([1, 2, 3], 4), 99)
+
+    # Aliasing check
+    bn = [[10, 20], [30, 40]]
+    s = get_boundary_nodes(bn, 1)
+    DT.set_boundary_node!(bn, (1, 2), 999)
+    @test s[2] == 999  # test that the reference s is updated too
+end