```{index}
single: procedural generation
single: terrain generation
```

# Capstone: a procedural world generator

A program that invents a landscape and draws it. Three stages turn noise
into terrain, and each is written three times: as plain Perl loops, as PDL
whole-array expressions, and as WGSL compute shaders. `--mesh all` runs all
three on every launch and prints how far they disagree, so the comparison is
known to be measuring one algorithm rather than three.

Everything after those stages is the plain Perl from this guide's earlier
chapters: the mesh from [meshes](meshes), the camera and rasteriser from
[software-rasteriser](software-rasteriser), the Lambert term from
[shading](shading), the writer from [model-formats](model-formats).

```text
world   512x512 grid, 6 octaves, 12 erosion passes, seed 1234
gpu     Quadro M2000M (Vulkan)

  backend   terrain   erosion   normals  readback     total  speedup
  perl      1.033s    4.945s    0.221s    0.000s    6.198s     1.0x
  pdl       0.253s    0.304s    0.017s    0.033s    0.607s    10.2x
  gpu       0.002s    0.004s    0.004s    0.034s    0.043s   145.2x

  agree   perl vs pdl   max dh = 0.00e+00   max dn.x = 0.00e+00
  agree   perl vs gpu   max dh = 2.95e-07   max dn.x = 1.26e-05
```

PDL and plain Perl agree to the last bit, because both compute in `double`.
The GPU differs in the seventh decimal, because WGSL storage buffers here
are `f32`. That is the expected difference, and it is worth printing rather
than assuming: a backend that has silently diverged shows up here as a large
number instead of a small one.

```{index}
single: compute shader
```

## Three stages, one shape

| Stage | Input | Output | Per cell |
|---|---|---|---|
| terrain | noise lattice | height | sum `N` octaves of bilinear samples |
| erosion | height | height | compare 8 neighbours |
| normals | height | 3 components | 4 neighbours, one `sqrt` |

Each produces one value per grid cell from a fixed neighbourhood, with no
cell depending on another cell's result within the same pass. That is the
shape a compute shader exists for, the shape PDL expresses without a loop,
and the shape the JIT compiles well. It is why one algorithm can have three
implementations that agree to the bit.

The renderer is deliberately not in that list. Rasterising is a scatter:
triangles write to overlapping pixels and the Z-buffer arbitrates. It stays
in Perl, and at 512 it costs `4.10s` against the field stages' `0.043s`.

```{index}
single: srand; reproducible noise
```

## One lattice, three backends

The noise lattice is built once, in Perl, and handed to whichever backend
runs:

```perl
srand $o{seed};
my @LAT = map { rand() } 0 .. $R * $R - 1;
```

**Invariant: every backend consumes byte-identical input.** A shader that
generated its own noise from an integer hash would produce a different field
from the Perl one, and "the three agree" would be untestable. This also
removes the portability trap in hash-based noise: `sin`-based hashes differ
between drivers, so terrain generated on one GPU would not reproduce on
another.

```{index}
single: value noise
single: fractal Brownian motion
single: ridged noise
single: smoothstep
single: compute shader; bounds check
```

## Terrain: octaves of value noise

Good terrain sums several **octaves**, each at twice the frequency and half
the amplitude of the last, which is **fractal Brownian motion**. Each octave
bilinearly samples the lattice with a smoothstep fade.

```wgsl
struct P { n: u32, r: u32, oct: u32, pad: u32, base: f32, p1: f32, p2: f32, p3: f32 };
@group(0) @binding(0) var<storage, read_write> lat: array<f32>;
@group(0) @binding(1) var<uniform> p: P;
@group(0) @binding(2) var<storage, read_write> out: array<f32>;

// NOT `smooth`: that is a reserved word in WGSL and the shader will not compile.
fn fade(t: f32) -> f32 { return t * t * (3.0 - 2.0 * t); }

fn cell(ix: i32, iy: i32, r: i32) -> f32 {
    return lat[u32((iy % r) * r + (ix % r))];
}

@compute @workgroup_size(8, 8)
fn main(@builtin(global_invocation_id) id: vec3<u32>) {
    if (id.x >= p.n || id.y >= p.n) { return; }   // grid need not divide by 8
    let r = i32(p.r);
    var sum: f32 = 0.0; var amp: f32 = 1.0; var freq: f32 = 1.0; var norm: f32 = 0.0;
    for (var o: u32 = 0u; o < p.oct; o = o + 1u) {
        // The per-octave offset decorrelates the octaves. Without it they
        // sample the same lattice cells and the sum has visible structure.
        let px = f32(id.x) / f32(p.n) * p.base * freq + f32(o) * 17.0;
        let py = f32(id.y) / f32(p.n) * p.base * freq + f32(o) * 17.0;
        let ix = i32(floor(px)); let iy = i32(floor(py));
        let a = cell(ix, iy, r);         let b = cell(ix + 1, iy, r);
        let c = cell(ix, iy + 1, r);     let d = cell(ix + 1, iy + 1, r);
        let ux = fade(px - floor(px));   let uy = fade(py - floor(py));
        let lo = a + (b - a) * ux;
        let hi = c + (d - c) * ux;
        // Ridged: fold EACH octave, never the sum. Squaring sharpens crests.
        let r = 1.0 - abs(2.0 * (lo + (hi - lo) * uy) - 1.0);
        sum = sum + amp * r * r;
        norm = norm + amp;
        amp = amp * 0.5; freq = freq * 2.0;
    }
    out[id.y * p.n + id.x] = sum / norm;
}
```

