Unconstrained Triangulations

Our first example considers constructing unconstrained Delaunay triangulations. To start, let us load in the packages we will need.

using DelaunayTriangulation
using CairoMakie # for plotting
using StableRNGs # for reproducibility
using LinearAlgebra # used for computing norms later

We consider just triangulating a random set of points. First, generating the points:

rng = StableRNG(123)
points = rand(rng, 2, 500) # just do rand(2, 500) if you are not concerned about the RNG
2×500 Matrix{Float64}:
 0.181026  0.669058   0.437991  0.763974  …  0.459817   0.312727  0.0634304
 0.36749   0.0427306  0.48329   0.924483     0.0794176  0.892313  0.606559

We now triangulate these points by using triangulate. We pass the rng as a keyword argument, but again if you are not concerned about the RNG (or set the seed using Random.seed!) then you can ignore this.

tri = triangulate(points; rng = rng)
Delaunay Triangulation.
   Number of vertices: 500
   Number of triangles: 980
   Number of edges: 1479
   Has boundary nodes: false
   Has ghost triangles: true
   Curve-bounded: false
   Weighted: false
   Constrained: false
fig, ax, sc = triplot(tri)
fig
Example block output

This tri object is our Triangulation, and we can interact with it in many ways.

Iterating over vertices

For example, we can iterate over the points in the triangulation using each_solid_vertex. Here we compute the centroid of the point cloud:

function compute_centroid(tri)
    s = [0.0, 0.0]
    for i in each_solid_vertex(tri)
        p = get_point(tri, i)
        s .+= p
    end
    s ./= num_solid_vertices(tri)
    return s
end
s = compute_centroid(tri)
2-element Vector{Float64}:
 0.5106629621450203
 0.477077379328767

We need to use the solid identifier because triangulations are made up of both solid and ghost vertices/edges/triangles, for reasons described in the manual. If we just use each_vertex(tri), then we also find a vertex -1 that corresponds to the boundary. For example, the points on the boundary can be obtained using:

get_neighbours(tri, -1)
Set{Int64} with 18 elements:
  456
  176
  425
  258
  325
  398
  197
  219
  319
  262
  112
  345
  371
  232
  245
  163
  468
  140

One reason to be careful with this especially is that get_point(tri, -1) does actually correspond to a coordinate,

get_point(tri, -1)
(0.500547641525688, 0.497711096515429)

(This is the pole of inaccessibility for the domain; see here for more details.) You can use each_vertex or each_ghost_vertex to consider all types of vertices or only the ghost vertices. If you just want the vertices, use each_vertex(tri), which will also include the ghost vertex.

each_vertex(tri)
Set{Int64} with 501 elements:
  319
  185
  420
  365
  422
  263
  242
  183
  224
  177
  77
  172
  103
  59
  211
  358
  366
  414
  116
  ⋮ 

Iterating over edges

We can also iterate over the edges of the triangulation using each_solid_edge, or each_edge for both solid and ghost edges and each_ghost_edge for only the ghost edges. To give an example, here we compute the average length of an edge.

function compute_mean_edge_length(tri)
    ℓ = 0.0
    for e in each_solid_edge(tri)
        u, v = edge_vertices(e)
        p, q = get_point(tri, u, v)
        ℓ += norm(p .- q)
    end
    ℓ /= num_solid_edges(tri)
    return ℓ
end
ℓ = compute_mean_edge_length(tri)
0.05645693808033534

By default, the triangulation has ghost edges, so each_edge and each_solid_edge are not the same.

each_edge(tri)
Set{Tuple{Int64, Int64}} with 1497 elements:
  (453, 388)
  (276, 329)
  (17, 12)
  (382, 127)
  (4, 197)
  (366, 352)
  (34, 313)
  (407, 106)
  (114, 56)
  (194, 90)
  (263, 136)
  (498, 388)
  (345, 446)
  (343, 47)
  (375, 33)
  (412, 433)
  (292, 304)
  (290, 127)
  (274, 334)
  ⋮ 

Note also that the edges are all given as unordered, so the set of edges only includes one of (i, j) and (j, i) for each edge (i, j).

Iterating over triangles

Similarly, we can iterate over the triangles using each_solid_triangle, each_ghost_triangle, or each_triangle. By default, ghost triangles are included in the output. Here we compute the area of the domain by getting the area of each triangle.

area(p, q, r) = 0.5 * ((getx(q) - getx(p)) * (gety(r) - gety(p)) - (gety(q) - gety(p)) * (getx(r) - getx(p)))
function compute_triangulation_area(tri)
    A = 0.0
    for T in each_solid_triangle(tri)
        i, j, k = triangle_vertices(T)
        p, q, r = get_point(tri, i, j, k)
        A += area(p, q, r)
    end
    return A
end
A = compute_triangulation_area(tri)
0.9683345161527972

(You can compute areas like this using get_area(tri).) You can access the set of triangles using get_triangles(tri):

get_triangles(tri)
Set{Tuple{Int64, Int64, Int64}} with 998 elements:
  (59, 426, 44)
  (163, 365, 295)
  (383, 473, 297)
  (499, 260, 216)
  (185, 134, 472)
  (239, 43, 274)
  (300, 31, 288)
  (113, 135, 221)
  (461, 7, 430)
  (230, 257, 316)
  (57, 203, 80)
  (277, 284, 233)
  (74, 144, 474)
  (398, 484, 44)
  (288, 31, 483)
  (361, 368, 111)
  (450, 77, 207)
  (220, 327, 315)
  (151, 329, 276)
  ⋮ 

