diff --git a/NEWS.md b/NEWS.md
index 0cdffce66..ce374e7bd 100644
--- a/NEWS.md
+++ b/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!()`.
diff --git a/Project.toml b/Project.toml
index 1d6340334..ec8eb0f58 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "DelaunayTriangulation"
 uuid = "927a84f5-c5f4-47a5-9785-b46e178433df"
 authors = ["Daniel VandenHeuvel <danj.vandenheuvel@gmail.com>"]
-version = "1.6.4"
+version = "1.6.5"
 
 [deps]
 AdaptivePredicates = "35492f91-a3bd-45ad-95db-fcad7dcfedb7"
diff --git a/docs/Project.toml b/docs/Project.toml
index 2b3a5fcf2..10db242cf 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -34,3 +34,6 @@ SimpleGraphs = "0.8"
 StableRNGs = "1.0"
 StatsBase = "0.34"
 Test = "1.11"
+
+[sources]
+DelaunayTriangulation = { path = ".." }
diff --git a/docs/src/softwarecomparisonc.png b/docs/src/softwarecomparisonc.png
index 79d1d3316..bb7bdcab5 100644
Binary files a/docs/src/softwarecomparisonc.png and b/docs/src/softwarecomparisonc.png differ
diff --git a/src/algorithms/triangulation/check_args.jl b/src/algorithms/triangulation/check_args.jl
index 5831cab04..8078a6f28 100644
--- a/src/algorithms/triangulation/check_args.jl
+++ b/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 @@ If `boundary_nodes` are provided, meaning [`has_boundary_nodes`](@ref), then the
    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 @@ end
 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 @@ function Base.showerror(io::IO, err::InconsistentOrientationError)
         # 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 @@ function Base.showerror(io::IO, err::InconsistentOrientationError)
 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 @@ function has_unique_points!(skip_points, points)
     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 @@ function has_consistent_connections_multiple_curves(boundary_nodes)
     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
diff --git a/src/algorithms/triangulation/main.jl b/src/algorithms/triangulation/main.jl
index be2dbc9d6..721d69b17 100644
--- a/src/algorithms/triangulation/main.jl
+++ b/src/algorithms/triangulation/main.jl
@@ -86,9 +86,10 @@ Here are some warnings to consider for some of the arguments.
 
 !!! 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 @@ Here are some warnings to consider for some of the arguments.
 
 - `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 @@ function triangulate(
     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
diff --git a/src/data_structures/mesh_refinement/boundary_enricher.jl b/src/data_structures/mesh_refinement/boundary_enricher.jl
index 1b45ab05d..6c4579a18 100644
--- a/src/data_structures/mesh_refinement/boundary_enricher.jl
+++ b/src/data_structures/mesh_refinement/boundary_enricher.jl
@@ -388,13 +388,13 @@ end
 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 @@ function check_args(enricher::BoundaryEnricher)
     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
 
 """
diff --git a/src/geometric_primitives/boundary_nodes.jl b/src/geometric_primitives/boundary_nodes.jl
index 76b2a5558..c80d17880 100644
--- a/src/geometric_primitives/boundary_nodes.jl
+++ b/src/geometric_primitives/boundary_nodes.jl
@@ -759,3 +759,25 @@ end
 @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
\ No newline at end of file
diff --git a/test/Project.toml b/test/Project.toml
index cbe42aef0..1c2d80d8b 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -37,3 +37,6 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 Aqua = "0.8.7"
+
+[sources]
+DelaunayTriangulation = {path = ".."}
\ No newline at end of file
diff --git a/test/helper_functions.jl b/test/helper_functions.jl
index abbe9994a..8bd18a401 100644
--- a/test/helper_functions.jl
+++ b/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
diff --git a/test/interfaces/boundary_nodes.jl b/test/interfaces/boundary_nodes.jl
index cd1f61d43..acfb252bc 100644
--- a/test/interfaces/boundary_nodes.jl
+++ b/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
diff --git a/test/runtests.jl b/test/runtests.jl
index 96fcc941c..be0bce873 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -22,13 +22,13 @@ function safe_include(filename; name=filename, push=true, verbose=true) # Workar
     end
 end
 
-@testset verbose = true "DelaunayTriangulation.jl" begin
-    @testset verbose = true "Aqua" begin
+@testset verbose = false "DelaunayTriangulation.jl" begin
+    @testset verbose = false "Aqua" begin
         Aqua.test_all(DelaunayTriangulation; ambiguities=false, project_extras=false, unbound_args=false) # don't care about julia < 1.2
         Aqua.test_ambiguities(DelaunayTriangulation) # don't pick up Base and Core...
     end
 
-    @testset verbose = true "Triangulation" begin
+    @testset verbose = false "Triangulation" begin
         safe_include("triangulation/rectangle.jl")
         safe_include("triangulation/bowyer_watson.jl")
         safe_include("triangulation/triangulate.jl")
@@ -38,14 +38,14 @@ end
         safe_include("triangulation/weighted.jl")
     end
 
-    @testset verbose = true "Interfaces" begin
+    @testset verbose = false "Interfaces" begin
         safe_include("interfaces/triangles.jl")
         safe_include("interfaces/edges.jl")
         safe_include("interfaces/points.jl")
         safe_include("interfaces/boundary_nodes.jl")
     end
 
-    @testset verbose = true "Data Structures" begin
+    @testset verbose = false "Data Structures" begin
         safe_include("data_structures/adjacent.jl")
         safe_include("data_structures/adjacent2vertex.jl")
         safe_include("data_structures/graph.jl")
@@ -63,20 +63,20 @@ end
         safe_include("data_structures/polygon_hierarchy.jl", verbose=false)
     end
 
-    @testset verbose = true "Predicates" begin
+    @testset verbose = false "Predicates" begin
         safe_include("predicates/certificate.jl")
         safe_include("predicates/boundaries_and_ghosts.jl")
         safe_include("predicates/general.jl")
         safe_include("predicates/index_and_ghost_handling.jl")
     end
 
-    @testset verbose = true "Utilities" begin
+    @testset verbose = false "Utilities" begin
         safe_include("utils.jl")
         safe_include("geo_utils.jl")
         safe_include("helper_function_tests.jl")
     end
 
-    @testset verbose = true "Point Location" begin
+    @testset verbose = false "Point Location" begin
         safe_include("point_location/brute_force.jl")
         safe_include("point_location/select_initial_point.jl")
         safe_include("point_location/select_initial_triangle_interior_node.jl")
@@ -86,7 +86,7 @@ end
         safe_include("point_location/find_polygon.jl")
     end
      
-    @testset verbose = true "Operations" begin
+    @testset verbose = false "Operations" begin
         safe_include("operations/add_triangle.jl")
         safe_include("operations/delete_triangle.jl")
         safe_include("operations/add_ghost_triangles.jl")
@@ -101,27 +101,27 @@ end
         safe_include("operations/delete_holes.jl")
     end
 
-    @testset verbose = true "Constrained Triangulation" begin
+    @testset verbose = false "Constrained Triangulation" begin
         safe_include("constrained_triangulation/segment_location.jl")
         safe_include("constrained_triangulation/segment_insertion.jl")
     end
 
-    @testset verbose = true "Refinement" begin
+    @testset verbose = false "Refinement" begin
         safe_include("refinement/refine.jl")
         safe_include("refinement/curve_bounded.jl")
     end
 
-    @testset verbose = true "Voronoi" begin
+    @testset verbose = false "Voronoi" begin
         safe_include("voronoi/voronoi.jl")
         safe_include("voronoi/power.jl")
     end
 
-    @testset verbose = true "Makie" begin
+    @testset verbose = false "Makie" begin
         safe_include("makie/makie.jl")
     end
 
-    @testset verbose = true "Run the documentation examples" begin
-        @testset verbose = true "Check that the applications in the docs run" begin
+    @testset verbose = false "Run the documentation examples" begin
+        @testset verbose = false "Check that the applications in the docs run" begin
             app_dir = joinpath(dirname(dirname(pathof(DelaunayTriangulation))), "docs", "src", "literate_applications")
             app_files = readdir(app_dir)
             for file in app_files
@@ -131,7 +131,7 @@ end
             isfile(mp4_path) && rm(mp4_path)
         end
 
-        @testset verbose = true "Test the tutorials" begin
+        @testset verbose = false "Test the tutorials" begin
             tut_dir = joinpath(dirname(dirname(pathof(DelaunayTriangulation))), "docs", "src", "literate_tutorials")
             tut_files = readdir(tut_dir)
             for file in tut_files
@@ -139,7 +139,7 @@ end
             end
         end
 
-        @testset verbose = true "Test the readme example" begin
+        @testset verbose = false "Test the readme example" begin
             safe_include("readme_example.jl")
         end
     end
diff --git a/test/triangulation/check_args.jl b/test/triangulation/check_args.jl
index abd53ae4d..f5842208e 100644
--- a/test/triangulation/check_args.jl
+++ b/test/triangulation/check_args.jl
@@ -5,49 +5,47 @@ _test_throws(e1, e2=e1) = @static VERSION ≥ v"1.9" ? e1 : e2
 
 @testset "check_dimension" begin
     points = rand(2, 50)
-    boundary_nodes = nothing
     hierarchy = DT.construct_polygon_hierarchy(points)
-    @test DT.check_args(points, boundary_nodes, hierarchy)
+    @test DT.check_args(points, nothing, nothing, hierarchy)
 
     points = rand(3, 50)
-    boundary_nodes = nothing
     hierarchy = DT.construct_polygon_hierarchy(points)
-    @test_logs (:warn, "The provided points are not in the plane. All but the first two coordinates of each point will be ignored.") DT.check_args(points, boundary_nodes, hierarchy)
+    @test_logs (:warn, "The provided points are not in the plane. All but the first two coordinates of each point will be ignored.") DT.check_args(points, nothing, nothing, hierarchy)
 end
 
 @testset "A simple case" begin
     points = rand(2, 50)
     boundary_nodes = nothing
     hierarchy = DT.construct_polygon_hierarchy(points)
-    @test DT.check_args(points, boundary_nodes, hierarchy)
+    @test DT.check_args(points, nothing, nothing, hierarchy)
 end
 
 @testset "Not enough points" begin
     points = rand(2, 2)
-    boundary_nodes = nothing
     hierarchy = DT.construct_polygon_hierarchy(points)
-    @test_throws _test_throws(DT.InsufficientPointsError) DT.check_args(points, boundary_nodes, hierarchy)
-    @test_throws _test_throws("InsufficientPointsError: The provided point set has 2 points, but triangulations require at least three points.", DT.InsufficientPointsError) DT.check_args(points, boundary_nodes, hierarchy)
+    @test_throws _test_throws(DT.InsufficientPointsError) DT.check_args(points, nothing, nothing, hierarchy)
+    @test_throws _test_throws("InsufficientPointsError: The provided point set has 2 points, but triangulations require at least three points.", DT.InsufficientPointsError) DT.check_args(points, nothing, nothing, hierarchy)
     @test_throws _test_throws(DT.InsufficientPointsError) triangulate(points, predicates=rt())
 end
 
 @testset "Duplicate points" begin
     points = [(1.0, 1.0), (2.0, 2.0), (5.5, 17.3), (1.0, 1.0), (0.0, 0.5), (1.7, 5.5), (2.0, 2.0), (2.0, 2.0), (25.5, 17.3), (5.5, 17.3)]
     boundary_nodes = nothing
+    segments = nothing
     hierarchy = DT.construct_polygon_hierarchy(points)
     skip_points = Int[]
     @test_logs (
         :warn,
-        "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:\n  (1.0, 1.0) at indices [1, 4]\n  (2.0, 2.0) at indices [2, 7, 8]\n  (5.5, 17.3) at indices [3, 10]\nTo suppress this warning, call `DelaunayTriangulation.toggle_warn_on_dupes!()`.\n") DT.check_args(points, boundary_nodes, hierarchy; skip_points)
+        "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:\n  (1.0, 1.0) at indices [1, 4]\n  (2.0, 2.0) at indices [2, 7, 8]\n  (5.5, 17.3) at indices [3, 10]\nTo suppress this warning, call `DelaunayTriangulation.toggle_warn_on_dupes!()`.\n") DT.check_args(points, boundary_nodes, segments, hierarchy; skip_points)
     @test skip_points == [4, 7, 8, 10]
     skip_points = [13, 15]
-    DT.check_args(points, boundary_nodes, hierarchy; skip_points)
+    DT.check_args(points, boundary_nodes, segments, hierarchy; skip_points)
     @test skip_points == [13, 15, 4, 7, 8, 10]
     @test_logs (
         :warn,
-        "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:\n  (1.0, 1.0) at indices [1, 4]\n  (2.0, 2.0) at indices [2, 7, 8]\n  (5.5, 17.3) at indices [3, 10]\nTo suppress this warning, call `DelaunayTriangulation.toggle_warn_on_dupes!()`.\n") triangulate(points, predicates=rt())
+        "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:\n  (1.0, 1.0) at indices [1, 4]\n  (2.0, 2.0) at indices [2, 7, 8]\n  (5.5, 17.3) at indices [3, 10]\nTo suppress this warning, call `DelaunayTriangulation.toggle_warn_on_dupes!()`.\n") triangulate(points, predicates=rt())
     DT.toggle_warn_on_dupes!()
-    @test_nowarn DT.check_args(points, boundary_nodes, hierarchy)
+    @test_nowarn DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test_nowarn triangulate(points; predicates=rt())
     points = [(1.0, 1.0), (2.0, 2.0), (5.5, 17.3), (1.0, 1.0), (0.0, 0.5), (1.7, 5.5), (2.0, 2.0), (2.0, 2.0), (25.5, 17.3), (5.5, 17.3)]
     tri = triangulate(points; predicates=rt())
@@ -56,7 +54,7 @@ end
     @test !DT.has_vertex(tri, 8)
     @test !DT.has_vertex(tri, 10)
     @test DT.validate_triangulation(tri)
-    tri = triangulate(points; predicates=rt(), skip_points = (1,))
+    tri = triangulate(points; predicates=rt(), skip_points=(1,))
     @test !DT.has_vertex(tri, 1)
     @test !DT.has_vertex(tri, 4)
     @test !DT.has_vertex(tri, 7)
@@ -68,92 +66,96 @@ end
 @testset "Orientation and connectivity of a single boundary curve" begin
     points = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)]
     boundary_nodes = [1, 2, 3, 4, 1]
+    segments = nothing
     hierarchy = DT.construct_polygon_hierarchy(points, boundary_nodes)
-    @test DT.check_args(points, boundary_nodes, hierarchy)
+    @test DT.check_args(points, boundary_nodes, segments, hierarchy)
     boundary_nodes[5] = 3
-    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
-    @test_throws _test_throws("InconsistentConnectionError: The boundary ends in vertex 3 but starts at vertex 1.", DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
+    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
+    @test_throws _test_throws("InconsistentConnectionError: The boundary ends in vertex 3 but starts at vertex 1.", DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test_throws _test_throws(DT.InconsistentConnectionError) triangulate(points; boundary_nodes, predicates=rt())
 
     boundary_nodes = [4, 3, 2, 1, 4]
     hierarchy = DT.construct_polygon_hierarchy(points, boundary_nodes)
-    @test_throws _test_throws(DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, hierarchy)
-    @test_throws _test_throws("InconsistentOrientationError: The orientation of the boundary curve with index 1 should be positive, but it is negative. You may be able to fix this by passing the curve as reverse(curve).", DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, hierarchy)
+    @test_throws _test_throws(DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, segments, hierarchy)
+    @test_throws _test_throws("InconsistentOrientationError: The orientation of the boundary curve with index 1 should be positive, but it is negative. You may be able to fix this by passing the curve as reverse(curve).", DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test_throws _test_throws(DT.InconsistentOrientationError) triangulate(points; boundary_nodes, predicates=rt())
 
     boundary_nodes = [[1, 2, 3, 4, 1]]
     hierarchy = DT.construct_polygon_hierarchy(points, boundary_nodes)
-    @test DT.check_args(points, boundary_nodes, hierarchy)
+    @test DT.check_args(points, boundary_nodes, segments, hierarchy)
     boundary_nodes[1][1] = 2
-    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
-    @test_throws _test_throws("InconsistentConnectionError: The boundary ends in vertex 1 but starts at vertex 2.", DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
+    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
+    @test_throws _test_throws("InconsistentConnectionError: The boundary ends in vertex 1 but starts at vertex 2.", DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test_throws _test_throws(DT.InconsistentConnectionError) triangulate(points; boundary_nodes, predicates=rt())
 
     boundary_nodes = [[4, 3, 2, 1, 4]]
     hierarchy = DT.construct_polygon_hierarchy(points, boundary_nodes)
-    @test_throws _test_throws(DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, hierarchy)
-    @test_throws _test_throws("InconsistentOrientationError: The orientation of the boundary curve with index 1 should be positive, but it is negative. You may be able to fix this by passing the curve as reverse(reverse.(curve)).", DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, hierarchy)
+    @test_throws _test_throws(DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, segments, hierarchy)
+    @test_throws _test_throws("InconsistentOrientationError: The orientation of the boundary curve with index 1 should be positive, but it is negative. You may be able to fix this by passing the curve as reverse(reverse.(curve)).", DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test_throws _test_throws(DT.InconsistentOrientationError) triangulate(points; boundary_nodes, predicates=rt())
 
     boundary_nodes = [[[1, 2, 3, 4, 1]]]
     hierarchy = DT.construct_polygon_hierarchy(points, boundary_nodes)
-    @test DT.check_args(points, boundary_nodes, hierarchy)
+    @test DT.check_args(points, boundary_nodes, segments, hierarchy)
     boundary_nodes[1][1][end] = 2
-    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
-    @test_throws _test_throws("InconsistentConnectionError: The boundary ends in vertex 2 but starts at vertex 1.", DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
+    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
+    @test_throws _test_throws("InconsistentConnectionError: The boundary ends in vertex 2 but starts at vertex 1.", DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test_throws _test_throws(DT.InconsistentConnectionError) triangulate(points; boundary_nodes, predicates=rt())
 
     boundary_nodes = [[[4, 3, 2, 1, 4]]]
     hierarchy = DT.construct_polygon_hierarchy(points, boundary_nodes)
-    @test_throws _test_throws(DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, hierarchy)
-    @test_throws _test_throws("InconsistentOrientationError: The orientation of the boundary curve with index 1 should be positive, but it is negative. You may be able to fix this by passing the curve as reverse(reverse.(curve)).", DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, hierarchy)
+    @test_throws _test_throws(DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, segments, hierarchy)
+    @test_throws _test_throws("InconsistentOrientationError: The orientation of the boundary curve with index 1 should be positive, but it is negative. You may be able to fix this by passing the curve as reverse(reverse.(curve)).", DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test_throws _test_throws(DT.InconsistentOrientationError) triangulate(points; boundary_nodes, predicates=rt())
 end
 
 @testset "Orientation and connectivity of a sectioned boundary curve" begin
     points = [(0.0, 0.0), (0.5, 0.0), (1.0, 0.0), (1.0, 1.0), (0.5, 1.0), (0.0, 1.0), (0.0, 0.5), (0.0, 0.25)]
     boundary_nodes = [[1, 2, 3], [3, 4], [4, 5, 6], [6, 7, 8, 1]]
+    segments = nothing
     hierarchy = DT.construct_polygon_hierarchy(points, boundary_nodes)
-    @test DT.check_args(points, boundary_nodes, hierarchy)
+    @test DT.check_args(points, boundary_nodes, segments, hierarchy)
     boundary_nodes[1][3] = 5
-    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
-    @test_throws _test_throws("InconsistentConnectionError: Segment 1 ends at vertex 5 but the next segment, segment 2, starts at vertex 3.", DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
+    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
+    @test_throws _test_throws("InconsistentConnectionError: Segment 1 ends at vertex 5 but the next segment, segment 2, starts at vertex 3.", DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test_throws _test_throws(DT.InconsistentConnectionError) triangulate(points; boundary_nodes, predicates=rt())
     boundary_nodes[1][3] = 3
     boundary_nodes[4][4] = 2
-    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
-    @test_throws _test_throws("InconsistentConnectionError: Segment 4 ends at vertex 2 but the next segment, segment 1, starts at vertex 1.", DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
+    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
+    @test_throws _test_throws("InconsistentConnectionError: Segment 4 ends at vertex 2 but the next segment, segment 1, starts at vertex 1.", DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test_throws _test_throws(DT.InconsistentConnectionError) triangulate(points; boundary_nodes, predicates=rt())
     boundary_nodes[4][4] = 1
     boundary_nodes[1][1] = 8
-    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
-    @test_throws _test_throws("InconsistentConnectionError: Segment 4 ends at vertex 1 but the next segment, segment 1, starts at vertex 8.", DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
+    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
+    @test_throws _test_throws("InconsistentConnectionError: Segment 4 ends at vertex 1 but the next segment, segment 1, starts at vertex 8.", DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test_throws _test_throws(DT.InconsistentConnectionError) triangulate(points; boundary_nodes, predicates=rt())
 
 
     points = [(0.0, 0.0), (0.5, 0.0), (1.0, 0.0), (1.0, 1.0), (0.5, 1.0), (0.0, 1.0), (0.0, 0.5), (0.0, 0.25)]
     boundary_nodes = [[1, 8, 7, 6], [6, 5, 4], [4, 3], [3, 2, 1]]
+    segments = nothing
     hierarchy = DT.construct_polygon_hierarchy(points, boundary_nodes)
-    @test_throws _test_throws(DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, hierarchy)
-    @test_throws _test_throws("InconsistentOrientationError: The orientation of the boundary curve with index 1 should be positive, but it is negative. You may be able to fix this by passing the curve as reverse(reverse.(curve)).", DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, hierarchy)
+    @test_throws _test_throws(DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, segments, hierarchy)
+    @test_throws _test_throws("InconsistentOrientationError: The orientation of the boundary curve with index 1 should be positive, but it is negative. You may be able to fix this by passing the curve as reverse(reverse.(curve)).", DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test_throws _test_throws(DT.InconsistentOrientationError) triangulate(points; boundary_nodes, predicates=rt())
 
     points = [(0.0, 0.0), (1.0, 0.0), (0.0, 1.0)]
     boundary_nodes = [[1, 2], [2, 3, 1]]
+    segments = nothing
     hierarchy = DT.construct_polygon_hierarchy(points, boundary_nodes)
-    @test DT.check_args(points, boundary_nodes, hierarchy)
+    @test DT.check_args(points, boundary_nodes, segments, hierarchy)
     boundary_nodes[1][2] = 3
-    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
-    @test_throws _test_throws("InconsistentConnectionError: Segment 1 ends at vertex 3 but the next segment, segment 2, starts at vertex 2.", DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
+    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
+    @test_throws _test_throws("InconsistentConnectionError: Segment 1 ends at vertex 3 but the next segment, segment 2, starts at vertex 2.", DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test_throws _test_throws(DT.InconsistentConnectionError) triangulate(points; boundary_nodes, predicates=rt())
     boundary_nodes[1][2] = 2
     boundary_nodes[2][3] = 5
-    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
-    @test_throws _test_throws("InconsistentConnectionError: Segment 2 ends at vertex 5 but the next segment, segment 1, starts at vertex 1.", DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
+    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
+    @test_throws _test_throws("InconsistentConnectionError: Segment 2 ends at vertex 5 but the next segment, segment 1, starts at vertex 1.", DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
     boundary_nodes[2][3] = 1
     boundary_nodes[1][1] = 3
-    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
-    @test_throws _test_throws("InconsistentConnectionError: Segment 2 ends at vertex 1 but the next segment, segment 1, starts at vertex 3.", DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
+    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
+    @test_throws _test_throws("InconsistentConnectionError: Segment 2 ends at vertex 1 but the next segment, segment 1, starts at vertex 3.", DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test_throws _test_throws(DT.InconsistentConnectionError) triangulate(points; boundary_nodes, predicates=rt())
 end
 
@@ -163,33 +165,35 @@ end
         (0.3, 0.3), (0.7, 0.3), (0.7, 0.7), (0.3, 0.7),
     ]
     boundary_nodes = [[[1, 2, 3, 4, 5, 6, 7, 1]], [[11, 10, 9, 8, 11]]]
+    segments = nothing
     hierarchy = DT.construct_polygon_hierarchy(points, boundary_nodes)
-    @test DT.check_args(points, boundary_nodes, hierarchy)
+    @test DT.check_args(points, boundary_nodes, segments, hierarchy)
     boundary_nodes[1][1][end] = 8
-    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
-    @test_throws _test_throws("InconsistentConnectionError: The boundary curve with index 1 ends in vertex 8 but starts at vertex 1.", DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
+    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
+    @test_throws _test_throws("InconsistentConnectionError: The boundary curve with index 1 ends in vertex 8 but starts at vertex 1.", DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test_throws _test_throws(DT.InconsistentConnectionError) triangulate(points; boundary_nodes, predicates=rt())
     boundary_nodes[1][1][end] = 1
     boundary_nodes[1][1][1] = 4
-    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
-    @test_throws _test_throws("InconsistentConnectionError: The boundary curve with index 1 ends in vertex 1 but starts at vertex 4.", DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
+    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
+    @test_throws _test_throws("InconsistentConnectionError: The boundary curve with index 1 ends in vertex 1 but starts at vertex 4.", DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test_throws _test_throws(DT.InconsistentConnectionError) triangulate(points; boundary_nodes, predicates=rt())
     boundary_nodes[1][1][1] = 1
     boundary_nodes[2][1][end] = 5
-    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
-    @test_throws _test_throws("InconsistentConnectionError: The boundary curve with index 2 ends in vertex 5 but starts at vertex 11.", DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
+    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
+    @test_throws _test_throws("InconsistentConnectionError: The boundary curve with index 2 ends in vertex 5 but starts at vertex 11.", DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test_throws _test_throws(DT.InconsistentConnectionError) triangulate(points; boundary_nodes, predicates=rt())
 
     boundary_nodes = [[[1, 7, 6, 5, 4, 3, 2, 1]], [[11, 10, 9, 8, 11]]]
+    segments = nothing
     hierarchy = DT.construct_polygon_hierarchy(points, boundary_nodes)
-    @test_throws _test_throws(DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, hierarchy)
-    @test_throws _test_throws("InconsistentOrientationError: The orientation of the boundary curve with index 1 should be positive, but it is negative. You may be able to fix this by passing the curve as reverse(reverse.(curve)).", DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, hierarchy)
+    @test_throws _test_throws(DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, segments, hierarchy)
+    @test_throws _test_throws("InconsistentOrientationError: The orientation of the boundary curve with index 1 should be positive, but it is negative. You may be able to fix this by passing the curve as reverse(reverse.(curve)).", DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test_throws _test_throws(DT.InconsistentOrientationError) triangulate(points; boundary_nodes, predicates=rt())
     boundary_nodes = [[[1, 2, 3, 4, 5, 6, 7, 1]], [[11, 8, 9, 10, 11]]]
     hierarchy = DT.construct_polygon_hierarchy(points, boundary_nodes)
-    @test_throws _test_throws(DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, hierarchy)
-    @test_throws _test_throws("InconsistentOrientationError: The orientation of the boundary curve with index 2 should be negative, but it is positive. You may be able to fix this by passing the curve as reverse(reverse.(curve)).", DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, hierarchy)
-    @test_throws _test_throws(DT.InconsistentOrientationError) triangulate(points; boundary_nodes, predicates=rt())
+    @test_throws _test_throws(DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, segments, hierarchy)
+    @test_throws _test_throws("InconsistentOrientationError: The orientation of the boundary curve with index 2 should be negative, but it is positive. You may be able to fix this by passing the curve as reverse(reverse.(curve)).", DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, segments, hierarchy)
+    @test_throws _test_throws(DT.InconsistentOrientationError) triangulate(points; boundary_nodes, segments, predicates=rt())
 end
 
 @testset "Orientation and connectivity of a disjoint boundary" begin
@@ -352,14 +356,15 @@ end
     curves = [J_curve, U_curve, L_curve, I_curve, A_curve_outline, A_curve_hole, dot_1, dot_2, dot_3, dot_4]
     boundary_nodes, points = convert_boundary_points_to_indices(curves)
     hierarchy = DT.construct_polygon_hierarchy(points, boundary_nodes)
-    @test DT.check_args(points, boundary_nodes, hierarchy)
+    segments = nothing
+    @test DT.check_args(points, boundary_nodes, segments, hierarchy)
 
     dot_1 = [[Z1, J2, I2, H2], [H2, G2, F2, E2, D2, C2, B2, A2, Z1]]
     curves = [J_curve, U_curve, L_curve, I_curve, A_curve_outline, A_curve_hole, dot_1, dot_2, dot_3, dot_4]
     boundary_nodes, points = convert_boundary_points_to_indices(curves)
     hierarchy = DT.construct_polygon_hierarchy(points, boundary_nodes)
-    @test_throws _test_throws(DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, hierarchy)
-    @test_throws _test_throws("InconsistentOrientationError: The orientation of the boundary curve with index 7 should be positive, but it is negative. You may be able to fix this by passing the curve as reverse(reverse.(curve)).", DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, hierarchy)
+    @test_throws _test_throws(DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, segments, hierarchy)
+    @test_throws _test_throws("InconsistentOrientationError: The orientation of the boundary curve with index 7 should be positive, but it is negative. You may be able to fix this by passing the curve as reverse(reverse.(curve)).", DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test_throws _test_throws(DT.InconsistentOrientationError) triangulate(points; boundary_nodes, predicates=rt())
 end
 
@@ -369,7 +374,8 @@ end
     enricher_I = DT.BoundaryEnricher(points_I, curve_I)
     points, boundary_nodes = get_points(enricher_I), get_boundary_nodes(enricher_I)
     hierarchy = DT.get_polygon_hierarchy(enricher_I)
-    @test DT.check_args(points, boundary_nodes, hierarchy, DT.get_boundary_curves(enricher_I))
+    segments = nothing
+    @test DT.check_args(points, boundary_nodes, segments, hierarchy, DT.get_boundary_curves(enricher_I))
     @test DT.check_args(enricher_I)
 
     curve_II = [[1, 2, 3, 4, 5], [5, 6, 7, 8, 9], [9, 10, 11, 1]]
@@ -381,7 +387,8 @@ end
     enricher_II = DT.BoundaryEnricher(points_II, curve_II)
     points, boundary_nodes = get_points(enricher_II), get_boundary_nodes(enricher_II)
     hierarchy = DT.get_polygon_hierarchy(enricher_II)
-    @test DT.check_args(points, boundary_nodes, hierarchy, DT.get_boundary_curves(enricher_II))
+    segments = nothing
+    @test DT.check_args(points, boundary_nodes, segments, hierarchy, DT.get_boundary_curves(enricher_II))
     @test DT.check_args(enricher_II)
 
     curve_III = [[[1, 2, 3, 4, 5], [5, 6, 7, 8, 9], [9, 10, 11, 1]], [[15, 14, 13, 12], [12, 15]]]
@@ -394,7 +401,8 @@ end
     enricher_III = DT.BoundaryEnricher(points_III, curve_III)
     points, boundary_nodes = get_points(enricher_III), get_boundary_nodes(enricher_III)
     hierarchy = DT.get_polygon_hierarchy(enricher_III)
-    @test DT.check_args(points, boundary_nodes, hierarchy, DT.get_boundary_curves(enricher_III))
+    segments = nothing
+    @test DT.check_args(points, boundary_nodes, segments, hierarchy, DT.get_boundary_curves(enricher_III))
     @test DT.check_args(enricher_III)
 
     curve_IV = [CircularArc((1.0, 0.0), (1.0, 0.0), (0.0, 0.0))]
@@ -402,16 +410,18 @@ end
     enricher_IV = DT.BoundaryEnricher(points_IV, curve_IV)
     points, boundary_nodes = get_points(enricher_IV), get_boundary_nodes(enricher_IV)
     hierarchy = DT.get_polygon_hierarchy(enricher_IV)
-    @test DT.check_args(points, boundary_nodes, hierarchy)
+    segments = nothing
+    @test DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test DT.check_args(enricher_IV)
     curve_IV = [CircularArc((1.0, 0.0), (1.0, 0.0), (0.0, 0.0), positive=false)]
     points_IV = NTuple{2,Float64}[]
     enricher_IV = DT.BoundaryEnricher(points_IV, curve_IV)
     points, boundary_nodes = get_points(enricher_IV), get_boundary_nodes(enricher_IV)
     hierarchy = DT.get_polygon_hierarchy(enricher_IV)
-    @test_throws _test_throws(DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, hierarchy, DT.get_boundary_curves(enricher_IV))
+    segments = nothing
+    @test_throws _test_throws(DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, segments, hierarchy, DT.get_boundary_curves(enricher_IV))
     str = "If 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."
-    @test_throws _test_throws("InconsistentOrientationError: The orientation of the boundary curve with index 1 should be positive, but it is negative. You may be able to fix this by passing the curve as reverse(curve).\n$str", DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, hierarchy, DT.get_boundary_curves(enricher_IV))
+    @test_throws _test_throws("InconsistentOrientationError: The orientation of the boundary curve with index 1 should be positive, but it is negative. You may be able to fix this by passing the curve as reverse(curve).\n$str", DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, segments, hierarchy, DT.get_boundary_curves(enricher_IV))
     @test_throws _test_throws(DT.InconsistentOrientationError) DT.check_args(enricher_IV)
     @test_throws _test_throws("InconsistentOrientationError: The orientation of the boundary curve with index 1 should be positive, but it is negative. You may be able to fix this by passing the curve as reverse(curve).\n$str", DT.InconsistentOrientationError) DT.triangulate(NTuple{2,Float64}[]; boundary_nodes=[CircularArc((1.0, 0.0), (1.0, 0.0), (0.0, 0.0), positive=false)], predicates=rt())
 
@@ -420,7 +430,8 @@ end
     enricher_V = DT.BoundaryEnricher(points_V, curve_V)
     points, boundary_nodes = get_points(enricher_V), get_boundary_nodes(enricher_V)
     hierarchy = DT.get_polygon_hierarchy(enricher_V)
-    @test DT.check_args(points, boundary_nodes, hierarchy)
+    segments = nothing
+    @test DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test DT.check_args(enricher_V)
 
     curve_VI = [
@@ -432,7 +443,8 @@ end
     enricher_VI = DT.BoundaryEnricher(points_VI, curve_VI)
     points, boundary_nodes = get_points(enricher_VI), get_boundary_nodes(enricher_VI)
     hierarchy = DT.get_polygon_hierarchy(enricher_VI)
-    @test DT.check_args(points, boundary_nodes, hierarchy)
+    segments = nothing
+    @test DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test DT.check_args(enricher_VI)
 
     curve_VII = [
@@ -443,7 +455,8 @@ end
     enricher_VII = DT.BoundaryEnricher(points_VII, curve_VII)
     points, boundary_nodes = get_points(enricher_VII), get_boundary_nodes(enricher_VII)
     hierarchy = DT.get_polygon_hierarchy(enricher_VII)
-    @test DT.check_args(points, boundary_nodes, hierarchy)
+    segments = nothing
+    @test DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test DT.check_args(enricher_VII)
 
     curve_VIII = [
@@ -459,7 +472,8 @@ end
     enricher_VIII = DT.BoundaryEnricher(points_VIII, curve_VIII)
     points, boundary_nodes = get_points(enricher_VIII), get_boundary_nodes(enricher_VIII)
     hierarchy = DT.get_polygon_hierarchy(enricher_VIII)
-    @test DT.check_args(points, boundary_nodes, hierarchy)
+    segments = nothing
+    @test DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test DT.check_args(enricher_VIII)
 
     curve_IX =
@@ -475,7 +489,8 @@ end
     enricher_IX = DT.BoundaryEnricher(points_IX, curve_IX)
     points, boundary_nodes = get_points(enricher_IX), get_boundary_nodes(enricher_IX)
     hierarchy = DT.get_polygon_hierarchy(enricher_IX)
-    @test DT.check_args(points, boundary_nodes, hierarchy)
+    segments = nothing
+    @test DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test DT.check_args(enricher_IX)
 
     curve_X = [
@@ -492,7 +507,8 @@ end
     enricher_X = DT.BoundaryEnricher(points_X, curve_X)
     points, boundary_nodes = get_points(enricher_X), get_boundary_nodes(enricher_X)
     hierarchy = DT.get_polygon_hierarchy(enricher_X)
-    @test DT.check_args(points, boundary_nodes, hierarchy)
+    segments = nothing
+    @test DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test DT.check_args(enricher_X)
 
     curve_XI = [
@@ -519,7 +535,8 @@ end
     enricher_XI = DT.BoundaryEnricher(points_XI, curve_XI)
     points, boundary_nodes = get_points(enricher_XI), get_boundary_nodes(enricher_XI)
     hierarchy = DT.get_polygon_hierarchy(enricher_XI)
-    @test DT.check_args(points, boundary_nodes, hierarchy)
+    segments = nothing
+    @test DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test DT.check_args(enricher_XI)
     curve_XI = [
         [
@@ -545,7 +562,8 @@ end
     enricher_XI = DT.BoundaryEnricher(points_XI, curve_XI)
     points, boundary_nodes = get_points(enricher_XI), get_boundary_nodes(enricher_XI)
     hierarchy = DT.get_polygon_hierarchy(enricher_XI)
-    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, hierarchy)
+    segments = nothing
+    @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test_throws _test_throws(DT.InconsistentConnectionError) DT.check_args(enricher_XI)
     curve_XI = [
         [
@@ -573,8 +591,9 @@ end
     enricher_XI = DT.BoundaryEnricher(points_XI, curve_XI)
     points, boundary_nodes = get_points(enricher_XI), get_boundary_nodes(enricher_XI)
     hierarchy = DT.get_polygon_hierarchy(enricher_XI)
-    @test_throws _test_throws(DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, hierarchy, DT.get_boundary_curves(enricher_XI))
-    @test_throws _test_throws("InconsistentOrientationError: The orientation of the boundary curve with index 4 should be positive, but it is negative. You may be able to fix this by passing the curve as reverse(reverse.(curve)).\n$str", DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, hierarchy, DT.get_boundary_curves(enricher_XI))
+    segments = nothing
+    @test_throws _test_throws(DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, segments, hierarchy, DT.get_boundary_curves(enricher_XI))
+    @test_throws _test_throws("InconsistentOrientationError: The orientation of the boundary curve with index 4 should be positive, but it is negative. You may be able to fix this by passing the curve as reverse(reverse.(curve)).\n$str", DT.InconsistentOrientationError) DT.check_args(points, boundary_nodes, segments, hierarchy, DT.get_boundary_curves(enricher_XI))
     @test_throws _test_throws(DT.InconsistentOrientationError) DT.check_args(enricher_XI)
     @test_throws _test_throws("InconsistentOrientationError: The orientation of the boundary curve with index 4 should be positive, but it is negative. You may be able to fix this by passing the curve as reverse(reverse.(curve)).\n$str", DT.InconsistentOrientationError) triangulate(_points_XI; boundary_nodes=_curve_XI, predicates=rt())
 
@@ -594,7 +613,8 @@ end
     enricher_XII = DT.BoundaryEnricher(points_XII, curve_XII)
     points, boundary_nodes = get_points(enricher_XII), get_boundary_nodes(enricher_XII)
     hierarchy = DT.get_polygon_hierarchy(enricher_XII)
-    @test DT.check_args(points, boundary_nodes, hierarchy)
+    segments = nothing
+    @test DT.check_args(points, boundary_nodes, segments, hierarchy)
     @test DT.check_args(enricher_XII)
 end
 
@@ -604,4 +624,24 @@ end
     @test !DT.WARN_ON_DUPES[]
     DT.toggle_warn_on_dupes!()
     @test DT.WARN_ON_DUPES[]
+end
+
+@testset "Duplicate points appearing in boundary nodes / segments" begin
+    points = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.0, 0.0)]
+    hierarchy = DT.construct_polygon_hierarchy(points)
+    boundary_nodes = [1, 2, 3, 4, 5]
+    segments = nothing
+    DT.check_args(points, boundary_nodes, segments, hierarchy)
+    @test boundary_nodes == [1, 2, 3, 4, 1]
+    boundary_nodes = [[1, 2], [2, 3], [3, 4], [4, 5]]
+    DT.check_args(points, boundary_nodes, segments, hierarchy)
+    @test boundary_nodes == [[1, 2], [2, 3], [3, 4], [4, 1]]
+    boundary_nodes = [[[1, 2, 3, 4, 5]]]
+    DT.check_args(points, boundary_nodes, segments, hierarchy)
+    @test boundary_nodes == [[[1, 2, 3, 4, 1]]]
+    boundary_nodes = [[1, 2, 3, 4, 5]]
+    segments = Set(((1, 2), (4, 5)))
+    DT.check_args(points, boundary_nodes, segments, hierarchy)
+    @test boundary_nodes == [[1, 2, 3, 4, 1]]
+    @test segments == Set(((1, 2), (4, 1)))
 end
\ No newline at end of file