Delaunay triangulation: Bowyer-Watson algorithm

In my last two posts, I did a bit of prep work leading to an implementation of the Bowyer-Watson algorithm. Now that we have the geometry building blocks we need, we can attack the core of the algorithm, and perform a Delaunay triangulation on a list of points.

Rather than attempt to explain what a Delaunay triangulation is, I will leave that out, and simply illustrate on an example. Starting from a collection of 20 random points, like so:

let points =
    let rng = Random 0
    List.init 20 (fun _ ->
        {
            X = rng.NextDouble() * 100.
            Y = rng.NextDouble() * 100.
        }
        )

… their Delaunay triangulation produces something like this: a mesh of triangles connecting all the points, without any edges crossing each other, and where the triangles have “reasonably even” angles:

The Bowyer-Watson algorithm is one way to generate such a triangulation. In this post, I’ll go over implementing it in F#.

Implementing the algorithm

In the same spirit as the previous posts, I will simply follow the Wikipedia pseudo-code, and convert it to F#. Let’s dive into the middle of the algorithm:

// add all the points one at a time to the triangulation
for each point in pointList do
    badTriangles := empty set
    // 1) first find all the triangles that are no longer valid due to the insertion
    for each triangle in triangulation do 
        if point is inside circumcircle of triangle
            add triangle to badTriangles
    // 2) find the boundary of the polygonal hole
    polygon := empty set
    for each triangle in badTriangles do
        for each edge in triangle do
            if edge is not shared by any other triangles in badTriangles
                add edge to polygon
    // 3) remove them from the data structure
    for each triangle in badTriangles do
        remove triangle from triangulation
    // 4) re-triangulate the polygonal hole
    for each edge in polygon do
        newTri := form a triangle from edge to point
        add newTri to triangulation

The algorithm starts with a valid triangulation (the initial “super triangle”, see [previous post][]), and adds points one-by-one, checking (1) which of the existing triangles are no longer valid, and updating them (2, 3, 4). This suggests a function which takes in a collection of triangles and a point, and returns an updated collection of triangles, something like this:

addPoint: Point -> Triangle array -> Triangle array

Step 1 is pretty direct. We re-use some of the tools we wrote in the previous post, to separate the current triangles in 2 groups, good and bad triangles:

// 1) first find all the triangles that are no longer valid due to the insertion
for each triangle in triangulation do 
    if point is inside circumcircle of triangle
        add triangle to badTriangles

let addPoint (point: Point) (triangles: Triangle []) =
    let badTriangles, goodTriangles =
        triangles
        |> Array.partition (fun triangle ->
            triangle
            |> circumCircle
            |> Circle.contains point
            )

Note: I changed Circle.isInside to Circle.contains, which reads much better I think.

Done. Now to Step 2:

// 2) find the boundary of the polygonal hole
for each triangle in badTriangles do
    for each edge in triangle do
        if edge is not shared by any other triangles in badTriangles
            add edge to polygon

This is a bit trickier. For each bad triangle, we want to extract the sides that are not shared, that is, unique edges. The tricky bit is that 2 edges (2 points) are equal if the 2 points are equal, regardless of their order. That is, edge A,B is equal to edge B,A.

Rather than implement equality, I will be lazy here, and use a trick:

let edge (pt1: Point, pt2: Point) =
    if pt1.X < pt2.X
    then pt1, pt2
    elif pt1.X > pt2.X
    then pt2, pt1
    else
        if pt1.Y < pt2.Y
        then pt1, pt2
        else pt2, pt1

The edge function takes a pair of points, and returns a tuple where the points are sorted by X coordinates and Y coordinates. Assuming I did not make a logical error here, this should result in edge(A,B)=edge(B,A).

Armed with that, we can then go over the bad triangles, extract all the edges, count them, and keep only the edges that are unique:

let uniqueEdges =
    badTriangles
    |> Array.collect (fun triangle ->
        [|
            triangle.A, triangle.B
            triangle.B, triangle.C
            triangle.C, triangle.A
        |]
        )
    |> Array.map edge
    |> Array.countBy id
    |> Array.filter (fun (_, count) -> count = 1)
    |> Array.map fst

Step 3 is already done for us. We don’t need to remove anything, because we already extracted the good triangles in Step 1. All that is left to do is Step 4: create a triangle for each unique edge, connecting it to the point we are adding to the triangulation:

for each edge in polygon do
    newTri := form a triangle from edge to point
    add newTri to triangulation

Pretty direct again:

let newTriangles =
    uniqueEdges
    |> Array.map (fun (a, b) ->
        { A = a; B = b; C = point }
        )

Array.append goodTriangles newTriangles

… and we are done. The complete addPoint function looks like this:

