Introduction

Look at a mountain from a plane, then from a hilltop, then at a rock in your hand: the jagged outline barely changes. Benoit Mandelbrot made that observation the heart of The Fractal Geometry of Nature (1982): natural shapes are statistically self-similar, and their roughness can be measured by a single number, a fractional dimension.

The same year, Alain Fournier, Don Fussell and Loren Carpenter published Computer rendering of stochastic models in Communications of the ACM. They showed how to generate such surfaces cheaply by recursive random midpoint displacement. The best-known descendant of that idea is the diamond-square algorithm, used for decades in games, films and map generators.

The recipe on a (2k+1)×(2k+1)(2^k+1)\times(2^k+1) grid of heights:

  1. Give the four corners random heights.
  2. Square step: set the centre of every square to the average of its four corners, plus a random offset.
  3. Diamond step: set the midpoint of every edge to the average of its (up to four) neighbours, plus a random offset.
  4. Shrink the offsets by a factor 2−H2^{-H}, halve the step, and repeat until every cell is filled.
step = N - 1; amp = 1
while step > 1:
    square step  (centres)   : avg(4 corners) + amp * rand(-1, 1)
    diamond step (edge mids) : avg(neighbours) + amp * rand(-1, 1)
    amp *= 2^(-H);  step /= 2

The parameter H∈(0,1]H\in(0,1] is the roughness exponent: small HH gives spiky, chaotic terrain, HH near 1 gives smooth rolling hills. A competing and more flexible approach, fractional Brownian motion (fBm), sums several layers of smooth noise, called octaves, each twice as detailed and a little weaker than the one before. Both aim at the same statistics; they differ in the artifacts they leave behind.

Grow a Terrain

Pick a generator, then move the sliders. The map is a shaded relief: a light from the north-west makes slopes bright or dark, and everything below the water level is blue. Under it, a thin strip shows the cross-section along the middle row.

<!-- Terrain generator controls and canvases. -->
<div class="modes">
  <button id="btn-ds" type="button" class="mode on">Diamond-square</button>
  <button id="btn-fbm" type="button" class="mode">fBm noise</button>
</div>
<canvas id="map" width="129" height="129" aria-label="Shaded relief map of the generated terrain"></canvas>
<canvas id="profile" width="129" height="40" aria-label="Cross-section of the terrain along the middle row"></canvas>
<div class="stat" id="stat"></div>
<div class="grid">
  <label>Seed <b id="seed-val">7</b>
    <input type="range" id="seed" min="1" max="999" value="7">
  </label>
  <label>Roughness exponent H (low = rough) <b id="h-val">0.80</b>
    <input type="range" id="hurst" min="5" max="100" value="80">
  </label>
  <label id="oct-label">Octaves (fBm) <b id="oct-val">6</b>
    <input type="range" id="octaves" min="1" max="8" value="6">
  </label>
  <label>Water level <b id="water-val">35%</b>
    <input type="range" id="water" min="0" max="100" value="35">
  </label>
</div>
<div class="checks">
  <label class="chk"><input type="checkbox" id="erode"> Thermal erosion</label>
  <label class="chk" id="warp-label"><input type="checkbox" id="warp"> Domain warp (fBm)</label>
  <button id="rand" type="button" class="ghost">Random seed</button>
</div>

