Physically based shading (microfacet BRDF)

Deep dive

How it works

Parts:

• BRDF f(l, v) — reciprocal (f(l,v) = f(v,l)) and energy-bounded (∫ f (n·l) dl ≤ 1)
• material mapping — base colour, metallic, roughness to c_diff, F0, α
• D (GGX) — density of facets with normal h
• G (Smith, height-correlated) — fraction of facets visible from both l and v
• F (Schlick) — reflectance of a facet as a function of v·h
• light loop — the rendering equation as a sum over lights

Operations, per pixel and light (rendering equation, Kajiya 1986: L_o = ∫ f L_i (n·l) dω, for punctual lights L_o = Σ f · E_i · (n·l)):

1. Map the material (glTF 2.0 Appendix B): c_diff = base·(1−metallic), F0 = lerp(0.04, base, metallic), α = roughness² (Disney's remap: more perceptually linear; Burley 2012).
2. Vectors: h = normalize(v + l); clamp n·v, n·l away from 0.
3. D (GGX / Trowbridge-Reitz, Walter et al. 2007): α² / (π ((n·h)²(α²−1) + 1)²); normalised so that ∫ D (n·h) dω = 1.
4. G / V: with Smith's Λ, height-correlated G₂ = 1/(1 + Λ(v) + Λ(l)) (Heitz 2014). Folded with the denominator: V = G₂/(4 n·l n·v) = 0.5 / (n·l √(n·v²(1−α²)+α²) + n·v √(n·l²(1−α²)+α²)).
5. F (Schlick 1994): F = F0 + (1 − F0)(1 − v·h)⁵, per RGB channel.
6. Combine (glTF): f = (1−F)·c_diff/π + F·D·V; multiply by light irradiance and n·l.
7. Sum over lights, scaled by shadow visibility; add IBL; output linear HDR for tone mapping.

          n   h
           \  |  /  l (light)       f = (1-F) c/π + F·D·V
      v     \ | /                   diffuse lobe: all directions
        \    \|/                    specular lobe: around the mirror direction
   ───────────────── surface

 base, metallic, roughness ──► c_diff, F0, α ──┬──► D(n·h, α) ───────┐
 v, l ──► h, n·v, n·l, v·h ────────────────────┼──► V(n·v, n·l, α) ──┼──► f ──► ·E_i·(n·l) ──► Σ lights
                                               └──► F(F0, v·h) ──────┘

Code

use std::f32::consts::PI;
use std::hint::black_box;
use std::time::Instant;

type V3 = [f32; 3];
fn dot(a: V3, b: V3) -> f32 { a[0] * b[0] + a[1] * b[1] + a[2] * b[2] }
fn norm(a: V3) -> V3 { let k = dot(a, a).sqrt().recip(); [a[0] * k, a[1] * k, a[2] * k] }

/// GGX / Trowbridge-Reitz normal distribution; alpha = roughness^2.
fn d_ggx(nh: f32, a: f32) -> f32 {
    let a2 = a * a;
    let d = nh * nh * (a2 - 1.0) + 1.0;
    a2 / (PI * d * d)
}
/// Height-correlated Smith visibility V = G2 / (4 n.l n.v) (Heitz 2014, form as in Filament).
fn v_smith(nv: f32, nl: f32, a: f32) -> f32 {
    let a2 = a * a;
    0.5 / (nl * (nv * nv * (1.0 - a2) + a2).sqrt() + nv * (nl * nl * (1.0 - a2) + a2).sqrt())
}
/// Schlick's Fresnel approximation.
fn f_schlick(f0: f32, vh: f32) -> f32 { f0 + (1.0 - f0) * (1.0 - vh).powi(5) }

/// Metal-roughness BRDF f(l, v) for unit n, v, l (all pointing away from the surface).
/// Multiply by radiance * n.l outside. Diffuse is glTF's (1 - F) * Lambert.
fn evaluate(n: V3, v: V3, l: V3, base: V3, metallic: f32, rough: f32) -> V3 {
    let (nv, nl) = (dot(n, v), dot(n, l));
    if nv <= 0.0 || nl <= 0.0 { return [0.0; 3]; }
    let h = norm([v[0] + l[0], v[1] + l[1], v[2] + l[2]]);
    let (nh, vh) = (dot(n, h).min(1.0), dot(v, h).min(1.0));
    let a = (rough * rough).max(1e-3);
    let ds = d_ggx(nh, a) * v_smith(nv, nl, a);
    std::array::from_fn(|i| {
        let f0 = 0.04 * (1.0 - metallic) + base[i] * metallic;
        let f = f_schlick(f0, vh);
        (1.0 - f) * (1.0 - metallic) * base[i] / PI + f * ds
    })
}

/// Directional albedo E(v) = integral f n.l dl, midpoint rule, cosine-weighted.
fn albedo(cos_v: f32, base: V3, metallic: f32, rough: f32, m: usize) -> f32 {
    let n = [0.0, 0.0, 1.0];
    let v = [(1.0 - cos_v * cos_v).sqrt(), 0.0, cos_v];
    let mut s = 0.0f64;
    for i in 0..m {
        let ct = (1.0 - (i as f32 + 0.5) / m as f32).sqrt(); // pdf = cos/pi
        let st = (1.0 - ct * ct).sqrt();
        for j in 0..m {
            let phi = 2.0 * PI * (j as f32 + 0.5) / m as f32;
            let f = evaluate(n, v, [st * phi.cos(), st * phi.sin(), ct], base, metallic, rough);
            s += (f[1] * PI) as f64;
        }
    }
    (s / (m * m) as f64) as f32
}

fn main() {
    assert_eq!(f_schlick(0.04, 1.0), 0.04);
    // D integrates to 1 when weighted by n.h
    for r in [0.3f32, 0.5, 0.8, 1.0] {
        let (a, k) = (r * r, 20000);
        let dt = 0.5 * PI / k as f32;
        let s: f64 = (0..k).map(|i| { let t = (i as f32 + 0.5) * dt; (d_ggx(t.cos(), a) * t.cos() * t.sin()) as f64 }).sum();
        let integral = 2.0 * PI as f64 * s * dt as f64;
        assert!((integral - 1.0).abs() < 2e-3, "D integral {integral} at roughness {r}");
    }
    // reciprocity
    let mut seed = 12345u32;
    let mut rnd = || { seed = seed.wrapping_mul(1664525).wrapping_add(1013904223); (seed >> 8) as f32 / 16777216.0 };
    let mut dirs = Vec::new();
    while dirs.len() < 4096 { // two random hemisphere directions per pair
        let mut d = [0.0; 6];
        for k in 0..2 {
            let p = norm([rnd() - 0.5, rnd() - 0.5, rnd() * 0.5 + 0.05]);
            d[3 * k..3 * k + 3].copy_from_slice(&p);
        }
        dirs.push(d);
    }
    let (n, base) = ([0.0, 0.0, 1.0], [0.9, 0.5, 0.2]);
    for d in &dirs {
        let (v, l) = ([d[0], d[1], d[2]], [d[3], d[4], d[5]]);
        let (a, b) = (evaluate(n, v, l, base, 0.5, 0.6), evaluate(n, l, v, base, 0.5, 0.6));
        for c in 0..3 { assert!((a[c] - b[c]).abs() <= 1e-4 * a[c].abs().max(1e-3)); }
    }
    // energy: a white metal (F = 1) may not return more than it receives
    let cosv = [1.0f32, 0.7, 0.4, 0.1];
    for r in [0.3f32, 0.5, 0.7, 1.0] {
        let e: Vec<f32> = cosv.iter().map(|&c| albedo(c, [1.0; 3], 1.0, r, 400)).collect();
        assert!(e.iter().all(|&x| x <= 1.001), "metal albedo > 1: {e:?}");
        let d: Vec<f32> = cosv.iter().map(|&c| albedo(c, [1.0; 3], 0.0, r, 400)).collect();
        println!("roughness {r}: white metal {e:.3?}, white dielectric {d:.3?} (cos v = 1, .7, .4, .1)");
    }
    println!("metal at roughness 1, normal view: {:.1}% energy lost", 100.0 * (1.0 - albedo(1.0, [1.0; 3], 1.0, 1.0, 400)));
    // benchmark: ns per evaluate()
    let t = Instant::now();
    let mut acc = 0.0f32;
    for _ in 0..2000 {
        for d in &dirs { acc += evaluate(n, [d[0], d[1], d[2]], [d[3], d[4], d[5]], black_box(base), 0.5, 0.6)[1]; }
    }
    black_box(acc);
    println!("{:.1} ns per evaluate()", t.elapsed().as_nanos() as f64 / (2000.0 * 4096.0));
}

Asserts: Fresnel at normal incidence equals F0; ∫ D (n·h) dω = 1 within 2e-3; reciprocity to 1e-4 relative on 4096 random pairs; white-metal albedo at most 1.001. Measured on this Mac (Apple M4 Pro), release build: a white metal at normal view returns 0.99 of the light at roughness 0.3 and 0.31 at roughness 1 (69.3 % lost). A white dielectric with the glTF split returns up to 1.28 at roughness 0.3, grazing view (cos v = 0.1), so that split is not strictly energy-conserving there. evaluate: 5.3 to 5.8 ns over several runs.

Complexity and constants

• Per light and pixel: 2 square roots, 1 division, a fifth power and 3 channels of F: about 5.6 ns here (measured, scalar f32, one core, loop included).
• Arithmetic from that: 1920x1080 with 10 lights is about 0.12 s on one core; light count drives cost.
• Hoist F0, c_diff, α and n·v out of the light loop.

Efficiency

• Cost: measured on this Mac (Apple M4 Pro): 5.3 to 5.8 ns per RGB BRDF evaluation, scalar Rust, one core, directions only varying. Bound: arithmetic latency (sqrt, division) on a CPU; GPU cost not measured.
• Levers:
  • Cull lights per tile or cluster — work per pixel falls from all lights to the few that reach it — clustered shading — many lights
  • Prefiltered environment plus a 2D DFG lookup for ambient light — one fetch instead of a hemisphere integral (Karis 2013; Filament) — IBL
  • V without square roots (Filament: 0.5/(n·l(n·v(1−α)+α)+n·v(n·l(1−α)+α))) — no gain here, 5.3 against 5.5 ns, measured — GPU effect not measured
  • Energy compensation from the same DFG table — negligible cost per Filament — rough metals
  • Clamp roughness (0.045 to 0.089) — avoids fp16 trouble and specular aliasing (Filament) — mobile, analytic lights
• Naive → standard → state of the art:
  • every light for every pixel, separable G, no energy fix · cost grows with lights x pixels
  • height-correlated V, Schlick F, culled lights, split-sum IBL · the glTF/Filament/Bevy shape, about 5 ns per light here
  • plus energy-compensation lobe and layered coat · the same shape plus one or two table fetches

Implementation notes

• Work in linear light: decode sRGB base colour first (colour conversion); roughness and metallic stay linear.
• evaluate returns the BRDF only; the caller multiplies by E_i and n·l. Store α, not roughness, if the material is static.
• Keep α above a floor (code: 1e-3; more in fp16): D peaks at 1/(π α²).

Pitfalls

• Dividing by n·l or n·v at the horizon: clamp (glTF uses abs(n·v)).
• sRGB colours treated as linear, or roughness passed where α is expected.
• Energy: rough metals darken (69 % loss at roughness 1 here); the glTF split can exceed 1 for smooth dielectrics at grazing view.
• A diffuse term on metals: glTF sets c_diff to black for metallic = 1.
• D is 0/0 at n·h = 1 when α = 0: clamp roughness.

Variants and alternatives

In production

• Real-time: Filament and Bevy implement GGX, height-correlated Smith and Schlick (Bevy's pbr_lighting.wesl credits the Filament doc); Epic documented GGX, Schlick Fresnel and α = roughness² for Unreal Engine 4 (Karis 2013).
• Film: Disney (Burley 2012), Sony Imageworks (Kulla and Conty 2017); pbrt-v4's TrowbridgeReitzDistribution uses 1/(1+Λ(wo)+Λ(wi)).
• Web 3D: glTF 2.0 Appendix B defines this model; the Khronos Sample Viewer renders it.

Implementations

Sources

• Cook, Torrance, A Reflectance Model for Computer Graphics, ACM TOG 1(1), 1982. DOI
• Kajiya, The Rendering Equation, SIGGRAPH 1986. DOI
• Schlick, An Inexpensive BRDF Model for Physically-based Rendering, Computer Graphics Forum 13(3), 1994. DOI
• Walter, Marschner, Li, Torrance, Microfacet Models for Refraction through Rough Surfaces, EGSR 2007. PDF
• Heitz, Understanding the Masking-Shadowing Function in Microfacet-Based BRDFs, JCGT 3(2), 2014. paper
• Burley, Physically-Based Shading at Disney, SIGGRAPH 2012 course notes. PDF
• Karis, Real Shading in Unreal Engine 4, SIGGRAPH 2013 course notes. PDF
• Heitz, Hanika, d'Eon, Dachsbacher, Multiple-scattering microfacet BSDFs with the Smith model, ACM TOG 35(4), 2016. DOI
• Kulla, Conty, Revisiting Physically Based Shading at Imageworks, SIGGRAPH 2017 course. slides
• Khronos, glTF 2.0 Specification, Appendix B: BRDF Implementation. spec
• Google, Physically Based Rendering in Filament. doc