A few years ago, I became interested in Voronoi tile generation, and built a system for generating spherical tile maps. Since then, I’ve revisited the project periodically to make optimizations and perform experiments, and this month I managed to completely eliminate the most expensive step while at the same time improving the efficiency of a second step.

The project’s inspiration is the tile maps used by Civilization and other 4X games. My goal is to replace the artificiality of a perfect hex-grid with something that feels more organic. ‘Blue Noise’ is a term that is sometimes used to describe a distribution that is random, but exhibits some uniformity (e.g. all the points being roughly evenly spaced). The basic objective of the system is to produce a distribution of tiles that guarantees that tiles observe a minimal spacing and are all close to the same size, but are still visibly irregular.

My original system broke down into four discrete steps:

1. Generate a random distribution of points on a sphere.
2. Relax the points using repulsion.
   1. Repeat until the point distribution reaches the desired level of uniformity.
3. Run the Bowyer-Watson algorithm to generate a Delaunay Triangulation from the points.
4. Generate a mesh from the Delaunay Triangulation.

This approach works, but the technical side of my brain has never liked the need for the iterative relaxation step. For one thing, it is computationally expensive for large numbers of points, coming in at O(n²/2) per iteration, and typically requiring several iterations to achieve the desired degree of uniformity.

Triangular Subdivision

Subdivision is not exactly a novel procedural generation technique, but for this use-case it enables a spectacular boost in efficiency. My original process uses the Bowyer-Watson Triangulation implementation from an excellent open source library created by Rafael Kübler da Silva.

The Bowyer-Watson algorithm has no information about where each new point is within the triangle graph, and must check the point against each existing triangle to see if it needs to be replaced. The number of triangles approaches 2x the total number of points over the course of the algorithm, for an efficiency of O(n²).

So rather than hand the Bowyer-Watson algorithm a complete array of pre-relaxed points, I generate new point after each Bowyer-Watson iteration by choosing a position within the current largest triangle (longest perimeter). The new generation process looks like this:

1. Generate an initial tetrahedron or octahedron of triangles.
2. Perform iterative subdivision of triangles (modified Bowyer-Watson).
3. Generate Voronoi mesh as before.

Because I know that the new point is within a specific triangle, the search for triangles in need of replacement can be confined to that triangle’s immediate neighbors. I’m not worried about obtuse triangles at a distance, since if such a triangle existed, it would likely have been chosen for subdivision ahead of the current one, or will be smoothed out soon after.

This brings the efficiency of my approach to O(n * log(n)). The log(n) being the cost of the heap used to track which triangle currently has the longest perimeter. It occurs to me that I don’t actually need to always subdivide the largest triangle, I just need to continue subdividing until all triangles sides are below a desired threshold. With further modification I might be able to achieve an efficiency of O(n), but I really can’t complain about my current progress.