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.
- 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.
- 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.
- 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:
- Rotate all the points so that the last point in the array is on the south pole. (see stackoverflow for math details)
- Run Delaunay without that point, because the stereographic projection will project that point out to infinity.
- 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:
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.