Correct form of the GGX geometry term

Question

I'm trying to implement a microfacet BRDF in my raytracer but I'm running into some issues. A lot of the papers and articles I've read define the partial geometry term as a function of the view and half vectors: G1(v, h). However, when implementing this I got the following result:

(Bottom row is dielectric with roughness 1.0 - 0.0, Top row is metallic with roughness 1.0 - 0.0)

There's a weird highlight around the edges and a cut-off around n.l == 0. I couldn't really figure out where this comes from. I'm using Unity as a reference to check my renders so I checked their shader source to see what they use and from what I can tell their geometry term is not parametrized by the half vector at all! So I tried the same code but used to macro surface normal instead of the half vector and got the following result:

To my untrained eye this seems way closer to the desired result. But I have the feeling this is not correct? The majority of the articles I read use the half vector but not all of them. Is there a reason for this difference?

I use the following code as my geometry term:

float RayTracer::GeometryGGX(const Vector3& v, const Vector3& l, const Vector3& n, const Vector3& h, float a)
{
    return G1GGX(v, h, a) * G1GGX(l, h, a);
}

float RayTracer::G1GGX(const Vector3& v, const Vector3& h, float a)
{
    float NoV = Util::Clamp01(cml::dot(v, h));
    float a2 = a * a;

    return (2.0f * NoV) / std::max(NoV + sqrt(a2 + (1.0f - a2) * NoV * NoV), 1e-7f);
}

And for reference, this is my normal distribution function:

float RayTracer::DistributionGGX(const Vector3& n, const Vector3& h, float alpha)
{
    float alpha2 = alpha * alpha;
    float NoH = Util::Clamp01(cml::dot(n, h));
    float denom = (NoH * NoH * (alpha2 - 1.0f)) + 1.0f;
    return alpha2 / std::max((float)PI * denom * denom, 1e-7f);
}

Answer 1

TL;DR: Your $G1$ formula is wrong.

---

Just to avoid confusion, I am assuming the isotropic version of the BRDF, the Smith microfacet model (as opposed to the V-cavity model), and the GGX microfacet distribution.

According to Heitz 2014, the masking/shadowing term $G1$ is

$$ \chi^{+}\left(\omega_{v}\cdot\omega_{m}\right) \frac{2}{1+\sqrt{1+\alpha_{o}^{2}\tan^{2}\theta_{v}}} $$

and according to Walter 2007, the formula is

$$ \chi^{+}\left(\frac{\omega_{v}\cdot\omega_{g}}{\omega_{v}\cdot\omega_{m}}\right)\frac{2}{1+\sqrt{1+\alpha^{2}\tan^{2}\theta_{v}}} $$

Where $\omega_{m}$ is the microfacet normal direction (halfway vector), $\omega_{g}$ is the main (geometric) normal direction (normal), $\omega_{v}$ is the incoming or outgoing direction, $\alpha$ is the isotropic roughness parameter, and $\chi^{+}\left(a\right)$ is the positive characteristic function, or the Heaviside step function (equals one if $a>0$ and zero otherwise).

As you can notice, the half vector $\omega_{m}$ is used only to make sure the $G1$ is zero if the geometrical configuration is forbidden. More precisely, it makes sure that the back surface of the microsurface is never visible from the $\omega_{v}$ direction on the front side of the macrosurface and vice versa (the latter case is meaningful only when also refractions are supported). If the calling code guarantees this, then you can obviously omit this parameter. That is probably the reason why they did so in Unity.

Your implementation, on the other hand, uses the half vector to compute the cosine of the $\omega_{v}$ direction with respect to the microfacet, which leads to computation of something else than the presented formulae.

If it is of any help, then this is my implementation of the $G1$ factor:

float SmithMaskingFunctionGgx(
    const Vec3f &aDir,  // the direction to compute masking for (either incoming or outgoing)
    const Vec3f &aMicrofacetNormal,
    const float  aRoughnessAlpha)
{
    PG3_ASSERT_VEC3F_NORMALIZED(aDir);
    PG3_ASSERT_VEC3F_NORMALIZED(aMicrofacetNormal);
    PG3_ASSERT_FLOAT_NONNEGATIVE(aRoughnessAlpha);

    if (aMicrofacetNormal.z <= 0)
        return 0.0f;

    const float cosThetaVM = Dot(aDir, aMicrofacetNormal);
    if ((aDir.z * cosThetaVM) < 0.0f)
        return 0.0f; // up direction is below microfacet or vice versa

    const float roughnessAlphaSqr = aRoughnessAlpha * aRoughnessAlpha;
    const float tanThetaSqr = Geom::TanThetaSqr(aDir);
    const float root = std::sqrt(1.0f + roughnessAlphaSqr * tanThetaSqr);

    const float result = 2.0f / (1.0f + root);

    PG3_ASSERT_FLOAT_IN_RANGE(result, 0.0f, 1.0f);

    return result;
}