Water and Terrain

A flat plane with a blue texture won't pass for an ocean, and a mountain range can't be modeled triangle by triangle. This chapter builds both from height fields: a height for every point on a horizontal plane. For water, the heights come from physics and statistics, moving waves with sizes measured at sea, summed by a fast Fourier transform. For land, they come from layered noise (Chapter 13). And for both, the hard part is drawing a surface that reaches the horizon with only as many triangles as the screen needs.

The simplest wave is a sine moving sideways: height A·sin(k·x − ω·t), with amplitude A, wavenumber k = 2π / wavelength and angular speed ω. Add a few with different directions, lengths and speeds and you get a passable pond. But real waves aren't sine-shaped: their crests are sharp and their troughs are broad.

The reason is that water doesn't just bob up and down. As a wave passes, each bit of water moves in a circle: forward at the crest, backward in the trough. So the water bunches up at the crests and spreads out in the troughs. The Gerstner wave (Franz Gerstner, 1802) models exactly this: move each point of the surface sideways as well as up, by a circle of radius A. The ratio Q = k·A is the wave's steepness; at Q = 1 the crest becomes a sharp point, and beyond it the surface loops over itself, which is where real waves break.

Water waves also obey a simple law, the dispersion relation for deep water: ω = √(g·k), where g is gravity. Long waves travel faster than short ones; a 100-meter swell moves at 45 km/h, a 1-meter ripple at 4.5 km/h. That's why a storm far away arrives as long, regular swell: the long waves outrun the short ones. A water renderer that gets this right looks right even to people who have never heard of it.

In code, one Gerstner wave is a few lines. It has a direction d, a wavelength and a steepness Q; the amplitude follows from the steepness, A = Q / k, and the speed from the dispersion relation. Then the surface point that started at xz moves around its circle:

let k = 2.0 * PI / wavelength;
let omega = sqrt(9.81 * k);
let amp = steepness / k;
let theta = k * dot(d, xz) - omega * t;
p += vec3f(d.x * amp * cos(theta), amp * sin(theta), d.y * amp * cos(theta));

A sea is a sum of these: start with p at the point's resting place and add every wave's circle. Where several waves' crests line up, their steepnesses add, so keep the total under 1 or those crests loop over.

The playground draws a sea of five Gerstner waves. The grid it displaces is shaped like the view: its rows bunch up near the camera and spread out toward the horizon, and each row is only as wide as the view at that distance. The normal comes from the derivatives of the displacement, how the surface moves when xz moves a little in x and in z. Those same derivatives say how squeezed the surface is, which paints foam on the sharpest crests; Part 3 comes back to that. Change the WAVES, or set CHOP to 0 for plain sine waves and to 1.5 to watch the crests fold.

A dozen Gerstner waves make a nice lake, but the open sea is a jumble of waves of every length and direction. Oceanographers have measured how much energy the sea holds at each wavelength for a given wind: a wave spectrum. Wind that blows long and steadily builds waves up to a peak wavelength that grows with the square of the wind speed; shorter waves carry less energy; waves running across or against the wind are weaker. The demo uses the Phillips spectrum, a simple shape that captures this; more measured ones (Pierson–Moskowitz, JONSWAP) are common in film and games.

The recipe, from Jerry Tessendorf's notes Simulating Ocean Water (2001), which nearly every ocean renderer since has followed:

1. Make a grid of wave vectors k: one wave per grid cell, 128 × 128 of them, with directions all around and lengths from the patch size down to two grid cells.
2. Give each a random amplitude and phase, with the amplitude's average size set by the spectrum. This happens once.
3. Every frame, advance each wave's phase by ω·t, with ω from the dispersion relation.
4. Add up all 16,384 waves at each of the 16,384 grid points.

