Ocean modelisation – Fourier

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].

Modeling terrains using Fourier is quite known. You get the Fouier Transform of a white noise (so that all frequencies appear), you filter it to select a subset of those frequencies and you use the inverse Fourier Transform on this subset. This gives you a heightmap representing your terrain ( http://paulbourke.net/fractals/noise/ — Article on « Frequency synthesis of Landscapes (and clouds) » ).

The same process can be applied to the oceans. When modeling the ocean using a Fourier spectra, the idea is to represent a set of waves as a sum of sinusoids with different frequencies. The trick is therefore to choose those frequencies carefully. Oceanographic research have actually an answer to this as several spectrum where defined to represent the frequency of waves. By applying a Fourier transform on a subset of those frequency, you can get an heightmap representing a very realistic ocean.

Here, we will used the formulas given by [3]. Those formulas and the corresponding code are already very well explained in the tutorial by Keith Lantz (http://www.keithlantz.net/2011/10/ocean-simulation-part-one-using-the-discrete-fourier-transform/) but i’ll reexplain it anyway.

The idea is to disturb a grid mesh of M*N points.


Those points will be displaced to form the ocean surface. The question is how to displace them. It all starts with an inverse Fourier Transform


where h(x, t) is the height for a point x of coordinates (x, y, z) at a time t. k are the frequencies of the complex exponentials and h~ is the Fourier Transform of h. Since we need as many frequencies as the number of points, we put k = ( (2πn / Lx ), (2πm / Lz) ) with n in [-N/2, N/2] and m in [-M/2, M/2].

To compute h with a discrete fourier transform, we need to sum the computed h~. We have


with ω(k)² = gk, g being the gravitational constant and k being the length of the vector k. This leaves the question of evaluating the complex number h0~, h0~* being its conjugate. Expressing h~ as a sum allows to maintain the complex conjugation property, h~(k, t) = h~*(-k, t).

We know that the waves follow a spatial spectrum similar as the Phillips spectrum for wind driven waves


where L = V²/g; V being the windspeed and w being the direction of the wind. The ^ symbol means that the vectors are normalized. A is a constant (0.0001 in our examples).

This formula is furthermore multiplied by


with l = 0.001 * L. This removes the waves too small to be seen.

From there, our h0~ takes samples from this spectrum as


where ξ are independent draws from a gaussian random number generator. Careful, h0 must be evaluated once at the initialisation of the ocean. Otherwise, at each time t, h will be evaluated from a different h~, and therefore will represent a different ocean.

We also have the following formula to compute the surface slopes.


To get the normal, we just compute : up_vector – slope

From there, we can code a heightmap generator to build our ocean.

Ocean heightmap. Grid of size 100x100. Computed in 77S.
Ocean heightmap. Grid of size 100×100. Computed in 77 seconds.

This gives an ocean surface quite realistic but it is not that satisfying yet. In real life, the ocean surface is not an heightmap

The ocean surface is not a heightmap !
The ocean surface is not a heightmap !

Therefore, Tessendorf presents the Choppy Waves algorithm. If it does not allow to create plunging waves as illustrated, it allow to creates sharp waves. This is very similar to Gerstner waves. Every point of the grid mesh is displaced on x and z to compress or dilate parts of the mesh; simply by adding D(x, t) to the x and z coordinates of the wave. This creates peak waves.


Heightmap ocean
Displaced grid ocean. Peak waves appear in the middle.


Leading to a final rendering

200×200 – 15 minutes to compute


Using a discrete Fourier Transform, we have the following code. However, as given in the caption of the pictures, computing the surface is very slow. This can be improved by using a Fast Fourier Transform, but this will be another story !

A small optimisation that can divide almost by 2 the computation times. For each h0~ that is computed, only part of them will impact strongly the final surface as some of them will lead to waves with amplitudes too small to be seen. What can be done in the initialisation therefore is to test the values of the h0~ and to keep only those above a thresold. Then, when computing h~, the sum will be done on much less k and will therefore run faster.

Grid of 100×100. No Thresholding (10 000 elements in the sum). Computed in 75 seconds


Grid of 100×100. Threshold at 0.01, leading to 5775 frequencies to compute the h~ sum. Computed in 44 seconds
class Ocean
    Ocean(float n, float m, 
          float lx, float lz, 
          Vector2 windir, float windspeed);
    void computeHeightmap(float time);
    Mesh* getMesh();
    int N, M; //Nombre de samples
    float Lx, Lz; //Distance entre deux samples
    Vector2 wind;
    Vector3* geometry;
    Vector3* normales;
    int nb_htilde;
    Complex* htilde0_tab;
    Complex* htilde0c_tab;
    Vector2* k_tab;
    unsigned int* indices;
    //h(x, t) + D(x, t) + n(x, t)
    void evaluate(Vector2 x, float t, 
                  float* h, Vector3* n, Vector2* disp);

    Complex htilde(Vector2 k, float t);    
    Complex htilde0(Vector2 k);
    double Phillips(Vector2 k);
    double dispersion(Vector2 k);
Ocean::Ocean(float n, float m, 
             float lx, float lz, 
             Vector2 windir, float windspeed)
    M = m;
    N = n;
    Lx = lx;
    Lz = lz;
    wind = windir.normalize() * windspeed;
    geometry = (Vector3*)malloc((M)*(N)*sizeof(Vector3));
    normales = (Vector3*)malloc((M)*(N)*sizeof(Vector3));
    htilde0_tab  = (Complex*)malloc((M)*(N)*sizeof(Complex));
    htilde0c_tab = (Complex*)malloc((M)*(N)*sizeof(Complex));
    k_tab = (Vector2*)malloc((M)*(N)*sizeof(Vector2));
    indices = (unsigned int*)malloc( 
                         (M-1)*(N-1)*sizeof(unsigned int)*6);

    nb_htilde = 0;    
    int index;

    for (int m = -M/2; m < M/2; m++)
    for (int n = -N/2; n < N/2; n++) {

        int n_prime = n + (N/2);
        int m_prime = m + (M/2);
        index = m_prime * (N) + n_prime;
        Vector2 k( 2 * M_PI * n / Lx,
        2 * M_PI * m / Lz);

        htilde0_tab[nb_htilde] = htilde0(k);
        if (htilde0_tab[nb_htilde].getModule() > 0.01)
            k_tab[nb_htilde] = k;
            htilde0c_tab[nb_htilde] = 

        geometry[index].x ( n * Lx / N );
        geometry[index].y ( 0.0f );
        geometry[index].z ( m * Lz / N );
        normales[index].x ( 0 );
        normales[index].y ( 1 );
        normales[index].z ( 0 );

    int indices_count = 0;
    for (int m_prime = 0; m_prime < M-1; m_prime++)
    for (int n_prime = 0; n_prime < N-1; n_prime++) {
        index = m_prime * (N) + n_prime;
        indices[indices_count++] = index; 
        indices[indices_count++] = index + 1;
        indices[indices_count++] = index + N;
        indices[indices_count++] = index + 1;
        indices[indices_count++] = index + N + 1;
        indices[indices_count++] = index + N;


/** To create and handle the meshes I use my own library that is not released **/
Mesh* Ocean::getMesh()
    Mesh* m = new Mesh();
    m->loadGeometry(geometry, (M)*(N));
    m->loadIndices(indices, (M-1)*(N-1)*6);
    m->addBuffer(M*N, 3, sizeof(Vector3), GL_FLOAT, normales);
    return m;

void Ocean::computeHeightmap(float time)
    int nb = 0;
    for(int n=-N/2; n<N/2; n+=1)
        //printf("%d / %d\n", nb, N);
        for(int m=-M/2; m<M/2; m+=1)
            int n_prime = n + (N/2);
            int m_prime = m + (M/2);
            int index = m_prime * (N) + n_prime;
            Vector2 x(n*Lx/N, m*Lz/M);
            float height;
            Vector3 normale;
            Vector2 disp;
            evaluate(x, time, &height, &normale, &disp);
            float lambda = 0.8;
                  ( n * Lx / N ) + lambda * disp.x() 
                  ( m * Lz / M ) + lambda * disp.y() 
            geometry[index].y( height );
            normales[index] = normale;
            normales[index].y( 1 );

void Ocean::evaluate(Vector2 x, float t, 
                     float* h, Vector3* norm_, Vector2* disp_)
    Complex height(0, 0);
    Vector3 normale;
    Vector2 disp;
    for(int i=0; i<nb_htilde; i++)
        Vector2 k = k_tab[i];
        if (k.length() < 0.000001)
        float kx = DotProduct(k, x);
        Complex h = htilde(k, t, i) * Complex(cos(kx), sin(kx));
        height += h;
        normale += Vector3(-k.x() * h.imaginary(), 
                           -k.y() * h.imaginary());
        disp += Vector2(-k.normalize().x() * h.imaginary(), 
                        -k.normalize().y() * h.imaginary());
    *h = height.real();
    *norm_ = normale;
    *disp_ = disp;

double Ocean::dispersion(Vector2 k)
    float w0 = 2.0f * M_PI / 200.0f;
    return floor(sqrt(9.8 * k.length()) / w0) * w0;

Complex Ocean::htilde(Vector2 k, float t, int idx)
    int n_prime = (k.x() * Lx) / (2* M_PI);
    n_prime += N/2;
    int m_prime = (k.y() * Lz) / (2* M_PI);
    m_prime += M/2;
    int index = idx;
    float omegat = dispersion(k) * t;
    float cos_ = cos(omegat);
    float sin_ = sin(omegat);
    Complex c0(cos_,  sin_);
    Complex c1(cos_, -sin_);
    return htilde0_tab[index] * c0 + htilde0c_tab[index] * c1;

Complex Ocean::htilde0(Vector2 k)
    double Ph = Phillips(k);
    Complex c = gaussianRandom();
    Complex r = c * sqrt(Ph/2.0);
    return r;

double Ocean::Phillips(Vector2 k)
    float V = wind.length();
    float L = (V*V) / 9.8;
    float A = 0.00001;

    float kl = k.length();
   if (kl < 0.000001)
        return 0.0;
    float kw = DotProduct(k.normalize(), wind.normalize());
    double ph = exp(-1 / (kl*L*kl*L) );
    ph /= kl*kl*kl*kl;
    ph *= kw*kw;
    ph *= A;
    float l = L * 0.001;
    ph *= exp(-kl*kl*l*l);
    return ph;

[1] A. Fournier and W. T. Reeves, “A simple model of ocean waves,” in ACM Siggraph Computer Graphics, 1986, vol. 20, pp. 75–84.
[2] G. A. Mastin, P. A. Watterberg, and J. F. Mareda, “Fourier synthesis of ocean scenes,” Computer Graphics and Applications, IEEE, vol. 7, no. 3, pp. 16–23, 1987.
[3] J. Tessendorf and others, “Simulating ocean water,” Simulating Nature: Realistic and Interactive Techniques. SIGGRAPH, vol. 1, 2001.

Laisser un commentaire

Entrez vos coordonnées ci-dessous ou cliquez sur une icône pour vous connecter:

Logo WordPress.com

Vous commentez à l'aide de votre compte WordPress.com. Déconnexion / Changer )

Image Twitter

Vous commentez à l'aide de votre compte Twitter. Déconnexion / Changer )

Photo Facebook

Vous commentez à l'aide de votre compte Facebook. Déconnexion / Changer )

Photo Google+

Vous commentez à l'aide de votre compte Google+. Déconnexion / Changer )

Connexion à %s