Sine, Cosine and the Unit Circle#

We very often use sine and cosine functions in shaders to be able to pass from a linear set of values to a periodic one (in $sin(t)$, both $t_0 = x$ and $t_1= x+ 2\pi$ will yield the same result) but how are they actually defined?
For this we need to take the unit circle, that is a circle in the cartesian plane whose radius is equal to one.

You can mathematically describe this circle by utilising Pythagoras’ Theorem and have that: $x^2 + y^2 = 1$.
You’ll have that you’ll be able to solve this equation only for values of x and y that place a point on the red line. Thus at any given angle $\alpha$ between the x axis and an arbitrary line going through the center, there will be a corresponding pair of x and y values that will place a point on the unit circle. That’s where sine and cosine come in, where they are the function that returns, respectively, the y coordinate and the x coordinate of the point for any given angle.

$cos(\alpha)^2 + sin(\alpha)^2 = x^2 + y^2 = 1$

Derivatives and Partial Derivatives#

Derivatives are a foundational block to be able to go from a function that only computes displacement of a mesh to computing the normal vector at any point on the mesh, thus unlocking the possibility to shade a quad with much more detail beyond the one available to us at a vertex level. Derivatives inform us on the difference between two points in a function. If you take any two points $x_0$ and $x_1$ for the function $f(x)$, the change between the two points is defined as:

$\Delta f(x_0,x_1) = f(x_1) - f(x_0)$

The derivative is what we get when we try to get two points infinitesimally close together.

$\frac{\delta f(x)}{\delta x} = f'(x) = \underset{h\to 0}{lim}\ \Delta f(x, x + h)$

For a constant value this would be 0, because whatever point you choose its value is always the same, for a linear function of the type $f(x) = k + mx$ it would be $m$, since the constant value $k$ yields a 0 change over two points while $m$ is the slope with which the values of the function change. Here is a quick recap of useful derivatives.

Partial derivatives are instead derivatives when it comes to multiple dimensions. Since the derivatives inform us about the rate of change of a function in only one dimension, partial derivatives allow us to isolate one dimension at a time. They can be computed by assuming all other dimensions remain constant (it’s as if you’re taking a very small slice of the function where this is true)
e.g. with $f(x,y,z) = Axy + Byz + Czx$, I won’t be rigorous and use $K$ to indicate that component is fully constant for that dimension.

Reigning in Complexity: Fractional Brownian Motion#

Being able to select a full list of values for the wave function definitely opens up a big space for customization of the look of the final wave, but authoring 10+ parameters and fine tuning them definitely becomes overwhelming quickly, but we still need enough waves to get good detail. We can take into consideration energy conservation from before to leverage a mathematical concept to add constraints on the way we generate the waves, without forgoing the ability to modify the final look.
Fractional Brownian Motion is the generalization of Brownian Motion, which, in few words, is the concept of a sequence of random values, where each one is not correlated to the next one. In the case of FBM, instead, we do allow for some correlation, as we will be seeing now.

If we redefine $a_i, \omega _i, \rho _i$ in the following way, we can much more easily scale the number of waves without having to handle $3n$ parameters each time.

So we ultimately define a sum of wave functions where at each step, starting from the base values of $A,W,P$ for ampltidue, frequency and phase:

• The amplitude decreases by a factor of $R_a$
• The phase and the frequency increase by a factor of $R_\rho$ and $R_\omega$

We will call the factors ramps from now on, as they ultimately describe how steep the values ramp down or up.

So ultimately now we have at any point a pseudo random number where each value still correlates with the others, given the parameters are all computed from the same starting values. The full mathematical definition of FBM is more complex and precise while in our case it indicates just the emulation of random behavior by layering multiple functions together with constants that are correlated across each layer.

So now from having $3n$ parameters for our function, we only have 6 (starting + ramp) and they still allow for a good amount of customization.
Here is a quick hlsl snippet summarizing the current wave generation process. And here is an example shader showing the result of the sum of sines with 32 iterations. The colour of the line also describes the slope of the normal line over the length of the wave.

1
#define WaveFunc sin
2
#define WaveFuncDerivative cos
3
#define WAVE_NUM 32
4

5
float GetWaveHeight(float x, float time, float A, float W, float P, float amplitudeRamp, float frequencyRamp, float phaseRamp){
6
    float h = 0;
7
    for (int i=0; i< WAVE_NUM; i++){
8
        h += A * pow(amplitudeRamp, i) * WaveFunc(
9
            W * pow(frequencyRamp,i) * x +
10
            P * pow(phaseRamp, i) * time);
11
    }
12
    return h;
13
}
14

15
float GetWaveNormalLine(float x, float time, float A, float W, float P, float amplitudeRamp, float frequencyRamp, float phaseRamp){
16
    float m = 0;
17
    for (int i=0; i< WAVE_NUM; i++){
18
        m += A * pow(amplitudeRamp, i) * W * pow(frequencyRamp,i) *
19
            WaveFuncDerivative(
20
                W * pow(frequencyRamp,i) * x +
21
                P * pow(phaseRamp, i) * time);
22
    }
23
    // This is the slope of the line normal to the function,
24
    // described as y = (-1/m)*x. If we wanted the direction
25
    // of the normal vector we could just get it as float2(-m, 1)
26
    return -1/m;
27
}