**The bounds check is not optional.** A dispatch covers whole workgroups, so
a 100-cell grid with `workgroup_size(8)` runs 104 invocations per axis and
the last four would write past the buffer.

**The fold is what makes ridges, and it must be applied per octave.**
`1 - |2v - 1|` reflects a noise value about its midpoint, turning a smooth
hill into a crease; squaring it sharpens the crest. Doing that to each
octave and then summing gives ridgelines at every scale.

Folding the **sum** instead is a trap worth naming, because the code looks
almost identical and the output looks plausible until you measure it. The
sum of six octaves, normalised by their total amplitude, is a weighted
average of six independent samples, so by the central
limit theorem it clusters tightly around 0.5, and 0.5 is exactly where the
fold peaks. The commonest value therefore maps to the *maximum* height. The
result is a plateau with dips punched in it rather than mountains, and the
height histogram gives it away at once:

```text
folding the sum      0:0%  1:2%  2:3%  3:6%  4:9% 5:12% 6:18% 7:22% 8:19% 9:8%
folding each octave  0:3% 1:10% 2:15% 3:17% 4:17% 5:15% 6:12%  7:7%  8:3% 9:1%
```

The second is what terrain looks like: much low ground, a thin tail of
peaks. `--histogram` prints this for any run.

```{index}
single: thermal erosion
single: PDL; slices for neighbour lookup
single: ping-pong buffers
```

## Erosion: the same pass three ways

Raw fractal noise looks synthetic because real landscapes have been worn.
**Thermal erosion** is the cheapest correction: material sitting at a steeper
angle than the rock holds slides to lower neighbours. The talus angle is that
steepest held slope.

Plain Perl states it directly:

```perl
for my $pass (1 .. $EROSION) {
    my @n = @h;
    for my $y (0 .. $N - 1) {
        my $ym = $y > 0 ? $y - 1 : 0;
        my $yp = $y < $N - 1 ? $y + 1 : $N - 1;
        for my $x (0 .. $N - 1) {
            my $xm = $x > 0 ? $x - 1 : 0;
            my $xp = $x < $N - 1 ? $x + 1 : $N - 1;
            my $c = $h[$y * $N + $x];
            my $moved = 0;
            for my $yy ($ym, $y, $yp) {
                for my $xx ($xm, $x, $xp) {
                    next if $xx == $x && $yy == $y;
                    my $d = $c - $h[$yy * $N + $xx];
                    $moved += ($d - $TALUS) * $RATE if $d > $TALUS;
                }
            }
            $n[$y * $N + $x] = $c - $moved;
        }
    }
    @h = @n;
}
```

PDL replaces the two inner loops with eight whole-array expressions. The
neighbour lookup is the part worth getting right: pad the field once with
edge replication, then read each neighbour as a **lazy slice** of the padded
copy. A slice is a view, so a pass costs one padded copy rather than one copy
per neighbour.

