Implementing Bowyer-Watson algorithm for delaunay triangulation

Question

I am trying to implement the following Bowyer-Watson algorithm to implement Delaunay Triangulation.

function BowyerWatson (pointList)
  // pointList is a set of coordinates defining the points to be triangulated
  triangulation := empty triangle mesh data structure
  add super-triangle to triangulation // must be large enough to completely contain all the points in pointList
  for each point in pointList do // add all the points one at a time to the triangulation
     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
     polygon := empty set
     for each triangle in badTriangles do // find the boundary of the polygonal hole
        for each edge in triangle do
           if edge is not shared by any other triangles in badTriangles
              add edge to polygon
     for each triangle in badTriangles do // remove them from the data structure
        remove triangle from triangulation
     for each edge in polygon do // re-triangulate the polygonal hole
        newTri := form a triangle from edge to point
        add newTri to triangulation
  for each triangle in triangulation // done inserting points, now clean up
     if triangle contains a vertex from original super-triangle
        remove triangle from triangulation
  return triangulation

The algorithm takes O(N log N) operations to triangulate N points as claimed. But is there any way to calculate it from the above algorithm? I mean which part of the above algorithm takes log N times, which results in (N log N) for n points? special degenerate cases exist where this goes up to O(N2) as written in wikipedia. What is the meaning of this degenerate case?

Answer 1

The pseudocode in the question, the Bowyer/Watson algorithm as described in the Wikipedia article, is O(n^2). The article does state this. The Wikipedia article is less clear -- and I would argue is in fact wrong -- about what one needs to do to make the algorithm O(n log n), doing so is not totally trivial.

The current Wikipedia article says the following

Efficiency can be improved in a number of ways. For example, the triangle connectivity can be used to locate the triangles which contain the new point in their circumcircle, without having to check all of the triangles - by doing so we can decrease time complexity to O(n log n).

I dispute that "triangle connectivity" alone can be used to achieve O(n log n) behavior.

What one needs to do at each iteration of the algorithm is given a new point P find the "bad triangles" that introducing P to the existing triangulation will create, the triangles which contain P in their circumcircles. The "triangle connectivity" property alluded to in the Wikipedia article is that these bad triangles will be a connected component in the adjacency graph of the triangulation so if you happen to know one such bad triangle you can find the others easily by performing a graph traversal constrained only to adjacent bad triangles. However, that is if you have one bad triangle to start from. An obvious choice would be the triangle that contains P.

The primary source of the O(n log n) time complexity claim, "Efficient Unstructured Mesh Generation by Means of Delaunay Triangulation and Bowyer-Watson Algorithm", S. Rebay, 1991, explains the situation as follows: (emphasis mine)

The algorithm efficiency depends on how quickly the search for triangles to be deleted at each point insertion is performed and this is made much easier by knowledge of the neighboring triangles to each triangle. In fact, since all the triangles to be deleted are always contiguous, a tree search among neighboring triangles can be used to find all the other triangles to be deleted after the first one. In the typical case, the number of triangles to be deleted at each point insertion does not depend on the number of all existing triangles. As a consequence, if the information pertaining to the the neighboring triangles is available and an O(log N) multidimensional search for the first triangle to be deleted is employed, the algorithm can compute the Delaunay triangulation of a set of N points in O(N log N) operations. In special cases, however, the number of triangles to be deleted at each point insertion can be very large. In the worst possible situation, when all existing triangles have to be deleted at each point insertion, the operation count of the Bowyer-Watson algorithm degrades to O(N^2).

To paraphrase, if you maintain the triangle adjacency graph and you have some way of finding the triangle containing an arbitrary point in O(log n) time then the Bowyer-Watson algorithm is O(n log n) because the number of triangles in a bad triangle component is going to be so small for typical triangulations that we can consider it to be lower than some small constant.

Wikipedia, Rebay, etc. do not discuss the "multi-dimensional search" you need to make the algorithm O(n log n) but I can offer three suggestions, from the literature and otherwise:

