Voronoi Diagram for Heightmap Generation
A quick way of thinking about the Voronoi diagram is that it contains a nearest neighbors graph of points with polygon construction defining lines whose construction are the closest points between two such neighboring points (called sites). The polygon is represented as a cell. Today it is probably more customary to use algorithms such as the Fortune algorithm which speeds up the process and offers greater flexibility relative to say grid/hashing computations. Though I wouldn’t go into further detail on the Fortune method, it does offer interesting geometric insights via quadratic functions, binary search trees, and sweep lines. Having any experience in Nearest Neighbor Search methods, one might be more familiar with KD trees which offer multidimensional search features and Binary search trees as an analog. The binary search tree is a backbone to the speed of the Fortune method, and thus I’d suggest deeper reading if interested as Nearest Neighbor Searches have much utility to offer beyond the scope of the Voronoi diagram. The Voronoi diagram in heightmap generation offers interesting control utility in visual graphs given the possibility of supply site points which can be then be used for tessellating polygons. The history of Voronoi diagrams, for instance, date back to its use in geographic information visualizations, for instance, by British physician John Snow in 1854 to illustrate geographic visual features of a cholera outbreak and its relation to water sources. I’d suggest further reading at https://en.wikipedia.org/wiki/Voronoi_diagram for other applications.
While constructing the Voronoi diagram today is as simple as providing cell site points in a two dimensional form (e.g. a point is {x: 1, y: 2} in a given array of such objects, rendering heightmap does have added work.
The problem here involves using such polygon to form the points in polygon which then can be used to render the heightmap. I’ve used in JavaScript a package called turf.js which provides some basic GIS functionality including computation of centroids, finding area of polygon, and point in polygon testing. Constructing a point in polygon, in this case, can be done creating a cell bounding box (min/max x and y positions of the polygon forming box), constructing all integer valued x,y tuple positions, and then point in polygon filtering these points.
Typically for each pixel on a heightmap, there will be an assigned grayscale value (most commonly it is 8 bits) or 256 different gray scale values ranging from black to white. Thus in the Voronoi diagram the problem is finding the two dimensional integer set of points in Voronoi cell, and finding its heightmap designated value. One could, of course, randomly assign, or not so randomly assign a value for all such pixels in Voronoi cell, but there are added gradations that could be provided using distance metrics. This can be done measuring the distance say from site position to point in cell. This type of visualization, for instance gives a sense of area relation, as smaller cells will vary closer to whiter intensity nearest the cell edges, while larger cells will have darker shading nearing cell edges. However there is more likely to be sensed discontinuities between cell edges. Other methods could involve distance metrics between site position to point in polygon relative to the minimum on the set of distances between site position and all cell edges nearest points.
Or in pseudo code:
distance(site to point_in_cell)/min(set of distances from site to nearest point on each polygon edge)
This renders a 2d projection of a sphere contained in such cell. A maximum on this set renders the same only that the 2d projection of a sphere is clipped by cell boundaries. Again edge discontinuities are to be appreciated in this case. Another possible distance metric involves use in constructing rays and forming a metric for normalizing the distance between site and point in polygon forming a ray where divisor in the maximum distance of the ray which is an intersecting point on cellular edge. In my code example, I use a ray distance metric and normalize the value. In this way, the edge boundaries should be equal for each ray projection, and thus eliminating edge boundary discontinuities. You could visualize these as polyhedral projections with varying degrees of white intensity falloff to site. The inverse of this would form visually a polyhedral landscape where cell edges are valley height minimums. For more aesthetic visual appeal, it is important noting that the sites of cells are not necessarily the centroid of the cell, and thus, I purposely used a distance metric from centroid of cells to the point in polygon and then applied the norming divisor mentioned above for such ray. The heightmap produced by such, does have a self similarity in field all minimums are the same and all maximums are between site point (or centroids if applicable) and cell edges.
If one wanted to vary locally a set of cells, so that relative another distinct set of cells, has translated height, how to do it is merely choice of cell set. Though again, avoiding height translation jump discontinuities, one might consider a clustering method or some manner of nearest neighbor computation with a given seed choice and then defining a metric defining the boundary of such clustering set. Perhaps, walking or node traversal (Djikstra’s method) could help if wanting to traverse cell edges and the cell points forming say the river or streams of the heightmap. In this way, one could defines a seed at river terminus from one boundary point, and then define a nodal headwater, all along adding increased complexity of branching along the way in selection. The highlands then represent, for instance, all cells that are increasingly distant from stream node points in such set construction. One could do basic interpolation for this.
Another method I’ve used, is area computation of cells. In this case formulating a parameter t for interpolation by using a norm divisor that is difference between maximum cell area and minimum cell area.
Thus a formula for area weight/rank with normalization might look like denoting minimum on the set of cell areas minCellSet, and the maximum on the set of cell areas maxCellSet:
(Input cell area- minCellset)/(maxCellSet-minCellSet)
This forms a normalized area rank for cell area distributions though it is worth noting the distribution of this area ranking form can vary, so it may be likely that when forming range sets by ranking, one might want to introduce additional statistical methods indicating more precisely rank distributions when t parameterizing interpolative methods on ranges and integrating this without additional height translating methods.
Specifically, one method could involves area weight ranking each specific edge node cumulatively by its cell edge point graph. Thus each cell forms an area rank that is cumulative summed for such edge point cell. In this case, a Voronoi diagram will include three such cells for every edge point. A simple hashmap of the string position value tuple makes tracking these points fairly easy, and thus dividing the cumulative sum by three constructs a method graduating the edge points of all cell by area ranking.
In my algorithmic example, I’ve allowed for manual adjustment of area ranking selection ranges. For instance, low ranking area clusters are linearly interpolated between smaller height value ranges relative to mid ranking area clusters, and the highest ranking area clusters (usually the edges of the diagram) are interpolated between smaller height value ranges. We can call this area rank_height mapping Arank_height_map.
Each point in cell thus will be determinative to an overall translated height not only by its relation to site (or centroid as preferred), but also its nearest edge point set. Parameterizing t for overall height translation is done by computing the linearly interpolated height between nearest cell edge ray intersection interpolated position, and then computing via Vector math the t value of t in the formula B +(A-B)t where A and B are the nearest edge points on cell. I’ve used a bezier computation for computing line to line intersections. In this case, ray from centroid (or site) to point in polygon and segment AB mentioned above. Supplying t to the linear interpolation between Arank_height_map of A and B respectively gives us an overall height translation for point in polygon.
One can proceed with all of the tools mentioned above with randomly assigned site positions to generate a heightmap. However, given means to measure cell area and apply this to height scaling also implies that site selection can be clustered, for instance, with higher density, or conversely which then lead to smaller or larger area cell distribution types that are used in simulating terrain feature types.
For rivers and streams, for instance, one could formulate a random walk. In this case selecting two random points. Computing the vector difference (A-B), computes the direction of the segment between such points, and then using random rotation angle, allows variance of the course. Step increments can be randomly assigned in so far as the magnitude length of this vector. Iteratively at each step, re evaluating the present position and destination is done until the random walk closes required distance. Branching subsidiary streams, can be performed, using random selection sets of points on the primary river walk. In this case repeating the process, however, using a probability of rotations within desired ranges and factoring the primary river’s overall direction. The overall branching start and end positions distance is fractional, for instance, to the primary river, and one can also assign probability distributions within range on these values. Feature considerations are segment crossings if segment collisions are to be avoided. For instance, segment to segment intersection testing is iterated in the procedure of branch start and end nodes. Another feature selection is branch segment to segment proximity. Testing proximity range involving, point to nearest line on point testing iteratively. Any tests here done before procedurally iterating the random walk of each branch, and subsequently repeating random selection of end node as required.
Branching habit is another matter, and one might find that the purely random walk given higher spatial density means likely that branched streams will inevitably collide. Of course, this gets into a bit more complication resolving anti collision methods on the iterative walk. In this case, I’ve used K-D Trees for more rapidly hashing nearest neighbor node sets. Important features worth consideration in such methodology are freedom of movement and avoiding crossings (if required). Bounding boxes and point in polygon tests can be one anti collider solution method here. Another possibility utilizes nearest neighbor nodes and computing the neighboring stream line for a given segment. A given walk choice selection could be tested for bezier line segment intersection. Iteratively rotating the randomly assigned walk vector direction is primarily weighted in decision, and then secondarily magnitude length of the direction vector is incrementally decreased, re iterating checks on step rotation angle. If using curve approximations to stream and river, one would need to employ curve to curve intersection testing. Distance testing is also another important feature relative to range of spacing between branched streams. In other cases, branch crossings may be desired. For instance, simulating river delta systems which in reality is a little less cumbersome problem.
Constructing a point site sample set of river and branch set allows us to generatively construct, for instance, random site points orthogonally relative sample set with, perhaps, some rotational variance. Noting that important features of scale distance length and distribution of points renders changes to visual articulation here. Measuring and controlling distribution of points relative to distance scale and randomly assigned is performed in a way differentiating say flood plain and channel.
Once having constructed a site sample point set, computing Voronoi diagram, generating individual height by spatial metrics described above, I’ve used a formulation which considers both local peak scale height and a base height by spatial factor.
let nh = getProfileheight(heightprofile,(p3d/scaleD));
var hf1 = lowaddhf*((1-min([1,nh]))+1)/2;
Here, lowaddhf is a area and position computed to height interpolated value ranging from [0,1] which is multiplied against the point relative peak in polygon height interpolated value nh relative to addition of lowaddhf averaged computation.
Texturing the heightmap can be additionally preformed using algorithmic overlays, for instance, using common noising features. In this case, I’ve used two distinct noising features (at different scales). For this, noising is set as the primary overlay value, and the computed voronoi field value (above), set as the secondary set, and then setting this return value as the primary for another overlay application relative finer textural noising features. The Heightmap shown above computed with 20,000 sites.
If further interested, code can be found at:
