Voroni Diagrams and Graphs IV: Efficient Computation of the Graph and it's simplexes.

Before we jump into algorithms, a few things. The algorithm I'm about to lay out is a version of the Bowyer-Watson algorithm I've referenced in a previous post. It's been pretty heavily modified since, and it was not my initial intent for it to be thins similar. Nonetheless, I think that my implementation has better asymptotic time than any implementation of the original I've yet seen, since (as we'll see) the most intensive parts of the algorithm (most of which were related to searches over the point set) have been offloaded to a simple graph search. We'll also be looking at things through the lens of the Voroni Graph, rather than the Delaunay Triangulation.



This brings us to a second key difference between what I'm about to propose and Bowyer-Watson, and it's not as nice as an improved asymptotic time. While Bowyer-Watson returns exactly the set of simplexes which are the Delaunay triangulation, what I'm about to lay out will find all simplexes in the Voroni Graph. These simplexes include the entire Delaunay triangulation, but it may also contain a few extras. While in large real world data sets, I suspect that the occurrence of these extra triangles will be fairly rare, it is pretty easy to cook up examples which will produce arbitrarily many extra triangles. The simplest ones we covered in the last post in this series.

I've looked for a way to naturally check for these triangles during the course of the algorithm, and I've found a few, but I haven't really explored what they do to average run time (I'm not optimistic, since I think the required circumsphere checks are fairly expensive.), so I'll hold off explaining them until a later post when I've implemented enough of this to actually do some testing. In the mean time, I think the simplest thing to do is just to check each simplex we get afterwards, since, again, I think that for a realistic data sets, the incidence of these false positives will be very low.

On to the actual algorithm. In the first post of this series, we laid out a fairly inefficient method for creating the Voroni graph. In essence, for each point, we thought about combining every possible Voroni graph the point could be in, then we eliminated connections that we found were redundant. If you tried to implement this algorithm naively, you probably noticed that it didn't scale very well. The reason for this is pretty simple. For each of n points, we found (n-1) simple graphs. For each of these new graphs, we compared it to (on average) a less than (n-1)/2 other graphs to decide if it was indeed part of the cell of our point. Therefore, we'd expect our runtime to be O(n³). There's a lot of checking going on under the hood there, and while this isn't a horrible asymptotic time, it's not a terribly practical one either. Using the properties of the Voroni graph itself, we can do much better.

There's also an even bigger problem with the naive approach. It doesn't help us find the Delaunay  Triangulation of the set at all. Even once we have the Voroni Graph, we will still need to search it for triangles. With a naive approach, this operation is O(n³).  While there are obviously improvements we can make to a naive guess and check approach, I claim this is also unnecessary.

I propose the following approach: the addition of a single point to set of points will only change that set's Voroni graph near the new point, and it will only change the graph in very specific ways. For instance, the addition of a new point will never cause the formation of an edge between two other points that previously did not have an edge. Since we can check to see if there will be an edge between any point in the graph and any new point ( We laid out the groundwork in part I of these posts. ), we can easily create both the Voroni Graph and the exhaustive list of (candidate) Delaunay triangles in this graph as follows:

1. Initialize an empty list of triangles (or similar) T, and an empty Graph V.
2. For each point p in our list of points to add to the graph...
1. Use a greedy search to find a vertex in V ( call it q ) which will have an edge with p. If there is no such q ( A.K.A. V is empty ) continue to the next point .
2. (Recursive) Add the edge between p and q. Check the neighbors of q to see if they have edges with p and or no longer form edges with q.
1. For each neighbor q' that no longer forms an edge with q
1. remove any triangles involving the edge qq'.
2. For each q' that does form an edge
1. check the neighbors of q' for all of the above, and so on.
3. Look for new triangles involving p.

That's it! If you've ever seen the Bowyer-Watson algorithm, you'll probably recognize what I've done above. I'd go so far as to say that what I've written down above IS Bowyer-Watson, just mathematically rephrased to take advantage of an associated graph.

It's worth noting that step 2, while, as we'll see, constant time, it's not as trivial as I've made it out to be, especially if we are trying to optimize our implementation. We'll devote the next post more or less exclusively to this topic.

Back on topic: What kind of runtimes do we get out of this scheme? Let's review what we know figured out in the second post of this series:

• Each node of a sufficiently large graph has on average 6 neighbors. ( AKA The number of edges is O(1) ).
• (Conjectured) Given a 2D point (and no other information), it takes O( n^1/2 ) time to find a node that will have an edge with that point.

Suppose we have a list of n points for which we wish to find T and V.

Clearly, step 1 is O( 1 ). The asymptotic runtime of step 2, can be found as follows: suppose we are adding a point when V has m elements (m is still large enough we can neglect special cases). We'll take O( m ^1/2 ) steps to find the first point with an edge, but because the number of neighbors of each node is constant, the rest should be O( 1 ) (assuming, of course, that are triangles are stored in an appropriately fast data structure for this sort of. If they are not, we may not be O( 1 )).

To approximate the total runtime of step two, let's use our definition for big-O notation to say that each run will take less than a * m^1/2 for some constant a. Summing this from m = 0 up until m = n,

we will get an expression for how long step two will complete. Clearly, we've got at worst O( n^3/2 ) time, since we can use that m^1/2 <= n^1/2 to get a total time less than a * n^3/2. This may seem like a big, rather unfavorable approximation, but I think it's just about the best we can do without mucking things up. I'll put a link to the Wolfram Alpha Query you can use to gain some insight here. Look at the expansion for infinite n in particular. Notice that the highest power is n^3/2. It's those negative terms that are problematic, but likely not significant.

Also, if we need to check each triangle in T to make sure it is Delaunay, however, we'll get an additional step that is O(n²), since there are O(n) triangles at the end of the day, and without additional cleverness, each circumcircle check is itself O(n). Honestly, this is probably the limiting step for this algorithm, and it's something that, as previously mentioned, needs a lot of work.

Regardless, even though O( n^3/2 ) to find the Voroni Graph is a lot better than O( n³ ), I claim that we can probably pare it down even further to O( n * log( n ) ) with some work. If all goes as planned, we'll talk about my reasons for that claim (and potentially an implementation) in the next post.

Also in the next post(s):

• I'll (hopefully) have some C++ code on github to demonstrate all this.
• I'll post some test data which I've collected from several sources, including the original scans from the Stanford Bunny, all in a common format we'll be using from here on out.
• We'll talk about extending the algorithm presented here to allow us to merge two (or more) preexisting Voroni Graphs together.
• We'll also start to talk about practical applications for the Delaunay triangulation and closely related ideas.