Image Based Real-time Environment Lighting

Image-based lighting refers to storing the radiance of environment light from each direction with a texture, which is then used for shading calculation. This texture is called environment texture.

The original reflection equation is

$$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.$$

Here we do not consider lighting other than environment light. And we assume that the environment light is the same for any shading point $\mathbf{p}$. Thus the environment texture can be denoted as $L_i(\boldsymbol{\omega}_i)$．Accordingly, the reflection equation can be written as

$$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(\boldsymbol{\omega}_i)\cos\theta_i\:d\Omega_i.$$

In this article, the geometric information at the point $\mathbf{p}$ that we need is only the normal vector $\mathbf{n}$. So the above formula can be rewritten as

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

Diffuse Reflection

Simplified Fresnel Term

If we only consider the diffuse reflection, the reflection equation is

$$\begin{aligned} L_o(\mathbf{n},\boldsymbol{\omega}_o)&=\int_{H^2}k_{d} f_{\text {lambert}}L_i(\boldsymbol{\omega}_i)\cos\theta_i\:d\Omega_i\\ &=\frac{\mathrm{base\_color}(1-\text {metalness})}{\pi} \int_{H^2}\left(1-F\left(\mathbf{h}, \boldsymbol{\omega}_o\right)\right)L_i(\boldsymbol{\omega}_i)\mathbf{n}\cdot\boldsymbol{\omega}_i \:d\Omega_i, \end{aligned}$$

where

$$\begin{aligned} F\left(\mathbf{h}, \boldsymbol{\omega}_o\right)&=F_{0}+\left(1-F_{0}\right)(1-(\mathbf{h} \cdot \boldsymbol{\omega}_o))^{5}, \\ F_{0}&=\mathrm{lerp}(0.16\:\mathrm{specular}^2, \mathrm{base\_color}, \mathrm{metalness}). \end{aligned}$$

In real-time rendering, a common approach is to precompute complex integrals. Generally, for any $\mathbf{n}$ and $\boldsymbol{\omega}_o$, we aim to obtain the corresponding $L_o(\mathbf{n}, \boldsymbol{\omega}_o)$. On one hand, there are infinitely many possible values for $\boldsymbol{\omega}_o$. On the other hand, considering the rotation of objects, $\mathbf{n}$ also has infinitely many values. It is clearly impossible to precompute $L_o(\mathbf{n}, \boldsymbol{\omega}_o)$ for all $\mathbf{n}$ and $\boldsymbol{\omega}_o$. The most naive idea is to precompute for a limited number of $\mathbf{n}$ and $\boldsymbol{\omega}_o$ and then interpolate the results. However, this requires creating a 4D LUT (lookup table), which consumes a vast amount of space and time. A more economical method is to simplify the integrand in the reflection equation:

$$\left(1-F\left(\mathbf{h}, \boldsymbol{\omega}_o\right)\right)L_i(\boldsymbol{\omega}_i)\mathbf{n}\cdot\boldsymbol{\omega}_i,$$

allowing the terms involving $\boldsymbol{\omega}_o$ to be moved outside the integral. Thus, we only need to precompute complex integrals for a limited number of $\mathbf{n}$ and store the results in a 2D LUT. The computations regarding $\boldsymbol{\omega}_o$ can be performed in real time.

Note that in the integrand, only the Fresnel term $F\left(\mathbf{h}, \boldsymbol{\omega}_o\right)$ includes $\boldsymbol{\omega}_o$. Since $\mathbf{h}$ contains $\boldsymbol{\omega}_i$, this term cannot be directly moved outside the integral. We can replace the original Fresnel term $F\left(\mathbf{h}, \boldsymbol{\omega}_o\right)$ with a simplified version $F^*\left(\mathbf{n}, \boldsymbol{\omega}_o\right)$, which is expressed as follows:

$$\begin{aligned} F^*\left(\mathbf{n},\boldsymbol{\omega}_o\right) &=F_{0}+\left(S-F_{0}\right)\left(1-\mathbf{n}\cdot \boldsymbol{\omega}_o\right)^{5},\\ S &=\max \left\{1-\text { roughness }, F_{0}\right\}. \end{aligned}$$

Thus, the reflection equation simplifies to:

$$\begin{aligned} L_o(\mathbf{n},\boldsymbol{\omega}_o) \approx\frac{\mathrm{base\_color}(1-\text {metalness})\left(1-F^*\left(\mathbf{n},\boldsymbol{\omega}_o\right)\right)}{\pi} \int_{H^2}L_i(\boldsymbol{\omega}_i)\mathbf{n}\cdot\boldsymbol{\omega}_i \:d\Omega_i. \end{aligned}$$

Afterwards, we only need to precompute the following integral for a finite number of $\mathbf{n}$:

$$I_{\mathrm{diff}}\left(\mathbf{n}\right)=\int_{H^2}L_i(\boldsymbol{\omega}_i)\mathbf{n}\cdot\boldsymbol{\omega}_i \:d\Omega_i.$$

Note that this integral represents the irradiance received by a shading point with normal $\mathbf{n}$, hence it is also referred to as an irradiance map.

If the environmental texture is a cubic box map, each of its texels corresponds to a direction vector. We can let $\mathbf{n}$ cover all these direction vectors and store the value of $I_{\mathrm{diff}}\left(\mathbf{n}\right)$ in the corresponding texel. During rendering, $L_o(\mathbf{n},\boldsymbol{\omega}_o)$ is calculated as follows:

$$L^*_o(\mathbf{n},\boldsymbol{\omega}_o)=\frac{\mathrm{base\_color}(1-\text {metalness})\left(1-F^*\left(\mathbf{n},\boldsymbol{\omega}_o\right)\right)}{\pi}I_{\mathrm{diff}}\left(\mathbf{n}\right).$$

Precomputation of the Irradiance Map

We can compute $I_{\mathrm{diff}}\left(\mathbf{n}\right)$ using the Monte Carlo integration method. Given $\mathbf{n}$, let $\eta_{\mathbf{n}}$ be a probability distribution over the hemisphere $H^2$ with the probability density function as

$$pdf_{\mathbf{n}}(\boldsymbol{\omega}_i) =\frac{\mathbf{n}\cdot\boldsymbol{\omega}_i}{\pi}.$$

We utilize $\eta_{\mathbf{n}}$ for importance sampling to obtain the statistic of $I_{\mathrm{diff}}\left(\mathbf{n}\right)$:

$$\widehat{I}_{\mathrm{diff}}\left(\mathbf{n}\right)=\frac{\pi}{N}\sum_{k=1}^NL_i(\boldsymbol{\omega}_i^{(k)}),\quad \boldsymbol{\omega}_i^{(k)}\sim\eta_{\mathbf{n}},$$

where $N$ is the number of samples. Below, we will specifically describe how to generate samples for $\eta_{\mathbf{n}}$. We will provide two methods: the coordinate transformation method and the quaternion rotation method.

Coordinate Transformation Method：
Define

$$\mathbf{u}= \begin{cases} (0,0,1) &\mathrm{if}\; \mathbf{n}\ne(0,0,\pm1), \\ (1,0,0) & \mathrm{if}\; \mathbf{n}=(0,0,\pm1). \end{cases}$$

The matrix $T_{\mathbf{n}}$ of the basis vectors of the tangent space in world space can be represented as

$$T_{\mathbf{n}}=\left(\mathbf{t},\mathbf{b},\mathbf{n}\right),$$

where

$$\begin{aligned} \mathbf{t}&=\frac{\mathbf{u}\times\mathbf{n}}{\Vert\mathbf{u}\times\mathbf{n}\Vert},\\ \mathbf{b}&=\mathbf{n}\times\mathbf{t}. \end{aligned}$$

Suppose $\xi_{1}^{(k)}, \xi_{2}^{(k)}$ are independent and uniformly distributed $U(0,1)$, then

$$T_{\mathbf{n}}\boldsymbol{\omega}_i^{(k)}=T_{\mathbf{n}} \begin{bmatrix} \sqrt{\xi_{1}} \cos \left(2 \pi \xi_{2}\right)\\ \sqrt{\xi_{1}} \sin \left(2 \pi \xi_{2}\right)\\ \sqrt{1-\xi_{1}} \end{bmatrix}$$

follows the distribution $\eta_{\mathbf{n}}$.

Quaternion Rotation Method：