Addressing Wave Peaks#

One other thing missing from the wave function in order to give it a more natural look is peaks. Beyond actual breaking waves (which are much more complex phenomenon that are closer to the realm of full fluid simulation), peaks are one of the characteristic features of waves. In order to have this, instead of playing around with the way we generate the parameters such as amplitude, phase etc… we can change the actual kernel of the sum. We are currently using the sine function, but we do have options that have a more peaky appearance: Gerstner Waves and Exponential Sines. This will change our functions a bit.

Gerstner Waves operate a bit differently, utilizing the input value to generate both coordinates in the x and y dimension, the function for the y dimension stays the same, utilizing $x_0$ as input of the sine function, while x becomes the following, wehre $q_i$ is another factor accounting for the steepness of the wave:

Here is the new functions when using the exponential sine instead:

For simplicity I’ll focus on the exponential sine, as Gerstner Waves are notoriously prone to easily breaking when parameters go out of range, thus needing more specific curation. Here is the same shader as before, but utilizing the exponential sine as kernel function instead.

Directing the waves#

In two dimensions, the waves can go into one of two directions: right or left, but this behavior scales quickly when going 3D, as now the direction of a wave can be any of the 360 degrees over the surface of the horizontal plane.

The way we can model this behavior is simpler than it looks. Looking at it from the top, you can imagine the sine wave as a function that will give values to $z$ based on some contribution from the $x$ and $y$ variables. We thus define factors $k_x,k_y$ which represent the contribution of $x$ and $y$ to the sine function computing $z$. We also see that $k_x$ and $k_y$ can describe a vector that represents the direction of the wave, rewriting $H(x,y)$ in a shorter form.

Let’s also rewrite this in a $y$-up system, given the final HLSL code will be written for Unity, which uses a $y$-up, left handed coordinate system. We also rewrite the derivative for each dimension, which we will need to compute the partial derivatives of the function.

This applies to any of the previous kernel functions, so now we can rewrite the previous 2D FBM algorithm in its three dimensional form, keeping in mind that we are going to randomize $D$ for each iteration of the FBM, yielding a different $Di$ at each step. This will allow us to sum together wave signals developing in different directions, replicating the effect of waves being generated from different sources (wind, previous waves bouncing off the shore etc…)

Computing the normal#

In the 2D case of wave generation we were able to obtain the line that is normal to any point of the function by first computing the slope of the tangent line (i.e. the derivative of the function at the point) and then finding the line perpendicular to it. We’ll operate similarly to find the normal vector at any point, starting from the partial derivatives over the two horizontal directions $x,z$. With the derivatives we’ll compute the two vectors $\vec{T},\vec{B}$, respectively the vector tangent and binormal of the displaced surface. Given the sine wave now develops over two dimensions, we’ll need two lines tangent to the function, instead of one.

The best way to visualize what is being done is to have a dividi et impera approach. We take a slice of the space so thin over the x and z axis,that the other variable will practically be constant and we can go back to considering the function a 2D function. (Thin slice over x for the z derivative, the opposize for the x derivative). We then compute the derivative for that faux 2D function, thus getting a vector that is necessarily tangent to it. Since the two vectors lie on perpendicular planes ($\vec{T}$ on $x$ plane, $\vec{B}$ on the $z$ plane), they are also perpendicular to each other.

From here, getting the normal vector is just a matter of computing the cross product between $B$ and $T$, as it’s by definition perpendicular to both.

So we now have:

• A wave-line signal we can use to displace a vertex along the vertical axis.
• A derivative of the signal in both horizontal dimensions, allowing us to compute the normal to the surface at any point.

Here is an hlsl snippet summarizing what we can compute so far.

1
#define WaveFunc ExpSine
2
#define WaveFuncDerivative ExpSineDer
3
#define WAVE_NUM 32
4

5
float ExpSine(float x){
6
    return exp(sin(x)-1.0);
7
}
8
float ExpSineDer(float x){
9
    return exp(sin(x)-1.0) * cos(x);
10
}
11

12

13
float Get3DWaveHeight(float x, float z, float2 waveDir, float time, float A, float W, float P, float amplitudeRamp, float frequencyRamp, float phaseRamp){
14
    float h = 0;
15
    for (int i=0; i< WAVE_NUM; i++){
16
        h += A * pow(amplitudeRamp, i) * WaveFunc(
17
            W * pow(frequencyRamp,i) * (waveDir.x * x + waveDir.y * z) +
18
            P * pow(phaseRamp, i) * time);
19
    }
20
    return h;
21
}
22

