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
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.