Linear-Time Approximate Spherical Gaussian Filtering Originally published on November 29, 2017

Physically based rendering using LTASG-filtered environment maps.

Blinn-Phong approximates a Gaussian distribution as the specular exponent increases [Lyon1993]. [Olano2010] has shown the following relation in terms of the angle $\theta$ between $n$ and $h$:

[Lyon1993] Lyron, R. 1993. Phong shading reformulation for hardware rendering. Tech. Rep. 43, Apple.

[Olano2010] Olano, M., & Baker, D. (2010, February). LEAN mapping. In Proceedings of the 2010 ACM SIGGRAPH symposium on Interactive 3D Graphics and Games (pp. 181-188). ACM.

$\cos^{s}(\theta)\approx\exp\left(-\frac{s}{2}\tan^{2}\theta\right)$

This makes the spherical Gaussian blur an excellent and appropriate choice for generating environmental maps. By extending the Gaussian blur's separable property, it is possible to implement it in $O(K)$ (where $K$ is a kernel radius) for reasonably small $\sigma$ values.

# Related Works

## Pre-filtered Mipmapped Radiance Environment Maps

AMD's CubeMapGen has been a popular choice to generate pre-filtered environment maps. However, it is designed for a offline generation and is too slow for a real-time use.

In three.js, PMREMGenerator is responsible for pre-filtering environment maps. It is implemented as a fragment shader that performs Monte Carlo sampling, and when the sample count is set as low as 32 it is capable of running at 60fps on Intel HD Graphics 400011. [Colbert2007] describes a practical implementation of GPU-based importance sampling for environment map pre-filtering

[Colbert2007] Colbert, M., & Krivanek, J. (2007). GPU-based importance sampling. GPU Gems, 3, 459-476.

## Mapping Gloss Values to Mip Levels

I wanted to have a constant kernel radius of $K$ (= 8) pixels for every mip level. $\sigma$ should be a somewhat smaller value than $K$ in order to fit the significant part of the Gaussian distribution within the kernel radius. I chose $\sigma=K/r$ where $r=4$.

Under this condition and given that the image size of the base mip level is $N$ pixels, the relationship between the specular exponent $s$ and the mip level $n$ is found as following:

\begin{aligned} \sigma=1/\sqrt{s} &=\frac{K}{r}\cdot\frac{1}{2^{-n}N} \\ s &=\left(\frac{2^{-n}Nr}{K}\right)^{2} \\ &=0.25(2^{-n}N)^{2} \\ n &=\frac{1}{2}\log_{2}4sN^{2} \end{aligned}

# Algorithm

## Separate Filtering

The basic idea of the separate Gaussian filter is decomposing a n-dimensional Gaussian filter into $n$ cascaded one-dimensional Gaussian filters as shown in the following example where $n=2$:

\begin{aligned} G(x,y) & =\frac{1}{2\pi\sigma^{2}}\exp\left(-\frac{x^{2}+y^{2}}{2\sigma^{2}}\right)\\ G_{x}(x,y) & =\begin{cases} \frac{1}{2\pi\sigma^{2}}\exp\left(-\frac{x^{2}}{2\sigma^{2}}\right) & y=0\\ 0 & y\ne0 \end{cases}\\ G_{y}(x,y) & =\text{ditto.}\\ G & =G_{x}\circ G_{y} \end{aligned}

This decomposition allows a $n$-dimensional Gaussian filter to be implemented with the time complexity $O(K)$ instead of $O(K^{n})$.

At cost of accuracy, this idea can be extended for a wider variety of filters that locally resemble a Gaussian filter, examples of which include a spatially varying anisotropic Gaussian filter [Zheng2011].

[Zheng2011] Zheng, Z., & Saito, S. (2011, August). Screen space anisotropic blurred soft shadows. In SIGGRAPH Posters (p. 75).

To apply this technique, one has to find the functions $A_{1}(\vec{x}),\ldots,A_{k}(\vec{x})$ each of which define the axis direction and the standard deviation of the corresponding one-dimensional Gaussian filter. Note that $\vec{x}$ represents a point in a $n$-manifold $\Gamma$ embedded in a Euclidean space, and $A_{i}(\vec{x})$ must be a tangent vector of $\Gamma$ at $\vec{x}$. The axis functions must fulfill the following condition in order for the resulting filter to locally resemble a $n$-dimensional Gaussian filter:

$\mathrm{rank}(A_{1}(\vec{x})\ \cdots\ A_{k}(\vec{x}))\ge n$

In addition, from a practical perspective, $A_{1}(\vec{x}),\ldots,A_{k}(\vec{x})$ must be as smooth as possible because abrupt changes in them lead to visual artifacts.

For a spherical Gaussian blur ($\Gamma=S^{2}$, $n=2$), there exists no $A_{1}(\vec{x}),A_{2}(\vec{x})$ that satisfies this condition on every $\vec{x}\in\Gamma$, which is obvious from the "hairy ball theorem" stating that there exists no nonvanishing continuous tangent vector field on even-dimensional $n$-spheres. Therefore, at least 3 axis functions are required to realize a spherical Gaussian blur using this technique.