• When $\mathbf{n}=(0,0,1)$, we can use inverse transform sampling. Suppose $\xi_{1}^{(k)}, \xi_{2}^{(k)}$ are independent and uniformly distributed $U(0,1)$, then

  $$\boldsymbol{\omega}_i^{(k)}=\left(\sqrt{\xi_{1}} \cos \left(2 \pi \xi_{2}\right), \sqrt{\xi_{1}} \sin \left(2 \pi \xi_{2}\right), \sqrt{1-\xi_{1}}\right)$$

  follows the distribution $\eta_{(0,0,1)}$.

• When $\mathbf{n}=(0,0,-1)$, $-\boldsymbol{\omega}_i^{(k)}$ follows the distribution $\eta_{(0,0,-1)}$.

• When $\mathbf{n}=(n_x,n_y,n_z)\ne(0,0,\pm 1)$, rotate the vector $(0,0,1)$ along the rotation axis

  $$(0,0,1)\times\mathbf{n}=(-n_y,n_x,0)$$

  to vector $\mathbf{n}$, the quaternion is

  $$q_\mathbf{n} = \cos\frac{\alpha}{2}+\sin\frac{\alpha}{2}\left(-\frac{n_y}{\sqrt{n_x^2+n_y^2}}i+\frac{n_x}{\sqrt{n_x^2+n_y^2}}j\right),$$

  where $\alpha$ is the rotation angle. Since

  $$\cos\alpha=(0,0,1)\cdot \mathbf{n}=n_z.$$

  we have

  $$\begin{aligned} q_\mathbf{n}&=\sqrt{\frac{1+\cos\alpha}{2}}+\sqrt{\frac{1-\cos\alpha}{2}}\left(-\frac{n_y}{\sqrt{n_x^2+n_y^2}}i+\frac{n_x}{\sqrt{n_x^2+n_y^2}}j\right)\\ &=\sqrt{\frac{1+n_z}{2}}+\sqrt{\frac{1-n_z}{2}}\left(-\frac{n_y}{\sqrt{n_x^2+n_y^2}}i+\frac{n_x}{\sqrt{n_x^2+n_y^2}}j\right). \end{aligned}$$

  From this, we can compute the rotation matrix for the quaternion rotation $\mathrm{rot}_q:p\mapsto qpq^{-1}$

  $$R_\mathbf{n}= \begin{pmatrix} \tfrac{n_x^2n_z+n_y^2}{n_x^2+n_y^2}&-n_xn_y(1+n_z)&n_x\\ -n_xn_y(1+n_z)&\tfrac{n_x^2+n_y^2n_z}{n_x^2+n_y^2}&n_y\\ -n_x&-n_y&n_z \end{pmatrix}.$$

  It is evident that $R_\mathbf{n}\boldsymbol{\omega}_i^{(k)}$ follows the distribution $\eta_{\mathbf{n}}$.

Specular Reflection

Now we consider only specular reflection. The reflection equation is given by

$$\begin{aligned} L_o(\mathbf{n},\boldsymbol{\omega}_o)&=\int_{H^2}f_s(\mathbf{n},\boldsymbol{\omega}_i, \boldsymbol{\omega}_o)L_i(\boldsymbol{\omega}_i)\cos\theta_i\:d\Omega_i\\ &=\int_{H^2}\dfrac{F\left(\mathbf{h}, \boldsymbol{\omega}_o\right) G(\mathbf{n},\boldsymbol{\omega}_i, \boldsymbol{\omega}_o) D(\mathbf{n},\mathbf{h})}{4(\mathbf{n} \cdot \boldsymbol{\omega}_o)(\mathbf{n} \cdot \boldsymbol{\omega}_i)}L_i(\boldsymbol{\omega}_i)\mathbf{n}\cdot\boldsymbol{\omega}_i \:d\Omega_i. \end{aligned}$$

To simplify the equation, we can bring $L_i(\boldsymbol{\omega}_i)$ out of the integral. Thus, we express $L_o(\mathbf{n},\boldsymbol{\omega}_o)$ as follows

