Constrained Triangulations

Outer Boundary

This tutorial now considers the case where, rather than only having constrained segments, we have a constrained outer boundary. This is especially useful as it allows us to, for example, have a non-convex boundary. To start, let us load in the packages we will need.

using DelaunayTriangulation
using CairoMakie

Now, we define some of the points we will be triangulating.

pts = [
    (-7.36, 12.55), (-9.32, 8.59), (-9.0, 3.0), (-6.32, -0.27),
    (-4.78, -1.53), (2.78, -1.41), (-5.42, 1.45), (7.86, 0.67),
    (10.92, 0.23), (9.9, 7.39), (8.14, 4.77), (13.4, 8.61),
    (7.4, 12.27), (2.2, 13.85), (-3.48, 10.21), (-4.56, 7.35),
    (3.44, 8.99), (3.74, 5.87), (-2.0, 8.0), (-2.52, 4.81),
    (1.34, 6.77), (1.24, 4.15),
]

22-element Vector{Tuple{Float64, Float64}}:
 (-7.36, 12.55)
 (-9.32, 8.59)
 (-9.0, 3.0)
 (-6.32, -0.27)
 (-4.78, -1.53)
 (2.78, -1.41)
 (-5.42, 1.45)
 (7.86, 0.67)
 (10.92, 0.23)
 (9.9, 7.39)
 ⋮
 (2.2, 13.85)
 (-3.48, 10.21)
 (-4.56, 7.35)
 (3.44, 8.99)
 (3.74, 5.87)
 (-2.0, 8.0)
 (-2.52, 4.81)
 (1.34, 6.77)
 (1.24, 4.15)

To define a boundary, we need to provide a counter-clockwise sequence of indices corresponding to the boundary points, and the first index must match the last index so the boundary is closed. While we could include in pts the boundary points that we want to include, and then write down the indices of the points within pts, this is cumbersome and often tedious to get correct. So, we instead provide the function convert_boundary_points_to_indices which takes in a vector of coordinates, and then returns the correct set of indices. Here is how we use it:

boundary_points = [
    (0.0, 0.0), (2.0, 1.0), (3.98, 2.85), (6.0, 5.0),
    (7.0, 7.0), (7.0, 9.0), (6.0, 11.0), (4.0, 12.0),
    (2.0, 12.0), (1.0, 11.0), (0.0, 9.13), (-1.0, 11.0),
    (-2.0, 12.0), (-4.0, 12.0), (-6.0, 11.0), (-7.0, 9.0),
    (-6.94, 7.13), (-6.0, 5.0), (-4.0, 3.0), (-2.0, 1.0), (0.0, 0.0),
]
boundary_nodes, pts = convert_boundary_points_to_indices(boundary_points; existing_points = pts);

The keyword argument existing_points is so that the points in boundary_points get appended (in-place) to pts, as we see:

pts

42-element Vector{Tuple{Float64, Float64}}:
 (-7.36, 12.55)
 (-9.32, 8.59)
 (-9.0, 3.0)
 (-6.32, -0.27)
 (-4.78, -1.53)
 (2.78, -1.41)
 (-5.42, 1.45)
 (7.86, 0.67)
 (10.92, 0.23)
 (9.9, 7.39)
 ⋮
 (-1.0, 11.0)
 (-2.0, 12.0)
 (-4.0, 12.0)
 (-6.0, 11.0)
 (-7.0, 9.0)
 (-6.94, 7.13)
 (-6.0, 5.0)
 (-4.0, 3.0)
 (-2.0, 1.0)

The boundary_nodes is then these indices:

boundary_nodes

21-element Vector{Int64}:
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
  ⋮
 35
 36
 37
 38
 39
 40
 41
 42
 23

To now triangulate, we use the boundary_nodes keyword argument. Like in the last tutorial, we also give a comparison to the unconstrained version.

tri = triangulate(pts)
cons_tri = triangulate(pts; boundary_nodes)

Delaunay Triangulation.
   Number of vertices: 28
   Number of triangles: 34
   Number of edges: 61
   Has boundary nodes: true
   Has ghost triangles: true
   Curve-bounded: false
   Weighted: false
   Constrained: true

fig = Figure()
ax1 = Axis(
    fig[1, 1], xlabel = "x", ylabel = L"y",
    title = "(a): Unconstrained", titlealign = :left,
    width = 300, height = 300,
)
ax2 = Axis(
    fig[1, 2], xlabel = "x", ylabel = L"y",
    title = "(b): Constrained", titlealign = :left,
    width = 300, height = 300,
)
triplot!(ax1, tri)
triplot!(ax2, cons_tri, show_constrained_edges = true, show_convex_hull = true)
resize_to_layout!(fig)
fig

Notice now that the boundary in (b) is not convex, as is clear from the convex hull shown in red. You can access the convex hull using get_convex_hull(cons_tri). We also note that the triangulation no longer contains every point in pts, as by default all triangles away from the boundary are deleted, so that we do actually have a boundary. If for some reason you do not want this behaviour, use delete_holes = false:

full_tri = triangulate(pts; boundary_nodes, delete_holes = false)
fig, ax, sc = triplot(full_tri, show_constrained_edges = true, show_convex_hull = true)

