Delaunay triangulation: hitting an impasse

During the recent weeks, I have been making slow but steady progress implementing Delaunay triangulation with the Bowyer-Watson algorithm. However, as I mentioned in the conclusion of my previous post, I spotted a bug, which I hoped would be an easy fix, but so far no such luck: it has me stumped. In this post, I will go over how I approached figuring out the problem, which is interesting in its own right. My hope is that by talking it through I might either get an idea, or perhaps someone else will spot what I am missing!

Anyways, let’s get into it. Per the Wikipedia entry:

a Delaunay triangulation […] of a set of points in the plane subdivides their convex hull into triangles whose circumcircles do not contain any of the points; that is, each circumcircle has its generating points on its circumference, but all other points in the set are outside of it. This maximizes the size of the smallest angle in any of the triangles, and tends to avoid sliver triangles.

The Bowyer-Watson algorithm seemed straightforward to implement, so that’s what I did. For those interested, my current code is here, and it mostly works. Starting from a list of 20 random points, I can generate a triangulation like so:

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

With a bit of tinkering, I can render the result using SVG:

This looks like what I would expect: all the points are connected in a mesh of triangles.

However, if I take a subset of the points, like so:

points
|> List.truncate 9
|> BowyerWatson.delaunay

… I then get the following result:

This is still pretty, but it can’t be right. From the Wikipedia entry:

a Delaunay triangulation […] of a set of points in the plane subdivides their convex hull into triangles

The result should be a convex hull, that is, the outer boundary should be convex. Or, stated differently, the outer boundary should never “bend inwards”, and yet, we can see such a bend on the figure above, in the West to South-East boundary. What is going on here?

My first assumption was of course Denial. Certainly, my code is correct, I must be hitting a weird edge case with that particular set of points. So I generated a list of Delaunay diagrams, taking more and more points:

This does unfortunately not look like a fluke. I can observe bends inwards occurring multiple times, for example in images 7, 8, 13.

So I switched from Denial to Acceptance: I must have made a bone-headed coding error somewhere. Fine, but… how do I even approach understanding what might be going wrong? The first problem occurs on image 7, and there are already 10 points involved. This is going to be a nightmare to work with.

My first step was to try and shrink the problem to the smallest failing case I could. With a bit of tinkering and intuition, I ended up reproducing the issue with a much smaller setup:

let points =
    [ 0; 2; 4; 8 ]
    |> List.map (fun i -> points.[i])
points
|> BowyerWatson.delaunay

This involves only 4 points now, and the result is most definitely not a convex polygon:

Great. Now that we have a small test case, let’s dig in. We’ll run the triangulation before and after the step where things go off the rails, including the outer super-triangle, which is relevant to how the triangulation is updated:

Before

After

Here we can see how, as we add a point in the South-East area, while the triangulation gets properly adjusted in the North-East area, the bottom section shows a bend inwards appearing. Some of the triangles that should be recomputed are not.

How does the algorithm decide if a triangle needs to be recomputed? Per the Bowyer-Watson pseudo code,

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

When a point is added, we take the existing triangles, and, if the point is inside the circumcircle of a triangle, we add it to the list of “bad triangles” that need to be recomputed. Let’s visualize these circles then:

And this is where I hit an impasse. On this diagram, we can see the point we are adding to the triangulation, in the South-East area. The triangle West of it must be recomputed to avoid creating the issue we observe, but its circumcircle, while coming close to our point, does not contain it. In other words, the algorithm concludes that this section of the triangulation is fine - when it is not.

Parting thoughts

So where does this leave me? I am genuinely stumped here. From what I can tell, my implementation is faithful to the pseudo-code. I can’t see an obvious error, and yet, the result is clearly incorrect.

My best guess at the moment is that either I missed a subtlety in the algorithm, or that perhaps there is an unspoken assumption / missing detail, perhaps around triangles that are connected to the corners of the outer super-triangle.

I suspect the issue has to do with the initial super-triangle. The algorithm begins by creating the so-called “super triangle”, which, per the pseudo-code outline, should just be

large enough to completely contain all the points

add super-triangle to triangulation
// must be large enough to completely contain all the points in pointList

However, that triangle is not unique, and, if I were to move a bit the position of the South corner of the triangle, I could move it so that it still contains all the points, and causes a re-computation of the problem area. Long story short, I suspect that the source of the issue might be around what conditions that initial triangle must satisfy, or perhaps around recomputations involving its corners.

At any rate, I am stuck. My next steps are first to sleep on it, and then possibly to see if the original source articles might help. In the meantime, hopefully you will have found something of interest in this bug hunt, with its unusual graphical debugging. And if you happen to know what I am doing wrong, I would love to hear it :)