How to Implement Fractional Brownian Motion in GLSL

Fractional Brownian Motion (fBm) is a procedural technique used in computer graphics to generate complex, natural-looking patterns such as terrain, clouds, and turbulent fluids. In GLSL, fBm is implemented by running a loop that accumulates multiple layers—or octaves—of a base noise function, scaling the frequency up and the amplitude down with each successive iteration. This article details the mathematical principles behind fBm, provides a modular GLSL implementation, and explores techniques to optimize and customize the resulting procedural textures.

The Concept of Octave Accumulation

Natural phenomena exhibit self-similarity across different scales: large structures define the broad silhouette, while progressively smaller structures add micro-detail. Fractional Brownian Motion models this behavior by summing scaled evaluations of a base noise function (such as Perlin, Simplex, or Value noise).

The accumulation relies on four primary parameters:

• Octaves: The total number of noise layers computed. Higher values yield more fine detail at the cost of performance.
• Amplitude: The weight or intensity of each octave.
• Frequency: The spatial density or scale of the noise features.
• Lacunarity: The factor by which the frequency increases each octave (commonly set to 2.0).
• Gain (or Persistence): The factor by which the amplitude decreases each octave (commonly set to 0.5).

Mathematically, the fBm function \(f(\mathbf{x})\) can be defined as:

\[f(\mathbf{x}) = \sum_{i=0}^{N-1} A \cdot g^i \cdot \text{noise}(\mathbf{x} \cdot f_0 \cdot l^i)\]

where \(N\) is the number of octaves, \(g\) is the gain, \(l\) is the lacunarity, \(A\) is the initial amplitude, and \(f_0\) is the initial frequency.

GLSL Implementation

To implement fBm in a GLSL fragment shader, you begin with a base 2D or 3D noise function and iteratively evaluate it within a loop.

1. The Base Noise Function

A standard 2D pseudo-random gradient or value noise function acts as the fundamental building block:

// Basic pseudo-random hash
float hash(vec2 p) {
    p = fract(p * vec2(123.34, 456.21));
    p += dot(p, p + 45.32);
    return fract(p.x * p.y);
}

// 2D Value Noise
float noise(vec2 p) {
    vec2 i = floor(p);
    vec2 f = fract(p);

    // Quintic interpolation curve
    vec2 u = f * f * f * (f * (f * 6.0 - 15.0) + 10.0);

    float a = hash(i);
    float b = hash(i + vec2(1.0, 0.0));
    float c = hash(i + vec2(0.0, 1.0));
    float d = hash(i + vec2(1.0, 1.0));

    return mix(mix(a, b, u.x), mix(c, d, u.x), u.y);
}

2. The Accumulation Loop

The fbm function accumulates the base noise across a set number of octaves:

#define NUM_OCTAVES 6

float fbm(vec2 st) {
    float value = 0.0;
    float amplitude = 0.5;
    float frequency = 1.0;
    float lacunarity = 2.0;
    float gain = 0.5;

    for (int i = 0; i < NUM_OCTAVES; i++) {
        value += amplitude * noise(st * frequency);
        frequency *= lacunarity;
        amplitude *= gain;
    }

    return value;
}

Eliminating Directional Artifacts with Rotation

Evaluating standard grid-aligned noise repeatedly with a lacunarity of 2.0 can cause visible grid artifacts or directional bias. A standard optimization involves applying a rotation matrix and a spatial shift to the sample coordinate between octaves.

float fbmRotated(vec2 st) {
    float value = 0.0;
    float amplitude = 0.5;
    
    // 2D rotation matrix (~27.5 degrees)
    mat2 rot = mat2(cos(0.5), sin(0.5), 
                   -sin(0.5), cos(0.5));

    for (int i = 0; i < NUM_OCTAVES; i++) {
        value += amplitude * noise(st);
        st = rot * st * 2.0 + vec2(100.0); // Rotate, scale, and shift
        amplitude *= 0.5;
    }

    return value;
}

This transforms the coordinate space non-linearly across layers, breaking up repetitive patterns and generating more organic results.

Normalization and Dynamic Range

When summing a geometric series where \(gain = 0.5\), the maximum theoretical sum converges to:

\[\sum_{i=0}^{\infty} (0.5)^{i+1} = 1.0\]

If the initial amplitude is \(0.5\) and the gain is \(0.5\), the output will naturally stay within the \([0.0, 1.0]\) range as long as the base noise outputs values in that range. If your base noise function returns values in the \([-1.0, 1.0]\) range (standard for gradient and simplex noise), you can normalize the final result after accumulation:

float normalizedFbm = fbm(st) * 0.5 + 0.5;

Variations: Turbulence and Ridged Multifractal

The accumulated value can be modified during iteration to produce distinct material characteristics:

• Turbulence: Taking the absolute value of the noise (abs(noise(st))) creates sharp creases resembling fire or marble veins.
• Ridged Noise: Inverting the absolute noise (1.0 - abs(noise(st))) produces sharp ridges, which are useful for realistic mountain ranges.