Recursive Gaussian Random Fields

Almost every procedural surface you have ever seen (terrain, clouds, marble, water, the bumps on a rock) is the same object underneath: a gaussian random field with a power law spectrum, built by adding random detail one scale at a time. The recipes have different names (fBm, octave noise, midpoint displacement, Diamond-Square, spectral synthesis) and they come from different decades and different fields, but they are all approximating the same thing, and they all end up with a single exponent doing all the work.

This post is the theory of that object. No engine code, just the math and three demos to poke at.

What a gaussian random field is

A random field is a random variable attached to every point of space. Instead of one number with a distribution, you have a function $F(x)$ where each $F(x)$ is random, and the interesting part is how the values at different points relate to each other.

It is called gaussian when any finite set of points you pick has a jointly gaussian distribution. That is a strong condition and it buys you a lot, because a multivariate gaussian is completely determined by two things:

$$\mu(x) = \mathbb{E}[F(x)] \qquad C(x, y) = \mathbb{E}\big[(F(x) - \mu(x))(F(y) - \mu(y))\big]$$

a mean and a covariance function. That is the whole model. No third moments, no skew, nothing else to specify. If you can write down $C$, you have described the field completely, and everything anybody generates in practice is some way of realizing a particular $C$.

Two extra assumptions make it tractable:

• Stationary (or homogeneous): $C(x,y)$ depends only on the offset $h = y - x$, so the statistics look the same everywhere. No special origin.
• Isotropic: $C$ depends only on $|h|$, so the statistics look the same in every direction too.

With both, the entire field is described by a single function of a single scalar: $C(|h|)$.

Why gaussian, and not something else

Two reasons, and they are both good ones.

The first is the central limit theorem. Sum up a lot of small independent random contributions (which is exactly what every construction below does) and the result converges to a gaussian whether you asked for it or not. You end up with a gaussian field by accident, so you may as well design for it.

The second is that among all distributions with a given covariance, the gaussian is the one with maximum entropy. It is the least committed choice you can make: it assumes the covariance and literally nothing else. If all you know about your terrain is "how correlated are two points at distance h", the gaussian field is the honest answer.

There is a third reason that is more of an aesthetic argument, but it shows up constantly in graphics code. The isotropic gaussian is the only distribution that both factors into independent per axis components and is rotationally symmetric, because

$$e^{-x^2} e^{-y^2} e^{-z^2} = e^{-(x^2 + y^2 + z^2)} = e^{-r^2}$$

depends only on the radius. That is why the standard way to get a uniformly random direction in 3D is to normalize three independent gaussians, and why the naive alternative (normalize a vector sampled inside a cube) is subtly wrong:

The cube version is denser along the eight cube diagonals, because the corners stick out further than the faces do. Gaussians have no corners to stick out.

The spectrum is the covariance

For a stationary field there is a second, usually more useful description. The Wiener-Khinchin theorem says the power spectral density is the Fourier transform of the covariance function:

$$S(k) = \hat{C}(k)$$

so specifying $S$ and specifying $C$ are the same act. And in the Fourier basis, a stationary gaussian field is trivial: the Fourier coefficients are independent complex gaussians whose variance is $S(k)$. All the correlation structure in space turns into independent random numbers with the right variances in frequency.

The family everybody uses in graphics is the power law:

$$S(k) \propto |k|^{-\beta}$$

Only one parameter. The reason this specific family dominates is scale invariance: a power law is the only spectrum with no built in feature size. Zoom in on it and, after rescaling the amplitude, you get statistically the same picture. Which is exactly the property mountains, coastlines and clouds have, and exactly what you want when the camera can get arbitrarily close.

Fractional Brownian motion and the Hurst exponent

The canonical scale invariant gaussian field is fractional Brownian motion. Its defining property is how the variance of increments grows with distance:

$$\mathbb{E}\big[|F(x + h) - F(x)|^2\big] \propto |h|^{2H}$$

$H \in (0, 1)$ is the Hurst exponent, and it is the only knob.

• $H = 0.5$ is ordinary Brownian motion, increments independent, variance growing linearly with distance.
• $H > 0.5$ means increments are positively correlated. Going up tends to keep going up, so you get smooth, persistent, rolling shapes.
• $H < 0.5$ means increments are anticorrelated. Going up tends to be followed by going down, so you get rough, jittery, noisy shapes.

fBm is self affine rather than self similar: scaling the domain by $\lambda$ requires scaling the amplitude by $\lambda^H$ to get a statistically identical field. Space and value scale differently, which is why the vertical exaggeration on procedural terrain is always a fudge factor you have to tune.

The exponent connects to the spectrum and to the fractal dimension of the graph:

$$\beta = 2H + d \qquad\qquad D = d + 1 - H$$

where $d$ is the dimension of the domain. A 2D height field with $H = 0.5$ has $\beta = 3$ and $D = 2.5$, sitting halfway between a smooth surface and something that fills its bounding volume. This is the formula of the whole subject, and every construction below is just a different way of hitting a target $H$.

Construction 1: octave summation (fBm noise)

The one everybody writes first, because it is five lines:

$$F(x) = \sum_{k=0}^{N} A_0 \, g^k \, noise(x \cdot L_0 \lambda^k)$$

Each octave $k$ evaluates a smooth band limited noise function at a frequency multiplied by the lacunarity $\lambda$ (usually 2) and an amplitude multiplied by the gain $g$ (usually 0.5). Frequency goes up geometrically, amplitude goes down geometrically, and the ratio between those two rates is what you are actually setting:

