Wrong output while implementing GGX importance sampling

I'm receiving a non-energy conserving output while trying to importance sample the GGX Distributionof the Microfacet model, which is generally 3-4 times bigger than the resulting PDF

I'm following this blog post: Sampling microfacet BRDF | Cao Jiayin Blog

And I'm sampling the $\theta$ angle with the following formula / code:

$$\theta = acos(\sqrt{ \frac{1 - \epsilon}{\epsilon(\alpha^2 - 1) +1} })$$

float sampleTheta(float roughness) {
float random = rnd();
return acos(
sqrt(              (1.0f - random) /
((roughness * roughness - 1.0f) * random + 1.0f)
)
);
}

The GGX Distribution function i'm using is:

$$\frac{ \alpha^2 }{\pi((\alpha^2 - 1) cos^2\theta + 1)^2}$$

float MicroFacetMaterial::Distribution(vec3 normal, vec3 halfvector, float alpha) {
float NoH = dot(normal, halfvector);
float NoH2 = NoH * NoH;
float alpha2 = alpha * alpha;

float a1 = (NoH2 * (alpha2 - 1.0f) + 1.0f);

return alpha2 / (M_PI * a1 * a1);
}

And finally, the blog post derived this PDF to use along with the BRDF:

$$p_{h}(\theta) = \frac{2\alpha^2 \ cos\theta \ sin\theta}{((\alpha^2 - 1)cos^2\theta + 1 )^2 } \\ \ \\ \ \\ p(\theta) = \frac{p_{h}(\theta)}{ 4 (\omega o \cdot \omega h)}$$

float get_pdf(vec3 wo, vec3 normal, vec3 halfvector, float roughness) {
float NoH = dot(normal, halfvector);
float NoH2 = NoH * NoH;
float sin = sqrt(1.0f - NoH2);
float alpha2 = roughness * roughness;

float a1 = (NoH2 * (alpha2 - 1.0f) + 1.0f);

float pdf_h = (2.0f * alpha2 * NoH * sin) / (a1 * a1);

return pdf_h / (4.0f * dot(wo, halfvector));
}

And as a final step I'm calculating incident radiance as follow for an RGB non-recursive path tracer

float cosT = clampDot(newdirection, normal);
vec3 halfVector = HalfwayVector(-ray.direction, newdirection);

float pdf = get_pdf(-ray.direction, normal, halfVector, _roughness);

float distribution = Distribution(normal, halfVector, _roughness);
vec3 fresnel = FresnelColor(newdirection, halfVector, clampDot(halfVector, -ray.direction));
float geometry = Geometry(-ray.direction, newdirection, normal, halfVector);

float denominator = 4.0f * clampDot(normal, -ray.direction) * clampDot(normal, newdirection);

vec3 brdf = (distribution * fresnel * geometry) / (denominator);

mask = (albedo * brdf * cosT) / pdf;

The output I'm getting is a model that behaves almost like a lightsource - while debugging, I saw the Distribution term normally is up to 4 times bigger than the respective PDF, I've also decided not to post the Geometry and Fresnel term since at least from the debugging I've done those are not contributing as much as the Distribution term (as expected)

If anyone can confirm the above looks alright, I'll try to look for the error elsewhere

• In the distribution term, you should take care that whenever cos(theta) is less or equal to 0, the distribution term should return 0 (Xi operator) graphics.cornell.edu/~bjw/microfacetbsdf.pdf – Nadir Apr 4 '18 at 8:44
• somehow I missed it, thanks for helping, though this doesn't really change the unfortunate output of my path tracer – Row Rebel Apr 4 '18 at 22:57
• I have gone a little bit deeper into your code, and have notice you are not taking into account energy conservation. When sampling microfacet normals, you need 2 terms, although the second one is usualy 2*PI*random, you need the microfacet normal to compute the fresnel term for the microfacet model. Then, you would use the half vector to compute another fresnel term and use it as 1 - newFresnelTerm to ensure energy conservation. The final formulae for the brdf should look like diffuse * (1- fresnelDiffuse) + specular * brdf, and then use it in the rendering equation – Nadir Apr 7 '18 at 11:17
• To be clear, the material I'm simulating has the diffuse component set to 0, I'm currently only simulating the specular part of it I thought the above sampling code would sample the incoming direction, and that the halfvector between wo and wi is the microfacet normal, is this correct? If not I may have to re-evaluate what I wrote – Row Rebel Apr 7 '18 at 14:32
• Yes, that was meant for the reflected ray. I have review my microfacet implementation to check for comparsion. I have noticed that when computing the final microfacet brdf, you should check whether the denominator is 0, in which case the brdf should return 0, otherwise you could run into a indetermination of the type 0 / 0. – Nadir Apr 7 '18 at 15:15