```perl
sub padded {
    my ($P) = @_;
    my $E = $N + 1;
    my $Q = zeroes($N + 2, $N + 2);
    $Q->slice("1:$N,1:$N")   .= $P;
    $Q->slice("0:0,1:$N")    .= $P->slice("0:0,:");        # left column
    $Q->slice("$E:$E,1:$N")  .= $P->slice("-1:-1,:");      # right column
    # Rows 1 and N already carry the correct left and right edges, so
    # copying them outward fills the top and bottom rows AND the corners.
    $Q->slice("0:$E,0:0")    .= $Q->slice("0:$E,1:1");
    $Q->slice("0:$E,$E:$E")  .= $Q->slice("0:$E,$N:$N");
    return $Q;
}

sub nbr {
    my ($Q, $dx, $dy) = @_;
    return $Q->slice(sprintf "%d:%d,%d:%d", 1+$dx, $N+$dx, 1+$dy, $N+$dy);
}

for my $pass (1 .. $EROSION) {
    my $Q = padded($Hp);
    my $moved = zeroes($N, $N);
    for my $dy (-1, 0, 1) {
        for my $dx (-1, 0, 1) {
            next if $dx == 0 && $dy == 0;
            # lclip(0) is the branchless `if (d > talus)`: everything at or
            # below the talus angle contributes exactly zero.
            $moved += (($Hp - nbr($Q, $dx, $dy)) - $TALUS)->lclip(0) * $RATE;
        }
    }
    $Hp = $Hp - $moved;
}
```

The two Perl-level loops that remain run 8 times per pass, not `8 * N * N`
times. That is the whole PDL technique: the loop nest moves from Perl into
one C pass per expression.

The shader keeps the neighbour loop, because on the GPU it runs inside a
single invocation that is already one of `N * N` running at once:

```wgsl
fn at(x: i32, y: i32, n: i32) -> f32 {
    return src[u32(clamp(y, 0, n - 1) * n + clamp(x, 0, n - 1))];
}

@compute @workgroup_size(8, 8)
fn main(@builtin(global_invocation_id) id: vec3<u32>) {
    let n = i32(p.n);
    if (id.x >= p.n || id.y >= p.n) { return; }
    let x = i32(id.x); let y = i32(id.y);
    let h = at(x, y, n);
    var moved: f32 = 0.0;
    for (var dy: i32 = -1; dy <= 1; dy = dy + 1) {
        for (var dx: i32 = -1; dx <= 1; dx = dx + 1) {
            if (dx == 0 && dy == 0) { continue; }
            let d = h - at(x + dx, y + dy, n);
            if (d > p.talus) { moved = moved + (d - p.talus) * p.rate; }
        }
    }
    dst[u32(y * n + x)] = h - moved;
}
```

**Invariant: `src` and `dst` must be different buffers.** Invocations run in
an undefined order, so a shader that updated the grid in place would have
each cell reading a mixture of this pass's and last pass's neighbours. The
result is not merely noisy, it is non-deterministic between runs. Note that
the Perl version obeys the same rule for the same reason, with `my @n = @h`.

The Perl side alternates them, which is the standard **ping-pong**:

```perl
for my $i (1 .. $EROSION) {
    my ($s, $d) = ($i % 2) ? ($bufA, $bufB) : ($bufB, $bufA);
    $gpu->compute($erode, 'main', [$s, $pe, $d], $g, $g);
}
my $final = ($EROSION % 2) ? $bufB : $bufA;      # parity decides which holds it
```

Buffers bind in the order given, so `[$s, $pe, $d]` supplies bindings 0, 1
and 2. The dispatch is two-dimensional because the workgroup is: `$g` is
`ceil(N/8)` on each axis.

**Symptom: the terrain differs slightly between identical runs.** Cause: a
compute shader reading and writing one buffer.

```{index}
single: heightfield
single: normal; heightfield
single: central differences
single: gather and scatter
```

## Normals: analytic, not accumulated

The [meshes chapter](meshes) computes vertex normals by accumulating face
normals into the vertices that share them. That is correct for arbitrary
geometry, but it is a **scatter**: several faces add into one vertex, so the
order matters and the operation does not parallelise without atomics.

A heightfield does not need it. The surface is a function of `x` and `z`, so
the normal follows from the slope by central differences:

```text
n = normalize(-dh/dx, 1, -dh/dz)
```

That is a **gather**: each vertex reads its four neighbours and writes only
itself. Choosing this formulation is what lets the third stage be the same
algorithm in all three backends.

```perl
my $dx = ($h[$y * $N + $xp] - $h[$y * $N + $xm]) / (($xp - $xm) * $S);
my $dz = ($h[$yp * $N + $x] - $h[$ym * $N + $x]) / (($yp - $ym) * $S);
my $l = sqrt($dx * $dx + 1 + $dz * $dz);
my $i = $y * $N + $x;
$nx[$i] = -$dx / $l; $ny[$i] = 1 / $l; $nz[$i] = -$dz / $l;
```

Interior cells span `2 * S` between neighbours; clamped edge cells span only
`S`. Rather than special-case the border, PDL builds the divisor as a field
so the expression stays branchless:

```perl
my $wx = zeroes($N, $N) + 2 * $S;
$wx->slice("0:0,:") .= $S; $wx->slice("-1:-1,:") .= $S;
```

```{index}
single: readback; cost
single: unpack; bulk transfer
single: PDL; get_dataref
```

## Getting the answer back

At 512 the GPU spends `0.010s` computing and `0.034s` returning the result.
**Readback, not compute, is the cost of using a GPU for this.** Three rules
follow, and they are what makes that `0.034s` rather than seconds.

**Transfer in bulk, never per element.** One `unpack` over the whole buffer,
never a loop of accessors:

```perl
my @nx = unpack 'f*', $nb[0]->read;
```

The same rule applies coming out of PDL. A piddle's data slab is reachable
as a string, and reading it in one go beats looping with `at` by 7.6x at
`512 * 512`:

```perl
sub unpiddle { return unpack 'd*', ${ $_[0]->sever->get_dataref } }
```

**Give each output its own buffer.** The shader writes `nx`, `ny` and `nz`
to three storage buffers rather than one interleaved `xyzxyz` buffer, so each
comes back as one `unpack`. De-interleaving in Perl is an `N * N` loop, and
at `1024 * 1024` it costs `0.232s` against a `0.002s` dispatch.

**The same trick works in the other direction.** This program's renderer
wants Perl arrays, so it stops at the `unpack`. When the destination is a
piddle instead, assign the raw buffer straight into its data slab rather
than building a list on the way; that moves 3 MB in `0.003s`:

```perl
my $G = zeroes(float, 3, $N, $N);
my $ref = $G->get_dataref;
$$ref = $raw;
$G->upd_data;
```

```{index}
single: terrain; materials
single: smoothstep; material blending
single: lighting; terrain relief
```

## Materials and light

Picking a colour by height alone produces a relief map, not a landscape.
Three corrections turn one into the other, and each is visible the moment it
is missing.

**Blend across a band, do not switch at a point.** A hard threshold draws a
contour line wherever two materials meet. Interpolating with a smoothstep
over a narrow band removes it:

```perl
my $t = smoothstep(($h - $b + $BLEND) / (2 * $BLEND));
@c = map { $MATC[$i][$_] + ($MATC[$i+1][$_] - $MATC[$i][$_]) * $t } 0..2;
```

**Let steepness override altitude.** A cliff is bare rock at any height.
Without this every band follows an exact contour and the terrain reads as
stacked cut-outs. The surface normal already carries the answer: `n.y` is 1
on flat ground and 0 on a vertical face.

```perl
my $steep = smoothstep((0.80 - $ny) / 0.35);
return map { $c[$_] + ($ROCK->[$_] - $c[$_]) * $steep } 0..2;
```

**Put the boundaries at percentiles, not at fixed heights.** Which heights
count as "high" depends on the terrain, and a fixed cut point that suits one
seed strands another entirely above or below it. Taking the 15th, 50th and
85th percentiles of the actual field keeps the palette balanced whatever the
generator produces:

```perl
my @bin = (0) x 256;
$bin[ int((($_ - $lo) / $span) * 255.999) ]++ for @h;
my ($acc, $j) = (0, 0);
for my $i (0 .. 255) {
    $acc += $bin[$i];
    while ($j <= $#PCT && $acc >= $PCT[$j] * @h) { $BOUND[$j++] = $i / 255 }
}
```

For light, **a low sun is what gives a landscape its relief.** Near the
zenith every upward-facing triangle returns almost the same Lambert term and
the terrain flattens into a white sheet; `--sun-elevation 80` demonstrates
this. A dim cool fill from the sky hemisphere then keeps faces turned away
from the sun as blue shadow rather than flat grey:

```perl
my $sky = 0.5 + 0.5 * $n->[1];
my @col = map {
    my $c = $base[$_] * (0.10 + 0.78 * $lam) + $SKY[$_] * 0.20 * $sky;
    $c > 255 ? 255 : int $c;
} 0 .. 2;
```

```{index}
single: GPU; when it pays
single: PDL; versus GPU
```

## Choosing a backend

| Situation | Use |
|---|---|
| the grid is small, or the code must run anywhere | `perl` |
| the work is array arithmetic and stays in the process | `pdl` |
| the same small operation runs over a very large grid | `gpu` |
| the result is consumed cell by cell in Perl | `pdl`, and keep it in the piddle |

The `145.2x` in the opening table is over plain Perl. Against PDL the GPU is
`14.1x`, and `0.034s` of its `0.043s` is transfer.

