Constrained Triangulations

Here we describe the algorithm used for computing a constrained triangulation. We assume that we have some point set $\mathcal P$, a set of edges $\mathcal E$ to be inserted, and some boundary edges $\mathcal B$. The algorithm we implement is given here and is built upon the basic idea of Chew's algorithm for triangulating convex polygons, namely solving the point location problem by deleting an associated cavity in a clever way.

First, suppose we have computed the unconstrained Delaunay triangulation $\mathcal D\mathcal T(\mathcal P)$ of our point set. The algorithm for then computing the constrained Delaunay triangulation $\mathcal D\mathcal T(\mathcal P, \mathcal E, \mathcal B)$ then works incrementally, inserting edges one at a time. Let us, then, describe the procedure for inserting some edge $e \in \mathcal C$, where $\mathcal C = \mathcal E \cup \mathcal B$. We break the discussion into sections.

Segment Location

Similar to how in adding points into a triangulation the first step is point location, here the first step is segment location. Here, our aim is to find all triangles that intersect the edge $e$. This is easy to do if we simply remember that the jump-and-march algorithm walks through all triangles that intersect an initial scan line until stopping: this is exactly what we want. So, what we have done is modify our jump-and-march code such that the history of triangles walked through is recorded. This information is recorded into a PointLocationHistory struct. With this, we also note that since the jump-and-march algorithm will require rotating around an initial point (one of the indices of $e$), we should want to minimise the number of triangles we may need to rate around. Thus, we rotate $e$ such that its initial vertex is the one with the least degree, hence we have less triangles to rotate around initially, thus reducing the time spent searching. From here on, we let $e$ be this rotated form. This segment location is handled via the following function:

DelaunayTriangulation.locate_intersecting_trianglesFunction
locate_intersecting_triangles(tri::Triangulation, e, rotate=Val(true); rng::AbstractRNG=Random.default_rng()) where {C}

Returns a list of triangles intersecting the segment e in tri. If is_true(rotate), then e will be sorted such that initial(e) has smaller degree in tri than terminal(e).

More precisely, the returned values are:

  • intersecting_triangles: The triangles intersecting e.
  • collinear_segments: Any segments collinear with e, giving in order of appearance.
  • left_vertices: All vertices of intersecting_triangles appearing to the left of e.
  • right_vertices: All vertices of intersecting_triangles appearing to the right of e.

