Bowyer-Watson algorithm: how to fill "holes" left by removing triangles with super triangle vertices

Question

I am implementing the Bowyer-Watson algorithm as presented at Wikipedia. In my implementation, everything works as I would expect up until the last part of the pseudocode:

for each triangle in triangulation // done inserting points, now clean up
     if triangle contains a vertex from original super-triangle
        remove triangle from triangulation

If I follow the pseudocode here literally, I can end up with missing triangles in my Delaunay triangulation.

As an example, please consider the images below. The sites I am triangulating are rendered as blue circles. The triangles are rendered with black lines (excluding the image borders) and connect sites or bounding/super triangle vertices. The circumcircles are rendered with gray and their centers are rendered with red circles. The Voronoi cells are each painted with a different color to (hopefully) make the problem more apparent.

This image shows the state of the triangulation right before performing the steps listed in the pseudocode above. Note that two of the super triangle's vertices are beyond the right and the bottom of the image.

This image shows the step after removing any triangles that contain super triangle vertices without any further considerations:

The top three vertices should have a new triangle with a circumcenter at the point where the greenish/brownish cells meet. The problem is that the corner vertex that was shown in the "before" image was inside this circumcircle, so the regular processing of the algorithm never generated this triangle.

How do I express this edge case in pseudocode so I can check for and solve it? I would like to avoid some horrific "try every combination of sites that shared a triangle with a super triangle vertex for valid circumcircle" loop.

I read the Bowyer and Watson papers a couple years back and will read them again for my answer if necessary. I was hoping that (1) somebody else might have the answer available and (2) I could use Stack Overflow to look the answer up if I ever run into this question again.

---

Edit

So I have found a relatively cheap but imperfect work-around. My super triangle is programmatically determined to surround the sites' bounding box without intersecting its sides. This idea was caused by all sorts of frustrating problems with Java considering some of my calculated circumcenter coordinates or distances between coordinates to be infinite. This caution led me to make my super triangle so small that its vertices sometimes fell in valid triangles' circumcenters. Increasing the size of the super triangle has made the problem seem to disappear. However, it is possible that a triangle on the convex hull can be so obtuse that one of the vertices still could fall inside a valid circumcircle.

I think this means that my initial question is still valid in the face of floating-point number limitations. Is there a cheap way to guarantee that the Bowyer-Watson algorithm generates a valid triangulation?

Answer 1

I encountered the same problem when i was implementing Bowyer-Watson algorithm described here: http://paulbourke.net/papers/triangulate/. I couldn't find anything helpful on internet and even asked at my university, but with no result. After a while I came up with a solution. I started with discovery that for the problem to disappear, the vertices of bounding triangle should ideally lie at infinity, which is not practical. So what does triangles circumcircle look like if triangle has one or two vertices at infinity? It is just line going through the other points. So testing if point lies in triangles circumcircle changes to testing if point lies left or right of line.

Algorithm then looks like this:

1. Check if any of triangles vertices lies at infinity. In other words: check if triangle is sharing some vertices with bounding triangle.

2. If it's sharing all three vertices: trivial.

3. If it's sharing zero vertices: classical approach - check if distance from point to circumcenter is shorter than circumradius.

4. If it's sharing one vertex: check if point lies to the left/right of line defined by the other two vertices. one vertex in infinity

5. If it's sharing two vertices: check if point lies to the left/right of line defined by these two vertices but shifted to the third point. In other words: you take only the slope vector from line between these shared vertices and shift it so that the line passes through the third point. two vertices in infinity

Testing whether point lies to the left or right of line depends on your triangles winding order.