How to calculate correct normals on a deformed plane with fractal brownian motion

Question

Using OpenGL, I am vertically displacing vertices in a flat horizontal plane by the function to simulate water waves. I then calculate the normal vectors to each vertex by using the partial derivatives of this function with respect to x and z, constructing a tangent and binormal vector and taking the cross product to get the resulting normal vector. This is so I can calculate diffuse and specular lighting. This works great with a single octave of fractal brownian motion, however with multiple octaves, black spots start to appear and other visual artifacts, such as the surface splitting into strange hexagon like splodges.

Here are some examples of what I'm seeing:

Lit with 1 octave
Lit with 32 octave

Normal map with 1 octave
Normal map with 32 octaves

The code is also on GitHub and can be compiled from source very easily if you follow the instructions in the README, however the important parts of code are here.

Initial plane vertex and index generation (src/Application.cpp on github):

struct Vertex{
    float x;
    float y;
    float z;
};

// Plane vertices
std::vector<Vertex> vertices;
float size = 10.0f;
int subdivisions = 100;
float tilewidth = size / static_cast<float>(subdivisions);

for (int z = 0; z < subdivisions + 1; z++) {
    for (int x = 0; x < subdivisions + 1; x++) {
        vertices.push_back(Vertex{x * tilewidth, 0.0f, z * tilewidth});
    }
}

// Plane indices
std::vector<unsigned int> indices;
for (int z = 0; z < subdivisions; z++) {
    for (int x = 0; x < subdivisions; x++) {
        int index = (z * (subdivisions + 1)) + x;   
        indices.insert(indices.end(), {
            index,
            index + (subdivisions + 1) + 1,
            index + (subdivisions + 1),

            index, 
            index + 1,
            index + (subdivisions + 1) + 1,
        });
    }
}

Vertex shader (res/ocean_simple.shader)

#version 330 core
layout (location = 0) in vec3 aPos;

uniform mat4 model;
uniform mat4 view;
uniform mat4 projection;
uniform float time;

out vec3 normal;
out vec3 FragPos;

const int NUM_OCTAVES = 32;

#define pi 3.141592653589793

// Rotate a vec2
vec2 rotate(vec2 vec, float rot)
{
    float s = sin(rot), c = cos(rot);
    return vec2(vec.x*c-vec.y*s, vec.x*s+vec.y*c);
}

void main()
{
    float y = 0.0;
    float xz = 0.0;
    vec2 direction = vec2(cos(0.5), sin(0.5));  
    vec2 derivative = vec2(0.0, 0.0);
    float frequency = 1.0;
    float amplitude = 1.0;
    float timeMultiplier = 1.0;

    // iterate over all octaves
    for (int i = 0; i < NUM_OCTAVES; i++) {
        // (x,z) position dotted with our wave direction gives us a scalar to feed into our e^(sin(x)-1), making the wave go in that direction
        xz = dot(aPos.xz, direction);   
        // increase y height 
        y += amplitude * exp(sin(xz * frequency + time * timeMultiplier) - 1.0);    
        // increase partial derivatives with respect to x and z (i think these are correct partial derivatives?)
        derivative += amplitude * frequency * direction * cos(xz * frequency + time * timeMultiplier) * exp(sin(xz * frequency + time * timeMultiplier) - 1.0);
        // change frequency and amplitude and wave direction for next octave - fractal brownian motion effect
        amplitude *= 0.5;
        frequency *= 2.0;
        direction = rotate(direction, 0.125*pi);
    }

    // Normal vector is cross of tangent and binomrial vector
    normal = normalize(cross(vec3(1.0, derivative.x, 0.0), vec3(0.0, derivative.y, -1.0)));

    gl_Position = projection * view * model * vec4(aPos.x, y, aPos.z, 1.0);
    FragPos = vec3(model * vec4(aPos.x, y, aPos.z, 1.0));
}

Fragment Shader (res/ocean_simple.shader):

#version 330 core
out vec4 FragColor;
in vec3 normal;
in vec3 FragPos;

uniform vec3 viewPos;

vec3 sunDirection = normalize(vec3(cos(0.5), 0.5, sin(0.5)));

void main()
{
    // diffuse is dot product of normal vector and sun direction - lamberts cosine law  
    float diffuse = max(dot(normal, sunDirection), 0.0);

    // specular highlights
    vec3 viewDir = normalize(viewPos - FragPos);
    vec3 reflectDir = reflect(-sunDirection, normal);
    float spec = pow(max(dot(viewDir, reflectDir), 0.0), 32);
    
    FragColor = vec4(65, 107, 223, 225) / 255.0 * (diffuse + spec);
}

Does anyone know what could be causing the issue with my normal vectors causing these hexagon like black splodges?

Answer 1

My GLSL is extremely rusty, but I'm quite confident that your code has 2 major issues:

• Inaccuracies due to floating point arithmetic, causing the runtime result of the program to deviate way too much from the exact theoretical result.
• doing a discrete approximation of high-frequency sin/cos curves without enough vertices per period.

I'll talk about a reduced case from your vertex shader code: no time parameter, constant direction for the waves.

// NOTE: untested code (take it as pseudo-code).
const int NUM_OCTAVES = 32;

void main()
{
    float y = 0.0;
    float xz = aPos.x;
    vec2 direction = vec2(1.0, 0.0);  
    vec2 derivative = vec2(0.0, 0.0);
    float frequency = 1.0;
    float amplitude = 1.0;

    // iterate over all octaves
    for (int i = 0; i < NUM_OCTAVES; i++) {
        y += amplitude * exp(sin(xz * frequency) - 1.0);
        derivative += amplitude * frequency * direction * cos(xz * frequency) * exp(sin(xz * frequency) - 1.0);
        amplitude *= 0.5;
        frequency *= 2.0;
    }

    // Normal vector is cross of tangent and binomrial vector
    normal = normalize(cross(vec3(1.0, derivative.x, 0.0), vec3(0.0, derivative.y, -1.0)));

    gl_Position = projection * view * model * vec4(aPos.x, y, aPos.z, 1.0);
    FragPos = vec3(model * vec4(aPos.x, y, aPos.z, 1.0));
}

Floating point inaccuracy issues

Looking at the C++ code, aPos.x has values in the [0;10] interval. Let's consider what happens in the loop when aPos.x == 5.0 and i == 31:

• frequency == 2^31.
• xz * frequency == 5*2^31. It can be exactly encoded in a 32-bit float. However the next smallest value that can be encoded in a 32-bit float is 5*2^31 + 1024 !
• sin(xz * frequency). Let's put aside that some GPUs have extremely inaccurate sin/cos when the argument is far from [-2.pi;2.pi]. It doesn't make any sense to take the sin of something that has been forcefully truncated/rounded to the nearest multiple of 1024 due to floating point precision. I strongly suspect that this is a major source of inaccuracy in your shader.
• amplitude == 2^-31
• y += amplitude * exp(sin(xz * frequency) - 1.0);. The image of the i=0 harmonic is [exp(-2);1]. The image of the i=31 harmonic is [2^-31.exp(-2); 2^-31] (assuming the GPU always returns something in the interval [-1;1] for sin even when the argument is way too great). So for the y-coordinate of the position, the high-i harmonics of your shader are totally negligible compared to the low-i harmonics. Even if the i=31 harmonic is wrongly computed, the error won't be noticeable for the vertex position returned by the shader.
• derivative += amplitude * frequency * direction * cos(xz * frequency) * exp(sin(xz * frequency) - 1.0). Note that amplitude * frequency == 1. The derivative of the high-i harmonics have the same image than the i=0 harmonic's derivative. So here any inaccuracy in the high-i harmonics will cause a significant error in your normals, even if the position ends up looking fine.

Approximating sin functions with discrete points.

Let's put aside floating point issues, and go back to fundamental mathematics and its sane Real numbers. Let's consider the following function. Given a real number a>0:

Function:
f(x) = sin(a*x)
Derivative:
f'(x) = a*cos(a*x)

Your program basically tries to plot the curve of this function on screen. The vertex shader creates 100 consecutive points {x,f(x)} with 0 <= x <= 10. Then the 3d pipeline's rasterizer does a linear interpolation between consecutive points. It also needs to evaluate the derivative to get a normal for lighting computations.

Let's assume that the sequence of x is x(i) = 0.1*i, and a=10*pi:

f(x(i)) = sin(10*pi*0.1*i) = sin(pi*i) = 0

f'(x(i)) = 10*pi*cos(10*pi*0.1*i) = 10*pi*cos(i*pi)
If i % 2 == 0, f'(x(i)) = 10*pi
If i % 2 == 1, f'(x(i)) = -10*pi

So in terms of vertex position, your shader would generate a perfect plane at y=const. However it would generate normals almost tangential to the y=const plane. They would also have almost opposite direction between each consecutive vertex. This will look totally wrong on screen.

With a=10*pi, this function has a period of 2*pi/a = 0.2. It's only generating 2 vertices to represent a single sin period. This is simply not enough. I would expect the lowest requirement to be around 10 vertices per period, possibly higher.

Now back to your original shader: for the i=31 harmonic, we have a=2^31. The period of the sin is around 3.10^-9, so we have 3.10^-8 vertex/period... It just won't work.

Of course, the "fix" certainly isn't to turn the 100x100 vertex grid into a 3'000'000'000² grid. The i>10 harmonics have 0.1% of the amplitude of the i=0 harmonic. They are negligible and should be removed. How many harmonics to keep is subjective, but be sure to have enough vertices per sin period for the harmonic of highest frequency.