$$\begin{aligned} L_o(\mathbf{n},\boldsymbol{\omega}_o)&=\int_{H^2}f_s(\mathbf{n},\boldsymbol{\omega}_i, \boldsymbol{\omega}_o) L_{i}(\boldsymbol{\omega}_i)(\mathbf{n} \cdot \boldsymbol{\omega}_i) d\Omega_i \\ &=\frac{\int_{H^2}f_s(\mathbf{n},\boldsymbol{\omega}_i, \boldsymbol{\omega}_o) L_{i}(\boldsymbol{\omega}_i)(\mathbf{n} \cdot \boldsymbol{\omega}_i) d\Omega_i}{\int_{H^2}f_s(\mathbf{n},\boldsymbol{\omega}_i, \boldsymbol{\omega}_o)\mathbf{n} \cdot \boldsymbol{\omega}_i d\Omega_i}\int_{H^2} f_s(\mathbf{n},\boldsymbol{\omega}_i, \boldsymbol{\omega}_o)\mathbf{n} \cdot \boldsymbol{\omega}_i d\Omega_i, \end{aligned}$$

We aim to compute the parts containing $L_i(\boldsymbol{\omega}_i)$ and those not containing $L_i(\boldsymbol{\omega}_i)$ separately.

Split Sum - Part 1

Formula Derivation

Direct calculation

$$\frac{\int_{H^2}f_s(\mathbf{n},\boldsymbol{\omega}_i, \boldsymbol{\omega}_o) L_{i}(\boldsymbol{\omega}_i)(\mathbf{n} \cdot \boldsymbol{\omega}_i) d\Omega_i}{\int_{H^2}f_s(\mathbf{n},\boldsymbol{\omega}_i, \boldsymbol{\omega}_o)\mathbf{n} \cdot \boldsymbol{\omega}_i d\Omega_i}$$

is relatively complex. To reduce the dimensionality of the precomputed parameter space, the assumption here is that the lobe shape remains unchanged, always being the lobe when the ray is incident along the normal, i.e., when calculating the integral, let the values of $\boldsymbol{\omega}_o, \mathbf{n}$ be equal to the value of the reflection vector $\mathbf{r}(\boldsymbol{\omega}_o)=2(\boldsymbol{\omega}_o\cdot\mathbf{n})\mathbf{n}-\boldsymbol{\omega}_o$ of $\boldsymbol{\omega}_o$. In other words, use $\mathbf{r}(\boldsymbol{\omega}_o)$ to replace $\boldsymbol{\omega}_o$ and $\mathbf{n}$ in the integrand. Let $\mathbf{n}^*=\mathbf{r}(\boldsymbol{\omega}_o)$,

$$\begin{aligned} &\;f^*_s\left(\boldsymbol{\omega}_i, \boldsymbol{\omega}_o\right)\\ =&\;f_s\left(\mathbf{n}^*,\boldsymbol{\omega}_i, \mathbf{n}^*\right)\\ =&\;\dfrac{F\left(\mathbf{h}, \boldsymbol{\omega}_o\right)}{4(\mathbf{n}^* \cdot \boldsymbol{\omega}_o)(\mathbf{n}^* \cdot \boldsymbol{\omega}_i)} G(\mathbf{n}^*,\boldsymbol{\omega}_i, \boldsymbol{\omega}_o) D(\mathbf{n}^*,\mathbf{h})\\ =&\;\dfrac{F_{0}+\left(1-F_{0}\right)(1-( \mathbf{n}^*\cdot\mathbf{h} ))^{5}}{4(\mathbf{n}^* \cdot \boldsymbol{\omega}_i)}\frac{2(\mathbf{n}^* \cdot \boldsymbol{\omega}_i)}{(\mathbf{n}^* \cdot \boldsymbol{\omega}_i)(2-\alpha)+\alpha}D(\mathbf{n}^*,\mathbf{h})\\ =&\;\dfrac{F_{0}+\left(1-F_{0}\right)(1-( \mathbf{n}^*\cdot\mathbf{h} ))^{5}}{2(\mathbf{n}^* \cdot \boldsymbol{\omega}_i)(2-\alpha)+2\alpha}D(\mathbf{n}^*,\mathbf{h})\;\\ =&\;C(\boldsymbol{\omega}_i,\mathbf{n}^*)D(\mathbf{n}^*,\mathbf{h}), \end{aligned}$$

where

$$\mathbf{h}=\frac{\boldsymbol{\omega}_i+\mathbf{n}^*}{\Vert\boldsymbol{\omega}_i+\mathbf{n}^*\Vert},$$

$C(\boldsymbol{\omega}_i,\mathbf{n}^*)$ is defined as

