Realtime Ocean rendering - I'm lost

Dear fellow OpenGL-ers

I’ve spent the last week or so rendering a simple ocean using gerstner waves but having issues with tiling, so I decided to start rendering them “properly” and dip my toes into the murky waters of rendering a heightfield using an iFFT.

There are plenty of papers explaining the basic gist -

  1. calculate a frequency spectrum
  2. use this to create a heightfield using ifft to convert from frequency domain to spatial domain - animating with time t

Since the beginning of this journey I have learned about things like the complex plane, the complex exponent equation, the FFT in more detail etc but after the initial steps of creating an initial spectrum (rendering a texture full of guassian numbers with mean 0 and sd of 1, filtered by the phillips spectrum) I am still totally lost.

My code for creating the initial data is here (GLSL):

float PhillipsSpectrum(vec2 k){
 
  //kLen is the length of the vector from the centre of the tex
  float kLen = length(k);
  float kSq = kLen * kLen;
 
  // amp is wave amplitude, passed in as a uniform
  float Amp = amplitude;
 
  //L = velocity * velcoity / gravity
  float L = (velocity*velocity)/9.81;
 
  float dir = dot(normalize(waveDir),normalize(k));
 
  return Amp * (dir*dir) * exp(-1.0/(kSq * L * L))/ (kSq * kSq) ;
}
 
void main(){
 
   vec3 sums;
 
   //get screenpos - center is 0.0 and ranges from -0.5 to 0.5 in both
   //directions
   vec2 screenPos = vec2(gl_FragCoord.x,gl_FragCoord.y)/texSize - vec2(0.5,0.5);
 
   //get random Guass number
   vec2 randomGuass = vec2(rand(screenPos),rand(screenPos.yx));
 
   //use phillips spectrum as a filter depending on position in freq domain
   float Phil = sqrt(PhillipsSpectrum(screenPos));
   float coeff = 1.0/sqrt(2.0);
 
   color = vec3(coeff *randomGuass.x * Phil,coeff * randomGuass.y * Phil,0.0);
}  

which creates a texture like this:

wave.png

Now I am totally lost as how to :
a) derive spectrums in three directions from the initial texture

