Examples

Breakdown

The Equation

These animations were achieved using the following equation:

float h = A * sin(k*p - w*t);

A = amplitude;
k = spatial_angular_frequency; // (a.k.a. wavenumber)
p = position in space;
w = temporal_angular_frequency;
t = time;

This equation is a sinusoidal function; it represents periodic oscillations or waves. We can use the result of this function as a height to some vertex we are drawing in our world.

In a 3D world where +X is right, +Y is up, and +Z is towards the viewer, I create an xz-plane with a bunch of vertices. Then in the vertex shader, I calculate the height of each vertex (y-coordinate) based on this equation.

It’s more intuitive to parametrize a wave like so:

amplitude:  Half-height of wave, i.e., distance from crest/trough to still-water line.
wavelength: Length of one wave cycle.
speed:      How fast the wave moves.

From these parameters, we can derive:

// Wavenumber / Spatial angular frequency.
k = 2 * PI / wavelength;

// Temporal angular frequency.
w = k * speed;

// We can use it directly.
A = amplitude;

// Time is a value we increase every frame in the CPU by delta_time and we pass it to GPU in a constant buffer.
t = time;

// Position of the current vertex in some dimension. We can do cool patterns by modifying this value as I'll show later.
p = vertex.position.x and/or vertex.position.z;

With this, we can calculate the height of each vertex (of the plane) using that equation and we’d get something that resembles a sine wave moving on the surface of the plane.

The cool part is that we can combine the result of multiple sinusoidal waves where each wave has different parameters:

So, in HLSL, we have something like the following:

struct Water_Wave
{
    float amplitude;
    float spatial_angular_frequency;  // TAU / wavelength;
    float temporal_angular_frequency; // spatial_angular_frequency * speed;
};

cbuffer Constants : register(b0)
{
    float time;
    int wave_count;

    // ...
}

// We pass a bunch of random Waves from the CPU.
StructuredBuffer<Water_Wave> waves : register(t0);

float3 sine(float3 v_pos, Water_Wave wave)
{
    float k = wave.spatial_angular_frequency;
    float w = wave.temporal_angular_frequency;
    float p = v_pos.x; // -/+v_pos.x: wave moves left/right... -/+v_pos.z: wave moves forward/backward.
    float h = wave.amplitude * sin(k * p - w * time);

    // Displace along +Y axis.
    return float3(0, h, 0);
}

PS_Input vs(float3 position : POSITION)
{
    // ...

    // Position perturbation.
    float3 offset = 0.0f;
    for (int wave_index = 0; wave_index < wave_count; wave_index++) {
        offset += sine(position, waves[wave_index]);
    }

    // Final pos = original position + position perturbation. 
    float3 final_pos = position + offset;
    
    // final_pos is in local space here, we can transform it now.

    // ...
}

Different Wave Patterns

Instead of directly using the vertex position in the sine() function, we can achieve different wave patterns by interpreting the vertex position differently:

Directional

float get_wave_coord(float3 v_pos, Water_Wave wave)
{
    // Project vertex position on desired direction.
    float result = 0.0f;
    result = dot(v_pos.xz, wave.direction);
    return result;
}

float3 sine(float3 v_pos, Water_Wave wave)
{
    float k = wave.spatial_angular_frequency;
    float w = wave.temporal_angular_frequency;
    float p = get_wave_coord(v_pos, wave);
    float h = wave.amplitude * sin(k * p - w * time);

    // Displace along +Y axis.
    return float3(0, h, 0);
}

Radial

{
    // Distance from origin.
    result = length(v_pos);
}

Oscillating Radial

{
    // Oscillating wave along a circular path.
    result = sin(length(v_pos));
}

Spiral

{
    float angle  = atan2(v_pos.x, v_pos.z);
    float radius = length(v_pos);
    result       = sin(radius - angle);
}

Elliptical

{
    // Elliptical wave with different axis lengths.
    result = length(v_pos.xz / float2(2.0, 1.0f));
}

Latticed

{
    result = sin(v_pos.x * 2.0) + sin(v_pos.z * 2.0);
}

Calculating Normals

Normals are direction vectors that point directly away from a surface. We need normals for light calculations that make our waves look nice.

The way we calculate a normal of a surface at some point ‘p’ is by doing a cross product between two other vectors: the tangent and bitangent. The tangent and bitangent are vectors that “just touch” the surface at the point ‘p’. They are perpendicular to one another and doing a cross product between them gives us a third vector that is perpendicular to both (the normal):

Note: There are infinite tangent and bitangent lines for a given surface, I just have to ensure that whichever ones I choose, it must be always true that normal = cross(tangent, bitangent), noting that I use a right-handed system.

To retrieve the tangent and bitangent of a surface at some point, we need to do some calculus, so let’s take a step back:
Imagine a hill that rises and falls as you walk in the x-direction. The height value varies based on ‘x’. The rate of change of height tells us how steeply the surface rises or falls in the x-direction. This is the same thing as saying:

  • The rate of change of ‘y’ with respect to ‘x’
  • The derivative of ‘y’ with respect to ‘x’
  • dy/dx
  • The slope of the height of the surface at point ‘p’ in direction ‘x’

