Iterative Thermal Erosion in GLSL Heightmaps

Thermal erosion simulates the natural degradation of steep terrain slopes into stable angles of repose, known as the talus angle. In OpenGL Shading Language (GLSL), this process is evaluated iteratively over a terrain heightmap by computing height differentials between adjacent texels, calculating sediment displacement beyond a critical threshold, and updating the surface across successive render passes. Using ping-pong framebuffers or compute shaders, the GPU accelerates this relaxation process, resolving millions of texel interactions in parallel.

The Thermal Erosion Model

Thermal erosion, or material slumping, occurs when the slope between a given point and its neighbor exceeds a critical angle of repose (the talus threshold, \(T\)). The height differential is defined as:

\[\Delta h_i = h_{\text{center}} - h_i\]

If \(\Delta h_i > T\), the slope is unstable. A fraction of the excess material, governed by a transfer rate constant \(c\), shears off and is deposited onto the lower neighbor:

\[\Delta m_i = c \cdot (\Delta h_i - T)\]

To maintain numerical stability and avoid over-eroding a cell past the height of its neighbors, the total displaced material must not exceed the maximum allowable drop divided by the number of participating neighbors.

The GPU Gather Architecture

Because GPUs rely on a parallel gather pattern—where a shader thread writes strictly to its designated output texel rather than scattering values to adjacent memory addresses—a cell must calculate both its outgoing material and the incoming material from its neighbors.

In a single-pass formulation, the new height \(h'\) of a cell at coordinate \(\mathbf{x}\) is calculated as:

\[h'(\mathbf{x}) = h(\mathbf{x}) - \sum \text{Outflow}(\mathbf{x} \to \text{neighbors}) + \sum \text{Inflow}(\text{neighbors} \to \mathbf{x})\]

Evaluating both inflow and outflow simultaneously requires sampling not only the direct neighbors of the cell, but also the neighborhood of each neighbor to determine how much material those neighbors distributed elsewhere. To avoid excessive texture lookups, thermal relaxation is often split into two distinct sub-passes:

  1. Sediment/Flux Pass: Calculates and stores the amount of material shifting along grid directions.
  2. Advection/Update Pass: Reads the flux data, subtracts outgoing material, adds incoming material, and writes the updated height.

Alternatively, a simplified symmetric relaxation can evaluate incoming and outgoing sediment in a single pass if the time step (\(\Delta t\)) is kept sufficiently small to prevent oscillation.

Ping-Pong Framebuffer Pipeline

GLSL shaders cannot read from and write to the same texture memory simultaneously without causing race conditions. The iterative relaxation loop resolves this by using two identical floating-point textures (typically GL_R32F or GL_RGBA32F formats) set up in a ping-pong configuration:

  1. Texture A is bound as an input sampler; Texture B is bound as the render target (Framebuffer Object or compute write image).
  2. The GLSL shader executes across the grid, evaluating one erosion step.
  3. The bindings swap: Texture B becomes the input sampler, and Texture A becomes the render target.
  4. The process repeats for \(N\) iterations until the terrain converges to the desired stability.

GLSL Implementation Example

The following fragment shader demonstrates a single-pass relaxation step using a 4-neighborhood (von Neumann) stencil:

#version 330 core

out vec4 FragColor;

uniform sampler2D u_Heightmap;
uniform vec2 u_TexelSize;     // 1.0 / TextureDimensions
uniform float u_CellSize;      // World-space distance between texels
uniform float u_TalusAngle;    // Minimum height difference to trigger erosion
uniform float u_ErosionRate;   // Rate coefficient (0.0 to 0.5 for stability)

void main()
{
    vec2 uv = gl_FragCoord.xy * u_TexelSize;
    float h = texture(u_Heightmap, uv).r;

    // Neighbor offsets: North, South, East, West
    vec2 offsets[4] = vec2[](
        vec2(0.0, 1.0),
        vec2(0.0, -1.0),
        vec2(1.0, 0.0),
        vec2(-1.0, 0.0)
    );

    float totalOutflow = 0.0;
    float totalInflow = 0.0;
    float maxDiff = 0.0;

    // Pass 1: Calculate material leaving this cell
    for (int i = 0; i < 4; ++i) {
        float h_neighbor = texture(u_Heightmap, uv + offsets[i] * u_TexelSize).r;
        float diff = h - h_neighbor;

        if (diff > u_TalusAngle) {
            float excess = diff - u_TalusAngle;
            totalOutflow += excess;
            maxDiff = max(maxDiff, diff);
        }
    }

    // Clamp outflow to avoid inversion artifacts
    if (totalOutflow > 0.0) {
        totalOutflow = min(totalOutflow * u_ErosionRate, maxDiff * 0.5);
    }

    // Pass 2: Calculate material arriving from neighbors
    for (int i = 0; i < 4; ++i) {
        vec2 n_uv = uv + offsets[i] * u_TexelSize;
        float h_neighbor = texture(u_Heightmap, n_uv).r;
        float diff = h_neighbor - h;

        if (diff > u_TalusAngle) {
            // Estimate neighbor's local outflow distribution
            float neighborExcess = diff - u_TalusAngle;
            totalInflow += (neighborExcess * u_ErosionRate) * 0.25;
        }
    }

    float newHeight = h - totalOutflow + totalInflow;
    FragColor = vec4(newHeight, 0.0, 0.0, 1.0);
}

Numerical Stability and Boundary Conditions

Iterative relaxation methods on the GPU are susceptible to numerical oscillation when transfer rates are set too high. If a cell transfers more material than is necessary to achieve the talus threshold, the slope will invert in the following iteration, causing checkerboard artifacts. Setting the global transfer rate \(c \le \frac{1}{2D}\) (where \(D\) is the dimensionality or neighborhood count) guarantees monotonic convergence.

Boundary texels along the edge of the heightmap must be handled explicitly by clamping texture sampling (GL_CLAMP_TO_EDGE) or mirroring coordinates. When simulating open-boundary systems, borders can act as sinks where sediment falls off the map, whereas closed boundaries preserve total mass across the entire heightmap.