This function also returns information about any segments that are collinear with e and all vertices of the intersecting triangles to the left and to the right of e. For the collinear segments, these are processed via the functions fix_segments!, connect_segments!, and extend_segments!, breaking e into a smaller set of segments so that no segments are collinear anymore. We will assume that there are no collinear segments for simplicity. For the vertices to the left and to the right of e, these are needed as they define the outline of points to be deleted on each side of e, thus giving a polygonal cavity (possibly self-intersecting, but this detail doesn't actually matter) that we can triangulate individually. For example, consider the triangulation:

a = (0.0, 0.0)
b = (0.0, 1.0)
c = (0.0, 2.5)
d = (2.0, 0.0)
e = (6.0, 0.0)
f = (8.0, 0.0)
g = (8.0, 0.5)
h = (7.5, 1.0)
i = (4.0, 1.0)
j = (4.0, 2.5)
k = (8.0, 2.5)
pts = [a, b, c, d, e, f, g, h, i, j, k]
tri = triangulate(pts; delete_ghosts=false, randomise=false)
fig, ax, sc = triplot(tri)
lines!(ax, [get_point(tri, 2, 7)...], color=:blue, linewidth=2)
An edge through a triangulation

When we perform segment location on this example on the highlighted segment (2, 7), we find:

julia> e = (2, 7)
(2, 7)

julia> intersecting_triangles, collinear_segments, left_vertices, right_vertices = DelaunayTriangulation.locate_intersecting_triangles(tri, e);

julia> intersecting_triangles
8-element Vector{Tuple{Int, Int, Int}}:
 (4, 3, 2)
 (3, 4, 10)
 (10, 4, 9)
 (9, 4, 5)
 (9, 5, 10)
 (10, 5, 8)
 (8, 5, 6)
 (8, 6, 7)

julia> collinear_segments
Tuple{Int, Int}[]

julia> left_vertices
7-element Vector{Int}:
  7
  8
 10
  9
 10
  3
  2

julia> right_vertices
5-element Vector{Int}:
 2
 4
 5
 6
 7

The intersecting triangles gives the triangles intersected by e in order of occurrence. For left_vertices, these are given in counter-clockwise order, with the left- and right-most elements being the indices of e. See that the vertex 10 is repeated in left_vertices. This is because the point 9 creates a sort of dangling edge in the polygonal cavity that we have to delete, and we need to somehow know to insert this point back into the triangulation. What we do, then, is to imagine an ant walking around the polygonal cavity. The ant will walk from 10 to 9, but then it has to come back down to 10, so we include it twice to represent these two visits. To see the cavity, let us delete these triangles:

DelaunayTriangulation.delete_intersected_triangles!(tri, intersecting_triangles)
An edge through a triangulation with excavated cavities

The polygonal cavities on each side of the blue segment are what we need to re-triangulate separately. Let us now describe this procedure.

Triangulating the Polygonal Cavities

Now we need to triangulate each cavity. We will describe this only for the cavity above (2, 7) in the figure below, but the procedure is exactly the same for the other cavity. Let $\mathcal V = (v_1, \ldots, v_m)$ be the sequence of vertices in counter-clockwise order around the cavity when we insert the segment $e=(v_1,v_m)$. In this case, $\mathcal V = (7, 8, 10, 9, 10, 3, 2)$. First, just like in Chew's algorithm, we need to build up a linked-list representing the cavity. This is done via the function

DelaunayTriangulation.prepare_vertex_linked_listFunction
prepare_vertex_linked_list(V::AbstractArray{I}) where {I}

Given a list of polygon vertices V, defines a linked list (prev, next) of polygon vertices so that (prev[i], i, next[i]) define a trio of polygon vertices in counter-clockwise order, and defines, and returns shuffled_indices which is currently unshuffled.

The first and last entries of the returned values (prev, next, shuffled_indices) will not be populated, instead being 0.

(Note: The algorithm given by Shewchuk and Brown linked above also allocates a distance array storing orient determinants for this preparation of the linked list, giving values proportional to the distance from $e$. This is not exactly robust for us, since we need to compute the sign of a difference of two determinants. In particular, let $o_1$ and $o_2$ be two robust estimates for an orient determinant. To determine if o_1 < o_2 is the same as defining a predicate for the sign of o_1 - o_2, but this is problematic as, while the computation of o_1 and o_2 may be reliable, their difference is not. Instead, we recompute this predicate in a robust manner each time, trading performance for robustness. This predicate is defined by point_closest_to_line.)

Once we have prepared the linked list, we need to delete vertices from it in a random order, corresponding to deleting vertices from the polygon in a random order. Like in Chew's algorithm, this is done in such a way that we can reverse the process and automatically get point location without ever needing the jump-and-march algorithm. This deletion is handled via delete_polygon_vertices_in_random_order! which simply loops over each vertex, calling select_random_vertex and update_vertex_linked_list! at each iteration:

DelaunayTriangulation.delete_polygon_vertices_in_random_order!Function
delete_polygon_vertices_in_random_order!(tri::Triangulation, V, shuffled_indices, prev, next, u, v, rng::AbstractRNG=Random.default_rng())

Given a triangulation tri, a vertex list V, a set of shuffled_indices, a linked list (prev, next) for the poylgon vertices, and a segment (u, v) that was inserted in order to define the polygon V, deletes vertices of V, via their representation in (prev, next, shuffled_indices), in a random order.

DelaunayTriangulation.select_random_vertexFunction
select_random_vertex(tri::Triangulation, V, shuffled_indices, prev, next, range, u, v, rng::AbstractRNG=Random.default_rng())

Given a triangulation tri, a line through points with indices u and v, a shuffled set of indices shuffled_indices, and a doubly-linked list (prev, next) of vertex indices, selects a random vertex j ∈ range that is not closer to the line than both of its neighbours. V is the original vertex list.

DelaunayTriangulation.update_vertex_linked_list!Function
update_vertex_linked_list!(shuffled_indices, prev, next, i, j)

Let π = shuffled_indices. This function replaces next[prev[π[j]]] with next[π[j]], prev[next[π[j]]] with prev[π[j]], and interchanges π[i] and π[j]. This has the act of deleting V[π[j]] from the polygon, where V is the list of polygon vertices of the polygon being evacuated during segment insertion for a constrained triangulation.

An important note is that, while this insertion algorithm works even with dangling edges, self-intersections, etc., it does not work when a point has an interior angle exceeding 360 degrees. We can detect this case by finding a point in the polygon that is closer to $e$ than its two neighbours, which is the only time such an interior angle is possible (see the paper for a proof). Thus, select_random_vertex actually keeps sampling vertices to delete in a random order until this is not the case, making use of vertex_is_closer_than_neighbours:

DelaunayTriangulation.vertex_is_closer_than_neighboursFunction
vertex_is_closer_than_neighbours(tri::Triangulation, u, v, jᵢ, jᵢ₋₁, jᵢ₊₁)
vertex_is_closer_than_neighbours(tri::Triangulation, V, u, v, j, shuffled_indices, prev, next)

Given a triangulation tri and a line through points with indices u and v, tests if the point with index jᵢ is closer to the line than those with index jᵢ₋₁ and jᵢ₊₁, assuming all these points are to the left of the line. The second method extracts these latter two indices using the linked list (prev, next) of vertices and a shuffled set of indices shuffled_indices together with the original vertex list V.

Note

This function is useful for constrained triangulations since the algorithm used will not work if a point being inserted on the cavity has interior angle of 360° or greater. This is possible only if a vertex is closer to the line than its neighbours on the polygon.

Once this is all done, we ready to start triangulating the cavity. After adding an initial triangle from the three remaining vertices, we add points in one at a time, making use of a function add_point_cavity_cdt!. This function is defined by:

function add_point_cavity_cdt!(tri::Triangulation, u, v, w)
    x = get_adjacent(tri, w, v)
    if !edge_exists(x)
        insert_flag = true
    else
        p, q, r, s = get_point(tri, w, v, x, u) 
        incircle_test = point_position_relative_to_circle(p, q, r, s)
        orient_test = triangle_orientation(tri, u, v, w)
        insert_flag = !is_inside(incircle_test) && is_positively_oriented(orient_test)
    end
    if insert_flag
        add_triangle!(tri, u, v, w; protect_boundary=true, update_ghost_edges=false)
    else
        delete_triangle!(tri, w, v, x; protect_boundary=true, update_ghost_edges=false)
        add_point_cavity_cdt!(tri, u, v, x)
        add_point_cavity_cdt!(tri, u, x, w)
    end
    return nothing
end

In particular, x = get_adjacent(tri, w, v) is used to find the triangle on the other side of the edge $vw$ from $u$, the point being inserted. The only way that $wv$ should not be deleted is if the triangle $wvx$ does not exist, as detected via !edge_exists(x), or if u is not inside the circumcircle of $wvx$ and $u$ is on the correct side of the edge $vw$. This is detected by the computation of insert_flag. If insert_flag, just add the triangle. Otherwise, we need to delete the triangle $wvx$ and $uvw$ as they are no longer constrained Delaunay. This is done by flipping $vw$ onto $ux$. Once we have done this for each point, we have successfuly triangulated the cavity.

This is all handled via the triangulate_cavity_cdt function:

DelaunayTriangulation.triangulate_cavity_cdtFunction
triangulate_cavity_cdt(points, V;
    IntegerType::Type{I}=Int,
    EdgeType::Type{E}=NTuple{2,IntegerType},
    TriangleType::Type{Vs}=NTuple{3,IntegerType},
    EdgesType::Type{Es}=Set{EdgeType},
    TrianglesType::Type{Ts}=Set{TriangleType},
    rng::AbstractRNG=Random.default_rng()) where {I,E,Vs,Es,Ts}
triangulate_cavity_cdt(tri, V; rng::AbstractRNG=Random.default_rng())

Triangulates the cavity, represented as a counter-clockwise list of vertices V with indices corresponding to those in points, left behind when deleting triangles intersected in a triangulation by an edge. If a triangulation is provided, the points are used from that.

For our example, what we find is (the triangulation tri was updated to put the missing triangles back in from the last piece of code, note):

julia> add_edge!(tri, 2, 7) # calls triangulate_cavity_cdt on each cavity
A constrained edge through a triangulation

Here we used add_edge!(tri, 2, 7), which does all this pre-processing for us. Similarly, for adding many edges, the method

is useful (triangulate calls this internally).

Excavating Exterior Faces

When we define boundary curves, we typically want to delete any points and triangles exterior to them. The logic of the method we use for this is simple. Basically, we "plant" a seed in an exterior face, and use it to infect other points in this exterior face, continuing this spread until all exterior faces are found. The function that performs this is delete_holes!, with relevant docstrings below.

DelaunayTriangulation.has_interiors_within_interiorsFunction
has_interiors_within_interiors(tri::Triangulation)

Returns true if the triangulation has multiple curves and the first curve has a positive area and all other curves have negative areas, meaning there are some interior curves that are inside other interior curves. Returns false otherwise.

The points that we find are then processed one at a time, checking all adjoining triangles to see if their centroid is in the interior or exterior. Special case is taken at the boundary nodes.