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, the camera and rasteriser from software-rasteriser, the Lambert term from shading, the writer from model-formats.
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.
Three stages, one shape#
Stage | Input | Output | Per cell |
|---|---|---|---|
terrain | noise lattice | height | sum |
erosion | height | height | compare 8 neighbours |
normals | height | 3 components | 4 neighbours, one |
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.
One lattice, three backends#
The noise lattice is built once, in Perl, and handed to whichever backend runs:
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.
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.
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:
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.
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:
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.
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:
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:
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.
Normals: analytic, not accumulated#
The meshes chapter 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:
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.
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:
my $wx = zeroes($N, $N) + 2 * $S;
$wx->slice("0:0,:") .= $S; $wx->slice("-1:-1,:") .= $S;
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:
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:
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:
my $G = zeroes(float, 3, $N, $N);
my $ref = $G->get_dataref;
$$ref = $raw;
$G->upd_data;
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:
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.
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:
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:
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;
Choosing a backend#
Situation | Use |
|---|---|
the grid is small, or the code must run anywhere |
|
the work is array arithmetic and stays in the process |
|
the same small operation runs over a very large grid |
|
the result is consumed cell by cell in Perl |
|
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 |
|
|
|
|---|---|---|---|
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 and the rasteriser from software-rasteriser, assemble into one file:
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 |
|---|---|
| grid edge in cells |
| selects the landscape |
| how many frequencies of noise to sum |
| how many passes to weather it; 12 keeps the three-way run quick, 40 gives noticeably smoother ridges |
|
|
| orbit camera; the terrain spans 6 units |
| field of view, and the height the camera aims at |
| sun position; low elevation gives the most relief |
| additionally write the mesh as Wavefront OBJ |
| display the result in an SDL2 window |
| print the height distribution and the material boundaries |
| disable back-face culling |
|
|
--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, so a run ends with the picture rather than a filename:
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.
Winding, and the one line that decides it#
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.