# 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_triangles`

— Function`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)
```

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)`

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_list`

— Function`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_vertex`

— Function`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_neighbours`

— Function```
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`

.

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_cdt`

— Function```
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`

Here we used `add_edge!(tri, 2, 7)`

, which does all this pre-processing for us. Similarly, for adding many edges, the method

`DelaunayTriangulation.triangulate_constrained!`

— Function`triangulate_constrained(tri::Triangulation; rng=Random.default_rng())`

Inserts all constrained edges and boundary edges into `tri`

.

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.delete_holes!`

— Function`delete_holes!(tri::Triangulation)`

Deletes all the exterior faces to the boundary nodes specified in the triangulation `tri`

.

`DelaunayTriangulation.has_interiors_within_interiors`

— Function`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.