I propose the following axis functions ($\left\{ \vec{a_{1}},\vec{a_{2}},\vec{a_{3}}\right\}$ is an orthonormal basis of $\mathbb{R}^{3}$):

\begin{aligned} A_{1}(\vec{x}) & =\sigma(\vec{a_{1}}-\vec{x}(\vec{x}\cdot\vec{a_{1}}))\\ A_{2}(\vec{x}) & =\sigma(\vec{a_{2}}-\vec{x}(\vec{x}\cdot\vec{a_{2}}))\\ A_{3}(\vec{x}) & =\sigma(\vec{a_{3}}-\vec{x}(\vec{x}\cdot\vec{a_{3}})) \end{aligned}

Each of them represents a tangent vector along the latitude, assuming the points $\pm\vec{a_{i}}$ are the north and south poles of the sphere. If $\left\{ \vec{a_{1}},\vec{a_{2}},\vec{a_{3}}\right\}$ is substituted with the standard basis, they can be written more neatly as:

\begin{aligned} A_{1}(\vec{x}) & =\sigma(\vec{e_{x}}-x_{x}\vec{x})\\ A_{2}(\vec{x}) & =\sigma(\vec{e_{y}}-x_{y}\vec{x})\\ A_{3}(\vec{x}) & =\sigma(\vec{e_{z}}-x_{z}\vec{x})\mathbf{} \end{aligned}

These axis functions are unambiguously derived by combining the following conditions: the tangential condition (they follow the surface of a sphere), the uniform blur condition, and the latitudinal condition (each of them follows a set of latitudinal lines surrounding a corresponding axis).

## Implementation on Cube Map

For each one-dimensional filter ($i\in\{1,2,3\}$) and each cube face, there are two cases to handle:

1. $\pm\vec{a_{i}}$ is inside the face — In this case, the filter is implemented as a radial blur oriented toward the pole $\pm\vec{a_{i}}$.

2. $\pm\vec{a_{i}}$ is outside the face — In this case, the filter is implemented as a directional blur along the U or V direction.

We will only consider the positive Z cube face in the following discussion.

Given a texture coordinate $(u,v)$, the corresponding point $\vec{x}\in S^{2}$ is found as:

$\vec{x}=\frac{1}{\sqrt{1+u^{2}+v^{2}}}\begin{pmatrix}u\\ v\\ 1 \end{pmatrix}$

In the first case where $\pm\vec{a_{i}}$ is inside the face (hence $\vec{a_{i}}=\vec{e_{z}}$)

$A_{i}(\vec{x})=\sigma\begin{pmatrix}-\frac{u}{\sqrt{1+u^{2}+v^{2}}}\\ -\frac{v}{\sqrt{1+u^{2}+v^{2}}}\\ 1-\frac{1}{1+u^{2}+v^{2}} \end{pmatrix}$

By projecting it on the plane $z=1$ we obtain:

$\left.\frac{d}{dt}\frac{\vec{x}+A_{i}(\vec{x})\cdot t}{\vec{e_{z}}\cdot\left(\vec{x}+A_{i}(\vec{x})\cdot t\right)}\right|_{t=0}=\begin{pmatrix}-u\sigma\sqrt{1+u^{2}+v^{2}}\\ -v\sigma\sqrt{1+u^{2}+v^{2}}\\ 0 \end{pmatrix}$

In the second case where $\pm\vec{a_{i}}$ is inside the face, assuming $\vec{a_{i}}=\vec{e_{x}}$

$A_{i}(\vec{x})=\sigma\begin{pmatrix}1-\frac{u^{2}}{1+u^{2}+v^{2}}\\ -\frac{uv}{1+u^{2}+v^{2}}\\ -\frac{u}{1+u^{2}+v^{2}} \end{pmatrix}$

By projecting it on the plane $z=1$ we obtain:

$\left.\frac{d}{dt}\frac{\vec{x}+A_{i}(\vec{x})\cdot t}{\vec{e_{z}}\cdot\left(\vec{x}+A_{i}(\vec{x})\cdot t\right)}\right|_{t=0}=\begin{pmatrix}\sigma\sqrt{1+u^{2}+v^{2}}\\ 0\\ 0 \end{pmatrix}$

# Demonstration

hyper3d-envmapgen implements the proposed algorithm using TypeScript, Rust, and WebAssembly. It uses a CPU to execute the algorithm, making it possible to generate a PMREM without hindering the main thread or a GPU. An interactive WebGL demo is available here.

# References

[Lyon1993] Lyron, R. 1993. Phong shading reformulation for hardware rendering. Tech. Rep. 43, Apple.

[Olano2010] Olano, M., & Baker, D. (2010, February). LEAN mapping. In Proceedings of the 2010 ACM SIGGRAPH symposium on Interactive 3D Graphics and Games (pp. 181-188). ACM.

[Zheng2011] Zheng, Z., & Saito, S. (2011, August). Screen space anisotropic blurred soft shadows. In SIGGRAPH Posters (p. 75).

[Colbert2007] Colbert, M., & Krivanek, J. (2007). GPU-based importance sampling. GPU Gems, 3, 459-476.