/* Light panel layout; the canvases scale to the container width. */
* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; margin: 0; padding: .4rem; background: #f5f7fa; color: #222; }
.modes { display: flex; gap: .4rem; margin-bottom: .5rem; }
button { font: 600 13px system-ui, sans-serif; padding: .4rem .8rem; border: 1px solid #3d6fa8;
         background: #fff; color: #3d6fa8; border-radius: 8px; cursor: pointer; }
button.mode { flex: 1; }
button.on { background: #3d6fa8; color: #fff; }
button.ghost { margin-left: auto; }
canvas { display: block; width: 100%; max-width: 387px; border-radius: 8px; background: #ccd; }
#map { aspect-ratio: 1 / 1; height: auto; }
#profile { height: 70px; margin-top: .35rem; }
.stat { font-size: .82rem; margin: .4rem 0; color: #334; min-height: 2.4em; line-height: 1.35; }
.grid { display: grid; grid-template-columns: repeat(auto-fit, minmax(150px, 1fr)); gap: .45rem .9rem; }
label { display: flex; flex-direction: column; font-size: .82rem; font-weight: 600; color: #334; gap: .15rem; }
input[type=range] { width: 100%; accent-color: #3d6fa8; cursor: pointer; }
.checks { display: flex; flex-wrap: wrap; align-items: center; gap: .5rem 1rem; margin-top: .5rem; }
label.chk { flex-direction: row; align-items: center; gap: .35rem; }
.disabled { opacity: .4; pointer-events: none; }

const i18nFmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
const N = 129;
const $ = (id) => document.getElementById(id);

// Small seeded random generator (mulberry32) so a seed always gives the same terrain.
function rng(seed) {
  let a = seed >>> 0;
  return function () {
    a = (a + 0x6D2B79F5) >>> 0;
    let t = a;
    t = Math.imul(t ^ (t >>> 15), t | 1);
    t ^= t + Math.imul(t ^ (t >>> 7), t | 61);
    return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
  };
}

// Diamond-square on a 129 x 129 grid. Each level halves the step and shrinks the random amplitude by 2^-H.
function diamondSquare(seed, H) {
  const rnd = rng(seed), h = new Float64Array(N * N);
  const rr = () => rnd() * 2 - 1;
  h[0] = rr(); h[N - 1] = rr(); h[(N - 1) * N] = rr(); h[N * N - 1] = rr();
  const k = Math.pow(2, -H);
  let step = N - 1, amp = 1;
  while (step > 1) {
    const half = step / 2;
    // Square step - the centre of each square is the mean of its four corners plus noise.
    for (let y = half; y < N; y += step) {
      for (let x = half; x < N; x += step) {
        const m = (h[(y - half) * N + x - half] + h[(y - half) * N + x + half] +
                   h[(y + half) * N + x - half] + h[(y + half) * N + x + half]) / 4;
        h[y * N + x] = m + rr() * amp;
      }
    }
    // Diamond step - each edge midpoint is the mean of its (up to four) neighbours plus noise.
    for (let y = 0; y < N; y += half) {
      for (let x = ((y / half) % 2 === 0) ? half : 0; x < N; x += step) {
        let s = 0, c = 0;
        if (x - half >= 0) { s += h[y * N + x - half]; c++; }
        if (x + half < N) { s += h[y * N + x + half]; c++; }
        if (y - half >= 0) { s += h[(y - half) * N + x]; c++; }
        if (y + half < N) { s += h[(y + half) * N + x]; c++; }
        h[y * N + x] = s / c + rr() * amp;
      }
    }
    amp *= k;
    step = half;
  }
  return h;
}

// Lattice value noise with a smooth quintic fade, hashed from the seed.
function lat(ix, iy, s) {
  let h = Math.imul(ix, 374761393) ^ Math.imul(iy, 668265263) ^ Math.imul(s, 1103515245);
  h = Math.imul(h ^ (h >>> 13), 1274126177);
  h ^= h >>> 16;
  return (h >>> 0) / 4294967296 * 2 - 1;
}
const fade = (t) => t * t * t * (t * (t * 6 - 15) + 10);
function vnoise(x, y, s) {
  const x0 = Math.floor(x), y0 = Math.floor(y), u = fade(x - x0), v = fade(y - y0);
  const a = lat(x0, y0, s), b = lat(x0 + 1, y0, s), c = lat(x0, y0 + 1, s), d = lat(x0 + 1, y0 + 1, s);
  return (a + (b - a) * u) * (1 - v) + (c + (d - c) * u) * v;
}

// Fractional Brownian motion - sum octaves; frequency doubles (lacunarity 2), amplitude is multiplied by the gain 2^-H.
function fbm(seed, H, oct, warp) {
  const h = new Float64Array(N * N), gain = Math.pow(2, -H), base = 4;
  for (let y = 0; y < N; y++) {
    for (let x = 0; x < N; x++) {
      let px = x / (N - 1) * base, py = y / (N - 1) * base;
      if (warp) {
        // Domain warping - shift the sampling position by another noise field.
        const qx = vnoise(px * 0.75 + 5.2, py * 0.75 + 1.3, seed + 101);
        const qy = vnoise(px * 0.75 + 8.7, py * 0.75 + 3.1, seed + 202);
        px += 1.6 * qx; py += 1.6 * qy;
      }
      let sum = 0, amp = 1, f = 1;
      for (let o = 0; o < oct; o++) {
        sum += amp * vnoise(px * f + o * 17.3, py * f + o * 31.7, seed * 7 + o);
        amp *= gain; f *= 2;
      }
      h[y * N + x] = sum;
    }
  }
  return h;
}

// Thermal erosion - slopes steeper than the talus angle shed material to lower neighbours; total mass is conserved.
function erodeThermal(src, iters) {
  let lo = Infinity, hi = -Infinity;
  for (const v of src) { if (v < lo) lo = v; if (v > hi) hi = v; }
  const T = (hi - lo) * 0.012;
  let g = Float64Array.from(src);
  const dx = [-1, 0, 1, -1, 1, -1, 0, 1], dy = [-1, -1, -1, 0, 0, 1, 1, 1];
  for (let it = 0; it < iters; it++) {
    const delta = new Float64Array(N * N);
    for (let y = 0; y < N; y++) {
      for (let x = 0; x < N; x++) {
        const i = y * N + x;
        let dmax = 0, total = 0;
        for (let k = 0; k < 8; k++) {
          const nx = x + dx[k], ny = y + dy[k];
          if (nx < 0 || ny < 0 || nx >= N || ny >= N) continue;
          const d = g[i] - g[ny * N + nx];
          if (d > T) { total += d; if (d > dmax) dmax = d; }
        }
        if (total === 0) continue;
        const move = 0.5 * (dmax - T);
        for (let k = 0; k < 8; k++) {
          const nx = x + dx[k], ny = y + dy[k];
          if (nx < 0 || ny < 0 || nx >= N || ny >= N) continue;
          const d = g[i] - g[ny * N + nx];
          if (d > T) { delta[ny * N + nx] += move * d / total; delta[i] -= move * d / total; }
        }
      }
    }
    for (let i = 0; i < N * N; i++) g[i] += delta[i];
  }
  return g;
}

// Rescale heights to the range 0 to 1.
function normalize(src) {
  let lo = Infinity, hi = -Infinity;
  for (const v of src) { if (v < lo) lo = v; if (v > hi) hi = v; }
  const out = new Float64Array(N * N), r = (hi - lo) || 1;
  for (let i = 0; i < N * N; i++) out[i] = (src[i] - lo) / r;
  return out;
}

// Estimate H from the slope of log(mean height difference) against log(distance).
function measureHurst(h) {
  const ds = [1, 2, 4, 8, 16], xs = [], ys = [];
  for (const d of ds) {
    let s = 0, n = 0;
    for (let y = 0; y < N; y++) for (let x = 0; x + d < N; x++) {
      s += Math.abs(h[y * N + x + d] - h[y * N + x]);
      s += Math.abs(h[x * N + y] - h[(x + d) * N + y]);
      n += 2;
    }
    xs.push(Math.log(d)); ys.push(Math.log(s / n + 1e-12));
  }
  const mx = xs.reduce((a, b) => a + b) / xs.length, my = ys.reduce((a, b) => a + b) / ys.length;
  let num = 0, den = 0;
  for (let i = 0; i < xs.length; i++) { num += (xs[i] - mx) * (ys[i] - my); den += (xs[i] - mx) ** 2; }
  return num / den;
}

// Interface and rendering.
const cv = $('map'), ctx = cv.getContext('2d');
const pv = $('profile'), pctx = pv.getContext('2d');
let mode = 'ds', raw = null, rawKey = '', hmap = null;

function lerpC(a, b, t) { return [a[0] + (b[0] - a[0]) * t, a[1] + (b[1] - a[1]) * t, a[2] + (b[2] - a[2]) * t]; }
// Colour ramp for land, from sand to snow.
const RAMP = [[0, [214, 200, 150]], [0.06, [96, 150, 70]], [0.35, [62, 118, 56]], [0.6, [128, 110, 82]], [0.8, [150, 146, 142]], [1, [250, 250, 252]]];
function landColor(t) {
  for (let i = 1; i < RAMP.length; i++) {
    if (t <= RAMP[i][0]) return lerpC(RAMP[i - 1][1], RAMP[i][1], (t - RAMP[i - 1][0]) / (RAMP[i][0] - RAMP[i - 1][0]));
  }
  return RAMP[RAMP.length - 1][1];
}

function render() {
  const wl = +$('water').value / 100;
  const img = ctx.createImageData(N, N), Z = 38;
  const L = [-1, -1, 1.2], ll = Math.hypot(L[0], L[1], L[2]);
  for (let y = 0; y < N; y++) {
    for (let x = 0; x < N; x++) {
      const i = y * N + x, e = hmap[i];
      const xl = Math.max(x - 1, 0), xr = Math.min(x + 1, N - 1), yu = Math.max(y - 1, 0), yd = Math.min(y + 1, N - 1);
      const gx = (hmap[y * N + xr] - hmap[y * N + xl]) / (xr - xl) * Z;
      const gy = (hmap[yd * N + x] - hmap[yu * N + x]) / (yd - yu) * Z;
      const nl = Math.hypot(gx, gy, 1);
      const s = Math.max(0, (-gx * L[0] - gy * L[1] + L[2]) / (nl * ll));
      let c;
      if (e < wl) {
        const depth = (wl - e) / (wl || 1);
        c = lerpC([96, 170, 210], [20, 60, 120], Math.min(1, depth * 1.3));
      } else {
        const t = (e - wl) / ((1 - wl) || 1);
        const f = 0.3 + 1.05 * s;
        c = landColor(t).map((v) => v * f);
      }
      img.data[i * 4] = c[0]; img.data[i * 4 + 1] = c[1]; img.data[i * 4 + 2] = c[2]; img.data[i * 4 + 3] = 255;
    }
  }
  ctx.putImageData(img, 0, 0);
  // Draw the cross-section along the middle row.
  pctx.clearRect(0, 0, pv.width, pv.height);
  pctx.fillStyle = '#dfe6ee'; pctx.fillRect(0, 0, pv.width, pv.height);
  const yw = pv.height - wl * (pv.height - 4) - 2;
  pctx.fillStyle = '#5c9fd0'; pctx.fillRect(0, yw, pv.width, pv.height - yw);
  pctx.beginPath(); pctx.moveTo(0, pv.height);
  const row = 64;
  for (let x = 0; x < N; x++) pctx.lineTo(x, pv.height - hmap[row * N + x] * (pv.height - 4) - 2);
  pctx.lineTo(N - 1, pv.height); pctx.closePath();
  pctx.fillStyle = '#6b6a52'; pctx.fill();
  $('water-val').textContent = Math.round(wl * 100) + '%';
}

function rebuild() {
  const seed = +$('seed').value, H = +$('hurst').value / 100, oct = +$('octaves').value, warp = $('warp').checked;
  const key = [mode, seed, H, oct, warp].join('|');
  if (key !== rawKey) {
    raw = mode === 'ds' ? diamondSquare(seed, H) : fbm(seed, H, oct, warp);
    rawKey = key;
  }
  let g = raw;
  if ($('erode').checked) g = erodeThermal(raw, 40);
  hmap = normalize(g);
  const Hm = measureHurst(hmap);
  $('seed-val').textContent = seed;
  $('h-val').textContent = H.toFixed(2);
  $('oct-val').textContent = oct;
  $('stat').textContent = i18nFmt("Measured Hurst exponent H ≈ %h% · surface fractal dimension D = 3 − H ≈ %d%", { h: Hm.toFixed(2), d: (3 - Hm).toFixed(2) });
  render();
}

function setMode(m) {
  mode = m;
  $('btn-ds').classList.toggle('on', m === 'ds');
  $('btn-fbm').classList.toggle('on', m === 'fbm');
  $('oct-label').classList.toggle('disabled', m === 'ds');
  $('warp-label').classList.toggle('disabled', m === 'ds');
  rebuild();
}

$('btn-ds').addEventListener('click', () => setMode('ds'));
$('btn-fbm').addEventListener('click', () => setMode('fbm'));
['seed', 'hurst', 'octaves', 'erode', 'warp'].forEach((id) => $(id).addEventListener('input', rebuild));
$('water').addEventListener('input', render);
$('rand').addEventListener('click', () => { $('seed').value = 1 + Math.floor(Math.random() * 999); rebuild(); });
setMode('ds');
Result
Console

Things to try:

  • Drag H from 1.00 down to 0.10 and watch the terrain go from rolling hills to jagged chaos. The panel estimates the exponent from the generated heights using a log-log slope, so you can check that the knob does what the theory says.
  • Switch to fBm noise and lower the octaves to 1: only the biggest bumps remain. Add octaves one at a time to see each layer of detail appear.
  • Look at the diamond-square map near the grid lines. Faint straight creases along the axes are the algorithm's known weakness, discussed below. Turn on domain warp in fBm mode to see how bending the sampling coordinates breaks up the regular look.
  • Tick thermal erosion: steep slopes slump into their lower neighbours (a simple talus model with total mass conserved), which softens the peaks and gives scree-like ramps. Note how the measured HH goes up, because erosion removes small-scale roughness first.

The Real Complexity

Diamond-square. Let n=(2k+1)2n=(2^k+1)^2 be the number of grid points. Every point is written exactly once, from a constant number of already known neighbours, so the time is Θ(n)\Theta(n) in the best, average and worst case (the algorithm is data-independent) and the space is Θ(n)\Theta(n) for the heightmap, computed in place with O(1)O(1) extra memory. The number of levels is k=log⁡2(n−1)k=\log_2(\sqrt{n}-1).

fBm with mm octaves. Each point sums mm noise evaluations, each O(1)O(1), so the time is Θ(m n)\Theta(m\,n) and the space is Θ(n)\Theta(n) (or O(1)O(1) if you evaluate pixels on demand, which is why noise suits GPUs and infinite worlds). Since octaves double in frequency, beyond log⁡2\log_2 of the grid size they fall below one cell and add nothing: m≤log⁡2nm\le\log_2 n is the useful range.

Thermal erosion with II iterations over a grid of nn cells costs Θ(I n)\Theta(I\,n) with 8 neighbour checks per cell.

What the exponent means. A fractional Brownian motion BHB_H with Hurst exponent H∈(0,1)H\in(0,1), introduced by Mandelbrot and Van Ness in 1968, satisfies

E[(BH(t+τ)−BH(t))2]=σ2 ∣τ∣2H.\mathbb{E}\big[(B_H(t+\tau)-B_H(t))^2\big]=\sigma^2\,|\tau|^{2H}.

  • H=12H=\tfrac12 is ordinary Brownian motion (a random walk). H<12H<\tfrac12 is anti-persistent and jagged, H>12H>\tfrac12 is persistent and smooth.
  • A profile (a slice) has fractal dimension D=2−HD=2-H and a surface has D=3−HD=3-H, between 2 (flat) and 3 (space-filling). That is the number the demo prints.
  • The power spectrum of the surface falls off as S(k)∝∣k∣−(2H+2)S(\mathbf{k})\propto|\mathbf{k}|^{-(2H+2)}. In fBm terrain generation, adding an octave multiplies the frequency by the lacunarity ℓ=2\ell=2 and the amplitude by the gain g=ℓ−H=2−Hg=\ell^{-H}=2^{-H}, which reproduces exactly this spectrum. Diamond-square uses the same factor per level, so the same HH describes both.

Honest limits. Diamond-square is only an approximation of fBm. Points created at coarse levels keep their value and are never perturbed again, so the increments are not stationary: the statistics depend on where you are on the grid. That produces visible axis-aligned crease artifacts, which were pointed out in the computer-graphics literature in the 1980s (Gavin Miller's 1986 SIGGRAPH paper on terrain rendering is the usual reference). fBm from lattice noise avoids the grid-aligned creases, at the price of a constant factor more work. Neither model produces real geology: rivers, valleys and drainage need erosion. Thermal erosion is the cheap one; hydraulic erosion, which simulates rain, flow and sediment transport, is what makes terrain look carved (Musgrave, Kolb and Mace, SIGGRAPH 1989). In this demo the measured HH also drifts from the slider at high HH because the fit is made over scales of 1 to 16 cells on a 129-cell map.

Where It Matters

  • Games and film. Procedural worlds store a seed instead of gigabytes of geometry. Heightmaps from fBm, domain warping and erosion are a standard starting point for terrain, which artists then edit by hand.
  • Clouds, textures and coastlines. The same noise with a different mapping gives marble, fire, clouds and ragged shorelines. The classic smooth-noise building block is Perlin noise.
  • Mathematics. The roughness exponent connects to the idea of a fractional dimension, the language of the Mandelbrot set and of self-similar shapes in general; and every generator here depends on a good seeded random number generator.
  • Turning heights into solids. A heightmap is a 2.5D surface; to get caves and overhangs you move to 3D density fields and extract surfaces with marching cubes.
  • Science. Real topography, financial series and river levels show Hurst-like scaling, so fBm is used to model and to stress-test them. Estimating HH from data (as the demo does with a variogram slope) is a standard first step.

Conclusion

Diamond-square and fBm show how a tiny amount of structure, a random offset that shrinks by 2−H2^{-H} at every scale, can mimic what nature does over millions of years. Both run in time linear in the number of points (times the octave count for fBm), and both are controlled by one meaningful number, the Hurst exponent, which sets the fractal dimension D=3−HD=3-H of the surface.

They are also models, not geology: diamond-square leaves grid creases, and credible valleys need erosion. Knowing where the model stops being true is part of using it well. Next, try Perlin noise for the smooth building block, or the Mandelbrot set for a different kind of self-similarity.

Share this article

Pick a channel — or use your device's native share sheet.

Comments

Loading comments...

https://www.kipuhub.com/en/article/diamond-square-terrain/Content licensed under CC BY-NC 4.0.