The crossover is lower than it looks, because the GPU's cost is almost all
fixed:

| grid | `perl` | `pdl` | `gpu` |
|---|---|---|---|
| 32 | 0.086s | **0.006s** | 0.014s |
| 64 | 0.170s | **0.009s** | 0.014s |
| 96 | 0.302s | 0.018s | **0.015s** |
| 128 | 0.474s | 0.040s | **0.018s** |
| 192 | 0.938s | 0.073s | **0.017s** |
| 512 | 6.198s | 0.607s | **0.043s** |

The GPU column barely moves from 32 to 192 while the grid grows 36-fold: up
to a few hundred thousand cells the dispatch is free and the time is latency
and transfer. PDL crosses it at about 80 cells on a side. Below that the card
is pure overhead.

The card is not a faster computer; it is a machine for doing the same small
thing to a great many values, and the skill is recognising which part of a
problem has that shape.

## Running it

The stages above, plus the mesh builder from [meshes](meshes) and the
rasteriser from [software-rasteriser](software-rasteriser), assemble into one
file:

```bash
pperl worldgen.pl --size 512 --octaves 6 --erosion 12 --seed 1234 \
                  --mesh all --width 800 --height 600 \
                  --out world.png --obj world.obj --show
```

| Option | Effect |
|---|---|
| `--size` | grid edge in cells |
| `--seed` | selects the landscape |
| `--octaves` | how many frequencies of noise to sum |
| `--erosion` | how many passes to weather it; 12 keeps the three-way run quick, 40 gives noticeably smoother ridges |
| `--mesh` | `perl`, `pdl`, `gpu`, `all`, or `auto` for the best available |
| `--azimuth`, `--elevation`, `--distance` | orbit camera; the terrain spans 6 units |
| `--fov`, `--look` | field of view, and the height the camera aims at |
| `--sun-azimuth`, `--sun-elevation` | sun position; low elevation gives the most relief |
| `--obj` | additionally write the mesh as Wavefront OBJ |
| `--show` | display the result in an SDL2 window |
| `--histogram` | print the height distribution and the material boundaries |
| `--no-cull` | disable back-face culling |
| `--power` | `high-performance` or `low-power` adapter preference |

`--mesh pdl` needs a binary built with `--features pdl`, and `--mesh gpu`
needs `--features webgpu` plus a usable adapter. `--mesh auto` picks the
fastest that is present.

`--show` opens a window with the finished PNG through the native `SDL2`
binding of [the 2D guide's toolkit chapter](../graphics/gui-toolkits), so a
run ends with the picture rather than a filename:

```perl
my $win = SDL2::Window->new(title => "worldgen $N x $N", w => $W, h => $H);
my $ren = $win->renderer;
my $tex = SDL2::Image::load_texture($ren, $o{out});
$ren->clear;
$ren->copy($tex);
$ren->present;
```

`--obj` writes ASCII Wavefront OBJ, one line per vertex and one per face, so
it grows fast: 4 MB at `256`, 17 MB at `512`. It is there to open the terrain
in a modelling tool, not as a transport format for large grids.

`--no-cull` should produce a byte-identical image, because the Z-buffer
already resolves what the cull merely avoids computing. If it does not, the
winding is wrong.

```{index}
single: winding order; screen y flip
```

## Winding, and the one line that decides it

```perl
next if $area == 0;                       # degenerate
next if $o{cull} && $area >= 0;           # back-face cull (screen y is down)
my $sgn = $area < 0 ? -1 : 1;             # inside test follows the winding
```

The projection flips `y`, because the screen's origin is top-left and the
mathematics has `y` going up. **A flip is a reflection, and a reflection
reverses orientation**, so every front-facing triangle arrives with
*negative* signed area. Culling on `$area <= 0` therefore discards exactly
the triangles it should keep, and an inside test written to accept only
positive weights then rejects every interior pixel of whatever survived.

Decide the sign once per triangle and let the coverage test follow it.

**Symptom: geometry is sparse or inside out, and inverting the cull changes
which half is missing rather than fixing it.** Cause: the cull and the
coverage test disagreeing about winding.

## What this demonstrates

The division of labour is the lesson. Three stages have the one shape a GPU
is built for, and on a 512 grid they drop from `6.2s` to `0.043s`. The parts
that do not have that shape stay in Perl, where they are readable, and the
renderer among them is now the program's bottleneck at `4.10s`.

Writing the same algorithm three times is not redundancy here. It is what
makes the number in the table trustworthy: the agreement line proves the
three are one algorithm, so the speedup column compares implementations
rather than problems.