This default behaviour does mean you need to be careful if you use DelaunayTriangulation.each_point or DelaunayTriangulation.each_point_index, as these iterators will contain all points, possibly iterating over points that aren't in the triangulation. For this reason, it is recommended that you use each_solid_vertex as a default.

There are multiple methods available for working directly with the boundary nodes. You can get the boundary nodes using get_boundary_nodes(tri):

get_boundary_nodes(cons_tri)

21-element Vector{Int64}:
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
  ⋮
 35
 36
 37
 38
 39
 40
 41
 42
 23

Later tutorials also consider other methods for working with the boundary where care needs to be taken with the boundary, or part of the boundary, being considered. For now, here is an example where we use get_right_boundary_node to iterate over the boundary in a counter-clockwise order, getting the area of the triangulation using

\[A = \dfrac{1}{2}\sum_{i=1}^n \left(y_i + y_{i+1}\right)\left(x_i - x_{i+1}\right).\]

Here is one implementation.

function shoelace_area(tri)
    bn = get_boundary_nodes(tri)
    n = num_boundary_edges(bn) # length(bn) - 1 in this case since bn[1] = bn[end]
    A = 0.0
    for i in 1:n
        vᵢ = get_boundary_nodes(bn, i)
        vᵢ₊₁ = get_boundary_nodes(bn, i + 1)
        pᵢ, pᵢ₊₁ = get_point(tri, vᵢ, vᵢ₊₁)
        xᵢ, yᵢ = getxy(pᵢ)
        xᵢ₊₁, yᵢ₊₁ = getxy(pᵢ₊₁)
        A += (yᵢ + yᵢ₊₁) * (xᵢ - xᵢ₊₁)
    end
    return A / 2
end
shoelace_area(cons_tri)

119.20499999999998

We also provide a map that contains the edges as the keys (not in order), and the values are Tuples (I, J) such that get_boundary_nodes(get_boundary_nodes(cons_tri, I), J) gives the corresponding edge. The first call, bn = get_boundary_nodes(cons_tri, I) is for obtaining the chain of boundary edges containing the boundary edge, and then get_boundary_nodes(bn, j) gets the actual edge.

get_boundary_edge_map(cons_tri)

