Growing an Island Out of Two Other Blog Posts

I ended the gaussian random fields post by pointing out that fBm, midpoint displacement and spectral synthesis are three spellings of the same sentence: fractal detail is a geometric series of decaying random perturbations, one per scale level. I ended the sunflower post with almost the exact same sentence, except the thing decaying was a branching random walk instead of a noise field.

Two posts, same closing paragraph. That is not a coincidence I can leave alone, so this post spends both of them on the same object: a heightmap. Branching walk for the mountains, inverse FFT for the fine grain, and one small new idea to shape a coastline. All in three.js, all live in your browser, ending in an actual island you can drag around.

That is the whole post, running. Here is how it is built, layer by layer.

The three ingredients

Macro structure. The sunflower branching walk places a root at the origin and spawns childCount children at distance $\alpha_d = \alpha_0 \beta^d$, recursing with the step shrinking every generation. I said in that post it is "a real tree, not a field", which is exactly the property I want for mountain ridges: a small number of wide, tall spines near the root, fanning out into a large number of thin, short ones near the leaves. That is what an eroded volcanic island actually looks like. Rivers and ridges are trees, not noise.

Fine detail. Spectral synthesis fills a complex buffer with gaussian noise shaped by $|k|^{-\beta/2}$, inverse FFTs it, and the real part is a stationary gaussian field with power spectrum $|k|^{-\beta}$, exact, one pass, no recursion. That is the rock texture, the bumps on top of the ridges, the stuff a branching tree of straight segments cannot give you on its own.

The coastline. Neither of the other two posts covers this, because it is a new small idea, the same one again in yet another domain. More on it below.

Shaping the coastline

Everything so far lives in a plane or a 3D volume. A coastline lives on a circle: it is a function of angle, $r(\theta)$, periodic by construction. A periodic 1D domain does not need an FFT, it is small enough to just write the Fourier series directly:

$$r(\theta) = R_0 \left(1 + A \sum_{m=1}^{M} \frac{\cos(m\theta + \varphi_m)}{m^{p}}\right)$$

$M$ harmonics, random phase $\varphi_m$ per harmonic, amplitude falling off as $m^{-p}$. That is the same amplitude rule as spectral synthesis, $|k|^{-\beta/2}$, with $k$ replaced by the harmonic index and the domain dropped from 2D down to 1D ($d=1$, $\beta = 2p$):

$$H = \frac{\beta - d}{2} = p - \frac12 \qquad\qquad D_{\text{curve}} = d + 1 - H = 2.5 - p$$

That $D_{\text{curve}}$ is the fractal dimension of the coastline itself, not of a height field, of the actual wiggly curve you'd get if you walked its perimeter with a ruler. This is Mandelbrot's coastline paradox by construction: real coastlines measured this way come out around $D \approx 1.2$ to $1.3$ (Britain is the famous number, about $1.25$), which back-solves to $p \approx 1.25$. Set the slider there and the jaggedness should look about right. Lower $p$ keeps the low harmonics loud and you get a smooth blob, basically a lens or an atoll. Higher $p$ kills them fast and you get fjords down to the harmonic cutoff.

Building the ridge skeleton

The branching walk from the sunflower post, flattened from 3D to the XZ plane, and instead of writing a particle it writes a line segment from parent to child:

function generateRidgeTree() {
  const segments = []
  let budget = NODE_BUDGET          // same guard as the original: if budget <= 0, stop
  const alpha0 = R0 * 0.5

  function recurse(x0, z0, stepAlpha, depth) {
    if (depth > ridgeDepth || budget <= 0) return
    for (let i = 0; i < ridgeBranches && budget > 0; i++) {
      const theta = randf() * Math.PI * 2
      const x1 = x0 + Math.cos(theta) * stepAlpha
      const z1 = z0 + Math.sin(theta) * stepAlpha
      segments.push({ x0, z0, x1, z1, depth })
      budget--
      recurse(x1, z1, stepAlpha * ridgeDecay, depth + 1)
    }
  }
  recurse(0, 0, alpha0, 0)
  return segments
}

NODE_BUDGET is the exact same idea as particleCount in the original C++, just renamed. Depth times branch factor grows fast (branch 4, depth 8 is $4^8 \approx 65000$ potential nodes), so the budget is what keeps this from either running forever or, worse, silently rasterizing sixty thousand ridge segments into a heightmap and grinding the tab to a halt.

Rasterizing ridges into a heightfield

A line segment is not a height. Each one gets turned into a ridge by a gaussian bump measured from the perpendicular distance to the segment, width and amplitude both shrinking with depth on the same $\beta^d$ schedule the walk already uses:

const width = Math.max(minWidth, RIDGE_BASE_WIDTH * Math.pow(ridgeDecay, depth * 0.6))
const amp   = RIDGE_MAX_HEIGHT * Math.pow(ridgeDecay, depth)
// ...
const contrib = amp * Math.exp(-perp2 / (width * width))
if (contrib > ridge[idx]) ridge[idx] = contrib

