Unconstrained Triangulations

It is simple to construct unconstrained triangulations. The method for this is through the triangulate function, shown below.

DelaunayTriangulation.triangulateFunction
triangulate(points::P; edges=nothing, boundary_nodes=nothing,
    IntegerType::Type{I}=Int,
    EdgeType::Type{E}=NTuple{2,IntegerType},
    TriangleType::Type{V}=NTuple{3,IntegerType},
    EdgesType::Type{Es}=Set{EdgeType},
    TrianglesType::Type{Ts}=Set{TriangleType},
    randomise=true,
    delete_ghosts=true,
    delete_empty_features=true,
    try_last_inserted_point=true,
    skip_points=Set{IntegerType}(),
    num_sample_rule::M=default_num_samples,
    rng::AbstractRNG=Random.default_rng(),
    point_order=get_point_order(points, randomise, skip_points, IntegerType, rng),
    recompute_representative_point=true,
    delete_holes=true,
    check_arguments=true
) where {P,I,E,V,Es,Ts,M}

Computes the unconstrained Delaunay triangulation of a set of points. If edges is provided, they will be inserted. If boundary_nodes is provided, a boundary will also be made from these nodes, with all triangles inside the boundaries deleted.

Arguments

  • points::P: The set of points to compute the triangulation of.

Keyword Arguments

  • edges=nothing: Any constrained edges to insert. If nothing, an unconstrained triangulation is built. The constrained edges should not intersect each other, and they should not cross over boundary edges.
  • boundary_nodes=nothing: Any boundaries to define. The specification of these boundary nodes is outlined in the boundary handling section of the documentation. All triangles away from a defined boundary are deleted if delete_holes.
  • IntegerType::Type{I}=Int: The integer type to use for indexing.
  • EdgeType::Type{E}=NTuple{2,IntegerType}: The type to use for representing edges.
  • TriangleType::Type{V}=NTuple{3,IntegerType}: The type to use for representing triangles.
  • EdgesType::Type{Es}=Set{EdgeType}: The type to use for representing collections of edges.
  • TrianglesType::Type{Ts}=Set{TriangleType}: The type to use for representing collections of triangles.
  • randomise=true: Whether to randomise the insertion order.
  • delete_ghosts=true: Whether to remove the ghost triangles at the end of the triangulation.
  • delete_empty_features=true: Whether to delete any empty neighbourhoods and adjacencies at the end of the triangulation.
  • try_last_inserted_point=true: When finding the next point, this decides if the previously inserted point should also be attempted.
  • skip_points=Set{IntegerType}(): Points to skip over when triangulationg, i.e. points to not include in the triangulation.
  • num_sample_rule::M=default_num_samples: A function of the form n -> Number, with n the number of points currently in the triangulation, that returns the number of points to sample during the point location steps.
  • rng::AbstractRNG=Random.default_rng(): The RNG to use.
  • point_order=get_point_order(points, randomise, skip_points, IntegerType, rng): The insertion order.
  • recompute_representative_point=true: At the end of the triangulation, will recompute the RepresentativePointList if true.
  • delete_holes=true: Whether to delete the exterior faces of all boundaries. There may be issues if you have boundary nodes but have this set to false - this kwarg is mostly for debugging.
  • check_arguments=true: Whether to check the arguments for validity.

Outputs

Returns a Triangulation.

In the code below, we give an example, and show how we can plot the result.

using DelaunayTriangulation, CairoMakie 
a = [1.5, 4.0]
b = [0.0, 3.5]
c = [2.0, 1.5]
d = [3.0, 2.5]
e = [2.5, 3.5]
f = [0.5, 3.0]
g = [2.5, -2.0]
h = [0.5, 1.5]
i = [0.0, 0.5]
j = [1.5, 3.0]
pts = [a, b, c, d, e, f, g, h, i, j]
tri = triangulate(pts)
fig, ax, sc = triplot(tri)
Triangulation

This object tri is a Triangulation.

julia> tri
Delaunay Triangulation.
    Constrained: false
    Has ghost triangles: false
    Number of points: 10
    Number of triangles: 12
    Number of edges: 27

As we describe in more detail in the data structures section in the sidebar, tri has several fields:

julia> propertynames(tri)
(:points, :triangles, :adjacent, :adjacent2vertex, :graph, :boundary_nodes, :boundary_edge_map, :boundary_map, :boundary_index_ranges, :constrained_edges, :all_constrained_edges, :convex_hull, :representative_point_list)

We explain each field below.

  • tri.points: This stores pts.

  • tri.triangles: This stores all the triangles. In this case,

julia> get_triangles(tri)
Set{Tuple{Int, Int, Int}} with 12 elements:
  (10, 5, 1)
  (9, 7, 3)
  (2, 6, 1)
  (3, 10, 8)
  (10, 4, 5)
  (9, 3, 8)
  (3, 4, 10)
  (8, 6, 2)
  (9, 8, 2)
  (10, 1, 6)
  (8, 10, 6)
  (3, 7, 4)

More generally, you can iterate over these triangles via each_triangle(tri). For example, the area of the triangulation could be computed as follows:

triangle_area(p, q, r) = 0.5 * (p[1] * q[2] + q[1] * r[2] + r[1] * p[2] - p[1] * r[2] - r[1] * q[2] - q[1] * p[2])
A = 0.0
for T in each_triangle(tri)
    i, j, k = indices(T)
    p, q, r = get_point(tri, i, j, k)
    A += triangle_area(p, q, r)
end 
  • tri.adjacent: This stores the adjacency relationships of the triangulation, mapping edges (u, v) to a vertex w so that (u, v, w) is a positively oriented triangle in tri. In this case, we have
