Phillips Spectrum Ocean Wave Displacement in GLSL

Simulating realistic ocean surfaces in real-time computer graphics relies on statistical oceanography and the fast Fourier transform. Synthesizing ocean wave displacement involves computing the Phillips spectrum in the frequency domain, propagating waves through time using the deep-water dispersion relation, and transforming the frequency amplitudes into spatial displacement maps using a two-dimensional Inverse Fast Fourier Transform (IFFT) executed across GLSL shaders.

1. The Phillips Spectrum Formulation

The Phillips spectrum models wind-driven equilibrium seas by calculating the statistical variance of wave heights for a given wave vector \(\mathbf{k} = (k_x, k_z)\). The mathematical formula for the Phillips spectrum \(P_h(\mathbf{k})\) is:

\[P_h(\mathbf{k}) = A \frac{e^{-\frac{1}{(k L)^2}}}{k^4} |\hat{\mathbf{k}} \cdot \hat{\mathbf{w}}|^2 e^{-k^2 l^2}\]

Where:

  • \(\mathbf{k}\) is the 2D wave vector, and \(k = |\mathbf{k}|\) is its magnitude.
  • \(\hat{\mathbf{k}} = \frac{\mathbf{k}}{k}\) is the unit directional vector of the wave.
  • \(\hat{\mathbf{w}}\) is the normalized wind direction vector.
  • \(V\) is the wind speed, and \(g\) is the gravitational constant (\(9.81 \, \text{m/s}^2\)).
  • \(L = \frac{V^2}{g}\) represents the largest possible wave generated by the wind speed.
  • \(A\) is a global scaling factor adjusting wave amplitudes.
  • \(l\) is an optional suppression factor to damp waves smaller than a specified capillary limit (\(l \ll L\)).

2. Generating the Initial Frequency Domain Amplitudes

To obtain real spatial heights after the transform, the frequency spectrum must satisfy Hermitian symmetry. Two independent Gaussian random variables \(\xi_r, \xi_i \sim \mathcal{N}(0, 1)\) generate the base complex amplitude \(\tilde{h}_0(\mathbf{k})\):

\[\tilde{h}_0(\mathbf{k}) = \frac{1}{\sqrt{2}} (\xi_r + i \xi_i) \sqrt{P_h(\mathbf{k})}\]

The conjugate spectrum for opposite wave vectors is computed as:

\[\tilde{h}_0^*(-\mathbf{k}) = \frac{1}{\sqrt{2}} (\xi_r - i \xi_i) \sqrt{P_h(-\mathbf{k})}\]

3. Time Evolution via the Dispersion Relation

Ocean waves evolve over time based on the gravitational dispersion relation for deep water:

\[\omega(\mathbf{k}) = \sqrt{g k}\]

The dynamic complex amplitude \(\tilde{h}(\mathbf{k}, t)\) at time \(t\) is calculated by advancing both the forward and conjugate frequencies:

\[\tilde{h}(\mathbf{k}, t) = \tilde{h}_0(\mathbf{k}) e^{i \omega(\mathbf{k}) t} + \tilde{h}_0^*(-\mathbf{k}) e^{-i \omega(\mathbf{k}) t}\]

Using Euler's identity (\(e^{i \theta} = \cos \theta + i \sin \theta\)), this evaluation is performed directly in a GLSL compute or fragment shader:

// GLSL snippet for time-dependent spectrum evaluation
vec2 complexMult(vec2 a, vec2 b) {
    return vec2(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
}

vec2 evaluateWaveSpectrum(vec2 h0, vec2 h0Conj, float kLen, float t) {
    float g = 9.81;
    float w = sqrt(g * kLen);
    
    vec2 exp_iwt = vec2(cos(w * t), sin(w * t));
    vec2 exp_inv_iwt = vec2(cos(w * t), -sin(w * t));
    
    return complexMult(h0, exp_iwt) + complexMult(h0Conj, exp_inv_iwt);
}

4. Inverse FFT and Spatial Displacement Synthesis

The physical surface displacement vector \(\mathbf{S}(\mathbf{x}, t) = (D_x, h, D_z)\) is synthesized by computing the 2D Inverse Discrete Fourier Transform (IDFT) of the complex frequency fields.

Vertical Height Field

The vertical displacement \(h(\mathbf{x}, t)\) at horizontal position \(\mathbf{x} = (x, z)\) is synthesized via:

\[h(\mathbf{x}, t) = \sum_{\mathbf{k}} \tilde{h}(\mathbf{k}, t) e^{i \mathbf{k} \cdot \mathbf{x}}\]

Horizontal Choppy Displacement

To create peaked wave crests and flattened troughs, horizontal displacement components \(\mathbf{D}(\mathbf{x}, t)\) are derived using the Hilbert transform of the wave height:

\[\mathbf{D}(\mathbf{x}, t) = \sum_{\mathbf{k}} -i \frac{\mathbf{k}}{k} \tilde{h}(\mathbf{k}, t) e^{i \mathbf{k} \cdot \mathbf{x}}\]

This splits into individual axes:

\[\tilde{D}_x(\mathbf{k}, t) = -i \frac{k_x}{k} \tilde{h}(\mathbf{k}, t)\] \[\tilde{D}_z(\mathbf{k}, t) = -i \frac{k_z}{k} \tilde{h}(\mathbf{k}, t)\]

5. GLSL GPU Implementation Pipeline

Running the summation directly requires \(\mathcal{O}(N^4)\) operations, which is too slow for real-time applications. Instead, an \(\mathcal{O}(N^2 \log_2 N)\) radix-2 Cooley-Tukey butterfly algorithm is implemented across textures:

  1. Spectrum Texture Generation: A compute shader populates an \(N \times N\) 32-bit floating-point texture containing \(\tilde{h}(\mathbf{k}, t)\), \(\tilde{D}_x(\mathbf{k}, t)\), and \(\tilde{D}_z(\mathbf{k}, t)\).
  2. Horizontal 1D FFT Passes: A fragment/compute pass executes \(\log_2(N)\) iterations along the rows, applying butterfly operations: \[X_{k} = X_{k, \text{even}} + e^{-2\pi i k / N} X_{k, \text{odd}}\]
  3. Vertical 1D FFT Passes: Another \(\log_2(N)\) iterations are performed along the columns of the transposed or ping-pong target textures.
  4. Vertex Shader Displacement: The real components of the resulting inverse transform map directly into the 3D vertex position: \[Vertex.xyz = \mathbf{x}_{original} + \lambda \begin{pmatrix} D_x(\mathbf{x}, t) \\ h(\mathbf{x}, t) \\ D_z(\mathbf{x}, t) \end{pmatrix}\] where \(\lambda\) is a choppiness factor preventing self-intersection when adjusted.