Dict{Tuple{Int64, Int64}, Tuple{Vector{Int64}, Int64}} with 20 entries:
  (23, 24) => ([23, 24, 25, 26, 27, 28, 29, 30, 31, 32  …  34, 35, 36, 37, 38, …
  (40, 41) => ([23, 24, 25, 26, 27, 28, 29, 30, 31, 32  …  34, 35, 36, 37, 38, …
  (25, 26) => ([23, 24, 25, 26, 27, 28, 29, 30, 31, 32  …  34, 35, 36, 37, 38, …
  (35, 36) => ([23, 24, 25, 26, 27, 28, 29, 30, 31, 32  …  34, 35, 36, 37, 38, …
  (30, 31) => ([23, 24, 25, 26, 27, 28, 29, 30, 31, 32  …  34, 35, 36, 37, 38, …
  (29, 30) => ([23, 24, 25, 26, 27, 28, 29, 30, 31, 32  …  34, 35, 36, 37, 38, …
  (42, 23) => ([23, 24, 25, 26, 27, 28, 29, 30, 31, 32  …  34, 35, 36, 37, 38, …
  (33, 34) => ([23, 24, 25, 26, 27, 28, 29, 30, 31, 32  …  34, 35, 36, 37, 38, …
  (38, 39) => ([23, 24, 25, 26, 27, 28, 29, 30, 31, 32  …  34, 35, 36, 37, 38, …
  (27, 28) => ([23, 24, 25, 26, 27, 28, 29, 30, 31, 32  …  34, 35, 36, 37, 38, …
  (26, 27) => ([23, 24, 25, 26, 27, 28, 29, 30, 31, 32  …  34, 35, 36, 37, 38, …
  (39, 40) => ([23, 24, 25, 26, 27, 28, 29, 30, 31, 32  …  34, 35, 36, 37, 38, …
  (24, 25) => ([23, 24, 25, 26, 27, 28, 29, 30, 31, 32  …  34, 35, 36, 37, 38, …
  (28, 29) => ([23, 24, 25, 26, 27, 28, 29, 30, 31, 32  …  34, 35, 36, 37, 38, …
  (41, 42) => ([23, 24, 25, 26, 27, 28, 29, 30, 31, 32  …  34, 35, 36, 37, 38, …
  (31, 32) => ([23, 24, 25, 26, 27, 28, 29, 30, 31, 32  …  34, 35, 36, 37, 38, …
  (32, 33) => ([23, 24, 25, 26, 27, 28, 29, 30, 31, 32  …  34, 35, 36, 37, 38, …
  (34, 35) => ([23, 24, 25, 26, 27, 28, 29, 30, 31, 32  …  34, 35, 36, 37, 38, …
  (37, 38) => ([23, 24, 25, 26, 27, 28, 29, 30, 31, 32  …  34, 35, 36, 37, 38, …
  (36, 37) => ([23, 24, 25, 26, 27, 28, 29, 30, 31, 32  …  34, 35, 36, 37, 38, …

In our case, the I is just boundary_nodes since we only have one contiguous boundary. To give an example, take

bem = get_boundary_edge_map(cons_tri)
e, (I, J) = first(bem)

(23, 24) => ([23, 24, 25, 26, 27, 28, 29, 30, 31, 32  …  34, 35, 36, 37, 38, 39, 40, 41, 42, 23], 1)

bn = get_boundary_nodes(cons_tri, I) # same as boundary_nodes for this problem; see the later tutorials
bn_j = get_boundary_nodes(bn, J)

This returns 23, which is the start of the edge e. The full edge is given by

get_boundary_nodes.(Ref(bn), (J, J + 1)) # Ref to not broadcast over bn

(23, 24)

To give an example, here's how we compute the perimeter of the triangulation. This only needs the edges, so we only consider the keys of the map.

function get_perimeter(tri)
    bem = get_boundary_edge_map(tri)
    ℓ = 0.0
    for e in keys(bem)
        u, v = edge_vertices(e)
        p, q = get_point(tri, u, v)
        ℓ += sqrt((getx(p) - getx(q))^2 + (gety(p) - gety(q))^2)
    end
    return ℓ
end
get_perimeter(cons_tri)

44.23794172896859

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

pts = [
    (-7.36, 12.55), (-9.32, 8.59), (-9.0, 3.0), (-6.32, -0.27),
    (-4.78, -1.53), (2.78, -1.41), (-5.42, 1.45), (7.86, 0.67),
    (10.92, 0.23), (9.9, 7.39), (8.14, 4.77), (13.4, 8.61),
    (7.4, 12.27), (2.2, 13.85), (-3.48, 10.21), (-4.56, 7.35),
    (3.44, 8.99), (3.74, 5.87), (-2.0, 8.0), (-2.52, 4.81),
    (1.34, 6.77), (1.24, 4.15),
]

boundary_points = [
    (0.0, 0.0), (2.0, 1.0), (3.98, 2.85), (6.0, 5.0),
    (7.0, 7.0), (7.0, 9.0), (6.0, 11.0), (4.0, 12.0),
    (2.0, 12.0), (1.0, 11.0), (0.0, 9.13), (-1.0, 11.0),
    (-2.0, 12.0), (-4.0, 12.0), (-6.0, 11.0), (-7.0, 9.0),
    (-6.94, 7.13), (-6.0, 5.0), (-4.0, 3.0), (-2.0, 1.0), (0.0, 0.0),
]
boundary_nodes, pts = convert_boundary_points_to_indices(boundary_points; existing_points = pts);

pts

boundary_nodes

tri = triangulate(pts)
cons_tri = triangulate(pts; boundary_nodes)

fig = Figure()
ax1 = Axis(
    fig[1, 1], xlabel = "x", ylabel = L"y",
    title = "(a): Unconstrained", titlealign = :left,
    width = 300, height = 300,
)
ax2 = Axis(
    fig[1, 2], xlabel = "x", ylabel = L"y",
    title = "(b): Constrained", titlealign = :left,
    width = 300, height = 300,
)
triplot!(ax1, tri)
triplot!(ax2, cons_tri, show_constrained_edges = true, show_convex_hull = true)
resize_to_layout!(fig)
fig

full_tri = triangulate(pts; boundary_nodes, delete_holes = false)
fig, ax, sc = triplot(full_tri, show_constrained_edges = true, show_convex_hull = true)

get_boundary_nodes(cons_tri)

function shoelace_area(tri)
    bn = get_boundary_nodes(tri)
    n = num_boundary_edges(bn) # length(bn) - 1 in this case since bn[1] = bn[end]
    A = 0.0
    for i in 1:n
        vᵢ = get_boundary_nodes(bn, i)
        vᵢ₊₁ = get_boundary_nodes(bn, i + 1)
        pᵢ, pᵢ₊₁ = get_point(tri, vᵢ, vᵢ₊₁)
        xᵢ, yᵢ = getxy(pᵢ)
        xᵢ₊₁, yᵢ₊₁ = getxy(pᵢ₊₁)
        A += (yᵢ + yᵢ₊₁) * (xᵢ - xᵢ₊₁)
    end
    return A / 2
end
shoelace_area(cons_tri)

get_boundary_edge_map(cons_tri)

bem = get_boundary_edge_map(cons_tri)
e, (I, J) = first(bem)

bn = get_boundary_nodes(cons_tri, I) # same as boundary_nodes for this problem; see the later tutorials
bn_j = get_boundary_nodes(bn, J)

get_boundary_nodes.(Ref(bn), (J, J + 1)) # Ref to not broadcast over bn

function get_perimeter(tri)
    bem = get_boundary_edge_map(tri)
    ℓ = 0.0
    for e in keys(bem)
        u, v = edge_vertices(e)
        p, q = get_point(tri, u, v)
        ℓ += sqrt((getx(p) - getx(q))^2 + (gety(p) - gety(q))^2)
    end
    return ℓ
end
get_perimeter(cons_tri)

This page was generated using Literate.jl.