# Ocean simulation part one: using the discrete Fourier transform

In this post we will implement the statistical wave model from the equations in Tessendorf's paper[1] on simulating ocean water. We will implement this model using a discrete Fourier transform. In part two we will begin with the same equations but provide a deeper analysis in order to implement our own fast Fourier transform. Using the fast Fourier transform, we will be able to achieve interactive frame rates. Below are two screen captures. The first is a rendering of the surface using some of the features from the post on lighting with GLSL with the addition of a simple fog factor. The second is a rendering of the surface mesh.

A rendering of the surface with lighting and fog applied.

A rendering of the surface mesh.

We begin with the equation for the wave height $$h(\vec{x},t)$$ at a location, $$\vec{x} = (x,z)$$, and time, $$t$$,

\begin{align*}
h(\vec{x},t) &= \sum_\vec{k} \tilde{h}(\vec{k},t)\exp(i\vec{k}\cdot\vec{x})
\end{align*}

We define $$\vec{k}=(k_x,k_z)$$ as,

\begin{align*}
k_x &= \frac{2\pi n}{L_x}\\
k_z &= \frac{2\pi m}{L_z}\\
\end{align*}

where,

\begin{align*}
-\frac{N}{2} \le &n < \frac{N}{2} \\ -\frac{M}{2} \le &m < \frac{M}{2} \\ \end{align*} In our implementation our indices, $$n'$$ and $$m'$$, will run from $$0,1,...,N-1$$ and $$0,1,...,M-1$$, respectively. So we have, \begin{align*} 0 \le &n' < N\\ 0 \le &m' < M\\ \end{align*} Thus, $$n$$ and $$m$$ can be written in terms of $$n'$$ and $$m'$$, \begin{align*} n &= n'-\frac{N}{2}\\ m &= m'-\frac{M}{2}\\ \end{align*} $$\vec{k}$$ in terms of $$n'$$ and $$m'$$, \begin{align*} k_x &= \frac{2\pi (n'-\frac{N}{2})}{L_x}\\ &= \frac{2\pi n'- \pi N}{L_x}\\ k_z &= \frac{2\pi (m'-\frac{M}{2})}{L_z}\\ &= \frac{2\pi m'- \pi M}{L_z}\\ \end{align*} We can now write our original equation for the wave height in terms of $$n'$$ and $$m'$$ as, \begin{align*} h'(x,z,t) &= \sum_{m'=0}^{M-1} \sum_{n'=0}^{N-1} \tilde{h}'(n',m',t)\exp\left(\frac{ix(2\pi n' - \pi N)}{L_x} + \frac{iz(2\pi m' - \pi M)}{L_z}\right) \end{align*} We still need to look at the function, $$\tilde{h}$$, as a function of $$\vec{k}$$ and $$t$$ and our equivalent function, $$\tilde{h}'$$, as a function of our indices, $$n'$$ and $$m'$$, and $$t$$. $$\tilde{h}(\vec{k},t)$$ is defined as, \begin{align*} \tilde{h}(\vec{k},t) = &\tilde{h}_0(\vec{k}) \exp(i\omega(k)t) +\\ &\tilde{h}_0^*(-\vec{k})\exp(-i\omega(k)t) \\ \tilde{h}'(n',m',t) = &\tilde{h}_0'(n',m') \exp(i\omega'(n',m')t) +\\ &{\tilde{h}_0'}^*(-n',-m')\exp(-i\omega'(n',m')t) \\ \end{align*} where $$\omega(k)$$ is the dispersion relation given by, \begin{align*} \omega(k) = \sqrt{gk} \\ \end{align*} and $$g$$ is the local acceleration due to gravity, nominally $$9.8\frac{m}{\sec^2}$$. In terms of our indices, \begin{align*} \omega'(n',m') = \sqrt{g\sqrt{\left(\frac{2\pi n'- \pi N}{L_x}\right)^2+\left(\frac{2\pi m'- \pi M}{L_z}\right)^2}} \\ \end{align*} We are almost there, but we need $$\tilde{h}_0'$$ and the complex conjugate, $${\tilde{h}_0'}^*$$, in terms of our indices, \begin{align*} \tilde{h}_0(\vec{k}) &= \frac{1}{\sqrt{2}}(\xi_r+i\xi{i})\sqrt{P_h(\vec{k})} \\ \tilde{h}_0'(n',m') &= \frac{1}{\sqrt{2}}(\xi_r+i\xi{i})\sqrt{P_h'(n',m')} \\ \end{align*} where $$\xi_r$$ and $$\xi_i$$ are independent gaussian random variables. $$P_h(\vec{k})$$ is the Phillips spectrum given by, \begin{align*} P_h(\vec{k}) &= A\frac{\exp(-1/(kL)^2)}{k^4}|\vec{k}\cdot \vec{w}|^2 \\ \end{align*} where $$\vec{w}$$ is the direction of the wind and $$L=V^2/g$$ is the largest possible waves arising from a continuous wind of speed $$V$$. In terms of the indices, \begin{align*} P_h'(n',m') &= A\frac{\exp(-1/(kL)^2)}{k^4}|\vec{k}\cdot \vec{w}|^2 \\ \vec{k} &= \left(\frac{2\pi n'- \pi N}{L_x}, \frac{2\pi m'- \pi M}{L_z}\right) \\ k &= |\vec{k}| \\ \end{align*} We should now have everything we need to render our wave height field, but we will evaluate a displacement vector for generating choppy waves and the normal vector for our lighting calculations at each vertex. The displacement vector is given by, \begin{align*} \vec{D}(\vec{x},t) &= \sum_{\vec{k}}-i\frac{\vec{k}}{k}\tilde{h}(\vec{k},t)\exp(i\vec{k}\cdot \vec{x}) \\ \end{align*} and the normal vector, $$\vec{N}$$, is given by, \begin{align*} \epsilon(\vec{x},t) &= \nabla h(\vec{x},t) = \sum_{\vec{k}}i\vec{k}\tilde{h}(\vec{k},t)\exp(i\vec{k}\cdot \vec{x}) \\ \vec{N}(\vec{x},t) &= (0,1,0) - \left(\epsilon_x(\vec{x},t), 0, \epsilon_z(\vec{x},t)\right) \\ &= \left(-\epsilon_x(\vec{x},t), 1, -\epsilon_z(\vec{x},t)\right) \\ \end{align*} all of which we can rewrite in terms of the indices for ease of implementation. We'll begin our implementation by looking at our vertex structure which will double to hold our values for $$\tilde{h}_0'(n',m')$$ and $${\tilde{h}_0'}^*(-n',-m')$$. We will apply our wave height and displacement vector to ox, oy, and oz to obtain the current position, x, y, and z.

struct vertex_ocean {
GLfloat   x,   y,   z; // vertex
GLfloat  nx,  ny,  nz; // normal
GLfloat   a,   b,   c; // htilde0
GLfloat  _a,  _b,  _c; // htilde0mk conjugate
GLfloat  ox,  oy,  oz; // original position
};


We have another structure to return the resulting wave height, displacement vector, and normal vector at a location, $$\vec{x}$$, and time, $$t$$,

struct complex_vector_normal {	// structure used with discrete fourier transform
complex h;		// wave height
vector2 D;		// displacement
vector3 n;		// normal
};


The following is the declaration of our cOcean object. There are references to the fast Fourier transform which we will cover in the next post. In this implementation we are assuming $$M=N$$ and $$L=L_x=L_z$$ from the equations above. We declare the property, Nplus1, for tiling purposes in order for our vertices to connect to the tile at an adjacent location.

class cOcean {
private:
bool geometry;							// flag to render geometry or surface

float g;								// gravity constant
int N, Nplus1;							// dimension -- N should be a power of 2
float A;								// phillips spectrum parameter -- affects heights of waves
vector2 w;								// wind parameter
float length;							// length parameter
complex *h_tilde,						// for fast fourier transform
*h_tilde_slopex, *h_tilde_slopez,
*h_tilde_dx, *h_tilde_dz;
cFFT *fft;								// fast fourier transform

vertex_ocean *vertices;					// vertices for vertex buffer object
unsigned int *indices;					// indicies for vertex buffer object
unsigned int indices_count;				// number of indices to render
GLuint vbo_vertices, vbo_indices;		// vertex buffer objects

GLint vertex, normal, texture, light_position, projection, view, model;	// attributes and uniforms

protected:
public:
cOcean(const int N, const float A, const vector2 w, const float length, bool geometry);
~cOcean();
void release();

float dispersion(int n_prime, int m_prime);		// deep water
float phillips(int n_prime, int m_prime);		// phillips spectrum
complex hTilde_0(int n_prime, int m_prime);
complex hTilde(float t, int n_prime, int m_prime);
complex_vector_normal h_D_and_n(vector2 x, float t);
void evaluateWaves(float t);
void evaluateWavesFFT(float t);
void render(float t, glm::vec3 light_pos, glm::mat4 Projection, glm::mat4 View, glm::mat4 Model, bool use_fft);
};


In our constructor we set up some resources, define our vertices and indices for the vertex buffer objects, create our shader program, and, lastly, generate our vertex buffer objects. The destructor frees our resources, and the release method frees our vertex buffer objects and shader program before our OpenGL context is destroyed.

cOcean::cOcean(const int N, const float A, const vector2 w, const float length, const bool geometry) :
g(9.81), geometry(geometry), N(N), Nplus1(N+1), A(A), w(w), length(length),
vertices(0), indices(0), h_tilde(0), h_tilde_slopex(0), h_tilde_slopez(0), h_tilde_dx(0), h_tilde_dz(0), fft(0)
{
h_tilde        = new complex[N*N];
h_tilde_slopex = new complex[N*N];
h_tilde_slopez = new complex[N*N];
h_tilde_dx     = new complex[N*N];
h_tilde_dz     = new complex[N*N];
fft            = new cFFT(N);
vertices       = new vertex_ocean[Nplus1*Nplus1];
indices        = new unsigned int[Nplus1*Nplus1*10];

int index;

complex htilde0, htilde0mk_conj;
for (int m_prime = 0; m_prime < Nplus1; m_prime++) {
for (int n_prime = 0; n_prime < Nplus1; n_prime++) {
index = m_prime * Nplus1 + n_prime;

htilde0        = hTilde_0( n_prime,  m_prime);
htilde0mk_conj = hTilde_0(-n_prime, -m_prime).conj();

vertices[index].a  = htilde0.a;
vertices[index].b  = htilde0.b;
vertices[index]._a = htilde0mk_conj.a;
vertices[index]._b = htilde0mk_conj.b;

vertices[index].ox = vertices[index].x =  (n_prime - N / 2.0f) * length / N;
vertices[index].oy = vertices[index].y =  0.0f;
vertices[index].oz = vertices[index].z =  (m_prime - N / 2.0f) * length / N;

vertices[index].nx = 0.0f;
vertices[index].ny = 1.0f;
vertices[index].nz = 0.0f;
}
}

indices_count = 0;
for (int m_prime = 0; m_prime < N; m_prime++) {
for (int n_prime = 0; n_prime < N; n_prime++) {
index = m_prime * Nplus1 + n_prime;

if (geometry) {
indices[indices_count++] = index;				// lines
indices[indices_count++] = index + 1;
indices[indices_count++] = index;
indices[indices_count++] = index + Nplus1;
indices[indices_count++] = index;
indices[indices_count++] = index + Nplus1 + 1;
if (n_prime == N - 1) {
indices[indices_count++] = index + 1;
indices[indices_count++] = index + Nplus1 + 1;
}
if (m_prime == N - 1) {
indices[indices_count++] = index + Nplus1;
indices[indices_count++] = index + Nplus1 + 1;
}
} else {
indices[indices_count++] = index;				// two triangles
indices[indices_count++] = index + Nplus1;
indices[indices_count++] = index + Nplus1 + 1;
indices[indices_count++] = index;
indices[indices_count++] = index + Nplus1 + 1;
indices[indices_count++] = index + 1;
}
}
}

vertex         = glGetAttribLocation(glProgram, "vertex");
normal         = glGetAttribLocation(glProgram, "normal");
texture        = glGetAttribLocation(glProgram, "texture");
light_position = glGetUniformLocation(glProgram, "light_position");
projection     = glGetUniformLocation(glProgram, "Projection");
view           = glGetUniformLocation(glProgram, "View");
model          = glGetUniformLocation(glProgram, "Model");

glGenBuffers(1, &vbo_vertices);
glBindBuffer(GL_ARRAY_BUFFER, vbo_vertices);
glBufferData(GL_ARRAY_BUFFER, sizeof(vertex_ocean)*(Nplus1)*(Nplus1), vertices, GL_DYNAMIC_DRAW);

glGenBuffers(1, &vbo_indices);
glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, vbo_indices);
glBufferData(GL_ELEMENT_ARRAY_BUFFER, indices_count*sizeof(unsigned int), indices, GL_STATIC_DRAW);
}

cOcean::~cOcean() {
if (h_tilde)		delete [] h_tilde;
if (h_tilde_slopex)	delete [] h_tilde_slopex;
if (h_tilde_slopez)	delete [] h_tilde_slopez;
if (h_tilde_dx)		delete [] h_tilde_dx;
if (h_tilde_dz)		delete [] h_tilde_dz;
if (fft)		delete fft;
if (vertices)		delete [] vertices;
if (indices)		delete [] indices;
}

void cOcean::release() {
glDeleteBuffers(1, &vbo_indices);
glDeleteBuffers(1, &vbo_vertices);
}


Below are our helper methods for evaluating the disperson relation and the Phillips spectrum in addition to $$\tilde{h}_0'$$ and $$\tilde{h}'$$.

float cOcean::dispersion(int n_prime, int m_prime) {
float w_0 = 2.0f * M_PI / 200.0f;
float kx = M_PI * (2 * n_prime - N) / length;
float kz = M_PI * (2 * m_prime - N) / length;
return floor(sqrt(g * sqrt(kx * kx + kz * kz)) / w_0) * w_0;
}

float cOcean::phillips(int n_prime, int m_prime) {
vector2 k(M_PI * (2 * n_prime - N) / length,
M_PI * (2 * m_prime - N) / length);
float k_length  = k.length();
if (k_length < 0.000001) return 0.0;

float k_length2 = k_length  * k_length;
float k_length4 = k_length2 * k_length2;

float k_dot_w   = k.unit() * w.unit();
float k_dot_w2  = k_dot_w * k_dot_w;

float w_length  = w.length();
float L         = w_length * w_length / g;
float L2        = L * L;

float damping   = 0.001;
float l2        = L2 * damping * damping;

return A * exp(-1.0f / (k_length2 * L2)) / k_length4 * k_dot_w2 * exp(-k_length2 * l2);
}

complex cOcean::hTilde_0(int n_prime, int m_prime) {
complex r = gaussianRandomVariable();
return r * sqrt(phillips(n_prime, m_prime) / 2.0f);
}

complex cOcean::hTilde(float t, int n_prime, int m_prime) {
int index = m_prime * Nplus1 + n_prime;

complex htilde0(vertices[index].a,  vertices[index].b);
complex htilde0mkconj(vertices[index]._a, vertices[index]._b);

float omegat = dispersion(n_prime, m_prime) * t;

float cos_ = cos(omegat);
float sin_ = sin(omegat);

complex c0(cos_,  sin_);
complex c1(cos_, -sin_);

complex res = htilde0 * c0 + htilde0mkconj * c1;

return htilde0 * c0 + htilde0mkconj*c1;
}


We have the method, h_D_and_n, to evaluate the wave height, displacement vector, and normal vector at a location, $$\vec{x}$$.

complex_vector_normal cOcean::h_D_and_n(vector2 x, float t) {
complex h(0.0f, 0.0f);
vector2 D(0.0f, 0.0f);
vector3 n(0.0f, 0.0f, 0.0f);

complex c, res, htilde_c;
vector2 k;
float kx, kz, k_length, k_dot_x;

for (int m_prime = 0; m_prime < N; m_prime++) {
kz = 2.0f * M_PI * (m_prime - N / 2.0f) / length;
for (int n_prime = 0; n_prime < N; n_prime++) {
kx = 2.0f * M_PI * (n_prime - N / 2.0f) / length;
k = vector2(kx, kz);

k_length = k.length();
k_dot_x = k * x;

c = complex(cos(k_dot_x), sin(k_dot_x));
htilde_c = hTilde(t, n_prime, m_prime) * c;

h = h + htilde_c;

n = n + vector3(-kx * htilde_c.b, 0.0f, -kz * htilde_c.b);

if (k_length < 0.000001) continue;
D = D + vector2(kx / k_length * htilde_c.b, kz / k_length * htilde_c.b);
}
}

n = (vector3(0.0f, 1.0f, 0.0f) - n).unit();

complex_vector_normal cvn;
cvn.h = h;
cvn.D = D;
cvn.n = n;
return cvn;
}


Below we have a method to evaluate our waves using the h_D_and_n method. Note the tiling under the conditions, $$n'=m'=0$$, $$n'=0$$, and $$m'=0$$.

void cOcean::evaluateWaves(float t) {
float lambda = -1.0;
int index;
vector2 x;
vector2 d;
complex_vector_normal h_d_and_n;
for (int m_prime = 0; m_prime < N; m_prime++) {
for (int n_prime = 0; n_prime < N; n_prime++) {
index = m_prime * Nplus1 + n_prime;

x = vector2(vertices[index].x, vertices[index].z);

h_d_and_n = h_D_and_n(x, t);

vertices[index].y = h_d_and_n.h.a;

vertices[index].x = vertices[index].ox + lambda*h_d_and_n.D.x;
vertices[index].z = vertices[index].oz + lambda*h_d_and_n.D.y;

vertices[index].nx = h_d_and_n.n.x;
vertices[index].ny = h_d_and_n.n.y;
vertices[index].nz = h_d_and_n.n.z;

if (n_prime == 0 && m_prime == 0) {
vertices[index + N + Nplus1 * N].y = h_d_and_n.h.a;

vertices[index + N + Nplus1 * N].x = vertices[index + N + Nplus1 * N].ox + lambda*h_d_and_n.D.x;
vertices[index + N + Nplus1 * N].z = vertices[index + N + Nplus1 * N].oz + lambda*h_d_and_n.D.y;

vertices[index + N + Nplus1 * N].nx = h_d_and_n.n.x;
vertices[index + N + Nplus1 * N].ny = h_d_and_n.n.y;
vertices[index + N + Nplus1 * N].nz = h_d_and_n.n.z;
}
if (n_prime == 0) {
vertices[index + N].y = h_d_and_n.h.a;

vertices[index + N].x = vertices[index + N].ox + lambda*h_d_and_n.D.x;
vertices[index + N].z = vertices[index + N].oz + lambda*h_d_and_n.D.y;

vertices[index + N].nx = h_d_and_n.n.x;
vertices[index + N].ny = h_d_and_n.n.y;
vertices[index + N].nz = h_d_and_n.n.z;
}
if (m_prime == 0) {
vertices[index + Nplus1 * N].y = h_d_and_n.h.a;

vertices[index + Nplus1 * N].x = vertices[index + Nplus1 * N].ox + lambda*h_d_and_n.D.x;
vertices[index + Nplus1 * N].z = vertices[index + Nplus1 * N].oz + lambda*h_d_and_n.D.y;

vertices[index + Nplus1 * N].nx = h_d_and_n.n.x;
vertices[index + Nplus1 * N].ny = h_d_and_n.n.y;
vertices[index + Nplus1 * N].nz = h_d_and_n.n.z;
}
}
}
}


Finally, we have our render method. We have a static variable, eval. Under the condition that we are not using the fast Fourier transform, we simply update our vertices and normals during the first frame. This will permit us to navigate the first frame of our wave height field in our main function. We then specify the uniforms and attribute offsets for our shader program and then render either the surface or the geometry.

void cOcean::render(float t, glm::vec3 light_pos, glm::mat4 Projection, glm::mat4 View, glm::mat4 Model, bool use_fft) {
static bool eval = false;
if (!use_fft && !eval) {
eval = true;
evaluateWaves(t);
} else if (use_fft) {
evaluateWavesFFT(t);
}

glUseProgram(glProgram);
glUniform3f(light_position, light_pos.x, light_pos.y, light_pos.z);
glUniformMatrix4fv(projection, 1, GL_FALSE, glm::value_ptr(Projection));
glUniformMatrix4fv(view,       1, GL_FALSE, glm::value_ptr(View));
glUniformMatrix4fv(model,      1, GL_FALSE, glm::value_ptr(Model));

glBindBuffer(GL_ARRAY_BUFFER, vbo_vertices);
glBufferSubData(GL_ARRAY_BUFFER, 0, sizeof(vertex_ocean) * Nplus1 * Nplus1, vertices);
glEnableVertexAttribArray(vertex);
glVertexAttribPointer(vertex, 3, GL_FLOAT, GL_FALSE, sizeof(vertex_ocean), 0);
glEnableVertexAttribArray(normal);
glVertexAttribPointer(normal, 3, GL_FLOAT, GL_FALSE, sizeof(vertex_ocean), (char *)NULL + 12);
glEnableVertexAttribArray(texture);
glVertexAttribPointer(texture, 3, GL_FLOAT, GL_FALSE, sizeof(vertex_ocean), (char *)NULL + 24);

glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, vbo_indices);
for (int j = 0; j < 10; j++) {
for (int i = 0; i < 10; i++) {
Model = glm::scale(glm::mat4(1.0f), glm::vec3(5.f, 5.f, 5.f));
Model = glm::translate(Model, glm::vec3(length * i, 0, length * -j));
glUniformMatrix4fv(model, 1, GL_FALSE, glm::value_ptr(Model));
glDrawElements(geometry ? GL_LINES : GL_TRIANGLES, indices_count, GL_UNSIGNED_INT, 0);
}
}
}


In the code above there are references to the complex, vector2, and vector3 objects. These are available in the download. Any library that supports complex variables in addition to 2-dimensional and 3-dimensional vectors would be appropriate. Additionally, we have used Euler's formula in our code above. Euler's formula states for any real $$x$$,

\begin{align*}
\exp(ix) &= \cos x + i\sin x \\
\end{align*}

This post is simply a brute force evaluation of our ocean surface using the discrete Fourier transform. In the following post we will provide a deeper analysis of the equations used here to implement our own fast Fourier transform algorithm. The download for this project includes a fast Fourier transform, so have a look at it and we'll analyze it in the next post.

I've added a mouse handler in a vain similar to the keyboard handler and joystick handler from previous posts. The source code expects the keyboard device node to be located at /dev/input/event0 and the device node for the mouse to be located at /dev/input/event1. You can navigate using the mouse and the WASD keys.

One last thing to mention is in regards to the fog factor. Below is the vertex shader which is reminiscent of our previous shaders with the addition of the fog_factor variable. We simply divide the z coordinate of our vertex position in view space to retrieve a fog factor in the interval $$[0,1]$$.

#version 330

in vec3 vertex;
in vec3 normal;
in vec3 texture;

uniform mat4 Projection;
uniform mat4 View;
uniform mat4 Model;
uniform vec3 light_position;

out vec3 light_vector;
out vec3 normal_vector;
out vec3 halfway_vector;
out float fog_factor;
out vec2 tex_coord;

void main() {
gl_Position = View * Model * vec4(vertex, 1.0);
fog_factor = min(-gl_Position.z/500.0, 1.0);
gl_Position = Projection * gl_Position;

vec4 v = View * Model * vec4(vertex, 1.0);
vec3 normal1 = normalize(normal);

light_vector = normalize((View * vec4(light_position, 1.0)).xyz - v.xyz);
normal_vector = (inverse(transpose(View * Model)) * vec4(normal1, 0.0)).xyz;
halfway_vector = light_vector + normalize(-v.xyz);

tex_coord = texture.xy;
}


In the fragment shader we simply apply the fog factor to the final fragment color.

#version 330

in vec3 normal_vector;
in vec3 light_vector;
in vec3 halfway_vector;
in vec2 tex_coord;
in float fog_factor;
uniform sampler2D water;
out vec4 fragColor;

void main (void) {
//fragColor = vec4(1.0, 1.0, 1.0, 1.0);

vec3 normal1         = normalize(normal_vector);
vec3 light_vector1   = normalize(light_vector);
vec3 halfway_vector1 = normalize(halfway_vector);

vec4 c = vec4(1,1,1,1);//texture(water, tex_coord);

vec4 emissive_color = vec4(1.0, 1.0, 1.0,  1.0);
vec4 ambient_color  = vec4(0.0, 0.65, 0.75, 1.0);
vec4 diffuse_color  = vec4(0.5, 0.65, 0.75, 1.0);
vec4 specular_color = vec4(1.0, 0.25, 0.0,  1.0);

float emissive_contribution = 0.00;
float ambient_contribution  = 0.30;
float diffuse_contribution  = 0.30;
float specular_contribution = 1.80;

float d = dot(normal1, light_vector1);
bool facing = d > 0.0;

fragColor = emissive_color * emissive_contribution +
ambient_color  * ambient_contribution  * c +
diffuse_color  * diffuse_contribution  * c * max(d, 0) +
(facing ?
specular_color * specular_contribution * c * max(pow(dot(normal1, halfway_vector1), 120.0), 0.0) :
vec4(0.0, 0.0, 0.0, 0.0));

fragColor = fragColor * (1.0-fog_factor) + vec4(0.25, 0.75, 0.65, 1.0) * (fog_factor);

fragColor.a = 1.0;
}


If you have any questions or comments, let me know. We will continue with this project in the next post.

References:
1. Tessendorf, Jerry. Simulating Ocean Water. In SIGGRAPH 2002 Course Notes #9 (Simulating Nature: Realistic and Interactive Techniques), ACM Press.

## 52 thoughts on “Ocean simulation part one: using the discrete Fourier transform”

1. Keith,

I am impressed with your ocean simulations.

I am looking for paid help to use/adopt your formulation/code in Unity3d.

If interested and available, please let me know.

Thanks,

DA

• Hey Dan,

That sounds like an interesting project, but I don't think I'll find the time to work on it.

Help yourself to the code, though. If you can find a place to use it, go for it.

Keith

2. Håvard says:

Your conjugate method is missing something?

Shouldn't

complex complex::conj() { return complex(this->a, this->b); }

be

complex complex::conj() { return complex(this->a, -this->b); }
?

• Hey Håvard..

Good eye.. I've updated the code for this project and the one on using the fft.

Thanks.

Where can I learn the math required to understand this? Thanks. I'm eager to try this out!

• Is there a particular section you are stuck on?

I just don't understand the equation for wave height.

• Hey.. I'm not sure what kind of background you have, but you may start with some research on the Fourier transform. Also, looking through the code and seeing how that correlates with the derivation might help. If you have a specific question, it might be easier for me to suggest a direction for you.

4. David says:

Hi

Thank you for a great tutorial.
I´m still puzzled about the tiling part with the last three "if statements"
what are they for exactly ?

• Hey David,

Because we can tile this seamlessly, we've defined the size of the grid as Nplus1. When we evaluate the wave heights, we only evaluate those heights on the N by N grid. The last row and column are reserved to equal the wave heights of the first row and column, respectively. If we were dealing only with points, we wouldn't need this extra row and column, but we want to fill that row and column with primitives, so those vertices have been added to draw the triangles. The first corner is a special case for adding the primitives to the last corner.

Since the height of the last vertices match those of the first, we can fit the second tile next to the first, and it will appear seamless.

• David says:

I see ... thank you !

5. Btrack says:

first Thanks for your post .
i compile it under windows but the code

there is no this func would you tell how can i deal with it? sorry for my poor english

• Hey..

The createProgram() function is defined in src/glhelper.cc. That source file is compiled and placed in lib/glhelper.o. This object file should be linked as the final part of the build process.

It should all be there. Do you have any output from the make process?

• Btrack says:

Thansk a lot for your replay . With the src/glhelper.cc ,the problem above is resolved。
But because I compile it in windows ,so the project miss a lot of header files 。then I install the ubuntu， and it sucessful compiled under ubuntu 。but my notebook compter don‘t have the independent video card，so when i run the main， the consequence is just a black frame 。Is it correct？

• That could be.. there isn't a lot of error checking. The shaders, src/oceanv.sh and src/oceanf.sh, were written for GLSL v3.3. The createProgram() function might be spitting out the info logs.

• Btrack says:

thanks for your replay but i still can't execute it,maybe i should buy a new computer.
now I have another question and I would appreciate if you can give me some advices.
first ,the simulation of the sea surface is composed of many small triangles ,now I number the triangles from 1 to N*N (the amount of triangle is N*N),then irradiating the sea ,if the angle of incidence of the light is very large，then some triangles will be blocked by others， now my question is how to know which triangles is blocked and output the number of the blocked triangles

• Hey..

I think the number of triangles should be $$2 \cdot N \cdot N$$. The number of quads would be $$N \cdot N$$. My first thought would be to cast a ray from each point on a triangle towards your light source and see if they intersect any other object on the way there, but that would only work if all three points intersect the same triangle on the way to the light source. You might look into shadow casting. That might give you what you're looking for.

6. Jon says:

after i compile and execute it ,it shows

error C0000: syntax error, unexpected $end at token "" (0) : error C0501: type name expected at token "" (1) : error C0000: syntax error, unexpected$undefined at token ""
(1) : error C0501: type name expected at token ""
(1) : error C1038: declaration of "E" conflicts with previous declaration at (1)

Vertex info
-----------
(0) : error C0000: syntax error, unexpected $end at token "" (0) : error C0501: type name expected at token "" Fragment info ------------- (1) : error C0000: syntax error, unexpected$undefined at token ""
(1) : error C0501: type name expected at token ""
(1) : error C1038: declaration of "E" conflicts with previous declaration at (1)

• Hey Jon..

What graphics card do you have? Have you checked if there is a latest driver for your card?

• Btrack says:

i install a new graphics card gtx650 and it still occurs error like

error C0000: syntax error, unexpected $end at token "" (0) : error C0501: type name expected at token "" (1) : error C0000: syntax error, unexpected$undefined at token ""
(1) : error C0501: type name expected at token ""
(1) : error C1038: declaration of "E" conflicts with previous declaration at (1)

Vertex info
-----------
(0) : error C0000: syntax error, unexpected $end at token "" (0) : error C0501: type name expected at token "" Fragment info ------------- (1) : error C0000: syntax error, unexpected$undefined at token ""
(1) : error C0501: type name expected at token ""
(1) : error C1038: declaration of "E" conflicts with previous declaration at (1)

and i feel like breaking down

another question,if i make it in wiondows,where should i rewrite sorry for my english
and thanks for your all replay

• Btrack says:

thanks a lot. I successfully compiled it and it runs OK in windows

7. Ferdin says:

Hi, good article!
I have a question: i adjust Length to 32;
hTilde_0 always return 0, why?
Thanks!

• Hey Ferdin..

I apologize for the delay. I've tried changing the length to 32, and it appears to evaluate properly on my side. Are you still having issues with the length?

8. Wallace says:

Fantastic article, are you able to explain more on how the displacement vector and the normal vector are derived?

Another question is on the conjugate of the value returned by the function: h'_0(-n',-m',t), shouldn't that function be without the t variable?

• Hey Wallace..

You are correct. The equation for $${\tilde{h}_0'}^*$$ should be independent of $$t$$. I've fixed this. It's been a while since I've worked with this project, so let me get back to you on your first question.

9. Sen says:

Hi Keith,

I was just wondering if the random numbers Zi and Zr are generated on every frame to calculate h(x,t)? I am doing a similar project and my waves are very jumpy at every frame!

Sen

• Hey Sen,

If I recall, the random numbers are generated once to evaluate $$\tilde{h}_0$$.

• Sen says:

Hi Keith,
Thanks so much for the reply! I must have misunderstood that part - I'll try and fix it with your advice in mind.

Thanks

P.s. I forgot to mention, great tutorial!

• Sen says:

Hi Keith,

I did mis-interpret the equation, now the random jumpy-ness has stopped! Thanks!

I still had a couple of small questions if thats ok:

1. When we work out h(x,t) = sum_over_k( h(k,t)*exp(ik.x) ), we are supposed to get a real number for the height - but the the equation for "choppiness" is almost identical to this one, only its multiplied by i*(k/k_magnitude) - surely this will yield a complex number multiplying a vector - how do we interpret this complex number for the shift in (x,z)? Or am i going wrong again!

2. what library are you using to deal with complex numbers? I will be using the "FFTW" library - which would you recommend?

Thanks again!

• Hey Sen..

1. You're right.. it will be a complex number. I just took the real part of that number and applied it to the displacement. See the h_D_and_n() where D is evaluated. It looks like, technically, I missed a multiplication by $$-1$$. However, that seems to be corrected in the evaluateWaves() method by letting lambda = -1.0.

2. I wrote a simple object to deal with complex numbers. In the second part of this project I implemented a FFT to use that complex number object. Beyond that, I'm not very familiar with libraries that implement a FFT. It looks like the FFTW library you mentioned would be a good choice.

10. Connor says:

hey i was just wondering if you could explain what a complex or send me into the right direction to researching it is I've tried googleing it but all i get are complex shapes also in the hTilde formula you have sin and cos but i cant see them in the equation can you explain why they're in there? cheers

• Hey Connor..

Here are a couple of links for you..

• Connor says:

thanks you are a huge help. didn't know why i just didn't just look on Wikipedia 😛

• Connor says:

hey i was also wondering in h_D_and_n why do you cos and sin k_dot_x? i am following your tutorial writing this all from scratch so i understand it more and so far Ive had no to little trouble. but i couldn't figure out why you did that

hope to hear from you soon 🙂

11. Andre says:

Dear Keith,

Would have a version of this simulator where you take account, for example, a wave directional spectrum and the resulting image is in gray levels?

André.

• Hey André..

I haven't done much with this simulation beyond what you see here. I'd be interested in checking out your project. Let me know if you come up with something.

12. Abraham says:

Hi, I ve the installed the drivers nVidia and visual studio 12, I can`t run any example.
How run Yuo project example?

13. Drato says:

Hey Keith,

first of all a big thanks for your tutorial on water sim.

I have a question concerning the FFT functions. Im trying to use kiss_fft instead of your functions, do you have any hints on how to do this, I have used the kiss_fft functions like this

for (int m_prime = 0; m_prime < tessendorfShader.N; m_prime++) {
}
for (int n_prime = 0; n_prime < tessendorfShader.N; n_prime++) {
}

so the offset is done beforehand(by getting the pointer &h_tilde[n_prime] or &h_tilde[m_prime * tessendorfShader.N] as the kiss functions do not provide an offset.

And another question if I wanted to map the sampled vertices to a texture when having just UVs (so going from 0.0 to 1.0) how would I map n_prime and m_prime is the following correct?

int n_prime = floor(u*(float)N+0.5);
int m_prime = floor(v*(float)N+0.5);
int index = m_prime * Nplus1 + n_prime;

Hope you can help me with that,

Best Regards!!!

Drato

14. gromit says:

Hi,

What is the difference between parameters L (L_x and L_y) and N (N, M)? You've set both to 32 in your main.h.

I usually see DFT equations where the denominator of the exponent is set to N and M.

Thanks 🙂

15. Alec says:

Hi,

thank you for this great tutorial. I'm currently trying to port your code for DFT to GPU. For this I calculated hTilde_0() and hTilde_0().conj on CPU in class constructor as you did and send the result to vertex shader, which workes fine up to this point.

In vertex shader hTilde() and h_D_and_n() are calculated. First thing is, that when I leave N to be 64, the program crashes, because of too many calculations. Setting N down works, but applying h_d_and_n.h.a to the height value of current position makes the result still look like the butterfly pattern and not at all like waves.
Am I missing some important calculation step here?

Alec

16. Matt says:

Hi, would you please tell me how can I compile this program on CodeBlocks ?
I'm Bachelor student and I just try to figure out how these program gonna be worked 😀
Thanks so much

17. Mihai Bairac says:

Hello Keith! Amazing tutorial! Now I am getting a grasp on how to use FFT in my own OpenGL ocean simulator. I have created an ocean simulator based on projected grid concept and simplex noise using tesselation shaders. Now I am on FFT. I want ot first create a version on CPU and then another one on GPU using compute shaders. I was wondering why do you use Nplus1 for the grid instead of just N? I know you are having some particular cases for tiling, still I don't really get the reson for that Nplus1. Cloud you please explain?

Thank you.

18. Mihai Bairac says:

Hello Keith! Amazing article!

I am trying to create an OpenGL ocean simulator myself. I did create one using the projected grid method, simplex noise & tesellation shaders. Now I am on FFT. I want to do it on CPU then on GPU.

I was wondering about the Nplus1 variable. Wouldnt it be enought to use only N to create a N x N grid.I know you have some particular tiling cases, still. I don't know exactly why one would need to add 1. I saw other projects where N & Nplus1 were used, but I just don't get it.

Thank you.

19. Mihai Bairac says:

Hello again,

I found the answer .... I can't believe I missed this. You used N & Nplus1, because it was easier to manuveur with these variables than N-1 and N. N-1 & N are usually present when using indexed drawing. I was accustomed to N and N-1 than Nplus1 and N. Probably it's just a matter of choise.

Thanks.

20. vineet says:

can explain me something about philips spectrum? physically what it is.

21. Rika says:

Hello Keith,

I'm not sure about whether your math is correct behind calculation h0*(-k). You basically just change it from -k to using -n' and -m'. Mathematically that is incorrect, I think...so -k = -(2 * pi* n)/ Lx and since n = n'-N/2 the result of -K would be
-k = - (2 * pi * (n' - N/2)) / Lx
-k = - (2 * pi * n' - pi*N) / Lx
-k = (- 2 * pi * n' + pi * N) / Lx
and in your code implementation it looks like this
-k = (-2 * pi * n' - pi *N)/ Lx
so in the end the sing before pi * N is different and should give a different result. I might be mistaking so please correct me if i'm wrong.

22. Pingback: Ocean Scene | kosmonaut games

23. Postrediori says:

Hi Keith!

Thanks a lot for your tutorial! This is one of the best description of Tessendorf's paper so far. I'm really impressed that after some down-porting of OpenGL to support 2.0 this model still runs smoothly on ancient PCs like uni-core Athlon and GeForce 4 MX. This implementation has a lot of potential!

Thanks)