$$C(\boldsymbol{\omega}_i,\mathbf{n}^*)=\dfrac{F_{0}+\left(1-F_{0}\right)(1-( \mathbf{n}^*\cdot\mathbf{h} ))^{5}}{2(\mathbf{n}^* \cdot \boldsymbol{\omega}_i)(2-\alpha)+2\alpha}.$$

Since

$$\int_{H^2_{\mathbf{n}^*}}D(\mathbf{n}^*,\mathbf{h})\cos\theta_{\mathbf{h}} \:d\Omega_{\mathbf{h}}=1,$$

where $\cos\theta_{\mathbf{h}}=\mathbf{n}^* \cdot \mathbf{h}$, we can use Monte Carlo integration to importance sample

$$pdf(\mathbf{h})=D(\mathbf{n}^*,\mathbf{h})\cos\theta_{\mathbf{h}}=\dfrac{\alpha^{2}\cos\theta_{\mathbf{h}}}{\pi\left(\left(\alpha^{2}-1\right)\cos^2\theta_{\mathbf{h}}+1\right)^{2}}$$

Noting that

$$4 \cos \theta_{\mathbf{h}}d\Omega_\mathbf{h}=d\Omega_i,$$

we have

$$\begin{aligned} &\hspace{1.35em}\frac{\int_{H^2}f_s(\mathbf{n},\boldsymbol{\omega}_i, \boldsymbol{\omega}_o) L_{i}(\boldsymbol{\omega}_i)(\mathbf{n} \cdot \boldsymbol{\omega}_i) d\Omega_i}{\int_{H^2}f_s(\mathbf{n},\boldsymbol{\omega}_i, \boldsymbol{\omega}_o)\mathbf{n}^* \cdot \boldsymbol{\omega}_i d\Omega_i}\\ &=\frac{\int_{H^2}f_s(\mathbf{n},\boldsymbol{\omega}_i, \boldsymbol{\omega}_o) L_{i}(\boldsymbol{\omega}_i)(\mathbf{n}\cdot \boldsymbol{\omega}_i) 4 \cos \theta_{\mathbf{h}}d\Omega_{\mathbf{h}}}{\int_{H^2}f_s(\mathbf{n},\boldsymbol{\omega}_i, \boldsymbol{\omega}_o)\mathbf{n}\cdot \boldsymbol{\omega}_i 4 \cos \theta_{\mathbf{h}}d\Omega_{\mathbf{h}}}\\ &\approx\frac{\sum_{k=1}^Nf_s^*(\boldsymbol{\omega}_i^{(k)}, \boldsymbol{\omega}_o) L_{i}(\boldsymbol{\omega}_i^{(k)})(\mathbf{n}^* \cdot \boldsymbol{\omega}_i^{(k)}) 4 \cos \theta_{\mathbf{h}^{(k)}}/pdf(\mathbf{h}^{(k)})}{\sum_{k=1}^Nf_s^*(\boldsymbol{\omega}_i^{(k)}, \boldsymbol{\omega}_o)\mathbf{n}^* \cdot \boldsymbol{\omega}_i^{(k)} 4 \cos \theta_{\mathbf{h}^{(k)}}/pdf(\mathbf{h}^{(k)})}\\ &=\frac{\sum_{k=1}^NC(\boldsymbol{\omega}_i^{(k)},\mathbf{n}^*)L_{i}(\boldsymbol{\omega}_i^{(k)})\mathbf{n}^* \cdot \boldsymbol{\omega}_i^{(k)}}{\sum_{k=1}^NC(\boldsymbol{\omega}_i^{(k)},\mathbf{n}^*)\mathbf{n}^* \cdot \boldsymbol{\omega}_i^{(k)}},\\ \end{aligned}$$

where

$$\mathbf{n}^* \cdot \boldsymbol{\omega}_i^{(k)}=2\cos^2\theta_{\mathbf{h}^{(k)}}-1=2(\mathbf{n}^* \cdot\mathbf{h}^{(k)})^2-1,$$

while $\mathbf{h}^{(k)}$ is a random variable on the hemisphere $H^2_{\mathbf{n}^*}$, following the probability distribution with probability density function

$$pdf(\mathbf{h})=\dfrac{\alpha^{2}\cos\theta_{\mathbf{h}}}{\pi\left(\left(\alpha^{2}-1\right)\cos^2\theta_{\mathbf{h}}+1\right)^{2}}$$

Sampling Scheme