Step 4 is the expensive one: 16,384 waves times 16,384 points is 268 million sine evaluations, per frame, for each quantity (height, and the two sideways displacements). The fast Fourier transform does the same sum in about N² log N steps instead of N⁴. It works by splitting the sum into its even and odd terms, two half-size sums that share most of their work, and splitting those again, down to pairs. Each combining step, a butterfly, takes two values and produces their sum and difference, one of them turned by a rotation.

On the GPU, a 2D FFT is a 1D FFT along every row, then along every column. The demo gives each row to one workgroup: its 64 threads load the 128 values into workgroup memory (Chapter 18), then run the seven stages, each thread doing one butterfly per stage, with a barrier between stages. A second dispatch does the columns. Three compute dispatches in all per frame: advance the spectrum, transform the rows, transform the columns. Heights and sideways displacements ride in the same pass, packed into one vec4f per point.

Here are the stages, after the row has been loaded into the workgroup array line in bit-reversed order. Thread t does butterfly number t of each stage, on two complex numbers stored as vec2f:

for (var s = 0u; s < 7u; s++) {           // 128 values, 7 stages
  let half = 1u << s;                     // pairs are 1, 2, 4, ... apart
  let j = t % half;
  let i0 = (t / half) * half * 2u + j;
  let i1 = i0 + half;
  let angle = PI * f32(j) / f32(half);
  let w = vec2f(cos(angle), sin(angle));  // the rotation
  let a = line[i0];
  let b = cmul(line[i1], w);              // complex multiply
  line[i0] = a + b;
  line[i1] = a - b;
  workgroupBarrier();                     // wait for the whole stage
}

