Bidirectional Reflectance Distribution Function

Given a reflection plane, the irradiance from the incident light at the point $\mathbf{p}$ can be written as

$$E_i(\mathbf{p})=\int_{H^2} L_i(\mathbf{p},\boldsymbol{\omega}_i)\cos\theta_i\:d\Omega_i,$$

where $d\Omega_i=\sin\theta_i\: d\theta_i\:d\phi_i$. The intuitive understanding of the above formula is that the incident light all over the hemisphere is divided into the incident light from many small cones, in which the irradiance contributed by the small cones in the direction $\boldsymbol{\omega}_i$ is $L_i(\mathbf{p},\boldsymbol{\omega}_i)\cos\theta_i\:d\Omega_i$. The total irradiance $E_i(\mathbf{p})$ from the incident light can be obtained by adding up the irradiance contributed by all the small cones.

The above formula is sometimes written in differential form

$$dE_i(\mathbf{p},\boldsymbol{\omega}_i)=L_i(\mathbf{p},\boldsymbol{\omega}_i)\cos\theta_i\:d\Omega_i.$$

If we want to consider the irradiance rate $L_o(\mathbf{p},\boldsymbol{\omega}_o)$ of the light emitted along the direction $\boldsymbol{\omega}_o$ at the point $\mathbf{p}$, the above idea of dividing and summing still works. Suppose that for every unit of irradiance added by a small cone in the direction $\boldsymbol{\omega}_i$, it can contribute more $f_{\mathrm{r}}\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)$ units of irradiance rate to the outgoing light in the direction $\boldsymbol{\omega}_o$, then the contribution of the small cone in the direction $\boldsymbol{\omega}_i$ to $L_o(\mathbf{p},\boldsymbol{\omega}_o)$ is

$$f_{\mathrm{r}}\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)dE_i(\mathbf{p},\boldsymbol{\omega}_i)=f_{\mathrm{r}}\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)L_i(\mathbf{p},\boldsymbol{\omega}_i)\cos\theta_i\:d\Omega_i.$$

Adding up the contributions of all the small cones gives the expression of $L_o(\mathbf{p},\boldsymbol{\omega}_o)$ as follows

$$L_o(\mathbf{p},\boldsymbol{\omega}_o)=\int_{H^2}f_{\mathrm{r}}\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)dE_i(\mathbf{p},\boldsymbol{\omega}_i)=\int_{H^2}f_{\mathrm{r}}\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)L_i(\mathbf{p},\boldsymbol{\omega}_i)\cos\theta_i\:d\Omega_i.$$

The above formula can also be written in the corresponding differential form

$$dL_o(\mathbf{p},\boldsymbol{\omega}_o,\boldsymbol{\omega}_i)=f_{\mathrm{r}}\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)dE_i(\mathbf{p},\boldsymbol{\omega}_i),$$

or equivalently

$$f_{\mathrm{r}}\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)=\frac{dL_o(\mathbf{p},\boldsymbol{\omega}_o,\boldsymbol{\omega}_i)}{dE_i(\mathbf{p},\boldsymbol{\omega}_i)}.$$

The function $f_{\mathrm{r}}\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)$ derived above is called bidirectional reflectance distribution function, or BRDF for short. The equation derived above

$$L_o(\mathbf{p},\boldsymbol{\omega}_o)=\int_{H^2}f_{\mathrm{r}}\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)L_i(\mathbf{p},\boldsymbol{\omega}_i)\cos\theta_i\:d\Omega_i.$$

is called reflection equation.

Properties

Physics based BRDF should have the following three properties:

1. Nonnegativity: $f_{\mathrm{r}}\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)\ge 0$;
2. Helmholtz reciprocity: $f_{\mathrm{r}}\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)=f_{\mathrm{r}}\left(\mathbf{p}, \boldsymbol{\omega}_i, \boldsymbol{\omega}_o\right)$;
3. Conservation of energy: $$\forall\boldsymbol{\omega}_i,\quad\int_{H^2} f_{\mathrm{r}}\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)\cos\theta_o\:d\Omega_o\le 1.$$

Calculation

There are three common ways to deal with integral