To apply the inverse transform sampling method, first parameterize $\mathbf{h}^{(k)}$ in the tangent space at $\mathbf{n}^*$ with spherical coordinates $(\theta_{\mathbf{h}^{(k)}},\varphi_{\mathbf{h}^{(k)}})\in[0,\pi/2]\times[0,2\pi]$. Let $\xi_{1}^{(k)}, \xi_{2}^{(k)}$ be independent and follow the uniform distribution $U(0,1)$. From

$$cdf_{\theta,\varphi}\left(\theta_{\mathbf{h}^{(k)}},\varphi_{\mathbf{h}^{(k)}}\right)=\int_{0}^{\varphi_{\mathbf{h}^{(k)}}} \mathrm{d} \varphi_{\mathbf{h}} \int_{0}^{\theta_{\mathbf{h}^{(k)}}} \frac{\alpha^{2} \cos \theta_{\mathbf{h}} \sin \theta_{\mathbf{h}} }{\pi\left(\left(\alpha^{2}-1\right) \cos ^{2} \theta_{\mathbf{h}} +1\right)^{2}} d\theta_{\mathbf{h}}$$

we know $\theta_{\mathbf{h}^{(k)}},\varphi_{\mathbf{h}^{(k)}}$ are independent. So we can derive the marginal distributions of $\theta_{\mathbf{h}^{(k)}},\varphi_{\mathbf{h}^{(k)}}$ respectively

$$\begin{aligned} cdf_{\theta}\left(\theta_{\mathbf{h}^{(k)}}\right)&=cdf_{\theta,\varphi}\left(\theta_{\mathbf{h}^{(k)}},2\pi\right)\\ &=\int_{0}^{2\pi} \mathrm{d} \varphi_{\mathbf{h}} \int_{0}^{\theta_{\mathbf{h}^{(k)}}} \frac{\alpha^{2} \cos \theta_{\mathbf{h}} \sin \theta_{\mathbf{h}} }{\pi\left(\left(\alpha^{2}-1\right) \cos ^{2} \theta_{\mathbf{h}} +1\right)^{2}} d\theta_{\mathbf{h}}\\ &=\int_{0}^{\theta_{\mathbf{h}^{(k)}}} \frac{2\alpha^{2} \cos \theta_{\mathbf{h}} \sin \theta_{\mathbf{h}} }{\left(\left(\alpha^{2}-1\right) \cos ^{2} \theta_{\mathbf{h}} +1\right)^{2}} d\theta_{\mathbf{h}}\\ &=\left.\frac{\alpha^{2}}{\alpha^{2}-1} \frac{1}{\left(\alpha^{2}-1\right) \cos^2\theta_{\mathbf{h}}+1}\right|_{0} ^{\theta_{\mathbf{h}^{(k)}}}\\ &=\frac{\alpha^{2}}{\alpha^{2}-1} \frac{1}{\left(\alpha^{2}-1\right)\cos^2\theta_{\mathbf{h}^{(k)}}+1}-\frac{1}{\alpha^{2}-1} \\ &=\frac{1}{\alpha^{2}-1} \frac{\alpha^2-1-\left(\alpha^{2}-1\right)\cos^2\theta_{\mathbf{h}^{(k)}}}{\left(\alpha^{2}-1\right)\cos^2\theta_{\mathbf{h}^{(k)}}+1} \\ &= \frac{1-\cos^2\theta_{\mathbf{h}^{(k)}}}{\left(\alpha^{2}-1\right)\cos^2\theta_{\mathbf{h}^{(k)}}+1}, \\ \end{aligned}$$ $$\begin{aligned} cdf_{\varphi}\left(\varphi_{\mathbf{h}^{(k)}}\right)&=cdf_{\theta,\varphi}\left(\frac{\pi}{2},\varphi_{\mathbf{h}^{(k)}}\right)\\ &=\int_{0}^{\varphi_{\mathbf{h}^{(k)}}} \mathrm{d} \varphi_{\mathbf{h}} \int_{0}^{\frac{\pi}{2}} \frac{\alpha^{2} \cos \theta_{\mathbf{h}} \sin \theta_{\mathbf{h}} }{\pi\left(\left(\alpha^{2}-1\right) \cos ^{2} \theta_{\mathbf{h}} +1\right)^{2}} d\theta_{\mathbf{h}}\\ &=\frac{1}{2\pi}\varphi_{\mathbf{h}^{(k)}}. \end{aligned}$$