// the helpers used above
// multiply two complex numbers stored as (real, imaginary)
fn cmul(a: vec2f, b: vec2f) -> vec2f {
  return vec2f(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
}

The barrier is what makes this work: no thread may start a stage until every thread has finished the one before, because each butterfly reads values that two other butterflies just wrote.

Choppy waves. The FFT gives heights, but that's the "sine" picture again, with round crests. The Gerstner idea carries over: compute a sideways displacement for each wave, pointing along its direction, and transform it too. Moving each point sideways by that much, times a choppiness factor, sharpens every crest at once.

Foam. Push choppiness too far and the surface folds over itself, exactly where real waves break. That's measurable: the Jacobian of the sideways displacement says how much a small patch of surface is stretched or squeezed. 1 means unchanged, smaller means squeezed together, below zero means folded inside out. Painting foam where it drops below a threshold puts whitecaps on the steepest crests, and keeping it for a few seconds afterward leaves the trails of foam seen behind breaking waves.

The shader already has what it needs. To find the normal, it measures how the displaced surface moves when the resting point moves a little along x and along z (one grid cell, e); the Jacobian is the area of the little parallelogram those two steps span, seen from above:

let du = (surface(xz + vec2f(e, 0.0)) - surface(xz - vec2f(e, 0.0))) / (2.0 * e);
let dv = (surface(xz + vec2f(0.0, e)) - surface(xz - vec2f(0.0, e))) / (2.0 * e);
let normal = normalize(cross(dv, du));
let jacobian = du.x * dv.z - du.z * dv.x;
let foam = smoothstep(0.8, 0.25, jacobian);

// the helpers used above
// The displaced point for rest position xz.
fn surface(xz: vec2f) -> vec3f {
  let s = sampleField(xz);
  let chop = u.params.x;                  // the choppiness setting
  return vec3f(xz.x - chop * s.y, s.x, xz.y - chop * s.z);
}
// The FFT result at xz: (height, dx, dz), blended between the four
// nearest grid points. The patch repeats every PATCH meters.
fn sampleField(xz: vec2f) -> vec3f {
  let g = xz / PATCH * f32(N) - 0.5;      // position in grid cells
  let i = vec2i(floor(g));
  let f = g - floor(g);
  let top = mix(texel(i), texel(i + vec2i(1, 0)), f.x);
  let bottom = mix(texel(i + vec2i(0, 1)), texel(i + vec2i(1, 1)), f.x);
  let v = mix(top, bottom, f.y);
  return vec3f(v.x, v.z, v.w);            // height, dx, dz
}
// One grid point, wrapping around at the patch's edges (N is a power of 2).
fn texel(i: vec2i) -> vec4f {
  let w = vec2u(i & vec2i(i32(N) - 1));
  return field[w.y * N + w.x];
}

Color. Water reflects and refracts (Chapter 26). At grazing angles the Fresnel term makes it a mirror of the sky; looking down, you see into it, and deep water returns a dark blue-green, the little light that scattered back up before being absorbed. Thin crests lit from behind glow green-blue, a cheap subsurface term. The sun's reflection on millions of tiny facets is the glitter path, a specular highlight on a very rough, very bumpy surface. And a shoreline adds depth-based color, refraction of the bottom, and caustics, which this chapter leaves out.

Tiling and distance. One patch, repeated, shows its pattern from above (try show patches). Real renderers add up several patches of different sizes, a cascade (a few kilometers down to a few meters), whose repeats don't line up. Far away, a single pixel covers many waves; sampling the surface there flickers, so detail is filtered away into an average roughness, as with mipmaps (Chapter 10) and roughness from normal maps (Chapter 23).

Drawing the sea is a level-of-detail problem too: the demo uses a grid of rings centered on the camera, dense at your feet and coarse at the horizon, whose vertices look up their heights from the FFT result. That's the ocean's version of the terrain techniques next.

A terrain is a height field too, stored as a heightmap image (from satellite surveys, an artist's sculpting, or an erosion simulation) or computed from noise. Chapter 13's fractal noise gives rolling hills. Mountains need more: ridged noise, 1 − |noise|, makes sharp crests; and damping each octave where the earlier ones are already steep leaves smooth valleys and rugged peaks, a quick imitation of erosion. The demo below, and the playground at the end, use these tricks.

Then the drawing problem. A 10-kilometer landscape at 2-meter resolution is 25 million vertices, but the screen only needs fine triangles near the camera; a triangle a few kilometers away should be hundreds of meters wide. The classic answer is a quadtree: start with a few big square patches; split any patch that's close to the camera into four, and split those again, until each patch is small enough for its distance. Every patch, whatever its size, is drawn with the same small grid mesh, say 16 × 16 cells, by instancing (Chapter 21): big patches have big cells.

Two problems remain, both at the boundary between a fine patch and a coarse one. Cracks: the fine patch has a vertex halfway along the shared edge, the coarse one doesn't, and if that vertex's height isn't exactly halfway between its neighbors, a gap opens. Popping: when a patch splits as the camera approaches, its shape jumps visibly.

CDLOD (continuous distance-dependent level of detail, Filip Strugar, 2009) fixes both with one trick: morphing. As a vertex approaches the distance where its patch would be replaced by a coarser one, every second vertex slides smoothly onto its even neighbor, so that by the boundary, the fine grid has become exactly the coarse grid. No cracks, because the seams match; no pops, because the change happens a little at a time. It all happens in the vertex shader, from the vertex's distance to the camera.

The other well-known scheme is the geometry clipmap (Frank Losasso and Hugues Hoppe, 2004): instead of a tree, a set of nested square rings of grid centered on the camera, each with twice the spacing of the one inside it, which slide along in whole grid steps as the camera moves. The heights live in a matching stack of textures, updated a strip at a time as the rings move. It's simple and very regular, made for streaming huge heightmaps, and it's how the ocean demo's rings work too. Both schemes are often combined with skirts (a strip hanging down from every patch edge, hiding any crack that's left) and GPU culling (Chapter 34) of patches outside the view.

A terrain is too big for one texture at close-up resolution, so it's painted with a few repeating materials (grass, dirt, rock, snow) and a splat map: a low-resolution texture whose channels say how much of each material goes where. Rules fill in what artists don't paint: rock where the slope is steep, snow where it's high and flat. Blending by a per-material height map, instead of by the weights alone, makes the transitions look natural: rock pokes through the snow in its cracks first.

On cliffs, the texture coordinates a heightmap gives (just x and z) stretch into long streaks. Triplanar mapping samples the texture three times, projected along each axis, and blends by the normal, so steep faces get an undistorted projection. Very large worlds go further, baking the blended result into a virtual texture (Chapter 34's streaming idea, for textures) that's only filled near the camera.

Height blending is a few lines in the shader. Each material has a height (its pebbles stand up, its cracks sink down), and each pixel adds that to the material's weight, keeps only the materials within a small depth of the highest, and normalizes:

let h = heights + weights;
let top = max(max(h.x, h.y), max(h.z, h.w)) - 0.2;
var b = max(h - top, vec4f(0.0));
b /= b.x + b.y + b.z + b.w;

Where one material clearly dominates, only it survives; in a transition, whichever is locally higher wins pixel by pixel, so dirt fills the gaps between grass tufts and snow fills the cracks in rock before it covers the slabs. A painted map one texel every 2 m ends up looking like it was painted at a few millimeters.

This fragment shader draws a landscape without any triangles: it marches each pixel's ray across the height field until it dips below the ground, stepping by a fraction of the ray's current height above the terrain (a cheap, slightly unsafe version of Chapter 15's sphere tracing). Change SCALE, the number of OCTAVES, or the WATER level; then edit height itself. The comments suggest ridges and terraces.

• Water and terrain are height fields, sums of waves of many sizes.
• In a Gerstner wave, water moves in circles, sharpening crests; deep-water waves obey ω = √(g·k), so long waves are faster.
• An ocean is a spectrum of random waves with measured energy per wavelength, summed every frame by an FFT: rows, then columns, in workgroup memory.
• Choppy sideways displacement sharpens crests; its Jacobian finds folds, where foam goes. Fresnel, deep-water scattering and sun glitter give the color; cascades of patch sizes hide tiling.
• Terrain level of detail draws the same small grid at different sizes, chosen by a quadtree; CDLOD morphing slides vertices onto the coarser grid by distance, avoiding cracks and pops. Geometry clipmaps use nested rings instead.
• Splat maps, slope and height rules, and triplanar mapping texture terrain.

Chapter 38 ends the book with the newest corner of graphics: images and scenes learned from data by neural networks and by optimization.

1. Folding crests. In the Gerstner playground, raise CHOP a little at a time. Around which value do the first crests start to fold, even though no single wave is steeper than 0.33?

   Hint

   Where all five crests happen to line up, their steepnesses add: 0.25 + 0.2 + 0.15 + 0.1 + 0.08 = 0.78. Times CHOP, that reaches 1 at about 1.3, and the rarer the lineup, the rarer the folds. Real ocean renderers scale choppiness down as more waves are added for the same reason.

2. Ridges. In the playground, replace noise(q) in the loop with the ridged version in the comment. What happens to the mountains? Then try it for the first octave only.

   Hint

   Ridged noise turns every valley of the noise into a sharp crest, so the mountains get knife-edge ridges. Using it only for the large octaves gives sharp ranges with smooth detail on top: select(noise(q), 1.0 - abs(noise(q) * 2.0 - 1.0), i < 2).

3. Swell speed. A swell with a wavelength of 200 meters: how fast does it travel? How long is its period, the time between crests passing a buoy?

   Hint

   k = 2π / 200 ≈ 0.0314, ω = √(9.81 × 0.0314) ≈ 0.555 per second. Speed is ω / k ≈ 17.7 m/s (64 km/h); the period is 2π / ω ≈ 11.3 seconds.

4. Why one level apart? In the quadtree figure, neighbors never differ by more than one level. Why would a fine patch next to a patch two levels coarser be a problem for morphing?

   Hint

   Morphing turns a patch into the next coarser grid at its far edge. Next to a patch two levels coarser, the shared edge would need to match a grid four times coarser, and the vertices in between would still be free to make cracks. Choosing each level's range large enough compared to its patch size guarantees the one-level rule.