23
float Get3DWaveNormal(float x, float z, float2 waveDir, float time, float A, float W, float P, float amplitudeRamp, float frequencyRamp, float phaseRamp){
24
    float mx = 0;
25
    float mz = 0;
26
    for (int i=0; i< WAVE_NUM; i++){
27
        mx += A * pow(amplitudeRamp, i) * W * pow(frequencyRamp,i) * waveDir.x *
28
            WaveFuncDerivative(
29
                W * pow(frequencyRamp,i) * (waveDir.x * x + waveDir.y * z) +
30
                P * pow(phaseRamp, i) * time);
31
        mz += A * pow(amplitudeRamp, i) * W * pow(frequencyRamp,i) * waveDir.z *
32
            WaveFuncDerivative(
33
                W * pow(frequencyRamp,i) * (waveDir.x * x + waveDir.y * z) +
34
                P * pow(phaseRamp, i) * time);
35
    }
36

37
    return normalize(-mx,1,-mz);
38
}

Cubemap Reflections#

Utilising simply the ambient colour when not reflecting the sun doesn’t actually give off the greatest visual effect, as in reality you rarely have a completely solid colour being reflected, due to atmospherical effects, clouds in the sky, objects in view etc… A cubemap greatly improves the visual quality of reflections in general, and here it is no exception. We already have the vector used to sample it in the reflected vector $\vec{R}$, so it will just be a matter of using the sampled cubemap as reflected colour.

1
    float3 FinalColor(
2
            ...
3
            float3 ambientColor
4
            samplerCUBE cubeMap
5
            ...){
6
        ...
7
        float3 R = reflect(V,N);
8
        float3 ambientColor = texCUBE(cubeMap, R);
9
        ...;
10
}

Improving Tessellation with Signal Processing#

Let’s start from the observation that there is a threshold (given by the amplitude of each fbm step) below which any new added wave won’t displace the surface enough to require higher vertex density. This will depend on some factors such as the starting amplitude and the ramp value for it, so it can change depending on the parameters. But we can generally say that there is a maximum frequency $\omega _i$ that will displace the wave enough to need tessellation. Starting from this $\omega _i$, we will need to “cut” a segment of length 1 enough times to accomodate for the displacements with this frequency. This is essentially finding the minimum viable sampling rate for a signal. We can compute it by using the theorem of Nyquist-Shannon.

We can then generalize this for any segment of a tri $\overline{AB}$.

So ultimately as long as our edge tessellation factor cuts the tri edge more than $N_{samples}$ times, we’re going to capture the wave signal with enough precision. This precision will still scale based on other factors and use cases.

There are also different cases where we don’t actually need higher vertex density, and can get away with just manipulating the normal at fragment stage.

• At a shallow angle and long distances, what matters is the outline of the surface, so we need to preserve visibility of high amplitude displacements.
• At closer distances, we want to have enough details to see the bigger waves actually cover what should be behind them.
• At closer distances we do also need enough “vertex resolution” to avoid artifacts given by the fact that the fragment will compute shading with higher precision.

We can thus reuse the Fresnel computation from before in order to make sure the tessellation is stronger when looking at the geometry at an angle and use distance factors that can be authored via the UI to manipulate the minimum and maximum range of the tessellation factor, as well as restrain the range of values for the fresnel effect on tessellation, as it might be undesirable in some cases (e.g. case 2 with a very close quad). This will allow us to have different configurations more prone to tessellation at the edges vs closer tessellation or a mix of both. Putting it all together yields the following computation to establish the edge tessellation factor:

Where:

• the normal $N$ and $V$ are computed at the midpoint between $p_1$ and $p_2$
• the $Max$ and $Min$ values are all UI controlled parameters.
• S is a UI controlled multiplier

As for the internal tessellation factor, we will compute an average across the three edges and scale it in order to have it grow together with the edge factor. Here is the HLSL code resulting from putting it all together.

1
    // p1 and p2 are in world space.
2
    float GetEdgeTessellationFactor(float3 p1, float3 p2, float freq, float3 Nmidpoint, float3 Vmidpoint, float2 fresnelTessRange, float2 distTessRange, float S){
3
        float3 pm = (p1 + p2) / 2;
4
        float sf = clamp(pow(1-dot(Nmidpoint, normalize(Vmidpoint)),5), fresnelTessRange.x, fresneltTessRange.y);
5
        float df = 1 - saturate( (length(viewDir) - distTessRange.x)/(distTessRange.y - distTessRange.x));
6
        float tessFactor = 2 * freq * length(p1-p2);
7
        tessFactor *= sf;
8
        tessFactor *= df;
9
        tessFactor *= S;
10
        return max(tessFactor,1);
11
    }
12

13
    float4 GetTessellationFactors(float3 p1,float3 p2,float3 p3, ...){
14
        float3 edgeFactors = float3(
15
            GetEdgeTessellationFactor(p1, p2, ...),
16
            GetEdgeTessellationFactor(p2, p3, ...),
17
            GetEdgeTessellationFactor(p3, p1, ...),
18
        );
19
        return float4(edgeFactors,
20
            (edgeFactors.x + edgeFactors.y + edgeFactors.z) * 0.2); // Here dividing by 5 yielded a better scaling of the internal factor compared to just the average.
21
    }

Here is the final tessellation result.