$$g = \lambda^{-H} \quad \Longleftrightarrow \quad H = -\frac{\ln g}{\ln \lambda}$$

The default gain = 0.5, lacunarity = 2 gives $H = 1$, which is smoother than most people expect. If your terrain looks like soft hills and you wanted crags, that is why: you want gain closer to 0.7.

Switch to "1D fBm" to see the octaves drawn separately above the sum. That view makes the trade obvious, since the last two or three octaves contribute nearly all the visible roughness while the first one contributes nearly all the shape.

A few things worth knowing about this construction:

• It is only an approximation of a power law spectrum. You are sampling the spectrum at a geometric series of frequencies instead of covering it continuously, so the spectrum comes out lumpy rather than smooth. Nobody notices, but it is why octave noise and true fBm are not identical objects.
• It is not exactly gaussian either, since each octave of gradient noise is bounded, not gaussian. With enough octaves the CLT drags it towards gaussian anyway.
• It costs $O(N)$ per evaluated point, but it is evaluable anywhere, with no grid and no precomputation, which is why it won in graphics. You can call it from a shader at an arbitrary position and get a consistent answer.

Construction 2: midpoint displacement

Older (Fournier, Fussell and Carpenter, 1982) and genuinely recursive. Start with the endpoints of an interval. Take the midpoint, set it to the average of its neighbours, add a random offset. Halve the offset scale, recurse into both halves. In 2D the same idea becomes Diamond-Square, alternating between diamond and square neighbourhoods on a grid.

The amplitude schedule is where the theory comes in. At recursion level $k$ the spacing is $2^{-k}$, and to hit a target Hurst exponent the displacement standard deviation has to be

$$\sigma_k \propto 2^{-kH}$$

which is the same geometric decay as octave summation with $\lambda = 2$ baked in permanently. The recursion halves the spacing, so lacunarity is not a free parameter here. The roughness constant people tune in Diamond-Square implementations is $2^{-H}$.

Try the "1D midpoint displacement" mode in the demo above and compare it to the 1D fBm mode at the same gain. Same character, different construction.

The reason this method is not a straight win, despite being $O(n)$ and dead simple:

• It is not stationary. Points created at early levels have their values fixed forever and never get perturbed again, while their neighbours keep receiving detail. The variance of the field ends up depending on where you are relative to the grid, and you get visible creases along the early subdivision lines. The standard fix is random successive additions (Voss): perturb every point at every level, not just the new midpoints.
• It is grid locked. You generate the whole array or nothing. There is no "evaluate at this arbitrary point", which is exactly the thing octave noise gives you.
• The interpolation is a heuristic. If you want it to be exact, the midpoint value should be drawn from the conditional gaussian distribution given the neighbours, which for fBm is not simply their average plus fixed noise.

That last point is the honest framing of the whole method: midpoint displacement is a cheap approximation to conditional sampling of a gaussian field. Which is why it is fast, and why it has artifacts.

Construction 3: spectral synthesis (the exact one)

If you want the actual gaussian field with the actual spectrum, do it in frequency space, where the coefficients are independent. Fill a complex buffer with independent gaussian noise scaled by $|k|^{-\beta/2}$, inverse FFT, take the real part:

$$F = \mathcal{F}^{-1}\Big[ |k|^{-\beta/2} \cdot \mathcal{N}(0,1) \Big]$$

That is it. One pass, no recursion, no octaves, and the result is exactly a stationary gaussian field with power spectrum $|k|^{-\beta}$, up to the discretization.

Left is the amplitude spectrum being generated (zero frequency in the middle), right is what comes out of the inverse transform. Slide $\beta$ from 2 to 5 and you can watch the spectrum's falloff steepen and the field go from sandpaper to hills. The readout tracks $H = (\beta - d)/2$ and $D = d + 1 - H$.

The trade offs versus the other two:

• Exact, in the sense that the resulting covariance is the one you asked for, with none of the lumpiness of octave sampling.
• Cost is $O(n^d \log n)$ for the whole domain, which is excellent per sample and useless if you only wanted one point.
• The output is periodic, since the DFT assumes it. This is either a bug (visible tiling on large terrain) or a feature (free seamless tiling textures), depending on your use case.
• You cannot evaluate it lazily. The FFT is all or nothing, which rules it out for infinite worlds.

There is a fourth method worth knowing exists, which is the fully general one: write down the covariance matrix $C$ over your sample points, Cholesky factorize it as $C = LL^T$, and multiply $L$ by a vector of independent gaussians. It works for any covariance, not just power laws, and it is $O(n^3)$, so it is only usable on small point sets. Circulant embedding is the trick that turns it back into an FFT for stationary fields on regular grids, which is really what spectral synthesis is doing under the hood.

Picking a method

In practice: octave noise if it has to live in a shader or an open world, spectral synthesis if you want a finite heightmap that is statistically correct (or a seamless tile), midpoint displacement if you want a fast heightmap and can live with the creases, and Cholesky if you are doing statistics rather than graphics.

The one idea

All three constructions are the same sentence with different nouns:

Fractal detail is a geometric series of decaying random perturbations, one per scale level.

Octave noise spends the decay on the amplitude of a continuous field. Midpoint displacement spends it on the size of a displacement during recursive subdivision. Spectral synthesis spends it on the variance of a Fourier coefficient. The exponent controlling the decay is $H$ every single time, and the fractal dimension is always $d + 1 - H$.

Once you see that, the parameter names in every noise library stop being arbitrary. gain, roughness, persistence, H, beta, fractal dimension are six spellings of one number, and converting between them is a couple of logarithms.