Computing Gerstner Waves in GLSL for Ocean Surfaces

Simulating realistic, non-linear ocean surfaces in real-time graphics is commonly achieved through analytical Gerstner waves. Unlike simple sine waves that only alter elevation, Gerstner waves displace surface vertices both vertically and horizontally toward wave crests, forming the sharp peaks and broad troughs characteristic of open water. In a GLSL rendering pipeline, this technique is implemented within the vertex shader, summing multiple wave equations to procedurally deform a flat mesh grid while analytically computing exact normals and tangents for lighting.

Mathematical Formulation

A single Gerstner wave displaces an initial flat point \(\mathbf{P}_0 = (x, 0, z)\) to a new position \(\mathbf{P} = (x', y', z')\). For a single wave \(i\), the displacement is defined as:

\[x' = x + \sum \left( Q_i A_i D_{i,x} \cos(k_i (\mathbf{D}_i \cdot \mathbf{P}_0) + \omega_i t) \right)\] \[z' = z + \sum \left( Q_i A_i D_{i,z} \cos(k_i (\mathbf{D}_i \cdot \mathbf{P}_0) + \omega_i t) \right)\] \[y' = \sum \left( A_i \sin(k_i (\mathbf{D}_i \cdot \mathbf{P}_0) + \omega_i t) \right)\]

Where:

\(\mathbf{D}_i = (D_{i,x}, D_{i,z})\): The normalized 2D horizontal direction vector of the wave.
\(A_i\): The amplitude (wave height).
\(L_i\): The wavelength, which determines the wave number \(k_i = \frac{2\pi}{L_i}\).
\(\omega_i\): The angular frequency, derived from the dispersion relation \(\omega_i = \sqrt{g \cdot k_i}\), where \(g\) represents gravity.
\(Q_i\): The steepness factor controlling the sharpness of the crest. To prevent the surface from self-intersecting and forming unnatural loops, the condition \(\sum Q_i \le 1\) must be maintained across all waves.

Analytical Normals and Tangents

Lighting calculations require accurate surface normals. Rather than sampling neighboring vertices using costly finite difference approximations, Gerstner waves allow for analytical derivation of the tangent (\(\mathbf{T}\)), binormal (\(\mathbf{B}\)), and normal (\(\mathbf{N}\)) vectors using partial derivatives:

\[\mathbf{T} = \frac{\partial \mathbf{P}}{\partial x} = \begin{pmatrix} 1 - \sum Q_i D_{i,x}^2 \sin(\theta_i) \\ \sum D_{i,x} \cos(\theta_i) \\ -\sum Q_i D_{i,x} D_{i,z} \sin(\theta_i) \end{pmatrix}\]

\[\mathbf{B} = \frac{\partial \mathbf{P}}{\partial z} = \begin{pmatrix} -\sum Q_i D_{i,x} D_{i,z} \sin(\theta_i) \\ \sum D_{i,z} \cos(\theta_i) \\ 1 - \sum Q_i D_{i,z}^2 \sin(\theta_i) \end{pmatrix}\]

The surface normal is calculated by taking the cross product of the binormal and tangent vectors:

\[\mathbf{N} = \mathbf{B} \times \mathbf{T}\]

GLSL Implementation

Below is a GLSL vertex shader implementation demonstrating how to sum multiple Gerstner waves to transform an input vertex and output the computed position and normal.

#version 330 core

layout(location = 0) in vec3 aPos;

uniform mat4 model;
uniform mat4 view;
uniform mat4 projection;
uniform float uTime;

out vec3 vWorldPos;
out vec3 vNormal;

struct GerstnerWave {
    vec2 direction;
    float amplitude;
    float steepness;
    float wavelength;
};

const int NUM_WAVES = 4;
uniform GerstnerWave uWaves[NUM_WAVES];

void evaluateGerstnerWave(GerstnerWave wave, vec3 gridPos, inout vec3 displacement, inout vec3 tangent, inout vec3 binormal) {
    float k = 2.0 * 3.14159265 / wave.wavelength;
    float c = sqrt(9.8 / k);
    vec2 d = normalize(wave.direction);
    float f = k * (dot(d, gridPos.xz) - c * uTime);
    float a = wave.amplitude;
    float q = wave.steepness / (k * a * float(NUM_WAVES));

    float sinF = sin(f);
    float cosF = cos(f);

    displacement.x += q * a * d.x * cosF;
    displacement.z += q * a * d.y * cosF;
    displacement.y += a * sinF;

    tangent += vec3(
        -q * d.x * d.x * sinF,
        d.x * cosF,
        -q * d.x * d.y * sinF
    );

    binormal += vec3(
        -q * d.x * d.y * sinF,
        d.y * cosF,
        -q * d.y * d.y * sinF
    );
}

void main() {
    vec3 gridPos = aPos;
    vec3 displacement = vec3(0.0);
    vec3 tangent = vec3(1.0, 0.0, 0.0);
    vec3 binormal = vec3(0.0, 0.0, 1.0);

    for (int i = 0; i < NUM_WAVES; ++i) {
        evaluateGerstnerWave(uWaves[i], gridPos, displacement, tangent, binormal);
    }

    vec3 finalPos = gridPos + displacement;
    vec3 calculatedNormal = normalize(cross(binormal, tangent));

    vWorldPos = vec3(model * vec4(finalPos, 1.0));
    vNormal = normalize(mat3(transpose(inverse(model))) * calculatedNormal);

    gl_Position = projection * view * vec4(vWorldPos, 1.0);
}

Performance and Precision Considerations

When computing wave sums in GLSL, evaluating four to eight waves yields sufficient visual chaos to break repetitive patterns while staying within tight frame budgets. To eliminate high-frequency precision jitter caused by large uTime values, time parameters should be wrapped modulo the least common multiple of wave periods on the CPU before binding uniforms to the shader program.