I have implemented a gradient noise function in GLSL in order to generate a heightmap for terrain. I am trying to get a normal vector for each pixel of the heightmap so I can calculate directional lighting from a “sun”. This normal vector will be stored in a texture along with the terrain height, so we can assume a flat surface will have a normal of vec3(0, 0, -1); (z axis).

The gradient itself is a 2D vector. I’m keeping it simple currently to one octave while I figure this out, and for that the noise and gradient for each pixel is generated like so:

vec3 random_pcg3d(uvec3 v)
{
    v = v * 1664525u + 1013904223u;
    v.x += v.y*v.z; v.y += v.z*v.x; v.z += v.x*v.y;
    v ^= v >> 16u;
    v.x += v.y*v.z; v.y += v.z*v.x; v.z += v.x*v.y;
    return vec3(v) * (1.0/float(0xffffffffu));
}

vec2 random_gradient(uvec3 p)
{
    vec3 uv = random_pcg3d(p);
    float r = sqrt(uv[0]);
    float phi = 2.0 * M_PI * uv[1];
    return vec2(r * cos(phi), r * sin(phi));
}

vec3 gradient_noise(vec2 pos, float gridSize)
{
    vec2 gridPos = pos * gridSize;
    uvec2 i = uvec2(gridPos);
    vec2 f = fract(gridPos);

    vec2 g11 = random_gradient(uvec3(i.x, i.y, 1));
    vec2 g12 = random_gradient(uvec3(i.x + 1u, i.y, 1));
    vec2 g21 = random_gradient(uvec3(i.x, i.y + 1u, 1));
    vec2 g22 = random_gradient(uvec3(i.x + 1u, i.y + 1u, 1));

    float d11 = dot(g11, f);
    float d12 = dot(g12, f - vec2(1.0, 0.0));
    float d21 = dot(g21, f - vec2(0.0, 1.0));
    float d22 = dot(g22, f - vec2(1.0, 1.0));

    f = smoothstep(0.0, 1.0, f);

    vec2 gradient_lerp_x0 = mix(g11, g12, f.x);
    vec2 gradient_lerp_x1 = mix(g21, g22, f.x);
    vec2 gradient  = mix(gradient_lerp_x0, gradient_lerp_x1, f.y);

    float q1 = mix(d11, d12, f.x);
    float q2 = mix(d21, d22, f.x);
    float noise = mix(q1, q2, f.y);

    return vec3(noise, gradient.x, gradient.y);
}

This is a tweak from a youtube tutorial that didn’t include gradient output, but I’ve found sources confirming the gradient can also be interpolated just as the noise value is.

Now from the gradient, I’m certain I should have all the information I need to calculate the surface normal in 3D space. I’ve managed to do this by sampling 2 points close to the pixel and estimating the gradient from the noise data alone which gives pretty good results down to a certain accuracy, but I would like to calculate it accurately as I already have the gradient data available.

This means I have the estimate version to compare my directly calculated results to, and I have tried many things to calculate the normal, but so far nothing I do even close to agrees with the estimate. I can’t be sure the estimate is definitely the more correct version, but it looks more promising when compared to the heightmap and gives the best visual results I’ve seen so far by a long way, so I’m assuming it might be the best I’ve managed.

Here are two attempts to calculate the 3D normal from the 2D gradient, the first of which has a few commented out attemtps to calculate the direction of the gradient in 3D space, and maybe this is the step I’m getting wrong:

vec3 calculate_surface_normal(float height, vec2 gradient)
{
    vec2 tangent_xy = normalize(gradient);
//    vec2 tangent_xy = vec2(1.0) / gradient;
    float tangent_z = length(gradient);
//    float tangent_z = 1;
//    float tangent_z = sqrt(1.0 - dot(tangent_xy, tangent_xy));

    vec3 gradient_forward = normalize(vec3(tangent_xy, tangent_z));
    vec3 gradient_right = vec3(rotate_by_phase(normalize(gradient), 0.25), 0);
    vec3 gradient_up = cross(gradient_right, gradient_forward);

    return gradient_up;
}

vec3 calculate_surface_normal(float height, vec2 gradient)
{
    // Compute the partial derivatives using finite differencing
    float dz_dx = gradient.x;
    float dz_dy = gradient.y;

    // Compute the surface normal using cross product
    vec3 tangent = normalize(vec3(1.0, 0.0, dz_dx));
    vec3 bitangent = normalize(vec3(0.0, 1.0, dz_dy));
    vec3 normal = normalize(cross(bitangent, tangent));

    // Negate z as up is simply -z
    normal.z *= -1.f;

    return normal;
}

I can compare the results by dispalying them together. These image are both showing (from left to right, then top to bottom):

• Noise value as heightmap
• Normal calculated from the gradient that was calculated while generating noise
• Normal estimated from comparing height value with two very close by regions
• Terrain coloured by height and lit using normal estimate (if above water level)

First normal calculation:

Second normal calculation:

And as it happens, with the way I’ve calculated the 3D tangent in the first function, it seems I get the same result between those functions… but nothing like the estimate method. There are similarities in the shape of details, but it looks like some of the regions are showing to point in different directions between the two methods.

Calculating the lighting using one of the gradient->normal functions produces results that don’t make sense given the detail I can see in the colour transitions, whereas the lighting in the images (calculated from the estimate) looks like it makes sense.