b) animate this according to time t like mentioned in this paper (https://developer.nvidia.com/sites/default/files/akamai/gamedev/files/sdk/11/OceanCS_Slides.pdf) on slide 5

I might be completely stupid and overlooking something really obvious - I’ve looked at a bunch of papers and just get lost in formulae even after acquainting myself with their meaning. Please help.

I’ve only scanned the slides briefly, but a few points:

Slide 8 gives the formula for the “height” as a function of wavenumber and time: H(k,t). This is the initial (complex) amplitude H[sub]0/sub multiplied by a rotation eiωt. One thing which isn’t made clear is that the angular frequency ω is related to the wavenumber k and the phase speed c[sub]p[/sub] by c[sub]p[/sub]=ω/k=λ/T. This wikipedia page may help with some of the background.

If you’re unfamiliar with the notation, a superscript asterisk indicates the complex conjugate: (a+bi)*=(a-bi).

If you’re unfamiliar with Euler’s forumla, the complex exponential eiωt is equivalent to cos(ωt)+i.sin(ωt). If you plot this expression as a function of t, you get a helix of unit radius and pitch 1/ω; the real part is cos(ωt) and the imaginary part sin(ωt), i.e. two sine waves of angular frequency ω, 90 degrees out of phase.

If you consider complex multiplication z=a*b in polar coordinates, the magnitude of z is the product of the magnitudes of a and b, and the argument (angle) of z is the sum of the arguments of a and b. Consequently, multiplying a complex number by eiθ rotates it by θ.

So for any wavenumber k (a function of the wavelength, k=2π/λ, where wavelength is a vector having both length and direction), the “height” of the wave as a function of time is the initial position H[sub]0/sub (a complex number representing the amplitude and phase) rotated by an angle proportional to time. The reason that H(k,t) has two terms is that you will have two waves travelling in opposite directions which combine to form a standing wave. The wavenumber for one is the negative of the other, and the rotation has the opposite direction (e-iωt).

This “height” is complex. The real part describes the vertical motion, while the imaginary part describes the horizontal motion in the direction of travel. The latter provides the definition of D(k,t). Multiplying by i swaps the real and imaginary components (so the real part of D is determined by the imaginary part of H). Multiplying by k/|k| converts the scalar value to a 2-D vector in the direction of k.

In short, the real part of H(k,t) gives you Z(k,t) while the real part of D(k,t) gives you X(k,t) and Y(k,t). Evaluating these for a specific t gives you three textures, X(k), Y(k) and Z(k), where k is the texture coordinate corresponding to a specific wavenumber (any number of waves with the same wavelength and direction but differing amplitudes and phases sum to a single wave).

Getting the final displacement map requires summing the values of all of theses waves at each point … which is exactly what an inverse FFT does.

Hi GClements,
thank you very much! I need to know more, as since I posted this and before your reply I found some code for the SDK sample that the slides represent. However it is in D3D and a little opaque to me. However from the code I can now “update the spectrum” for x,z and h values using the following code.

The two components of each (the real and imaginary) are stored in two textures, with ht stored in one and the X(k,t) and Y(k,t) stored in two components each in a four component texture, as below:


//waveUpdate.fp
#version 400

uniform float t;
uniform vec2  texSize;
uniform sampler2D spectrumTex;

layout(location = 1) out vec4 htOut;
layout(location = 2) out vec4 dXdTdYdTOut;

void main(){
	//sample the original frequency domain texture at two points - 
	//the corresponding texCoords and the inverse
	vec2 texCoord = gl_FragCoord.xy/texSize;
	vec2 otherTexCoord = vec2(1.0,1.0) - texCoord;


	vec3 h0_k   = texture(spectrumTex,texCoord).xyz;
	vec2 h0_mk  = texture(spectrumTex,otherTexCoord).xy;

	//get sin and cos values for current time t
	float sinV = sin(t);
	float cosV = cos(t);

	//get the real and imaginary parts into a vector
	vec2 ht;
	ht.x = (h0_k.x + h0_mk.x) * cosV - (h0_k.y + h0_mk.y) * sinV;
	ht.y = (h0_k.x - h0_mk.x) * sinV + (h0_k.y - h0_mk.y) * cosV;

	//get length from centre of freq domain
	//H(t) -> Dx(t), Dy(t)
	vec2 k = gl_FragCoord.xy - texSize * 0.5;
	float sqK = dot(k,k);

	float rsqr_k = 0.0;
	
	if(sqK > 1e-12){
	  rsqr_k = 1.0 /sqrt(sqK);
	}

	k *= rsqr_k;

	//find dX/dt and dY/dt
	vec2 dtX = vec2(ht.y * k.x, -ht.x * k.x);
	vec2 dtY = vec2(ht.y * k.y, -ht.x * k.y);

	//two render targets
	htOut            = vec4(ht,0.0,1.0);
	dXdTdYdTOut = vec4(dtX,dtY);
}

I will need to perform an IFFT on the results right? However, doesn’t that also get me a complex number made of a real and imaginary number? An example from the compute shader in the D3D sample goes like so:


[numthreads(COHERENCY_GRANULARITY, 1, 1)]
void Radix008A_CS2(uint3 thread_id : SV_DispatchThreadID)
{
	if(thread_id.x >= thread_count)
		return;

	// Fetch 8 complex numbers
	uint i;
	float2 D[8];
	uint iaddr = thread_id << 3;
	for (i = 0; i < 8; i++)
		D[i] = g_SrcData[iaddr + i];

	// Math
	FFT_forward_8(D);

	// Store the result
	uint omod = thread_id & (ostride - 1);
	uint oaddr = ((thread_id - omod) << 3) + omod;
	g_DstData[oaddr + 0 * ostride] = D[0];
	g_DstData[oaddr + 1 * ostride] = D[4];
	g_DstData[oaddr + 2 * ostride] = D[2];
	g_DstData[oaddr + 3 * ostride] = D[6];
	g_DstData[oaddr + 4 * ostride] = D[1];
	g_DstData[oaddr + 5 * ostride] = D[5];
	g_DstData[oaddr + 6 * ostride] = D[3];
	g_DstData[oaddr + 7 * ostride] = D[7];
}

In the above code ostride changes from iteration to iteration to do a correct butterfly operation I believe. I still don’t really see how these complex vals are then used to generate the heightmaps - do you use the real value only?
Also since I’m using results stored in textures I will have to ping pong the textures for output and storing intermediary values right? This level of math is very exhausting to me…

Yes. You’ll only want the real part. Actually, the complex part may end up being zero due to the definition of H(k,t).

Right.

There are FFT algorithms which can perform the transform in-place, but GPU-oriented implementations (e.g. Stockham auto-sort) tend to use distinct input and output textures, partly due to better cache coherence and partly because older hardware didn’t support reading and writing a texture simultaneously (many such implementations were implemented as fragment shaders which rendered a quad).

Thank you again. I’m going to try progressing with this in a few days as all the new concepts have made my head hurt! I think I’m being hard on myself as I feel an idiot for not understanding much despite reading a bunch of papers and looking at source code but it is a fairly complex problem - the papers published more than likely took months of research and I’m hard on myself because I can’t do a full implementation in 1 1/2 weeks.

At least I’ve learned about the transformation from frequency domain -> spatial domain and what that actually means, as well as complex numbers, euler’s formula, gerstner waves etc. I’m good with maths like basic kinematics and calculus but when people throw imaginary numbers and complex planes in the mix my brain decides to exit stage right.

I’m rambling, thanks again. Expect to see another post in a few days.

Probably the best book on complex variables

As a math subject complex variables can be quite beautiful, the various symmetries and their manifestations!