Comparing different sampling methods for GGX IBL

Question

I have made a demo that shades a sphere with an environment map. I wanted to compare 3 different ways of sampling the environment.

1. Uniform sampling.
2. Importance sampling with the NDF as suggested in this article: Real Shading in Unreal Engine 4 (2013)
3. Importance sampling with the VNDF (distribution of visible normals) as explained in this article: Sampling the GGX Distribution of Visible Normals (2018)

The program I wrote allows me to compare all those the sampling methods.

All sampling methods seem to work and look okay to me. But they look different to each other, so I was hoping someone can help me figure out why (could be that I made some mistake or it's just how it works).

First I'm going to show a comparison of how they look. After that, I will explain how I implemented each of them. Also, I have made a small project so you try it yourselves.

The program is implemented as an iterative path-tracer, so it takes about 20 seconds for it to start to look great (specially for uniform sampling).

With low roughness, it's alright, they all look the same. In the following image it's using roughness = 0.1. On the left, uniform sampling. On the right, NDF importance sampling.

Things start to look different with higher roughness.

Roughness = 0.35

Roughness = 0.6

The one that differs the most is the one that samples according to the NDF (Unreal article). It looks brighter at the center and darker on the borders.

Uniform sampling vs importance sampling the VNDF look more similar, but still not the same. This comparison really bothers me because, in the article, they claim that the sampling method should be exact. In order to better appreciate the difference we can flip the images.

Implementation

You can find all the code in this repo. I put all the relevant code in the main.cpp.

Now I'm going to explain how I implemented each method. The code is written in GLSL.

Uniform sampling

For generating uniform samples I have the following function. It takes a couple of pseudo-random numbers between 0 and 1 (seed), and the normal of the surface ($\vec{N}$). The samples that it generates are inside the hemisphere in the direction of $\vec N$.

vec3 uniformSample(vec2 seed, vec3 N)
{
    float phi = 2 * PI * seed.x;
    float r = sqrt(1 - seed.y*seed.y);
    vec3 v = vec3(r*cos(phi), seed.y, r*sin(phi));
    v = normalize(v);
    vec3 up = vec3(0, 1, 0);
    if(dot(up, N) > 0.99)
        up = vec3(1, 0, 0);
    vec3 X = normalize(cross(N, up));
    vec3 Z = cross(X, N);
    return X * v.x + N * v.y + Z * v.z;
}

And this is how I integrate all those uniform samples applying the Cook-Torrance formula. I'm reusing the same samples for computing the diffuse term. In this case the PDF is just the constant value $\frac{1}{2 \pi}$. So we just multiply by $2\pi$.

vec3 diffuse = vec3(0, 0, 0);
// uniform sampling
{
    vec3 specular = vec3(0.0, 0.0, 0.0);
    for(uint iSample = 0u; iSample < u_numSamples; iSample++)
    {
        uint sampleId = iSample + u_numSamples * u_numFramesWithoutChanging;
        uvec3 seedUInt = pcg_uvec3_uvec3(uvec3(gl_FragCoord.x, gl_FragCoord.y, sampleId));
        vec2 seed2 = vec2(makeFloat01(seedUInt.x), makeFloat01(seedUInt.y));
        vec3 L = uniformSample(seed2, N);
        vec3 H = normalize(V + L);
        vec3 env = textureLod(u_convolutedEnv, L, 0.0).rgb;

        float NoL = max(0.0001, dot(N, L));
        float NoH = max(0.0001, dot(N, H));
        float NoH2 = NoH * NoH;
        float rough4 = u_rough2*u_rough2;
        float q = NoH2 * (rough4 - 1.0) + 1.0;
        float D = rough4 / (PI * q*q);
        float G = ggx_G_smith(NoV, NoL, u_rough2);
        vec3 fr = F * G * D / (4 * NoV * NoL);
        diffuse += env * NoL;
        specular += env * fr * NoL;
    }
    diffuse *= 2 * k_diff * u_albedo / float(u_numSamples);
    specular *= 2 * PI / float(u_numSamples);
    o_uniform = vec4(diffuse + specular, 1);
}

Importance sampling the NDF

In this case I've just copied the code in the article (Unreal 2013) and adapted it to the OpenGL coordinate system. Unlike in the previous implementation, the sampling function doesn't output the sampling direction but the half vector. We will have to compute the sampling direction $\vec{L}$ by reflecting the view direction $\vec V$.

vec3 importanceSampleGgxD(vec2 seed, float rough2, vec3 N)
{
    float phi = 2.0 * PI * seed.x;
    float cosTheta = sqrt((1.0 - seed.y) / (1 + (rough2*rough2 - 1) * seed.y));
    float sinTheta = sqrt(1.0 - cosTheta*cosTheta);
    vec3 h;
    h.x = sinTheta * cos(phi);
    h.y = cosTheta;
    h.z = sinTheta * sin(phi);
    vec3 up = abs(N.y) < 0.999 ? vec3(0, 1, 0) : vec3(1, 0, 0);
    vec3 tangentX = normalize(cross(up, N));
    vec3 tangentZ = cross(tangentX, N);
    return h.x * tangentX + h.y * N + h.z * tangentZ;
}

For the integration, since we are using importance sampling, we use a different normalization factor according to the new PDF $ \frac{D (\vec N \cdot \vec H)} {4 (\vec V \cdot \vec H)}$.

// importance sampling NDF
{
    vec3 specular = vec3(0.0, 0.0, 0.0);
    for(uint iSample = 0u; iSample < u_numSamples; iSample++)
    {
        uint sampleId = iSample + u_numSamples * u_numFramesWithoutChanging;
        uvec3 seedUInt = pcg_uvec3_uvec3(uvec3(gl_FragCoord.x, gl_FragCoord.y, sampleId));
        seedUInt.xy += seedUInt.x;
        vec2 seed2 = vec2(makeFloat01(seedUInt.x), makeFloat01(seedUInt.y));
        vec3 H = importanceSampleGgxD(seed2, u_rough2, N);
        vec3 L = reflect(-V, H);
        if(dot(N, L) > 0) {
            vec3 env = textureLod(u_convolutedEnv, L, 0.0).rgb;
            float NoL = max(0.0001, dot(N, L));
            float NoH = max(0.0001, dot(N, H));
            float G = ggx_G_smith(NoV, NoL, u_rough2);
            specular += env * F*G* dot(L, H) / (NoL * NoH);
        }
    }
    specular /= float(u_numSamples);
    o_importanceNDF = vec4(diffuse + specular, 1);
}

Importance sampling the VNDF

This one attempts to sample considering the visibility of micro-normals, that means the geometry term $G$ of the Cook-Torrance equation is also cosidered. In the article, the code they provide assumes that the surface normal is aligned with the vertical axis (Z in the paper, and Y in my code). So I had to rotate the vectors accordingly.

vec3 importanceSampleGgxVD(vec2 seed, float rough2, vec3 N, vec3 V)
{
    mat3 orthoBasis; // build an ortho normal basis with N as the up direction
    orthoBasis[1] = N;
    orthoBasis[0] = abs(dot(N, vec3(0, 0, 1))) < 0.001 ?
                        cross(N, vec3(1, 0, 0)) :
                        cross(N, vec3(0, 0, 1));
    orthoBasis[0] = normalize(orthoBasis[0]);
    orthoBasis[2] = cross(orthoBasis[0], N);

    mat3 invOrthoBasis = transpose(orthoBasis);
    V = invOrthoBasis * V;

    // stretch view
    V = normalize(vec3(rough2*V.x, V.y, rough2*V.z));
    // orthonormal basis
    vec3 T1 = (V.y < 0.999) ? normalize(cross(V, vec3(0,0,1))) : vec3(1,0,0);
    vec3 T2 = cross(T1, V);
    // sample point with polar coordinates (r, phi)
    float a = 1.0 / (1.0 + V.y);
    float r = sqrt(seed.x);
    float phi = (seed.y<a) ? seed.y/a * PI : PI + (seed.y-a)/(1.0-a) * PI;
    float P1 = r*cos(phi);
    float P2 = r*sin(phi)*((seed.y<a) ? 1.0 : V.y);
    // compute normal
    vec3 h = P1*T1 + P2*T2 + sqrt(max(0.0, 1.0 - P1*P1 - P2*P2))*V;
    // unstretch
    h = normalize(vec3(rough2*h.x, max(0.0, h.y), rough2*h.z));
    return orthoBasis * h;
}