The direction of the tangent is along the surface of the hill.
So, the tangent at a point where the height is zero is: float3(1, 0, 0)
Otherwise the tangent will be tilted based on the rate of change of height with respect to the x-direction.

We know the function we’re using to calculate the height (assume simple sine wave propagating in x-direction):

float h = A * sin(k*pos.x - w*t);

The derivative of that is:

float dh_dx = A * k * cos(k*pos.x - w*t);

The tangent is then:

float3 tangent = normalize(float3(1, dh_dx, 0));

Now let’s do the same in 3D:

We’re dealing with the xz-plane; the normal of that plane is the y-axis/up-vector. There are infinite possibilities to choose from when it comes to the direction of the tangent and bitangent vectors. But remember, I just have to ensure that N = cross(T, B).

Let’s choose (T = z-axis) and (B = x-axis).
Now doing cross(T, B) correctly gives the y-axis (use the right-hand rule to see).

So:

float dh_dx = A * k * cos(k*pos.x - w*t);
float dh_dz = A * k * cos(k*pos.z - w*t);

float3 tangent   = normalize(float3(0, dh_dz, 1)); // Along z-axis
float3 bitangent = normalize(float3(1, dh_dz, 0)); // Along x-axis

float3 normal = cross(tangent, bitangent) = float3(-dh_dx, 1, -dh_dz);

We can omit the cross product calculation by directly doing float3 normal = float3(-dh_dx, 1.0, -dh_dz);
However, it’s super important to note that we have multiple waves, so we want to accumulate the “tilt” of the normal from all the waves:

float4 ps(PS_Input input) : SV_TARGET
{
    // ...

    // Normal perturbation.
    float3 offset  = 0.0f;
    for (int wave_index = 0; wave_index < wave_count; wave_index++) {
        offset += sine_normal(input.local_pos, waves[wave_index]);
    }
    
    // Final nor = original normal + normal perturbation.
    float3 final_nor = float3(0,1,0) + offset;

    // final_nor is in local space here, we can transform and normalize it now.

    // ... 
}

Basically, sine_normal() should return only float3(-dh_dx, 0, -dh_dz);

Different wave patterns will have different derivatives.

Directional

float3 sine_normal(float3 v_pos, Water_Wave wave)
{
    float k  = wave.spatial_angular_frequency;
    float w  = wave.temporal_angular_frequency;
    float A  = wave.amplitude;

    // This is the height rate of change in the wave's propagation direction.
    float dHeight_dr = A * k * cos(k * get_wave_coord(v_pos, wave) - w * time);

    // Directional wave pattern.
    float dHeight_dx = dHeight_dr * wave.direction.x;
    float dHeight_dz = dHeight_dr * wave.direction.y;

    return float3(-dHeight_dx, 0.0f, -dHeight_dz);
}

Radial

    // Radial waves - derivative with respect to radial distance.
    float2 grad_r = normalize(v_pos.xz);
    dHeight_dx    = dHeight_dr * grad_r.x;
    dHeight_dz    = dHeight_dr * grad_r.y;

Oscillating Radial

    float radius  = length(v_pos);
    float2 grad_r = normalize(v_pos.xz);
    dHeight_dx    = dHeight_dr * grad_r.x * cos(radius);
    dHeight_dz    = dHeight_dr * grad_r.y * cos(radius);

Spiral

    // Spiral wave - polar coordinates (angle and radius).
    float radius  = length(v_pos);
    float angle   = atan2(v_pos.x, v_pos.z);
    float2 grad_r = normalize(v_pos.xz);
    dHeight_dx    = dHeight_dr * grad_r.x * cos(radius - angle);
    dHeight_dz    = dHeight_dr * grad_r.y * cos(radius - angle);

Elliptical

    // Elliptical wave - scaled coordinates.
    float2 scaled_pos  = v_pos.xz / float2(2.0, 1.0);
    float2 grad_scaled = normalize(scaled_pos) / float2(2.0, 1.0);
    dHeight_dx         = dHeight_dr * grad_scaled.x;
    dHeight_dz         = dHeight_dr * grad_scaled.y;

Latticed

    // Lattice wave - derivatives in both x and z directions.
    dHeight_dx = dHeight_dr * 2.0 * cos(2.0 * v_pos.x);
    dHeight_dz = dHeight_dr * 2.0 * cos(2.0 * v_pos.z);

The main resource I can recommend is Acerola’s awesome videos:
https://www.youtube.com/watch?v=PH9q0HNBjT4

Note: There’s a small mistake in the video where the spatial angular frequency is shown to be 2 / wavelength when in reality it’s meant to be 2 * PI / wavelength

Wave simulation can be improved with Fractional Brownian Motion (FBM) and Fast Fourier Transforms (FFT), which is what I’m planning to look at next.

Thanks for reading ^^