Taking the max instead of summing is the important bit, it is what keeps two nearby ridges from adding up into an unrealistic spike where they cross. Real ridgelines meet at a shared peak, they don't stack.

The other important bit is cheap: because width shrinks with depth, deep segments (and there are a lot of them, the tree is bushiest at the leaves) only touch a small bounding box of grid cells each, while shallow segments (wide, but there are only a handful of them) touch a large box. The work stays roughly balanced across depths instead of blowing up, which is the same shape of argument as "why octree image compression is fast": most of your budget is spent where there is detail to spend it on.

The fine detail pass

This part is copy-pasted, more or less verbatim, from the spectral synthesis demo:

function synthDetailField(n, beta) {
  const re = new Float64Array(n * n), im = new Float64Array(n * n)
  for (let y = 0; y < n; y++) {
    const ky = y <= n / 2 ? y : y - n
    for (let x = 0; x < n; x++) {
      const kx = x <= n / 2 ? x : x - n
      if (kx === 0 && ky === 0) continue
      const k = Math.hypot(kx, ky)
      const amp = Math.pow(k, -beta / 2)
      re[y * n + x] = gauss() * amp
      im[y * n + x] = gauss() * amp
    }
  }
  fft2(re, im, n, true)     // inverse, real part is the field
  return normalize(re)
}

n has to be a power of two, radix-2 FFT demands it. I set it to 128 and reused the same 128 for the mesh resolution, so the terrain has 128×128 vertices and there is no resampling step. That is a simplification, not a requirement: the FFT resolution and the mesh resolution are genuinely independent, you could synthesize detail at 128×128 and bilinearly resample it onto a mesh of any size. I just didn't need to here.

The detail field gets added on top of the ridge field, scaled by the "Detail" slider, so you can dial it from bare fractal ridgelines up to full sandpapered rock.

Combining land and sea

Every grid cell gets a distance and an angle from the island center. Compare distance against the coastline radius at that angle, and squash the transition with a smoothstep instead of a hard cutoff:

const t = dist / coastRadius(theta)
const landMask = 1 - smoothstep(0.55, 1.02, t)
const h = landHeight * landMask + seabed * (1 - landMask)

landMask is 1 well inside the coastline, 0 well outside it, and eases between the two across a band instead of a cliff. seabed slopes gently down and away from shore instead of just being flat, so the ocean floor doesn't look like it was stamped out with a cookie cutter.

Coloring it like an island

No textures, this is all per-vertex color, computed from height and slope right before the geometry gets uploaded:

• underwater: dark to lighter blue with depth
• near sea level: sand
• low-to-mid elevation: grass, unless the slope is steep, then rock (that is what makes cliffs look like cliffs instead of green walls)
• high elevation: rock, then snow near the peak

A little per-vertex hash-based jitter breaks up the color bands so they don't look like contour lines painted on. It is not physically based, it's a handful of if statements on height and a dot product, and it is enough to read as "island" at a glance, which was the actual goal.

What's missing

No erosion. Real terrain gets its final shape from water moving material downhill for a very long time, which rounds gaussian bumps into the smoother, more asymmetric shapes real ridgelines have. This heightmap is the pre-erosion shape, the tectonic skeleton before millions of years of hydraulic sanding. Simulating that is a whole separate post (thermal and hydraulic erosion on a heightmap grid, basically a cellular automaton).

No infinite world. I flagged this exact limitation in the gaussian fields post's comparison table: spectral synthesis is $O(n \log n)$ for the whole domain and cannot be evaluated lazily, which "rules it out for infinite worlds". This demo hits that limitation on purpose, it generates one fixed island, not a chunk-streamed planet. An open-world version would need to swap construction 3 for construction 1 (plain octave noise) in the detail pass, since that one really is evaluable at an arbitrary point with no whole-grid step, at the cost of an approximate rather than exact spectrum. The ridge skeleton wouldn't need to change, since a single island's mountain tree is already a bounded, finite structure by construction (that boundedness is the whole point of the $\beta < 1$ convergence argument from the sunflower post).

No GPU. All three passes (FFT, ridge rasterization, coloring) run once per parameter change on the CPU, in plain JS, and it is fast enough for one island at 128×128 (well under a hundred milliseconds end to end). None of it would survive scaling up to game-quality terrain sizes without moving the ridge rasterization pass to a compute shader, but that is an engineering problem, not a math one, and the math is the part these two posts were actually about.

The one idea, a third time

Branching walk: decay spent on the step length of a discrete process in space, one tree. Spectral synthesis: decay spent on the amplitude of a Fourier coefficient in frequency, one field. Coastline: decay spent on the amplitude of a Fourier coefficient too, just on a circle instead of a plane, one curve.

Three domains, three data structures, one exponent controlling roughness in every single case, and if you set $\beta = 2p$ and stop caring whether the domain is a tree, a plane or a circle, it is genuinely the same formula three times. That is either a deep fact about scale-invariant randomness or a sign that I keep reaching for the same tool regardless of the problem. Probably both.