3D Ocean Shader Using Gerstner Waves

When I started making my own Pirate Game, I needed to create a good ocean for my game since 90% of the time, we will be fighting, sailing & looking at the ocean. After searching for quite some time, I came across Gerstner waves. Some other techniques includes Fast Fourier Transform (FFT), a simple sine-wave based ocean and so on. However, using Gerstner-wave equation shader is the balance between visuals & performance (for me).

To start, we need to understand the Gerstner wave equation. Later, we’ll apply it to our plain mesh & compute the vertex displacement for every vertex. Later, we’ll come towards fragment shader in which normal computations, foam, refraction & other effects will be created.

What is a Gerstner Wave?

I am assuming that you are aware of simple sine wave. Sine wave displaces vertices up and down. Gerstner wave, in addition to this sine-based y-axis displacement, also does horizontal-axes displacement of vertices using same sine-wave. Imagine a surface made up of many points. Suppose we apply sine function like this:

y = sin(x)

What happens? Each point at ‘x’ displaces in y axis based on its x value die to sine function. Now imagine if you also add another component to this motion in x-axis, merely another sine-wave. It will cause combined movement in y (up & down) and in x (left & right). If you see both motion in sync, you’ll see that it looks something like the below image (see green line):

trochoidal wave graph

Here the x-component is displaced by:

x = x + cos(x)

Why it is “x + cos(x)” and not directly “cos(x)”? Because we first set them to their real position which is x & then add cos(x) to it to add motion on top of its original position. But in y-direction, all of our vertices are at equal height for a plane, so doing y = y + sin(x) will not cause any effect (they are at y=0 already hehe).

So at this point our overall equation for a 2D side-view Trochoidal wave is:

wave(x, y, time) = (x + cos(x + time), y + sin(x + time))

But note that it works only in xy plane (works well for side-scrolling games). How to make it work in 3D environment? Since in 3D environment, we have another axis in horizontal along with x: the z-axis. If we keep making waves that only take x in input, the wave will not look much different than it did in 2D place. So we need to take into account both xz. Also, we have to find a way to define direction of wave in xz plane.

Gerstner wave in 3D

The new z axis is very similar to the x, but we need to determine the direction of our wave in xz. Direction can be determined by a 2D xz-vector. So now we simply multiply the x component of direction with x part of our equation & z component with z part of the equation & we also take dot product of direction, with, actual xz-position of a vertex in order to determine the displacement produced by sine to be in accordance with the direction. Thus our equation becomes:

wave(x, y, z, time) = (
x + direction.x * cos(dot(direction, vec2(x, z)) + time),
y + sin(dot(direction, vec2(x, z) + time),
z + direction.z * cos(dot(direction, vec2(x, z)) + time)
)

Here, the bigger the value of x-component of direction-vector, the more the wave goes in that x-direction and vice versa. It should be in range (-1,1) for descent output. Notice that dot product between direction & vec2(x,z) simply means that how much the direction is towards any of these axis. This is how it looks so far:

gerstner wave on 3d mesh without normals

Our code so far:

// direction.y represents z axis (since it is 2D vec)
vec3 gerstner(vec3 vertex, vec2 direction, float time){
	float displaced_x = vertex.x + direction.x * cos(dot(direction, vertex.xz) + time);
	float displaced_z = vertex.z + direction.y * cos(dot(direction, vertex.xz) + time);
	float displaced_y = vertex.y + sin(dot(direction, vertex.xz) + time);
	return vec3(displaced_x, displaced_y, displaced_z);
}

So far, we haven’t considered amplitude, steepness & speed. But it is simple to add them to our above equation. The trochoidal wave equation we have achieved above is the complete template for us to add some more properties on top of it to make it artist-tweakable. Also, keep in mind that we need to calculate normal vector for this function as well. And we’ll do it after we have completed the displacement function.

Now we will add properties such as amplitude, steepness, speed, wavelength to our wave function.

1. Wave speed: It is a simple float value that you can multiply with time in our equation.

2. Steepness: It is how sharp our wave can be. This must be between 0 & 1 otherwise weird behavior occurs. Multiply a simple float value outside of cosine in above equation. This works well as long as wavelength is not introduced in our equation. After we add wavelength, you’ll see weird behaviour when wavelength is not sync with steepness (and loops start to appear). To cope with this, instead of directly multiplying with steepness, multiply with (steepness/wavelength).

2. Amplitude: It is how high our wave can reach. In steepness above, we did not multiply it with y-axis calculation. So amplitude is only multiplied with y-axis (in my implementation here).

So our wave function is fully complete, I am writing complete vertex displacement code here. After that, we’ll proceed to calculating normals for it:

// direction.y represents z axis (since it is 2D vec)
vec3 gerstner(vec3 vertex, vec2 direction, float time, float speed, float steepness, float amplitude, float wavelength){
	float displaced_x = vertex.x + (steepness/wavelength) * direction.x * cos(wavelength * dot(direction, vertex.xz) + speed * time);
	float displaced_z = vertex.z + (steepness/wavelength) * direction.y * cos(wavelength * dot(direction, vertex.xz) + speed * time);
	float displaced_y = vertex.y + amplitude * sin(wavelength * dot(direction, vertex.xz) + speed * time);
	return vec3(displaced_x, displaced_y, displaced_z);
}

Calculating Normals

Normals are vectors that point away from a surface. In graphics, they are used for light calculations. For a surface deformed by a function (such as our above wave function), the normals, too will be calculated based on it.

For a simple visualization, think of a 2D wave. Now find a line tangent to the wave at some point, by taking the difference of function at x & x + a small offset from it. Then calculate normal vector simply as a line perpendicular to this tangent.

However, for a 3D plane, a single tangent vector cannot represent slope of surface for both of the directions (x & z). So we need two tangent vectors for each of xz directions. We call the other one birnormal vector. Then we calculate the cross product of tangent & binormal to get our normal vector.

For now, I have used a simplified version of normal vector equation & I am posting code for complete implementation of normal vector for our wave function:

vec3 gerstner_normal(vec3 vertex, vec2 direction, float time, float speed, float steepness, float amplitude, float wavelength) {
    float cosfactor = cos(wavelength * dot(direction, vertex.xz + speed * time));
	float sinfactor = sin(wavelength * dot(direction, vertex.xz + speed * time));
	float x_normal = -direction.x * wavelength * amplitude * cosfactor;
	float z_normal = -direction.y * wavelength * amplitude * cosfactor;
	float y_normal = 1.0 - (steepness/wavelength) * wavelength * amplitude * sinfactor;
	return vec3(x_normal, y_normal, z_normal);
}

At this stage, our ocean looks like this:

ocean shader using gerstner wave, with normals
Single wave, with a basic normal map, metallic, etc applied

For a more realistic ocean, you’ll need to make 2 or 3 such waves with varying properties & then add them together. The resulting ocean will look much better. I am concluding this post here.

The complete code (Godot engine)

shader_type spatial;

uniform vec2 testdirection;
uniform float speeed;
uniform float steeepness : hint_range(0.0, 1.0, 0.05) = 1.0;
uniform float aamplitude = 1.0;
uniform float waavelength = 1.0;

// direction.y represents z axis (since it is 2D vec)
vec3 gerstner(vec3 vertex, vec2 direction, float time, float speed, float steepness, float amplitude, float wavelength){
	float displaced_x = vertex.x + (steepness/wavelength) * direction.x * cos(wavelength * dot(direction, vertex.xz) + speed * time);
	float displaced_z = vertex.z + (steepness/wavelength) * direction.y * cos(wavelength * dot(direction, vertex.xz) + speed * time);
	float displaced_y = vertex.y + amplitude * sin(wavelength * dot(direction, vertex.xz) + speed * time);
	return vec3(displaced_x, displaced_y, displaced_z);
}

vec3 gerstner_normal(vec3 vertex, vec2 direction, float time, float speed, float steepness, float amplitude, float wavelength) {
    float cosfactor = cos(wavelength * dot(direction, vertex.xz + speed * time));
	float sinfactor = sin(wavelength * dot(direction, vertex.xz + speed * time));
	float x_normal = -direction.x * wavelength * amplitude * cosfactor;
	float z_normal = -direction.y * wavelength * amplitude * cosfactor;
	float y_normal = 1.0 - (steepness/wavelength) * wavelength * amplitude * sinfactor;
	return vec3(x_normal, y_normal, z_normal);
}

void vertex(){
	VERTEX = gerstner(VERTEX, normalize(testdirection), TIME, speeed, steeepness, aamplitude, waavelength);
	NORMAL = gerstner_normal(VERTEX, normalize(testdirection), TIME, speeed, steeepness, aamplitude, waavelength);
}

To add more

Water is not done yet. We need to add refraction for objects below the surface & underwater fog to make objects fade with depth. Also, we need to add foam around object boundaries. However, for each of these, I have written a separate post that you can read one by one.

Continue reading 3D refraction shader for our water surface.

More Reading

  1. https://developer.nvidia.com/gpugems/gpugems/part-i-natural-effects/chapter-1-effective-water-simulation-physical-models
  2. https://catlikecoding.com/unity/tutorials/flow/waves/
  3. https://fire-face.com/water/