Using the inverse transform of the cdf, we can obtain the sampling scheme of $\mathbf{h}^{(k)}$ in spherical coordinates

$$\begin{aligned} \theta_{\mathbf{h}^{(k)}}&=\arccos\sqrt{\frac{1-\xi_{1}^{(k)}}{1-(1-\alpha^2)\xi_{1}^{(k)}}},\\ \varphi_{\mathbf{h}^{(k)}}&=2\pi\xi_{2}^{(k)}. \end{aligned}$$

From this, we can derive the sampling scheme of $\boldsymbol{\omega}_i^{(k)}$ in spherical coordinates

$$\begin{aligned} \theta_{i}^{(k)}&=\arccos\frac{1-(1+\alpha^2)\xi_{1}^{(k)}}{1-(1-\alpha^2)\xi_{1}^{(k)}},\\ \varphi_{i}^{(k)}&=2\pi\xi_{2}^{(k)}. \end{aligned}$$

Since the above spherical coordinates are generated in the tangent space, we need to first convert them to Cartesian coordinates, and then use coordinate transformation to convert them to world space coordinates. Using the previously defined $T_{\mathbf{n}}$, the sampling scheme of $\boldsymbol{\omega}_i^{(k)}$ in world space can be written as

$$\boldsymbol{\omega}_i^{(k)}=T_{\mathbf{n}} \begin{bmatrix} \sin\theta_{i}^{(k)}\cos\varphi_{i}^{(k)}\\\\ \sin\theta_{i}^{(k)}\sin\varphi_{i}^{(k)}\\\\ \cos\theta_{i}^{(k)} \end{bmatrix}=T_{\mathbf{n}} \begin{bmatrix} \dfrac{2\alpha\cos2\pi\xi_{2}^{(k)}\sqrt{\xi_{1}^{(k)}\left(1-\xi_{1}^{(k)}\right)}}{1-(1-\alpha^2)\xi_{1}^{(k)}}\\ \\ \dfrac{2\alpha\sin2\pi\xi_{2}^{(k)}\sqrt{\xi_{1}^{(k)}\left(1-\xi_{1}^{(k)}\right)}}{1-(1-\alpha^2)\xi_{1}^{(k)}}\\ \\ \dfrac{1-(1+\alpha^2)\xi_{1}^{(k)}}{1-(1-\alpha^2)\xi_{1}^{(k)}} \end{bmatrix}.$$

Finally, we can derive the estimate for the first part of the split integral

$$S(\mathbf{n}^*,\alpha)=\frac{\sum_{k=1}^NC(\boldsymbol{\omega}_i^{(k)},\mathbf{n}^*)L_{i}(\boldsymbol{\omega}_i^{(k)})\mathbf{n}^* \cdot \boldsymbol{\omega}_i^{(k)}}{\sum_{k=1}^NC(\boldsymbol{\omega}_i^{(k)},\mathbf{n}^*)\mathbf{n}^*\cdot \boldsymbol{\omega}_i^{(k)}},$$

where

$$\begin{aligned} \mathbf{n}^*&=2(\boldsymbol{\omega}_o\cdot\mathbf{n})\mathbf{n}-\boldsymbol{\omega}_o,\\ C(\boldsymbol{\omega}_i^{(k)},\mathbf{n}^*)&=\dfrac{F_{0}+\left(1-F_{0}\right)(1-( \mathbf{n}^*\cdot\mathbf{h} ))^{5}}{2(\mathbf{n}^* \cdot \boldsymbol{\omega}_i^{(k)})(2-\alpha)+2\alpha}. \end{aligned}$$

We can precompute $S(\mathbf{n}^{*},\alpha)$ for some $\alpha$, generating multiple cubemap mip levels like a mipmap. For other $\alpha$, the value of $S(\mathbf{n}^*,\alpha)$ can be obtained by trilinear interpolation. During rendering, the first part of the split integral can be estimated as

