Phong Distribution

The normals of the microfacets are distributed according to the Phong Distribution:

D⁡(ω)=ap+22πcos⁡apθD⁡(θ,ϕ)=ap+22πcos⁡apθsin⁡θ\begin{aligned} \operatorname{D}(\omega) &= \frac{a_p +2}{2\pi}\cos^{a_p}\theta\\ \operatorname{D}(\theta,\phi) &= \frac{a_p +2}{2\pi}\cos^{a_p}\theta\sin\theta \end{aligned}

Walter et al. recommend a relation between the Beckmann aba_b and the phong apa_p:

ap=2ab−2−2\begin{aligned} a_p &= 2a_b^{-2} - 2 \end{aligned}

Here is the result:

function pow_cos_hemisphere_cap(theta_max,n){

    const u_1 = Math.random();
    const u_2 = Math.random();

    const theta = Math.acos(Math.pow(
        1.0 - u_1 * (1.0 - Math.pow(Math.cos(theta_max),n+1)), 
        1.0/(n+1)));
    const phi = u_2 * 2.0 * Math.PI;

    return [theta,phi];
}

function beckmann_width_to_phong(a) {
    return 2.0 / (a*a) - 2.0;
}

function sample_phong_normal(a) {

    return pow_cos_hemisphere_cap(Math.PI / 2.0, a + 1);
}

Algorithm

Sample microfacet normal according to the Phong distribution:

  1. Choose uniform numbers u1,u2u_1, u_2 in [0,1][0,1]
  2. Calculate θ=cos⁡−1(1−u1ap+2)\theta = \cos^{-1}(\sqrt[a_p+2]{1 - u_1} ) or θ=cos⁡−1(u1ap+2)\theta = \cos^{-1}(\sqrt[a_p+2]{u_1} )
  3. Calculate ϕ=2πu2\phi = 2\pi u_2

Optionally convert to cartesian coordinates with:

x=sin⁡θcos⁡ϕx = \sin\theta\cos\phi

y=sin⁡θsin⁡ϕy = \sin\theta\sin\phi

z=cos⁡θz = \cos\theta

PDF Used In Monte Carlo Integration

pm(ω)=ap+22πcos⁡apθcos⁡θ=ap+22πcos⁡ap+1θp_m(\omega) = \frac{a_p +2}{2\pi}\cos^{a_p}\theta\cos\theta = \frac{a_p +2}{2\pi}\cos^{a_p+1}\theta

pm(θ,ϕ)=ap+22πcos⁡apθcos⁡θsin⁡θ=ap+22πcos⁡ap+1θsin⁡θp_m(\theta,\phi) = \frac{a_p +2}{2\pi}\cos^{a_p}\theta\cos\theta\sin\theta = \frac{a_p +2}{2\pi}\cos^{a_p +1}\theta\sin\theta

Derivation

This distribution is just a special case of the power cosine hemisphere cap above, but with

θmax=π2,n=ap+1\theta_{max}=\frac{\pi}{2}, n = a_p+1

p(ω)=n+12π(1−cos⁡n+1θmax)cos⁡nθ=n+12πcos⁡nθ=ap+22πcos⁡ap+1θp(\omega) = \frac{n+1}{2\pi(1 - \cos^{n+1}\theta_{max}) }\cos^n\theta = \frac{n+1}{2\pi }\cos^n\theta=\frac{a_p+2}{2\pi }\cos^{a_p + 1}\theta

The formula given by Walter et al. comes from replacing 1−u11-u_1 by u1u_1, since both are uniformly distributed within [0,1][0,1].