Implementing Cook-Torrance BRDF in GLSL

This article provides a direct guide to formulating and implementing the Cook-Torrance microfacet specular reflection model inside a modern GLSL fragment shader. It outlines the core mathematical equation, defines the three fundamental microfacet functions—the Normal Distribution Function (NDF), the Geometry Shadowing Function, and the Fresnel factor—and concludes with a complete, production-ready GLSL fragment shader implementation.

The Cook-Torrance Specular Equation

The Cook-Torrance microfacet model calculates specular reflection by treating rough surfaces as microscopic mirrors (microfacets). The specular Bidirectional Reflectance Distribution Function (BRDF) is defined mathematically as:

\[f_{\text{specular}} = \frac{D \cdot F \cdot G}{4 (\mathbf{N} \cdot \mathbf{V})(\mathbf{N} \cdot \mathbf{L})}\]

Where:

\(\mathbf{N}\) is the surface normal vector.
\(\mathbf{V}\) is the view (camera) direction vector.
\(\mathbf{L}\) is the light source direction vector.
\(\mathbf{H}\) is the halfway vector, calculated as \(\text{normalize}(\mathbf{L} + \mathbf{V})\).
\(D\) is the Normal Distribution Function.
\(F\) is the Fresnel reflection coefficient.
\(G\) is the Geometric attenuation/shadowing-masking factor.

Component Formulations

Modern physically based rendering (PBR) pipelines standardly use the GGX/Trowbridge-Reitz variants for these terms.

1. Normal Distribution Function (\(D\)): GGX / Trowbridge-Reitz

Describes the alignment of microfacets relative to the halfway vector \(\mathbf{H}\), parameterized by roughness \(\alpha\):

\[D_{\text{GGX}}(\mathbf{N}, \mathbf{H}, \alpha) = \frac{\alpha^2}{\pi \left((\mathbf{N} \cdot \mathbf{H})^2 (\alpha^2 - 1) + 1\right)^2}\]

2. Fresnel Function (\(F\)): Schlick's Approximation

Determines the ratio of light reflected versus refracted based on the viewing angle, using the base surface reflectivity \(F_0\):

\[F_{\text{Schlick}}(\mathbf{V}, \mathbf{H}, F_0) = F_0 + (1.0 - F_0) (1.0 - (\mathbf{V} \cdot \mathbf{H}))^5\]

3. Geometric Shadowing Function (\(G\)): Smith's Model with Schlick-GGX

Accounts for microfacets occluding or casting shadows on one another. It combines geometry terms for both the view and light directions:

\[G(\mathbf{N}, \mathbf{V}, \mathbf{L}, k) = G_1(\mathbf{N}, \mathbf{V}, k) \cdot G_1(\mathbf{N}, \mathbf{L}, k)\]

Where \(G_1\) is defined as:

\[G_1(\mathbf{N}, \mathbf{X}, k) = \frac{\mathbf{N} \cdot \mathbf{X}}{(\mathbf{N} \cdot \mathbf{X})(1 - k) + k}\]

For direct analytic lights, \(k = \frac{(\alpha + 1)^2}{8}\), where \(\alpha = \text{roughness}^2\).

GLSL Fragment Shader Formulation

Below is the GLSL formulation implementing these functions and combining them into the final specular contribution.

#version 330 core

in vec3 vWorldPos;
in vec3 vNormal;

out vec4 FragColor;

uniform vec3 uCameraPos;
uniform vec3 uLightPos;
uniform vec3 uLightColor;
uniform vec3 uAlbedo;
uniform float uRoughness;
uniform float uMetallic;

const float PI = 3.14159265359;

// 1. Normal Distribution Function (Trowbridge-Reitz GGX)
float DistributionGGX(vec3 N, vec3 H, float roughness) {
    float a = roughness * roughness;
    float a2 = a * a;
    float NdotH = max(dot(N, H), 0.0);
    float NdotH2 = NdotH * NdotH;

    float denom = (NdotH2 * (a2 - 1.0) + 1.0);
    denom = PI * denom * denom;

    return a2 / max(denom, 0.0000001);
}

// 2. Geometry Function helper (Schlick-GGX)
float GeometrySchlickGGX(float NdotV, float roughness) {
    float r = (roughness + 1.0);
    float k = (r * r) / 8.0;

    float denom = NdotV * (1.0 - k) + k;
    return NdotV / max(denom, 0.0000001);
}

// Geometry Function (Smith's method)
float GeometrySmith(vec3 N, vec3 V, vec3 L, float roughness) {
    float NdotV = max(dot(N, V), 0.0);
    float NdotL = max(dot(N, L), 0.0);
    float ggx1 = GeometrySchlickGGX(NdotV, roughness);
    float ggx2 = GeometrySchlickGGX(NdotL, roughness);

    return ggx1 * ggx2;
}

// 3. Fresnel Equation (Fresnel-Schlick approximation)
vec3 FresnelSchlick(float cosTheta, vec3 F0) {
    return F0 + (1.0 - F0) * pow(clamp(1.0 - cosTheta, 0.0, 1.0), 5.0);
}

void main() {
    vec3 N = normalize(vNormal);
    vec3 V = normalize(uCameraPos - vWorldPos);
    vec3 L = normalize(uLightPos - vWorldPos);
    vec3 H = normalize(V + L);

    // Calculate base reflectivity F0 (0.04 for dielectrics, albedo for metals)
    vec3 F0 = vec3(0.04);
    F0 = mix(F0, uAlbedo, uMetallic);

    // Compute Cook-Torrance BRDF components
    float D = DistributionGGX(N, H, uRoughness);
    float G = GeometrySmith(N, V, L, uRoughness);
    vec3 F  = FresnelSchlick(max(dot(H, V), 0.0), F0);

    // Specular Cook-Torrance BRDF
    vec3 numerator = D * G * F;
    float denominator = 4.0 * max(dot(N, V), 0.0) * max(dot(N, L), 0.0) + 0.0001;
    vec3 specular = numerator / denominator;

    // Energy conservation: diffuse vs specular fraction
    vec3 kS = F;
    vec3 kD = vec3(1.0) - kS;
    kD *= 1.0 - uMetallic;

    // Outgoing radiance calculation
    float NdotL = max(dot(N, L), 0.0);
    vec3 diffuse = (kD * uAlbedo / PI);
    
    vec3 radiance = (diffuse + specular) * uLightColor * NdotL;

    FragColor = vec4(radiance, 1.0);
}

Key Considerations for Implementation

Epsilon Clamping: Small float additions (e.g., 0.0001) in denominators prevent hardware division-by-zero errors when surface angles approach grazing values.
Roughness Parameterization: Squaring perceptual roughness (\(\alpha = \text{roughness}^2\)) before passing it to GGX produces more visually linear specular highlights.
Energy Conservation: \(k_D\) and \(k_S\) ensure that the total reflected light (diffuse plus specular) does not exceed the incoming energy. Metallic surfaces absorb diffuse components entirely, which is represented by kD *= 1.0 - uMetallic.