The triangles are all positively oriented, meaning the triangles are given such that the corresponding points are traversed in counter-clockwise order.

Querying neighbour information

We can query neighbour information at several levels.

Points

For a given point, there are two type of neighbours: The neighbouring vertices, and the neighbouring triangles. The neighbours can be obtained using get_neighbours. For example, the set of vertices that share an edge with the fifth vertex is:

get_neighbours(tri, 5)
Set{Int64} with 4 elements:
  93
  117
  301
  214

The set of triangles that share an edge with the fifth vertex is obtained using get_adjacent2vertex. This returns a set of edges (v, w) such that, for a given vertex u, (u, v, w) is a positively oriented triangle in the triangulation. For example,

get_adjacent2vertex(tri, 5)
Set{Tuple{Int64, Int64}} with 4 elements:
  (93, 117)
  (301, 214)
  (117, 301)
  (214, 93)

means that the triangles that contain 5 as a vertex are (5, 93, 117), (5, 117, 301), (5, 301, 214), and (5, 214, 93). We can verify this:

filter(T -> 5 ∈ triangle_vertices(T), get_triangles(tri))
Set{Tuple{Int64, Int64, Int64}} with 4 elements:
  (93, 5, 214)
  (301, 5, 117)
  (93, 117, 5)
  (301, 214, 5)

These queries can also be applied to the ghost vertices, in which information about the boundary is provided.

get_neighbours(tri, -1)
Set{Int64} with 18 elements:
  456
  176
  425
  258
  325
  398
  197
  219
  319
  262
  112
  345
  371
  232
  245
  163
  468
  140
get_adjacent2vertex(tri, -1)
Set{Tuple{Int64, Int64}} with 18 elements:
  (468, 232)
  (398, 258)
  (176, 371)
  (425, 456)
  (232, 163)
  (112, 262)
  (258, 140)
  (319, 176)
  (262, 219)
  (163, 325)
  (371, 425)
  (219, 197)
  (140, 245)
  (325, 398)
  (456, 112)
  (197, 345)
  (245, 319)
  (345, 468)

Edges

For a given edge (u, v), the relevant neighbours are the vertices that are next to it so that a triangle is formed. We can find the vertex w such that (u, v, w) is a positively oriented triangle in the triangulation using get_adjacent(tri, u, v). For example,

get_adjacent(tri, 163, 365)
295

means that (163, 365, 295) is a positively oriented triangle, as we can verify:

DelaunayTriangulation.triangle_orientation(tri, 163, 365, 295)
Certificate.PositivelyOriented = 6

(The representation of this predicate using a DelaunayTriangulation.Certificate is described in more detail in the manual.) The other triangle adjoining the unordered edge (u, v), meaning the oriented edge (v, u), is obtained similarly:

get_adjacent(tri, 365, 163)
232

If an edge (u, v) is on the boundary, oriented so that there is no solid vertex w such that (u, v, w) is a triangle in the triangulation, then get_adjacent(tri, u, v) returns the ghost vertex. For example,

get_adjacent(tri, 398, 258)
-1

means that (398, 258) is a boundary edge and (398, 258, -1) is a ghost triangle. You can test for this case using DelaunayTriangulation.is_boundary_edge:

DelaunayTriangulation.is_boundary_edge(tri, 258, 398)
true

Just the code

An uncommented version of this example is given below. You can view the source code for this file here.

using DelaunayTriangulation
using CairoMakie # for plotting
using StableRNGs # for reproducibility
using LinearAlgebra # used for computing norms later

rng = StableRNG(123)
points = rand(rng, 2, 500) # just do rand(2, 500) if you are not concerned about the RNG

tri = triangulate(points; rng = rng)

fig, ax, sc = triplot(tri)
fig

function compute_centroid(tri)
    s = [0.0, 0.0]
    for i in each_solid_vertex(tri)
        p = get_point(tri, i)
        s .+= p
    end
    s ./= num_solid_vertices(tri)
    return s
end
s = compute_centroid(tri)

get_neighbours(tri, -1)

get_point(tri, -1)

each_vertex(tri)

function compute_mean_edge_length(tri)
    ℓ = 0.0
    for e in each_solid_edge(tri)
        u, v = edge_vertices(e)
        p, q = get_point(tri, u, v)
        ℓ += norm(p .- q)
    end
    ℓ /= num_solid_edges(tri)
    return ℓ
end
ℓ = compute_mean_edge_length(tri)

each_edge(tri)

area(p, q, r) = 0.5 * ((getx(q) - getx(p)) * (gety(r) - gety(p)) - (gety(q) - gety(p)) * (getx(r) - getx(p)))
function compute_triangulation_area(tri)
    A = 0.0
    for T in each_solid_triangle(tri)
        i, j, k = triangle_vertices(T)
        p, q, r = get_point(tri, i, j, k)
        A += area(p, q, r)
    end
    return A
end
A = compute_triangulation_area(tri)

get_triangles(tri)

get_neighbours(tri, 5)

get_adjacent2vertex(tri, 5)

filter(T -> 5 ∈ triangle_vertices(T), get_triangles(tri))

get_neighbours(tri, -1)

get_adjacent2vertex(tri, -1)

get_adjacent(tri, 163, 365)

DelaunayTriangulation.triangle_orientation(tri, 163, 365, 295)

get_adjacent(tri, 365, 163)

get_adjacent(tri, 398, 258)

DelaunayTriangulation.is_boundary_edge(tri, 258, 398)

This page was generated using Literate.jl.