Triangulations are now stored in their own data structures (preparing for refinement)

Author: DanielVandH

Project.toml

@@ -1,7 +1,7 @@
 name = "DelaunayTriangulation"
 uuid = "927a84f5-c5f4-47a5-9785-b46e178433df"
 authors = ["Daniel VandenHeuvel <[email protected]>"]
-version = "0.2.0"
+version = "0.3.0"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"

README.md

@@ -38,7 +38,8 @@ p4 = [-1.0, 2.0]
 p5 = [4.0, 2.0]
 p6 = [-2.0, -1.0]
 pts = [p1, p2, p3, p4, p5, p6]
-T, adj, adj2v, DG, HG = triangulate_berg(pts)
+tri, HG = triangulate_berg(pts)
+T, adj, adj2v, DG = tri # tri is a Triangulation struct, see ?Triangulation
 ```
 This function `triangulate_berg` creates the triangulation for the given `pts`, returning:
 - `T`: This is the set of triangles, of default type `Set{NTuple{3, Int64}}`. By default, triangles are represented as tuples of indices `(i, j, k)` so that `(pts[i], pts[j], pts[k])` are the points representing the triangle. All triangles `T` are positively oriented.
@@ -166,7 +167,9 @@ y5 = [0.8815, 0.8056, 0.80268, 0.73258, 0.6, 0.598, 0.5777, 0.525, 0.4346, 0.364
 x = [x1, x2, x3, x4, x5]
 y = [y1, y2, y3, y4, y5]
 ## Mesh 
-T, adj, adj2v, DG, pts, BN = generate_mesh(x, y, 0.05; gmsh_path=GMSH_PATH)
+tri, BN = generate_mesh(x, y, 0.05; gmsh_path=GMSH_PATH)
+T = tri.triangles 
+pts = tri.points
 ## Visualise 
 fig = Figure()
 ax = Axis(fig[1, 1])
@@ -358,7 +361,7 @@ This all applies equally as well to `triangulate_berg`.
 
 ```julia
 Random.seed!(_seed) # Keep the same insertion order for later
-T, adj, adj2v, DG = triangulate_berg(pts; 
+(T, adj, adj2v, DG), _ = triangulate_berg(pts; 
     IntegerType = Int32,
     EdgeType = CustomEdge, 
     TriangleType = CustomTriangle, 

src/DelaunayTriangulation.jl

@@ -44,6 +44,8 @@ export triangulate_bowyer
 export triangulate_berg
 export triangulate_structured
 export triangulate
+export Triangulation 
+
 #=
 export voronoi 
 export VoronoiTessellation

src/bowyerwatson.jl

@@ -56,6 +56,7 @@ end
         randomise=true,
         trim=true,
         try_last_inserted_point=true,
+        trim_empty_features=true,
         skip_pts=Set{Int64}())
 
 Triangulates the set of points `pts` using the Bowyer-Watson algorithm.
@@ -72,13 +73,17 @@ Triangulates the set of points `pts` using the Bowyer-Watson algorithm.
 - `randomise=true`: Whether to randomise the insertion order.
 - `trim=true`: Whether to remove the ghost triangles at the end.
 - `try_last_inserted_point=true`: At each stage, this should be `true` if the point that was last inserted should also be tried for initialising [`jump_and_march`](@ref).
+- `trim_empty_features = true`: Whether to remove keys from the `Adjacent` map that refer to deleted edges, and edges from the `DelaunayGraph` that refer to deleted edges.
 - `skip_pts`: There may be points you want to avoid adding into the triangulation, but they still exist in `pts`. If this is the case, add the indices for the points into a `Set` and set this keyword to this set.
 
 # Outputs 
+The output is a `Triangulation` struct, containing
+
 - `T`: The set of triangles. 
 - `adj`: The [`Adjacent`](@ref) map.
 - `adj2v`: The [`Adjacent2Vertex`](@ref) map.
 - `DG`: The `[DelaunayGraph`](@ref).
+- `pts`: The provided point set.
 """
 function triangulate_bowyer(pts;
     IntegerType::Type{I}=Int64,
@@ -124,7 +129,7 @@ function triangulate_bowyer(pts;
     trim_empty_features && clear_empty_keys!(adj)
     trim_empty_features && clear_empty_points!(DG)
     trim && remove_ghost_triangles!(T, adj, adj2v, DG)
-    return T, adj, adj2v, DG
+    return Triangulation(T, adj, adj2v, DG, pts)
 end
 
 function lazy_triangulate_bowyer(pts;

src/data_structures.jl

@@ -375,4 +375,50 @@ function add_edge!(G::HistoryGraph, T, V::Vararg{Tri,N}) where {Tri,N}
         add_edge!(G, T, V)
     end
     return nothing
+end
+
+###################################################
+#/
+#/
+#/ Triangulation 
+#/
+#/
+################################################### 
+"""
+    Triangulation{T,I,E,Es}
+
+Data structure for a triangulation.
+
+# Fields 
+- `triangles::T`: Structure containing the collection of triangles.
+- `adjacent::Adjacent{I,E}`: The adjacent map.
+- `adjacent2vertex::Adjacent2Vertex{I,Es,E}`: The adjacent-to-vertex map. 
+- `graph::DelaunayGraph{I}`: The Delaunay graph. 
+- `points::P`: The point set.
+"""
+struct Triangulation{T,I,E,Es,P}
+    triangles::T
+    adjacent::Adjacent{I,E}
+    adjacent2vertex::Adjacent2Vertex{I,Es,E}
+    graph::DelaunayGraph{I}
+    points::P
+end
+get_triangles(tri::Triangulation) = tri.triangles
+get_adjacent(tri::Triangulation) = tri.adjacent
+get_adjacent2vertex(tri::Triangulation) = tri.adjacent2vertex
+get_graph(tri::Triangulation) = tri.graph
+get_points(tri::Triangulation) = tri.points
+
+function Base.iterate(tri::Triangulation, state=1)
+    if state == 1
+        return (get_triangles(tri), 2)
+    elseif state == 2
+        return (get_adjacent(tri), 3)
+    elseif state == 3
+        return (get_adjacent2vertex(tri), 4)
+    elseif state == 4
+        return (get_graph(tri), 5)
+    elseif state == 5
+        return (get_points(tri), nothing)
+    end
 end

src/deberg.jl

@@ -23,6 +23,45 @@ function add_point_berg!(T, adj::Adjacent{I,E}, adj2v, DG, HG::HistoryGraph{Tri}
     return nothing
 end
 
+"""
+    triangulate_berg(pts;
+        IntegerType::Type{I}=Int64,
+        EdgeType::Type{E}=NTuple{2,IntegerType},
+        TriangleType::Type{V}=NTuple{3,IntegerType},
+        EdgesType::Type{Es}=Set{EdgeType},
+        TrianglesType::Type{Ts}=Set{TriangleType},
+        randomise=true,
+        trim=true,
+        trim_empty_features=true,
+        skip_pts=Set{Int64}()) where {I,E,V,Es,Ts}
+
+Triangulates the set of points `pts` using the de Berg's algorithm.
+
+# Inputs 
+- `pts`: The point set.
+
+# Keyword Arguments 
+- `IntegerType::Type{I}=Int64`: Type used to represent integers. 
+- `EdgeType::Type{E}=NTuple{2,IntegerType}`: Type used to represent edges. 
+- `TriangleType::Type{V}=NTuple{3,IntegerType}`: Type used to represent triangles. 
+- `EdgesType::Type{Es}=Set{EdgeType}`: Type used to represent a collection of edges.
+- `TrianglesType::Type{Ts}=Set{TriangleType}`: Type used to represent a collection of triangles. 
+- `randomise=true`: Whether to randomise the insertion order.
+- `trim=true`: Whether to remove the ghost triangles at the end.
+- `trim_empty_features = true`: Whether to remove keys from the `Adjacent` map that refer to deleted edges, and edges from the `DelaunayGraph` that refer to deleted edges.
+- `skip_pts`: There may be points you want to avoid adding into the triangulation, but they still exist in `pts`. If this is the case, add the indices for the points into a `Set` and set this keyword to this set.
+
+# Outputs 
+The output is a `Triangulation` struct, containing
+
+- `T`: The set of triangles. 
+- `adj`: The [`Adjacent`](@ref) map.
+- `adj2v`: The [`Adjacent2Vertex`](@ref) map.
+- `DG`: The `[DelaunayGraph`](@ref).
+- `pts`: The provided point set.
+
+The second output is `HG`, a [`HistoryGraph`](@ref).
+"""
 function triangulate_berg(pts;
     IntegerType::Type{I}=Int64,
     EdgeType::Type{E}=NTuple{2,IntegerType},
@@ -48,7 +87,7 @@ function triangulate_berg(pts;
     trim_empty_features && clear_empty_keys!(adj)
     trim_empty_features && clear_empty_points!(DG)
     trim && remove_bounding_triangle!(T, adj, adj2v, DG)
-    return T, adj, adj2v, DG, HG
+    return Triangulation(T, adj, adj2v, DG, pts), HG
 end
 
 function remove_bounding_triangle!(T::Ts, adj::Adjacent{I,E}, adj2v, DG) where {I,E,Ts}

src/gmsh.jl

@@ -142,7 +142,7 @@ function generate_mesh end
         end
     end
     T, adj, adj2v, DG = triangulate(elements, nodes, boundary_nodes)
-    return T, adj, adj2v, DG, nodes, unique.(boundary_nodes)
+    return Triangulation(T, adj, adj2v, DG, nodes), unique.(boundary_nodes)
 end
 @inline function generate_mesh(x::Vector{LinRange{Float64,Int64}}, y::Vector{LinRange{Float64,Int64}}, ref; mesh_algorithm=6, verbosity=0, gmsh_path="./gmsh-4.9.4-Windows64/gmsh.exe", num_threads=Base.Threads.nthreads())
     return generate_mesh(collect.(x), collect.(y), ref; verbosity, mesh_algorithm, gmsh_path, num_threads)

src/structured.jl

@@ -69,7 +69,7 @@ function triangulate_structured(a, b, c, d, Nˣ, Nʸ; return_boundary_types=fals
         BNy = Vector{Vector{Float64}}([[], [], [], []])
         if single_boundary
             BN = ch_idx
-            return T, adj, adj2v, DG, points, BN
+            return Triangulation(T, adj, adj2v, DG, points), BN
         end
         for i in ch_idx
             if points[1, i] == a
@@ -104,8 +104,8 @@ function triangulate_structured(a, b, c, d, Nˣ, Nʸ; return_boundary_types=fals
         BN[2] .= BN[2][idx₂]
         BN[3] .= BN[3][idx₃]
         BN[4] .= BN[4][idx₄]
-        return T, adj, adj2v, DG, points, BN
+        return Triangulation(T, adj, adj2v, DG, points), BN
     else
-        return T, adj, adj2v, DG, points
+        return Triangulation(T, adj, adj2v, DG, points)
     end
 end

src/utils.jl

@@ -467,7 +467,7 @@ set `pts` to those from de Berg's method. If both are identical, return `true`,
 out-of-place.
 """
 function compare_deberg_to_bowyerwatson(T, adj, adj2v, DG, pts)
