Cook-Torrance shader cuts off really weird when NdotL <= 0

Question

So i've been trying to implement the Cook-Torrance shader model in a toy project I'm working on and it looks quite good when looking at the right angle: But when you're looking from a shallow angle the visuals you get bright artifacts and around the cutoff.

The cutoff happens because I'm checking if NdotL > 0, but if I remove it things start to become even stranger: Colors get inverted, some kind of line emerges where NdotL == 0 and every fragment where NdotH < 0 becomes black, making it have an egg shape.

Here is the shader code:

#version 330 core
in vec3 Normal;
in vec3 FragPos;
in vec2 TexCoord;
in vec3 camPos;
in vec3 lightDir;

out vec4 color;

uniform sampler2D diffuseTexture;
uniform sampler2D glossTexture;
uniform sampler2D metalTexture;
uniform samplerCube cubemapTexture;
uniform vec3 lightPos;

float F(float ior, vec3 view, vec3 halfV) {
    float F0 = abs((1.0 - ior) / (1.0 + ior));
    F0 *= F0;
    float VoH = dot(view,halfV);
    float fresnel = F0+(1-F0) * pow(1 - VoH,5);
    return fresnel;
}

float chiGGX(float v) {
    return v > 0 ? 1 : 0;
}
float G(vec3 view, vec3 norm, vec3 halfV, float alpha) {
    float VoH2 = clamp(dot(view,halfV),0.0,1.0);
    float chi = chiGGX( VoH2 / clamp(dot(view,norm),0.0,1.0));
    VoH2 = VoH2 * VoH2;
    float tan2 = (1-VoH2) / VoH2;
    return (chi*2)/(1+sqrt(1+alpha*alpha*tan2));

}
float D(float roughness, vec3 norm, vec3 halfV) {
    float NdotH = max(dot(norm, halfV), 0.0);
    float r1 = 1.0 / ( 4.0 * roughness * roughness * pow(NdotH, 4.0));
    float r2 = (NdotH * NdotH - 1.0) / (roughness * roughness * NdotH * NdotH);
    return r1 * exp(r2);
}
void main()
{
    float gamma = 2.2f;
    float roughnessValue = texture(glossTexture, TexCoord).r;
    vec3 lightColor = vec3(1.0f, 0.8f, 1.0f)*4.0;

    vec3 norm = normalize(Normal);

    vec3 viewDir = normalize(camPos-FragPos);
    vec3 halfVector = normalize(lightDir + viewDir);
    float diff = max(dot(norm, lightDir), 0.0);

    float NdotL = dot(norm, lightDir);

    float spec = 0;
    if(NdotL > 0.0) {
        spec = ( F(1.45, viewDir, halfVector) * G(viewDir,norm,halfVector,roughnessValue) * D(roughnessValue,norm,halfVector)) / (3.14151828 * dot(norm, viewDir) * dot(norm, lightDir));
    }
    vec3 specular = spec * lightColor;
    vec3 ambient = vec3(0.05);
    vec3 diffuse = (1 - texture(metalTexture, TexCoord).r) * diff * lightColor + ambient;
    vec3 finalcolor = (diffuse * pow(texture(diffuseTexture, TexCoord).rgb, vec3(gamma))) + specular;

    color = vec4(finalcolor, 1.0f);
    color.rgb = pow(color.rgb, vec3(1.0/gamma));

}

I know there are some unused values but that's because the shader isn't completed yet.

Answer 1

Here is my independent implementation of Cook-Torrance based on Beckmann distribution:

layout(location = 0) in PerVertex
{
    special3 pos; // tangent to view
    vec2 texcoord;
    vec4 diffuse;
} IN;

layout(location = 0) out vec4 OUT;

layout(binding = 0) uniform sampler2D u_bump;
layout(binding = 1) uniform sampler2D u_raughness;

struct PerLight
{
    vec3 position;
    vec3 color;
};

layout(binding = 1) uniform lights_block
{
    int nlights;
    PerLight lights[4];
} LIGHTS;



float D(float m, float c) {
    float c2 = c*c, m2 = m*m, c2m2 = c2*m2;
    return exp((c2 - 1)/c2m2)/(3.14159*c2m2*c2);
}
float F(float R0, float NV) { return R0 + (1 - R0)*pow(1 - NV, 5); }
float G(float NL, float NV, float NH, float HV) { return min(1, 2*NH*min(NV, NL)/HV); }

void accumulate_light(special3 tangent, vec3 viewDir, float roughness, PerLight light, inout vec3 diffuse, inout vec3 specular)
{
    vec3 lightDir = quat_apply(tangent.q, light.position);

    if(lightDir.z > 0)
    {
        lightDir = normalize(lightDir);
        float NL = lightDir.z;
        diffuse += NL * light.color;

        float NV = viewDir.z;
        if(NV > 0)
        {
            vec3 halfDir = normalize(lightDir + viewDir);
            float NH = halfDir.z;
            float HV = dot(halfDir, viewDir);
            specular += D(roughness, NH)*F(0.034, NV)*G(NL, NV, NH, HV)/(4*NV*NL) * light.color;
        }
    }
}

void main()
{
    special3 tangent = {
        vec3(0),
        texture(u_bump, IN.texcoord.xy)
    };

    tangent = special_mul(IN.pos, tangent);
    tangent = special_inverse(tangent);

    const vec3 viewDir = normalize(tangent.v);

    vec3 diffuse = vec3(0.05);
    vec3 specular = vec3(0);
    float raughness = max(0.215 + texture(u_raughness, IN.texcoord.xy).r - .5, 0.001);
    for(int i = 0; i < LIGHTS.nlights; ++i)
        accumulate_light(tangent, viewDir, raughness, LIGHTS.lights[i], diffuse, specular);

    OUT = vec4(diffuse*IN.diffuse.rgb + specular, IN.diffuse.a);
}

Here are the images I get:

At a grazing angle you see a hard cutoff, but apparently this is how it is supposed to be for high roughness values:

If I add a bump map then the effect isn't visible anymore: