Delaunay triangulation with Bowyer-Watson: initial super triangle, revisited

In our last installment, I hit a roadblock. I attempted to implement Delaunay triangulations using the Bowyer-Watson algorithm, followed this pseudo-code from Wikipedia, and ended up with a mostly working F# implementation. Given a list of points, the code produces a triangulation, but occasionally the outer boundary of the triangulation is not convex, displaying bends towards the inside, something that should never supposed happen for a proper Delaunay triangulation.

While I could not figure out the exact issue, by elimination I narrowed it down a bit. My guess was that the issue was probably a missing unstated condition, probably related to the initial super-triangle. As it turns out, my guess was correct.

The reason I know is, a kind stranger on the internet reached out with a couple of helpful links (thank you!):

Bowyer-Watson algorithm: how to fill “holes” left by removing triangles with super triangle vertices
Bowyer-Watson algorithm for Delaunay triangulation fails, when three vertices approach a line

The second link in particular mentions that the Wikipedia page is indeed missing conditions, and suggests that the initial super triangle should verify the following property to be valid:

it seems that one should rather demand that the vertices of the super triangle have to be outside all circumcircles of any three given points to begin with (which is hard when any three points are almost collinear)

That doesn’t look overly complicated, let’s modify our code accordingly, and check if this fixes our problem!

First off, what does our original code do?

What we are after with the so-called super-triangle is a triangle that contains all the points we are triangulating. The strategy we followed is simple:

1) compute a rectangular box that contains all the points we want to triangulate,
2) compute a triangle that contains that rectangular box.

Below is how (1) looks in code; we omitted (2) because we are simply going to re-use it as-is, but you can look at this earlier post for details.

let superTriangle (points: seq<Point>): Triangle =

    // find the left and right edges 
    // of the containing box
    let xs =
        points
        |> Seq.map (fun pt -> pt.X)
        |> Array.ofSeq
    let xMin = xs |> Array.min
    let xMax = xs |> Array.max
    // find the top and bottom edges
    // of the containing box
    let ys =
        points
        |> Seq.map (fun pt -> pt.Y)
        |> Array.ofSeq
    let yMin = ys |> Array.min
    let yMax = ys |> Array.max

    // omitted: given a box, 
    // we can compute a triangle around it.

The missing condition we need to incorporate says that the corners of the super triangle

have to be outside all circumcircles of any three given points

Stated differently, the super triangle needs to be large enough to contain not just the points we want to triangulate, but also, for every single possible triangle they can form, the corresponding circumcircle.

As a reminder, the circumcircle of a triangle is the circle that passes through all the triangle corners:

Fortunately, we already have the code we need to compute that, because it is used in the algorithm itself. The details of the calculations are not particularly important, what matters here is that we have a function that, given a Triangle, will give us back a Circle like so:

let circumCircle (triangle: Triangle) =
    // Calculations details omitted, see link:
    // https://en.wikipedia.org/wiki/Circumcircle#Cartesian_coordinates_2
    { Center = center; Radius = radius }

With these tools in hand, we should be able to fix our code. Let’s try to reduce our problem to a (mostly) known problem, following this strategy:

1) enumerate every triangle we can form with the points we want to triangulate,
2) for each triangle, compute the circumcircle
3) compute a rectangular box that contains all the circles,
4) compute a triangle that contains that rectangular box.

(2) and (4) we know how to do already, all we need to focus on is (1) and (3).

So how can we enumerate the triangles we can form using a collection of points? Triangles that have the same corners are identical, regardless of the order (ABC is the same triangle as CBA or any other combination), so all we need is to enumerate all distinct combinations of 3 points we can form.

One way to go about it would be along these lines:

let triangles (points: seq<Point>) =
    let points = points |> Array.ofSeq
    let len = points.Length
    seq {
        for p1 in 0 .. len - 3 do
            for p2 in (p1 + 1) .. len - 2 do
                for p3 in (p2 + 1) .. len - 1 ->
                    {
                        A = points.[p1]
                        B = points.[p2]
                        C = points.[p3]
                    }
        }

That takes care of (1), now let’s consider (3). How can we find a rectangular box that contains a collection of circles?

In our earlier version, we proceeded by finding the left-most and right-most points in our collection, which define the left and right boundary of the rectangular box, and did the same for the top and bottom-most points.

We can reduce our new problem to that same problem, by converting each circle into 4 points that box it from the left, right, top and bottom, like so:

let xs =
    circles
    |> Seq.collect (fun circle ->
        seq {
            circle.Center.X - circle.Radius
            circle.Center.X + circle.Radius }
        )
    |> Array.ofSeq
let ys =
    circles
    |> Seq.collect (fun circle ->
        seq {
            circle.Center.Y - circle.Radius
            circle.Center.Y + circle.Radius
            }
        )
    |> Array.ofSeq

let xMin = xs |> Array.min
let xMax = xs |> Array.max

let yMin = ys |> Array.min
let yMax = ys |> Array.max

Essentially, we create a box around the circles first, and then find a box that contains all these boxes.

And… we are pretty much done at that point, all we need to do is bolt these 4 steps together:

let superTriangle (points: seq<Point>) =
    let points = points |> Array.ofSeq
    let len = points.Length
    let circles =
        // 1) enumerate all triangles
        seq {
            for p1 in 0 .. len - 3 do
                for p2 in (p1 + 1) .. len - 2 do
                    for p3 in (p2 + 1) .. len - 1 ->
                        {
                            A = points.[p1]
                            B = points.[p2]
                            C = points.[p3]
                        }
            }
        // 2) compute their circumcircles
        |> Seq.map circumCircle
    // 3) convert the circles into points that
    // box them left, right, top and bottom
    let xs =
        circles
        |> Seq.collect (fun circle ->
            seq {
                circle.Center.X - circle.Radius
                circle.Center.X + circle.Radius }
            )
        |> Array.ofSeq
    let ys =
        circles
        |> Seq.collect (fun circle ->
            seq {
                circle.Center.Y - circle.Radius
                circle.Center.Y + circle.Radius
                }
            )
        |> Array.ofSeq

    let xMin = xs |> Array.min
    let xMax = xs |> Array.max

    let yMin = ys |> Array.min
    let yMax = ys |> Array.max

    // 4) back to original algorithm: 
    // find triangle around the bounding boxs
    // Omitted, same as before

And we are done!

Parting thoughts

Does it work? Yes, it appears so. It’s not a proof, but when running the modified version, the issues we had previously are gone. Now we have a nicely convex triangulation!

The obvious drawback here is that the initial super-triangle computation went from O(n) to O(n^3). Instead of computing the bounding box in a single pass over the points, we need to compute every possible triangle first. That being said, correct and slow beats fast and wrong.

Another approach that was suggested involved using a super triangle with infinite coordinates. I am still wrapping my head around what that means practically, and might look into it later. There is also still the issue of points that could be aligned. But… these questions can wait, this is where we’ll stop for today! And thank you again, Stranger on the Internet, for helping me get un-stuck!