1. Store the bounding rectangles of all the triangles in the 2D range search data structure of your choice. For example, an R tree or an R* tree.

2. Use the hierarchy of the Delaunay Triangulation itself. Maintain the history of the triangulation, a tree of the "levels" of triangles. I believe that this is what CGAL's Delaunay implementation does.

3. Randomly sample k triangles, find the closest triangle t to the target point in the sample, perform an A* search of the adjacency graph from t to the triangle that contains the target point. This would probably perform well but is not formally O(log n).

Answer 2

Maybe it is a bit late to answer, but it might be useful to someone else.

A Delaunay triangulation have a circumcircle property, hence no point of the Delaunay triangulation can lie within the circumscribed circle of any triangle. The Bowyer-Watson algorithm add a point that does not verify this property. It can be shown that all triangles which circumscribed circle contain the new point are contiguous.

To obtain the theoretical NlogN you must use the contiguous fact. (I use a connectivity table)

=> You need to search for a first triangle in your triangulation where the circumcirle property is not respected( complexity log(N) )

=> Once you have it, delete contigued triangles (using connectivity table) (independent of total number of nodes)

=> Build new triangles (and update the table) (independent of total number of nodes)

And you need to do it N times, for each nodes. This result in a theoretical Nlog(N) complexity.

The algorithm given on wikipedia looks like a loop of loop of N size. Hence, it should automatically gave you a N^2 complexity.

This complexity should be the worst case scenario were all elements are deleted and no connectivity is available, as said Ripi2. In such case, you would have no choice but to search in all triangles. See https://doi.org/10.1006/jcph.1993.1097 for more informations

===========================================================================

You need to found a initial triangle where the point P is ( hence it is a triangle to be remeshed):

Choose a triangle T of your triangulation (if possible close to your point. See below for more). Compute G, barycenter of T. We denote t1,t2,t3 the edges of triangle T. L1 L2 and L3 are the triangles linked to T by edges t1, t2 and t3 (adjacency)

if [GP] cross t1, t2 or t3, define T=L1,L2,L3 respectively => Loop with new T if not : => You have found the triangle containing P

It is also possible to choose a "good initial guess" for T, using a initial coarse regular subdivision, and storing nodes of your triangles. This is similar to the idea 1 of Jwezorek below, but probably cost less:

Because you already have a good "box" around all your points, you can create a multidimensional box subdivided log(N)-th time. We will supposed that you have all the nodes of your Delaunay triangulation inside that subdivided box. When you search for a T close to the new point P, insert P in your subdivided box, and then locally search around it. Using one of the first node you found, use your connectiviy table to get a first good triangle T.

I do recommand that book, for those who wants more infos, and have time for reading: "Delaunay triangulation and meshing, Paul-Louis George, Houman Borouchaki." Edited answer due to Jwezorek questions

Answer 3

The algorithm described here is a naive incremental delaunay triangulation. The complexity is clearly O(N^2). "for each point in pointList do // add all the points one at a time to the triangulation" this main loop is O(N) and "for each triangle in triangulation do // first find all the triangles that are no longer valid due to the insertion" this inner loop is O(N).

And this algorithm is O(N^2) in all scenarios. However, there are other algorithms that do give O(N * Log(N)) or better complexity. See "Divide and Conquer" from here: https://en.wikipedia.org/wiki/Delaunay_triangulation

Answer 4

There are two main steps in the algorithm:

1. for each triangle in triangulation, if point is inside circumcircle ...
2. for each triangle in badTriangles, remove and add new tringles ...

The first step is a bad approach. There are better ways to locate a triangle where the new point to add lies on. Search for "walking" from an old triangle towards new point.

The second step can involve all triangles for a degenerated case where all points lie in a common big circle.

Answer 5

Anyway, thanks to the initial pseudo-code, very clear, I had fun implementing it with Scratch : https://scratch.mit.edu/projects/351857925

So thank you :-)