$$\frac{\int_{H^2}f_s(\mathbf{n},\boldsymbol{\omega}_i, \boldsymbol{\omega}_o) L_{i}(\boldsymbol{\omega}_i)(\mathbf{n}^* \cdot \boldsymbol{\omega}_i) d\Omega_i}{\int_{H^2}f_s(\mathbf{n},\boldsymbol{\omega}_i, \boldsymbol{\omega}_o)\mathbf{n}^* \cdot \boldsymbol{\omega}_i d\Omega_i}\approx S(\mathbf{r}(\boldsymbol{\omega}_o),\alpha)=S(2(\boldsymbol{\omega}_o\cdot\mathbf{n})\mathbf{n}-\boldsymbol{\omega}_o$$

Split Sum - Part 2

Formula Derivation

The right-hand side of the split integral can be simplified to

$$\begin{aligned} &\; \int_{H^2} f_s(\mathbf{n},\boldsymbol{\omega}_i, \boldsymbol{\omega}_o)\mathbf{n} \cdot \boldsymbol{\omega}_i d\Omega_i\\ =&\;\int_{H^2}\left(F_{0}+\left(1-F_{0}\right)(1-( \mathbf{n}\cdot\mathbf{h} ))^{5}\right) \frac{f_s(\mathbf{n},\boldsymbol{\omega}_i, \boldsymbol{\omega}_o)}{F\left(\mathbf{h}, \boldsymbol{\omega}_o\right)}\mathbf{n} \cdot \boldsymbol{\omega}_i d\Omega_i\\ =&\; F_{0}\int_{H^2} \frac{f_s(\mathbf{n},\boldsymbol{\omega}_i, \boldsymbol{\omega}_o)}{F\left(\mathbf{h}, \boldsymbol{\omega}_o\right)}\mathbf{n} \cdot \boldsymbol{\omega}_i \left(1-(1-( \mathbf{n}\cdot\mathbf{h} ))^{5}\right)d\Omega_i\\&+\int_{H^2} \frac{f_s(\mathbf{n},\boldsymbol{\omega}_i, \boldsymbol{\omega}_o)}{F\left(\mathbf{h}, \boldsymbol{\omega}_o\right)}\mathbf{n} \cdot \boldsymbol{\omega}_i(1-( \mathbf{n}\cdot\mathbf{h} ))^{5} d\Omega_i. \end{aligned}$$

Importance sampling is given by

$$\begin{aligned} A &\approx \frac{1}{N} \sum_{k=1}^{N} \frac{\frac{D\left(\boldsymbol{\omega}_{h}^{(k)}\right) G\left(\boldsymbol{\omega}_{o}, \boldsymbol{\omega}_{i}^{(k)}\right)}{4\left(\boldsymbol{\omega}_{o} \cdot \mathbf{n}\right)\left(\boldsymbol{\omega}_{i}^{(k)} \cdot \mathbf{n}\right)}\left(1-\left(1-\boldsymbol{\omega}_{o} \cdot \boldsymbol{\omega}_{h}^{(k)}\right)^{5}\right)\left(\boldsymbol{\omega}_{i}^{(k)} \cdot \mathbf{n}\right)}{\mathrm{pdf}\left(\boldsymbol{\omega}_{i}^{(k)}\right)}\\ &=\frac{1}{N} \sum_{k=1}^{N} \frac{G\left(\boldsymbol{\omega}_{o}, \boldsymbol{\omega}_{i}^{(k)}\right)\left(\boldsymbol{\omega}_{o} \cdot \boldsymbol{\omega}_{h}^{(k)}\right)\left(1-\left(1-\boldsymbol{\omega}_{o} \cdot \boldsymbol{\omega}_{h}^{(k)}\right)^{5}\right)}{\left(\boldsymbol{\omega}_{o} \cdot \mathbf{n}\right)\left(\boldsymbol{\omega}_{h}^{(k)} \cdot \mathbf{n}\right)}, \end{aligned}$$ $$B \approx \frac{1}{N} \sum_{k=1}^{N} \frac{G\left(\boldsymbol{\omega}_{o}, \boldsymbol{\omega}_{i}^{(k)}\right)\left(\boldsymbol{\omega}_{o} \cdot \boldsymbol{\omega}_{h}^{(k)}\right)\left(1-\boldsymbol{\omega}_{o} \cdot \boldsymbol{\omega}_{h}^{(k)}\right)^{5}}{\left(\boldsymbol{\omega}_{o} \cdot \mathbf{n}\right)\left(\boldsymbol{\omega}_{h}^{(k)} \cdot \mathbf{n}\right)}.$$