$$L_o(\mathbf{p},\boldsymbol{\omega}_o)=\int_{H^2}f_{\mathrm{r}}\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)L_i(\mathbf{p},\boldsymbol{\omega}_i)\cos\theta_i\:d\Omega_i.$$

1. The first kind of surface integral on the unit hemisphere, which can be parameterized by spherical coordinates

   $$L_o(\mathbf{p},\boldsymbol{\omega}_o)=\int_{0}^{2\pi}\int_{0}^{\frac{\pi}{2}}f_{\mathrm{r}}\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i(\theta_i,\phi_i)\right)L_{i}(\mathbf{p},\boldsymbol{\omega}_i(\theta_i,\phi_i))\cos\theta_i\sin\theta_i\:d\theta_i\:d\phi_i.$$

2. Convert to the integral on luminous surface $A$

   $$L_o(\mathbf{p},\boldsymbol{\omega}_o)=\int_{A} f_{\mathrm{r}}\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right)L_i(\mathbf{p},\boldsymbol{\omega}_i) \cos \theta_{i} \frac{\cos \theta^\prime \mathrm{d} A}{r^{2}},$$

where $\theta^\prime$ is the angle between the normal vector at some point on surface $A$ and $-\boldsymbol{\omega}_i$.

Instance: the BRDF applied by UE4

UE4 engine uses the following BRDF

$$f_{r}\left(\omega_{i}, \omega_{o}\right)=k_{d} f_{\text {lambert}}+f_{s}$$

where

$$\begin{array}{l} f_{\text {lambert}}=\dfrac{\text{base_color}}{\pi} \\ f_{s}=\dfrac{F\left(\mathbf{h}, \mathbf{v}\right) G(\mathbf{l}, \mathbf{v}) D(\mathbf{h})}{4(\mathbf{n} \cdot \mathbf{v})(\mathbf{n} \cdot \mathbf{l})} \\ k_{s}=F\left(\mathbf{h}, \mathbf{v}\right) \\ k_{d}=\left(1-k_{s}\right)(1-\text {metalness}). \end{array}$$

In the above formula, $\mathbf{n}$ is the unit normal vector, $\mathbf{l}$ is the unit incident vector, $\mathbf{v}$ is the unit reflection vector, $\mathbf{h}=(\mathbf{l}+\mathbf{v})/\Vert\mathbf{l}+\mathbf{v}\Vert$ is the microplane unit normal vector. Next, we will describe the specific forms of functions $F,D,G$.

Normal Distribution Function $D$

GGX (Trowbridge-Reitz) Distribution

$$\begin{align} \alpha&=\text {roughness}^{2}\\ D(\mathbf{h})&=\dfrac{\alpha^{2}}{\pi\left((\mathbf{n} \cdot \mathbf{h})^{2}\left(\alpha^{2}-1\right)+1\right)^{2}} \end{align}$$

Geometry Function $G$

GGX-Smith correlated joint approximate formula

$$\begin{align} \alpha&=\text {roughness}^{2} \\ \Lambda(\mathbf{v})&=(\mathbf{n} \cdot \mathbf{l})((\mathbf{n} \cdot \mathbf{v})(1-\alpha)+\alpha) \\ \Lambda(\mathbf{l})&=(\mathbf{n} \cdot \mathbf{v})((\mathbf{n} \cdot \mathbf{l})(1-\alpha)+\alpha) \\ G(\mathbf{l}, \mathbf{v})&=\frac{2(\mathbf{n} \cdot \mathbf{v})(\mathbf{n} \cdot \mathbf{l})}{\Lambda(\mathbf{v})+\Lambda(\mathbf{l})} \end{align}$$

Fresnel Term $F$

Fresnel-Schlick approximate formula

$$\begin{align} F\left(\mathbf{h}, \mathbf{v}\right)&=F_{0}+\left(1-F_{0}\right)(1-(\mathbf{h} \cdot \mathbf{v}))^{5} \\ F_{0}&=\operatorname{lerp}(0.16\:\text{specular}^2, \text{base_color}, \text{metalness}) \end{align}$$

where $\operatorname{lerp}(a,b,t)=(1-t)a+tb$ is the linear interpolation function.