Beckmann Distribution

The normals of the microfacets are distributed according to the Beckmann distribution:

$$\begin{aligned} \operatorname{D}(\omega) &= \frac{1}{\pi a_b^2\cos^4\theta}\exp{\frac{-\tan^{2}\theta}{a_b^2}}\\ \operatorname{D}(\theta,\phi) &= \frac{\sin\theta}{\pi a_b^2\cos^4\theta}\exp{\frac{-\tan^{2}\theta}{a_b^2}} \end{aligned}$$

Here is the result:

function sample_beckmann(a){

    const u_1 = Math.random();
    const u_2 = Math.random();
    
    // Math.log(u_1) works the same way
    const theta = Math.atan(Math.sqrt(-a * a * Math.log(1.0 - u_1)));
    const phi = u_2 * 2.0 * Math.PI;

    return [theta,phi];
}

Algorithm

Sample microfacet normal according to the Beckmann distribution:

1. Choose uniform numbers $u_1, u_2$ in $[0,1]$
2. Calculate $\theta = \tan^{-1}( \sqrt{-a_b^2\log (1-u_1)})$ or $\theta = \tan^{-1}( \sqrt{-a_b^2\log (u_1)})$
3. Calculate $\phi = 2\pi u_2$

Optionally convert to cartesian coordinates with:

$$x = \sin\theta\cos\phi$$

$$y = \sin\theta\sin\phi$$

$$z = \cos\theta$$

PDF Used In Monte Carlo Integration

$$p_m(\omega) = \frac{\cos\theta}{\pi a_b^2\cos^4\theta}\exp{\frac{-\tan^{2}\theta}{a_b^2}} = \frac{1}{\pi a_b^2\cos^3\theta}\exp{\frac{-\tan^{2}\theta}{a_b^2}}$$

$$p_m(\theta,\phi) = \frac{\cos\theta\sin\theta}{\pi a_b^2\cos^4\theta}\exp{\frac{-\tan^{2}\theta}{a_b^2}} = \frac{\tan\theta}{\pi a_b^2\cos^2\theta}\exp{\frac{-\tan^{2}\theta}{a_b^2}}$$

Derivation

Due to normalization, the densities are given as follows:

$$\begin{aligned} p(\omega) &= D(\omega)\cos\theta\\ p(\theta,\phi) &= p(\omega)\sin\theta\\ &= \frac{\sin\theta\cos\theta}{\pi a_b^2\cos^4\theta}\exp{\frac{-\tan^{2}\theta}{a_b^2}} \\ &= \frac{\tan\theta}{\pi a_b^2\cos^2\theta}\exp{\frac{-\tan^{2}\theta}{a_b^2}} \end{aligned}$$

Marginal And Conditional Densities:

$$\begin{aligned} p_\theta(\theta) &= \int_0^{2\pi}p(\theta,\phi)d\phi \\ &=\int_0^{2\pi}\frac{\tan\theta}{\pi a_b^2\cos^2\theta}\exp{\frac{-\tan^{2}\theta}{a_b^2}} d\phi \\ &= 2\pi\frac{\tan\theta}{\pi a_b^2\cos^2\theta}\exp{\frac{-\tan^{2}\theta}{a_b^2}} \\ &= 2\frac{\tan\theta}{ a_b^2\cos^2\theta}\exp{\frac{-\tan^{2}\theta}{a_b^2}} \\ &= 2\pi p(\theta, \phi) \end{aligned}$$

$$p(\phi \vert \theta) = \frac{p(\theta,\phi)}{p_\theta(\theta)} = \frac{p(\theta,\phi)}{2\pi p(\theta,\phi)} = \frac{1}{2\pi}$$

Compute CDFs

$$P_\theta(\theta) = \int_{0}^{\theta}2\frac{\tan\theta'}{ a_b^2\cos^2\theta'}\exp{\frac{-\tan^{2}\theta'}{a_b^2}}d\theta'$$

To solve this integral, we change variables and introduce $u = \tan\theta'$. That gives us $du = (\tan\theta')'d\theta' = \frac{1}{\cos^2\theta'}d\theta'$. Solving for $du$ (with the usual slight abuse of notation): $d\theta' = \cos^2\theta' du$.

As we will resubstitute later, for now let's disregard the limits:

$$\begin{aligned} P_\theta(\theta) &= \int 2\frac{\tan\theta'}{ a_b^2\cos^2\theta'}\exp{(\frac{-\tan^{2}\theta'}{a_b^2})} d\theta' \\ &= \frac{2}{a_b^2} \int \frac{u}{ \cos^2\theta'}\exp{(\frac{-u^{2}}{a_b^2})}\cos^2\theta' du\\ &= \frac{2}{a_b^2} \int u\exp{(\frac{-u^{2}}{a_b^2})} du \end{aligned}$$

The integral can now be solved by basic means by considering the derivative of the exponential.

$$\begin{aligned} P_\theta(\theta) &= \frac{2}{a_b^2} \int u\exp{(\frac{-u^{2}}{a_b^2})} du\\ &= \frac{2}{a_b^2} [-\frac{1}{2}a_b^2\exp(\frac{-u^2}{a_b^2})] \\ &= -[\exp(\frac{-u^2}{a_b^2})] \end{aligned}$$

Resubstituting $u = \tan\theta'$:

$$\begin{aligned} P_\theta(\theta) &= \frac{2}{a_b^2} \int u\exp{(\frac{-u^{2}}{a_b^2})} du\\ &= -[\exp(\frac{-\tan^2\theta'}{a_b^2})]_0^{\theta}\\ &= -(\exp(\frac{-\tan^2\theta}{a_b^2}) - \exp(\frac{-\tan^20}{a_b^2}) ) \\ &= -(\exp(\frac{-\tan^2\theta}{a_b^2}) - 1 ) \\ &= 1- \exp(\frac{-\tan^2\theta}{a_b^2}) \end{aligned}$$

$$P(\phi \vert \theta) = \int_{0}^{\phi}p(\phi ' \vert \theta)d\phi' = \int_{0}^{\phi} \frac{1}{2\pi }d\phi' =\frac{\phi}{2\pi}$$

Invert CDFs

$$\begin{aligned} u_1 &= P_\theta(\theta)\\ &= 1- \exp(\frac{-\tan^2\theta}{a_b^2}) \\ u_1 - 1 &=- \exp(\frac{-\tan^2\theta}{a_b^2}) \\ 1- u_1 &= \exp(\frac{-\tan^2\theta}{a_b^2}) \\ \log (1-u_1) &= \frac{-\tan^2\theta}{a_b^2} \\ -a_b^2\log (1-u_1)&= \tan^2\theta \\ \sqrt{-a_b^2\log (1-u_1)}&= \tan\theta \\ \tan^{-1}( \sqrt{-a_b^2\log (1-u_1)}) &= \theta\\ \theta &=\tan^{-1}( \sqrt{-a_b^2\log (1-u_1)}) \end{aligned}$$

Since $u_1$ is uniformly distributed in $[0,1]$, so is $1-u_1$, so we can use either of them for the generation.

$$\begin{aligned} u_2 &= P(\phi \vert \theta) \\ &= \frac{\phi}{2\pi}\\ \phi &= u_22\pi \end{aligned}$$