-    _T, _adj, _adj2v, _DG, _ = triangulate_berg(pts)
+    (_T, _adj, _adj2v, _DG), _ = triangulate_berg(pts)
     return compare_deberg_to_bowyerwatson(T, adj, adj2v, DG, _T, _adj, _adj2v, _DG)
 end
 

test/data_structures.jl

@@ -210,7 +210,7 @@ end
     @test collect(DT.edges(adj)) == [(2, 3), (5, 7)] || collect(DT.edges(adj)) == [(5, 7), (2, 3)]
 
     @inferred DT.get_edge(adj, 2, 3)
-    @inferred DT.get_edge(adj, (2,3))
+    @inferred DT.get_edge(adj, (2, 3))
     @inferred DT.get_edge(adj, -3, 50)
 end
 
@@ -306,7 +306,7 @@ end
     @test DT.get_neighbour(DG, 3) == Set([1])
     @test edges(DG) == DG.graph.E
 
-    pts = rand(SVector{2, Float64}, 250)
+    pts = rand(SVector{2,Float64}, 250)
     T, adj, adj2v, DG = DT.triangulate_bowyer(pts)
     for i in eachindex(pts)
         @test deg(graph(DG), i) == DT.num_neighbours(DG, i)
@@ -377,4 +377,71 @@ end
     DT.add_edge!(HH, T₁, T₂)
     DT.add_edge!(HH, T₁, T₃)
     @test DT.graph(H) == DT.graph(HH)
+end
+
+###################################################
+#/
+#/
+#/ Storing a Delaunay triangulation 
+#/
+#/
+###################################################
+@testset "Testing from the Bowyer-Watson algorithm" begin
+    pts = rand(2, 500)
+    Random.seed!(2929292)
+    T, adj, adj2v, DG = triangulate_bowyer(pts)
+    tri = Triangulation(T, adj, adj2v, DG, pts)
+
+    @test DT.get_triangles(tri) === T
+    @test DT.get_adjacent(tri) === adj
+    @test DT.get_adjacent2vertex(tri) === adj2v
+    @test DT.get_graph(tri) === DG
+    @test DT.get_points(tri) === pts
+
+    _T, _adj, _adj2v, _DG, _pts = tri
+    @test _T === T
+    @test _adj === adj
+    @test _adj2v === adj2v
+    @test _DG === DG
+    @test _pts === pts
+
+    Random.seed!(2929292)
+    _tri = triangulate_bowyer(pts)
+    @test _tri isa Triangulation
+    @test _tri.triangles == T
+    @test _tri.adjacent.adjacent == adj.adjacent
+    @test _tri.adjacent2vertex.adjacent2vertex == adj2v.adjacent2vertex
+    @test _tri.graph.graph == DG.graph
+    @test _tri.points == pts
+end
+
+@testset "Testing from de Berg's method" begin
+    pts = rand(2, 500)
+    Random.seed!(29292222292)
+    (T, adj, adj2v, DG), HG = triangulate_berg(pts)
+    @test HG isa DT.HistoryGraph
+    tri = Triangulation(T, adj, adj2v, DG, pts)
+
+    @test DT.get_triangles(tri) === T
+    @test DT.get_adjacent(tri) === adj
+    @test DT.get_adjacent2vertex(tri) === adj2v
+    @test DT.get_graph(tri) === DG
+    @test DT.get_points(tri) === pts
+
+    _T, _adj, _adj2v, _DG, _pts = tri
+    @test _T === T
+    @test _adj === adj
+    @test _adj2v === adj2v
+    @test _DG === DG
+    @test _pts === pts
+
+    Random.seed!(29292222292)
+    _tri, _HG = triangulate_berg(pts)
+    @test _tri isa Triangulation
+    @test _tri.triangles == T
+    @test _tri.adjacent.adjacent == adj.adjacent
+    @test _tri.adjacent2vertex.adjacent2vertex == adj2v.adjacent2vertex
+    @test _tri.graph.graph == DG.graph
+    @test _tri.points == pts
+    @test HG.graph == _HG.graph
 end