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 laterWe 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 RNG2×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.606559We 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: falsefig, ax, sc = triplot(tri)
fig
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.477077379328767We 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
140One 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.05645693808033534By 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
214The 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
140get_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)295means 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)232If 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)-1means 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)trueJust 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.