My latest video is live on YouTube! In this beautifully animated tour of Voronoi diagrams, Delaunay triangulation, and natural neighbor interpolation, I explain the algorithms used in my mb-geometry rubygem for generating Delaunay triangulations and Voronoi partitions. Then interpolation is given as a practical application of these algorithms.
This video is part of my ongoing video series, Code, Sound & Surround.
You can keep reading for a partial transcript of the video, and I’ve written more about Voronoi partitions, Delaunay triangulations, and mb-geometry previously.
Video transcript
This is a partial transcript of the video, excluding the intro and outro sections. This is taken from my original script, so may differ from the video where I improvised when recording, or rearranged things when editing.
Voronoi partitioning overview
Let’s use animation to get an intuitive feel for how Voronoi partitions and Delaunay triangulations behave. Feel free to use the video timestamps to jump to sections that interest you, but I do recommend watching everything in order.
Both of these algorithms are given a set of X/Y coordinates as inputs, and produce polygons based on those input points as outputs.
We’ll start with Voronoi diagrams. A Voronoi partition or Voronoi diagram splits the plane into convex polygons showing the regions of the plane closest to each input point. These regions are called cells. We often give each cell a different color to emphasize the boundaries between polygons.
With just one point, the entire plane is a solid color, showing that all points are closer to this input point than any other (because there is no other).
With two points, the plane is split in half. We can continue the pattern to three points, and so on.
Let’s go back to three points, and see what happens when we go from a line to a triangle. A square, a pentagon.
We can also do grids, both rectangular and hexagonal.
And we can get more creative, too.
Delaunay triangulation overview
Now let’s look at Delaunay triangulation. Delaunay triangulation gives us a set of triangles connecting all of our points together, without overlap, and avoiding really thin triangles as much as possible.
Voronoi diagrams and Delaunay triangulations are closely related to one another. Every edge of a Delaunay triangle is perpendicular to an edge of a Voronoi polygon. And if we draw a circle through the corners of each triangle, called a circumcircle, the center of each Delaunay circle, or circumcenter, corresponds to a vertex of a Voronoi polygon.
Here are some of the previous animations, but with the Delaunay triangulation displayed on top of the Voronoi diagram. Notice that triangulation isn’t possible until our points are no longer in a single line.
Delaunay triangulation algorithm details
Now that we’ve started building our visual intuition for how these algorithms behave, let’s dive into the details of Delaunay triangulation. There are multiple algorithms that can produce a Delaunay triangulation, some faster than others. But they all follow the same basic rule – no input point may lie within the circumcircle of any triangle, which prevents overlapping triangles and avoids very thin triangles.
I chose to use Ruby to implement a 1980 algorithm by Lee and Schachter, filling in the gaps with some small modifications. And, you have my apologies if I mispronounce any names in this video.
This is called a divide-and-conquer algorithm because it splits the input points into smaller and smaller groups (divide), then builds those groups back up into a full triangulation (conquer). Each subgroup forms what’s called a convex hull, with the outer boundary points connected in a counter-clockwise loop.
Let’s go step by step through a simple triangulation of 5 points. First the points are passed into the triangulation algorithm and sorted.
Then, we split the points in half along the X axis, using the Y axis to break any ties. The dashed vertical line shows where our code decided to split the points.
After this we look at the left half. If it has more than three points, we continue splitting. If it has three or fewer points, we can create the counterclockwise convex hull directly.
As we have two points in this case, we just connect them left-to-right and then right-to-left.
Then we look at the right half of the split.
Since we have a triangle, we connect the points in counterclockwise order. If we had three collinear points, we would instead connect them left-to-right and right-to-left, as we did with two points.
Now we have two convex hulls to merge.
We first use variations on the methods described in Lee and Schachter’s paper to find the lower and upper connecting tangents between the hulls.
This algorithm starts by testing the line between the rightmost and leftmost points of the two hulls, then walks step-by-step downward to find the lower tangent.
The same is done for the upper tangent, but walking upward instead of downward.
Now the lower and upper tangents are highlighted with dashed lines.
Having identified the bottom and top connecting segments, we use the merging algorithm from the paper to delete and create edges so that the Delaunay condition, that is no points within any triangle’s circumcircle, is satisfied for the merged points.
The algorithm works its way up between the two hulls, starting at the lower tangent. If the circumcircle of test triangle P1,P2,P3 contains the query point Q, an edge is deleted. This is indicated here by a yellow circle.
Otherwise, if all of the query points are outside the test circumcircle, an edge might be created. This is shown by a green circle.
This continues until the upper tangent is reached, at which point the merged convex hulls contain a valid Delaunay triangulation. The two hulls are then merged into one.
Now, for this relatively simple triangulation requiring only one subdivision, the algorithm is complete.
So let’s look at something a little more complex.
With a larger number of points, we have to split multiple times to get to just 2 or 3 points in a subgroup.
The algorithm runs the same way it did for 5 points, with the same walk up and down to find connecting tangents, and the same merging steps are repeated to combine subgroups into ever larger convex hulls.
When we run into a subgroup of three collinear points we can see that our algorithm links them in a line, just like it would have for two points.
When merging these collinear hulls, we only have to link the two closest points. The merging algorithm handles this implicitly.
Vertically collinear points are handled the same way.
Subgroups continue to merge until we reach the final merging step.
And the algorithm is complete.
Let’s watch one more triangulation, at full speed.
There you have it. That’s the Lee and Schachter divide-and-conquer algorithm for Delaunay triangulation, as implemented in my mb-geometry rubygem. You can view this code on GitHub.
Voronoi partitioning algorithm details
Earlier in the video we saw that there is a close relationship between Voronoi diagrams and Delaunay triangulations. It turns out that you can calculate one from the other. In my code, I opted to calculate the Delaunay triangulation first, and then derive the Voronoi partition from that. There are other approaches that might be better or worse depending on how the code will be used.
The first step in generating a Voronoi diagram this way is constructing the Delaunay triangulation. To avoid infinite polygons at the edges, we pick a bounding box to use as our display area, and reflect the input points across the walls of this bounding box. Because Voronoi polygon edges always occur half way between two input points, this guarantees we will have polygon edges at our bounding box.
Now with this expanded triangulation, we collect every triangle’s circumcenter, as these mark the vertices of our Voronoi polygons. Due to limits on floating point precision, sometimes vertices that should be in the same place aren’t. So we have to find all vertices that are close enough together that they should probably be in the same place, and merge them.
Next, we look at all the triangles connected to a given input point. We place their circumcenters into a counterclockwise list for that input point. That list of circumcenters forms the boundary for that input point’s Voronoi polygon.
Do this for every input point, and we have our Voronoi diagram.
Natural neighbor interpolation
Voronoi diagrams look nice, but they also have practical applications. As one example they are sometimes used in geospatial information systems. Suppose you take a bunch of measurements at random locations, let’s say elevation. A Voronoi diagram can tell you what the measurement was at the closest sample to any given point.
This is called nearest neighbor interpolation, and if you’ve worked with image scaling, that term “nearest neighbor” might sound familiar. Basically, if we have a bunch of samples, and we fill in any gaps between the samples with whatever sample was the closest, that’s nearest neighbor interpolation.
What if we want to blend between the sample points to approximate the value for a place we didn’t measure directly? In general, this is called interpolation, and there are lots of ways to do it. Some methods of interpolation only work on perfect rectangular grids, while others can work with random points. We’ll be looking at methods for interpolating random points.
In computer graphics, a linear blend of colors between the corners of a polygon is called Gouraud shading. In other disciplines this might be called linear barycentric interpolation.
So we can use our graphics system to perform simple linear interpolation on our dataset – that is, connecting our samples together with straight lines. We’ll draw the Delaunay triangulation using Gouraud shading.
This runs fast and obviously gives us a smoother result than nearest neighbor, but still leaves hard edges, and the interpolation can only ever blend between three input points, even if more of the input samples could have influenced the true value at our interpolated point. Also, the interpolation is fairly sensitive to changes in the input points.
We can do a little bit better for some use cases by interpolating with Voronoi polygons instead of Delaunay triangles. This is still limited in the number of input points that can influence an interpolated measurement, still has sharp edges, and still jumps around if input points are moved slightly.
We can do better still. Let’s jump back to our Voronoi diagram. If we add a new point somewhere to our graph, that point gets its own cell and that cell takes area away from the existing cells. We can use the amount of area taken from each existing cell to generate a weighted average. The lines shown here represent the weight assigned to each averaged cell.
This is called Natural Neighbor Interpolation.
We can move our new cell anywhere within our Voronoi diagram, and get an interpolated estimate for whatever value we are measuring – color, height, etc.
If we repeat this process along a grid of samples, we can plot a smooth, interpolated representation of our input points. This is very slow, but produces very nice results.
Looking at height again, natural neighbor interpolation gives us a very smooth interpolated curve. The sampled value at an input point is exactly that input value, so known inputs are preserved. And all neighboring input points are included in the average, so no relevant data is left out.
If we compare the different interpolation methods so far, natural neighbor interpolation provides the smoothest curves without changing known measurement points, but comes with a significant speed disadvantage.
Because it’s slow, we can choose to sample it at sparser intervals, and then use simpler linear interpolation or Gouraud shading from there, trading accuracy for speed.
Haha, that looks really freaking cool!
And that’s natural neighbor interpolation! There are many other types of interpolation I haven’t mentioned, so you can search online for spatial interpolation if you want to learn more, and I’ll put some links in the description.
By the way I’ll be applying natural neighbor interpolation to sound in a future video, so stay tuned for that.