Hydraulic & Thermal Erosion on Infinite Procedural Terrain

Overview

The goal of this project was to create a method for generating "infinite" procedural eroded terrain which would I would ultimately use for generating procedural planets. The most common method of generating infinite procedural terrain is to combine multiple octaves of noise generated by algorithms that take a 2d or 3d coordinate as an input and generate a value that can be interpreted as elevation in the case of heightmap terrain. This approach is popular becasue its easily parallelizable and thus fast and only requires the seed used to initialize the pseudorandom noise functions as an input. While various noise functions and other math can be composed in a much more complex function that is capable of generate interesting looking terrain features it is difficult to reproduce terrain features created by hydraulic, coastal and thermal erosion in a realistic, coherent manner; even attaining decent results requires an extensive amount of experimentation and tweaking of the procedural terrain function.

While the planetary terrain shown in the screenshots above may look interesting and somewhat resemble realistic terrain features, it is too repetitive and unrealistic to pass for real terrain; features such as tectonic uplift, water drainage channels, angle-of-repose slopes, etc. are missing and are difficult to reproduce without actually simulating plate tectonics and erosion. Although simulating hydraulic and thermal erosion on a heightmap terrain is fairly simple, doing so on an "infinite" or at least planetary-scale terrain seems intractable due to the fact that patches of terrain will be influenced by distant patches which may require a global simulation of the entire terrain at full resolution. One thought that might immediately come to mind is that simply limited the range at which remote patches of terrain could influence each other could allow one to generate an infinite terrain in independent patches as long as a sufficiently large "buffer" or "border" area is included in the erosion simulation - and it turns out this works and is exactly what I implemented. However, I first investigated alternative methods of generating eroded terrain including an interesting "erosion noise" function, a river network algorithm and ML methods such as CNNs or GANs to either "erode" an input terrain or directly generate an eroded terrain, respectively, and ultimately found that these methods would be some combination of not sufficiently realistic, too computationally expensive for real-time terrain generation or not clearly applicable to infinite terrain.

Method

The algorithm I settled on is fairly simple and takes a base terrain heightmap as input and runs a very simplistic simulation of hydraulic and thermal erosion on the heightmap, iteratively on the GPU. As stated previously, the only conceptual trick here is that as long as the base heightmap includes a border area of terrain present in adjacent patches of terrain, then the interior, usable part of the output heightmap will be consistent with neighboring patches of terrain - the exact size of the border depends on the maximum range of erosion, if the border is too small then seams and other artifacts may appear at the boundaries of terrain patches.

While this approach sounds limiting as we want to ideally keep the erosion range as small as possible to keep the unused parts of these heightmaps as small as possible, thus causing the erosion features to be limited in scale; by leveraging the fact that most terrain systems use a quadtree-LOD we can generate mutli-scale erosion features for free by simply reusing a patch's parent heightmap as the input for the next depth level of terrain patches - this is necessary for use on planetary-scale terrain (this is not implemented in the Unity project).

While the algorithm is simple in concept, the implementation was somewhat tricky as doing this on the GPU has to be done carefully to avoid race conditions and other issues that would lead to non-deterministic results. I adapted Sebastian Lague's code here for a simple droplet-based hydraulic erosion simulation on the GPU and noticed that it yields non-deterministic results as it simulates each droplet's entire lifetime, as it makes changes to the heightmap, in parallel - clearly this results in a race condition as droplets literally race across and carve the terrain in a highly non-deterministic way. Alleviating this race condition requires modifying the algorithm in a few ways: we must track the droplets across multiple compute disptaches, we have to update the droplets simultaneously at each step, we have to defer making changes to the heightmap until every droplet has been updated and we have to make sure that we use atomic operations to accumulate the changes that will be made to heightmap after every step. The final algorithm which includes a simple thermal erosion computed after the hydralic erosion step looks something like this (all steps are done in compute shaders):

1. generate procedural input terrain patch via noise algorithms
2. generate droplets for terrain patch
3. hydralic erosion for several iterations:
   1. update droplet position, velocity, sediment amount, etc., accumulating height delta to an integer buffer (necessary to use interlocked ops to avoid race conditions)
   2. copy accumulated height delta to heightmap
4. thermal erosion for several iterations
   1. transfer some fraction of the height difference between the neighboring pixel with the largest height differece, accumulating height delta same as above
   2. copy accumulated height delta to heightmap

I initially implemented this algorithm in Unity and later in the Nebula engine for use on procedural planets and ran into a bug/oversight in the Vulkan implementation of both engines that was not an in issue in Unity's DX12 mode. Since the droplets compute buffer is accessed and modified repeatedly in the hydraulic erosion compute kernel it is important to make sure that its memory property flag is set to VK_MEMORY_PROPERTY_DEVICE_LOCAL_BIT and not VK_MEMORY_PROPERTY_HOST_COHERENT_BIT as per the note on this page of Khronos Vulkan reference otherwise the shader will be unusably slow. I fixed this issue in my engine but at the time of this writing using this shader in Unity's Vulkan mode will be extremely slow.