Delaunay+Voronoi on a sphere


The second step is to construct a Delaunay Triangulation on these points on a sphere. A standard Delaunay library works on points in a 2D plane. Surprisingly, we can use existing 2D Delaunay libraries to run on points on a sphere. The key idea is to transform the data before running the algorithm.

  1. Project the points from the sphere onto an infinite plane using a stereographic projection. This maps the northern hemisphere onto points inside the unit circle, and the southern hemisphere onto points outside the unit circle.
  2. Run the unmodified Delaunay triangulation library on the points on the infinite plane. This general idea is useful — instead of modifying an algorithm, you can often modify the input data, run an unmodified algorithm, and then modify the output data.
  3. Wrap the results from the infinite plane back onto the sphere.

There will be a hole left over. The outer boundary of the Delaunay Triangulation is the convex hull of those points. That convex hull ends up wrapping into a small polygon near the south pole of the sphere. See this animation I made (also available in animated GIF and Quicktime movie formats).

We need to stitch this hole back together. It turns out I already do this kind of wrapping for my dual-mesh library, for different reasons, and I adapted that code to work here. The best way to do this is:

  1. Rotate all the points so that the last point in the array is on the south pole. (see stackoverflow for math details)
  2. Run Delaunay without that point, because the stereographic projection will project that point out to infinity.
  3. Stitch that point back in: for each half-edge on the convex hull, construct a triangle that connects the two points in that half-edge to the chosen south pole point. If there are N half-edges on the hull, this will produce N new triangles and 3N new half-edges.

However, rotating the points is not strictly necessary. I decided to run on the unrotated sphere, and add a new point at the south pole instead of using an existing one. It’s Good Enough For Now. I’m thinking that I’ll have half a million regions for my map generator, so one extra region at the pole won’t be noticeable. That said, the sphere rotate approach doesn’t look that hard, so I should probably implement it at some point.

Here’s the resulting Delaunay Triangulation:

Delaunay Triangulation on the sphere

Note: there’s a library that does all of this for you! See d3-geo-voronoi. I didn’t end up using it because my mapgen4 code runs on the data structures from Delaunator, and it was going to be more work to adapt it to use d3-geo-voronoi than to implement the projection and stitching myself. The projection is around 10 lines of code, and the stitching I already had implemented. My stitching approach seems to be a little different from d3-geo-voronoi’s; he adds four infinity points and I don’t add any. His is probably more robust.