julia> get_adjacent(tri)
Adjacent{Int, Tuple{Int, Int}}, with map:
DataStructures.DefaultDict{Tuple{Int, Int}, Int, Int} with 43 entries:
  (9, 3)  => 8
  (8, 9)  => 3
  (4, 7)  => -1
  (2, 1)  => -1
  (10, 1) => 6
  (2, 8)  => 6
  (10, 8) => 3
  (3, 9)  => 7
  (4, 5)  => 10
  (8, 3)  => 10
  (9, 8)  => 2
  ⋮       => ⋮

This object is iterable, allowing for you to do e.g.

for (uv, w) in get_adjacent(tri)
    u = initial(uv)
    v = terminal(uv)
    ...
end
  • tri.adjacent2vertex: This is a map that returns, given an index i, all other edges (j, k) such that (i, j, k) is a positively oriented triangle in the triangulation. In this case, we have
julia> get_adjacent2vertex(tri)
Adjacent{Int, Set{Tuple{Int, Int}}, Tuple{Int, Int}}, with map:
Dict{Int, Set{Tuple{Int, Int}}} with 11 entries:
  5  => Set([(10, 4), (1, 10)])
  8  => Set([(9, 3), (2, 9), (6, 2), (3, 10), (10, 6)])
  1  => Set([(6, 10), (10, 5), (2, 6)])
  6  => Set([(1, 2), (2, 8), (10, 1), (8, 10)])
  -1 => Set([(7, 9), (4, 7), (2, 1), (9, 2), (5, 4), (1, 5)])
  9  => Set([(7, 3), (3, 8), (8, 2)])
  3  => Set([(7, 4), (4, 10), (10, 8), (8, 9), (9, 7)])
  7  => Set([(3, 9), (4, 3)])
  4  => Set([(5, 10), (3, 7), (10, 3)])
  2  => Set([(8, 6), (9, 8), (6, 1)])
  10 => Set([(4, 5), (6, 8), (8, 3), (5, 1), (1, 6), (3, 4)])

This object is iterable, allowing for you to do e.g.

for (w, S) in get_adjacent2vertex(w)
    for (u, v) in S 
        ...
    end
end 
  • tri.graph: This is a graph that returns, given an index i, all other indices j such that (i, j) is an edge in the triangulation. In this case, we have
julia> get_graph(tri)
Graph
    Number of edges: 27
    Number of vertices: 11

julia> get_edges(tri)
Set{Tuple{Int, Int}} with 27 elements:
  (2, 9)
  (4, 5)
  (1, 2)
  (6, 8)
  (6, 10)
  (3, 7)
  (-1, 2)
  (4, 7)
  (3, 4)
  (1, 5)
  (-1, 9)
  (4, 10)
  (2, 8)
  (-1, 5)
  (1, 6)
  (3, 9)
  (7, 9)
  ⋮

julia> get_neighbours(tri)
Dict{Int, Set{Int}} with 11 entries:
  5  => Set([4, -1, 10, 1])
  8  => Set([6, 2, 10, 9, 3])
  1  => Set([5, 6, 2, 10, -1])
  6  => Set([2, 10, 8, 1])
  -1 => Set([5, 4, 7, 2, 9, 1])
  9  => Set([7, 2, -1, 8, 3])
  3  => Set([4, 7, 10, 9, 8])
  7  => Set([4, -1, 9, 3])
  4  => Set([5, 7, -1, 10, 3])
  2  => Set([6, -1, 9, 8, 1])
  10 => Set([5, 4, 6, 8, 3, 1])
  • tri.boundary_nodes: This is a list of all fixed boundary nodes in the triangulation. In our case, we have none. See the Gmsh section for an example. The actual nodes on the boundary in this case can be obtained via tri.convex_hull.

  • tri.boundary_edge_map: This is a Dict that maps all boundary edges to their position in tri.boundary_nodes. See the Gmsh section for an example.

  • tri.boundary_map: This would be a list mapping boundary indices to all the fixed boundary nodes in tri.boundary_nodes corresponding to that index. This map is empty in this case as we have no fixed boundary nodes,, but see the Gmsh section for an example.

  • tri.boundary_index_ranges: This is be a list mapping indices of boundary curves to all boundary indices belonging to that curve. In this case, we have

julia> get_boundary_index_ranges(tri)
OrderedCollections.OrderedDict{Int, UnitRange{Int}} with 1 entry:
  -1 => -1:-1

This tells us that whenever we see a -1 as a vertex, we have a ghost vertex corresponding to the outer curve, so e.g. if get_adjacent(tri, u, v) == -1, then (u, v) is an edge on the boundary. A better example is in the Gmsh section.

  • tri.constrained_edges: This would be the collection of constrained edges if we had any. See the constrained trangulation section.

  • tri.all_constrained_edges: This is a collection of all constrained edges currently in the triangulation, including the boundary edges. We have none here, but see the constrained triangulation section.

  • tri.convex_hull: This is the ConvexHull of tri.points. In this case, we have

julia> get_convex_hull(tri)
Convex hull.
    Indices:
7-element Vector{Int}:
 7
 4
 5
 1
 2
 9
 7
  • tri.representative_point_list: This is the Dict that maps curve indices to the coordinate used for representing corresponding boundary indices. Typically, these points are near the centroid of the curve; see the pole_of_inaccessibility function. In our case,
julia> tri.representative_point_list
Dict{Int, DelaunayTriangulation.RepresentativeCoordinates{Int, Float64}} with 1 entry:
  1 => RepresentativeCoordinates{Int, Float64}(1.51155, 1.43234, 0)