From Chaos to Continents: Voronoi-Based Procedural Terrain Generation with Hydrology and 3D Visualization

To create an irregular mesh, we partition the plane using a Voronoi diagram [8]. Given a set of points $S$ as sites, the Voronoi region $V_p$ corresponding to a site $p\in S$ contains all points in the plane closer to $p$ than to any other site in $S$.

We construct this mesh via its dual graph, the Delaunay triangulation. We utilize the delaunator library [12] to efficiently compute the triangulation. The circumcenters of the resulting Delaunay triangles serve as the vertices of the Voronoi polygons, while the edges between them define the adjacency relationships between cells.

Voronoi sites are generated uniformly at random, though naive sampling produces uneven distributions with high variance in cell area and edge length. To regularize the mesh, we apply Lloyd’s Relaxation [11], which iteratively converges to a Centroidal Voronoi Tessellation [9]. This redistributes sites more evenly, producing cells that are relatively uniform in area while retaining irregular boundaries; see Figure 1.

To generate elevation in our terrain, we use Perlin noise [14] implemented in [5], which provides spatially coherent values unlike uncorrelated sampling. We introduce a Land Anchor system to control terrain shape. Rather than a single radial mask for island formation, our method generalizes to $N$ arbitrary anchor points. The final elevation $E$ of a cell with Voronoi site at $(x,y)$ is defined as $E_{final}=\text{Perlin}(x,y)-{\left(d_{min}/R\right)}^2$ where $d_{min}$ is the Euclidean distance to the nearest land anchor, and $R$ is the map radius. This approach concentrates high elevation values around specific focal points while preserving the irregular, jagged nature of the noise along the coastlines; see Figure 2. Subsequently, a piecewise power curve is applied to the elevation data to compress low-lying areas into flat beaches and exponentially scale higher values into steep mountain peaks.

To texture the terrain, we simulate moisture using a secondary Perlin noise map modified by proximity to water. Cells are classified into biomes (e.g., Desert, Forest, Snow) using a Whittaker Biome Diagram [17], which maps $(Elevation,Moisture)$ pairs to terrain types.

Discrete biome classification can result in abrupt transitions between distinct zones; see Figure 3 (left). To mitigate this, we implement a smoothing pipeline:

1. 1.

   Generate a $256\times 256$ lookup table representing the Whittaker classification.

2. 2.

   Apply SciPy’s [16] Gaussian Blur ($\sigma =8$) to the table to generate boundary gradients.

3. 3.

   Assign cell colors using bilinear interpolation on this blurred table.

This results in gradient transition zones, such as a grassland buffer separating a forest from a beach; see Figure 3 (right).

We simulate hydrology by utilizing the graph structure of the mesh, defining river channels along Voronoi edges. A priority-flood algorithm [1] is implemented to model flow dynamics:

1. 1.

   Flow Routing: We construct a directed graph where each vertex connects to its lowest neighbor, creating a continuous path to the ocean.

2. 2.

   Flow Accumulation: We traverse the graph from local maxima to the coast, aggregating flow values at each step.

The resulting accumulation values determine the width of the river segments, generating a river network structure; see Figure 4.