let addPoint (point: Point) (triangles: Triangle []) =
    let badTriangles, goodTriangles =
        triangles
        |> Array.partition (fun triangle ->
            triangle
            |> circumCircle
            |> Circle.contains point
            )
    let edge (pt1: Point, pt2: Point) =
        if pt1.X < pt2.X
        then pt1, pt2
        elif pt1.X > pt2.X
        then pt2, pt1
        else
            if pt1.Y < pt2.Y
            then pt1, pt2
            else pt2, pt1

    let uniqueEdges =
        badTriangles
        |> Array.collect (fun triangle ->
            [|
                triangle.A, triangle.B
                triangle.B, triangle.C
                triangle.C, triangle.A
            |]
            )
        |> Array.map edge
        |> Array.countBy id
        |> Array.filter (fun (_, count) -> count = 1)
        |> Array.map fst

    let newTriangles =
        uniqueEdges
        |> Array.map (fun (a, b) ->
            { A = a; B = b; C = point }
            )

    Array.append goodTriangles newTriangles

Adding all the points

Now that we have a function to add a point to an existing triangulation, all we need is to chain that together, adding the points one by one.

Let’s start manually, for illustration. First, we need an initial super triangle that will contain all the points (see the first post in this series):

let points =
    let rng = Random 0
    List.init 20 (fun _ ->
        {
            X = rng.NextDouble() * 100.
            Y = rng.NextDouble() * 100.
        }
        )

let triangle = BowyerWatson.superTriangle points

This gives us the points we want to triangulate, and an initial triangle. We can now start adding the first point, which produces an array of 3 triangles:

let update1 =
    BowyerWatson.addPoint points.[0] [| triangle |]

And we can keep going:

let update2 =
    BowyerWatson.addPoint points.[1] update1

let update3 =
    BowyerWatson.addPoint points.[2] update2

After update 3, this produces the following triangulation:

Note that the triangulation still includes the initial, outer super-triangle. We will remove that in a minute. First, let’s see how we can write that loop and add all the points in one go.

Updates 1, 2 and 3 follow a clear pattern:

• Take the current collection of triangles, the State of the triangulation,
• Take the next Point in the collection of Points,
• Add the Point, which gives us a new collection of triangles,
• Repeat, using the new collection of triangles as a new State.

This is exactly what a fold does: starting from an initial state and a collection, we iterate over the collection, updating the state each time until there isn’t anything left to iterate over.

val fold:
   folder: ('State -> 'T -> 'State) ->
   state : 'State ->
   source: seq<'T>
        -> 'State

So we can iterate over all the points like so:

([| triangle |], points)
||> Seq.fold (fun triangulation point ->
    BowyerWatson.addPoint point triangulation
    )

Wrapping it all up

We are more or less done at that point. The only thing we need to do is remove any triangle that has an edge connected to one of the initial 3 points. Let’s do that, and wrap the whole thing into a function:

let delaunay (points: seq<Point>) =
    let corners (triangle: Triangle) =
        set [ triangle.A; triangle.B; triangle.C ]
    let initial = superTriangle points
    let initialCorners = corners initial
    ([| initial |], points)
    ||> Seq.fold (fun triangulation point ->
        addPoint point triangulation
        )
    |> Array.filter (fun triangle ->
        Set.intersect (corners triangle) initialCorners
        |> Set.isEmpty
        )

We write a small utility function, corners, which takes a triangle and extracts its 3 corners into a Set. This is useful, because we can use set functions, like Set.intersect, to verify if 2 triangles have a common corner, and, if so, remove them from the list. And… we are done, and can run a full Delaunay triangulation, and visualize it:

points
|> BowyerWatson.delaunay
|> List.ofArray
|> List.map (fun triangle ->
    polygon
        [
            Style.Fill (Fill.Color "lightyellow")
            Style.Stroke (Stroke.Color "black")
            Style.Stroke (Stroke.Width 0.2)
        ]
        [ triangle.A; triangle.B; triangle.C ]
    )
|> Plot.plot (400, 400)

Note: I am still working on that SVG generation DSL, once it’s in a presentable state I plan to blog about it.

Parting thoughts

I was expecting the core algorithm to give me more trouble than it did. As it turns out, the hardest part was getting the geometry functions done - wiring up the main loop was more or less a direct transcription of the pseudo code to F#.

That being said, I know I am not fully done, because there is a bug, which I spotted by accident. If I run the same code, using only the first 9 points instead of the full list, I get this:

This looks mostly correct, which is also commonly described as “wrong”. The problem is visible on the bottom edge, connecting the west and south-east corners. The triangulation looks pretty, but, if it is done right, the result should always be a convex hull. Stated differently: the outer polygon should never bend inwards.

There is a clear inwards bend on that edge, which tells me that something is wrong in my implementation somewhere. If you run the algorithm for more points, the problem fixes itself, but reappears a few iterations later. The algorithm as implemented seems to occasionally produce obviously incorrect triangulations. This might be a more general problem, but it seems related to adding points that are added close to the outside boundary.

The problem here is that I am not entirely sure yet how to approach fixing that issue. Fortunately, I have an example that reproduces the bug, but it involves adding 1 point to an 8 points mesh. Isolating what step goes wrong or even debugging will be a nightmare.

Anyways, that’s probably what I will be tackling next! In the meantime, if you can spot the bug… let me know :)