Modeling ocean waves can be a difficult task. Waves are generated from the interaction between the wind and the ocean surface and this interaction is usually overly simplified during modelisation. There are two main classes of solution to generate waves, either using the Fourier spectra of the ocean [2], or representing the waves as a sum of trochoids [1].

This second representation is called the Gerstner model. When Fournier and Reeves presented this model, they wanted to be able to modelise waves breaking along the shores as well as those in the middle of the ocean. Here, we focus on deep ocean wave so we will put aside the changes induced by the shores.

Their model present several interesting caracteristics

- Generate realistic waves
- Generate realistic motion
- Generate a wide range of waves, from ripples to high storm waves
- Allow some user control

Although this model is realistic, it does NOT give any physical answer. This model might be useless to an oceanographer but it looks like the real ocean. In computer graphics and rendering, that is all we ask 🙂

A small flaw with this paper would be that it only represent a single wavetrain. The rendered geometry is therefore lacking the real complexity of the ocean (which might be added by normal mapping though). However, it has the advantage that the whole ocean can be parametered only by a wind speed.

They use throchoïds to build their wave. This means that for each x point on the surface, this point is displaced on a circular orbit. This simulates the compression that happens in water motion, generating the waves.

In 1D with the y axis pointing up, the equations arewhere x is the position of the point and t the time. From this, any wavetrain can be generated depending on the parameters

- h : the height of the wave
- k : the spatial frequency of the wave
- ω: the temporal frequency of the wave

But if you chose those parameters freely and randomly, you might see this kind of artefact appear

This is the wave crossing itself. This kind of self intersection artefact appears when k*h > 1.

The **dispersion relation** of the waves gives us the relation between k and ω

In their paper, [1] propose several formula linking all the parameters to a wind speed V; allowing to generate a correct wavetrain with a single parameter.

- h = 0.5 * 0.001 * 7.065 * V
^{2.5} - ω = (g/V) * √ (2/3)

- k = ω² / g

where g is the gravitational constant.

Which gives as a primary result for deep ocean waves

float hfactor(float V) { return 2 * 0.001 * 7.065 * pow(V, 2.5); } float wfactor(float V) { return (9.8 / V) * sqrt(2.0/3.0); } float kfactor(float V) { float w = wfactor(V); return (w*w)/9.8; } void main() { vec3 dpos = position; float h = hfactor(windspeed); float k = kfactor(windspeed); float w = wfactor(windspeed); dpos.y = -h * cos(k*dpos.x - w*time); dpos.x += h * sin(k*dpos.x - w*time); vertex_position = vec4(dpos.xyz, 1); gl_Position = mvp * vertex_position; }

This is nice but if we stick to [1], we are still limited to a single wavetrain. Also, in [1], the question of wind direction is handled as a rotation of the wavetrain but they do not get into the details.

Therefore, we will move to [4]. This paper from Tessendorf extends the model of [1] to handle several wavetrains. Furthermore, he handles the wind direction by simple dot products in the equations. The issue is that the relations between the variables presented above do not hold anymore as the equations in which they are to be applied changed. The only relation that hold is the dispersion relation.

The new equations are

Where

- x is the displaced point, and x0 its position at rest.
- y is its coordinate in the up axis.
- k is now a vector, called the wavevector. Its norm gives the wavelength and its direction gives the direction of the wave.
- ω is the temporal frequency of the wave.
- h is the amplitude of the wave
- Φ is the phase of the wave

As you may note on those equations, what Tessendorf did was to sum several gerstner wavetrains with different parameters. The main issue I had while implementing this was to find a set of parameters that would give correct results.

In the end I went with generating random parameters. I used the following values :

- ω = random number in [1:6]
- |
**k**| = ω² / 9.8 - h = 1 / |k| * 0.05 (value chosen empirically).

By using this and a fixed k, you can get a nice moving plane but nothing realistic. To add some randomness to it, I switched to the polar coordinates;

**k** (x, y) becomes **k** (|**k**|, θ). For x > 0 and y > 0, θ = arctan(y/x) ( http://en.wikipedia.org/wiki/Polar_coordinate_system ).

Then, I added a random angle from [-0.75, 0.75] to θ.

How does it translates in code;

What I did was to load a grid of points on the GPU and then, in a vertex shader, I displaced all the vertices according to the equations of Tessendorf. I generate the parameters for the waves in a pre process, and moves it to the GPU in a 1D Texture. The parameters to write are ω, h (name A in the code), |**k**|sin(θ) and |**k**|cos(θ).

CAREFUL, the texture must NOT be filtered 🙂 ! You do not want the parameters from the first wavetrain to impact the second ! To avoid this use

glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_NEAREST_MIPMAP_NEAREST); glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);

and fill the texture with

unsigned int nbtrains = 128; float* waves = (float*)malloc(sizeof(float)*4*nbtrains); for(unsigned int i=0; i<nbtrains; i++) { float w = (float)( rand()%30 + 7 ); // w in [7:36] w /= 7.0; //w in [1:36/7] float k_norm = (w*w)/9.8; float A = 1.0 / k_norm; A *= 0.05;//( (float)( rand()%100 ) / 100.0 ); waves[i*4] = A; waves[i*4 + 1] = w; float theta_off = ( (float)(rand()%10 - 5)/7.0 ); float x = winddir.x(); float y = winddir.y(); float theta = atan(y / x); theta += theta_off; float sin_theta = sin(theta); float cos_theta = cos(theta); waves[i*4 + 2] = k_norm * sin_theta; waves[i*4 + 3] = k_norm * cos_theta; }

Then, all is left is to read this in the vertex shader.

vec2 xz = dpos.xz; float h = 0; for(int i=0; i<nb_trains; i++) { vec4 param = textureLod(waves_param, (float(i)/128.0) + 0.0001, 0); float A = param.x; float w = param.y; float sin_t = param.z; float cos_t = param.w; vec2 k = vec2(cos_t, sin_t); float phi = dot(k,dpos.xz) - w*time + i; xz -= normalize(k)*A*sin(phi); h += A * cos(phi); } vertex_position = vec4(xz.x, h, xz.y, 1);

To get a nice ocean mesh 🙂 !

Of course, with this solution the computations are looooooooooooooooooooooooooong …

There is a solution to increase the speed by using the FFT (Fast Fourier Transform) to compute the sum of sin/cos faster by I did not investigate it (yet)

*ACM Siggraph Computer Graphics*, 1986, vol. 20, pp. 75–84.

*Computer Graphics and Applications, IEEE*, vol. 7, no. 3, pp. 16–23, 1987.

*Computer Graphics Forum*, 2010, vol. 29, pp. 487–496.

*Simulating Nature: Realistic and Interactive Techniques. SIGGRAPH*, vol. 1, 2001.