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 nn and hh:

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

coss(θ)exp(s2tan2θ)\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)O(K) (where KK 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 KK (= 8) pixels for every mip level. σ\sigma should be a somewhat smaller value than KK in order to fit the significant part of the Gaussian distribution within the kernel radius. I chose σ=K/r\sigma=K/r where r=4r=4.

Under this condition and given that the image size of the base mip level is NN pixels, the relationship between the specular exponent ss and the mip level nn is found as following:

σ=1/s=Kr12nNs=(2nNrK)2=0.25(2nN)2n=12log24sN2\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 nn cascaded one-dimensional Gaussian filters as shown in the following example where n=2n=2:

G(x,y)=12πσ2exp(x2+y22σ2)Gx(x,y)={12πσ2exp(x22σ2)y=00y0Gy(x,y)=ditto.G=GxGy\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 nn-dimensional Gaussian filter to be implemented with the time complexity O(K)O(K) instead of O(Kn)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 A1(x),,Ak(x)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 x\vec{x} represents a point in a nn-manifold Γ\Gamma embedded in a Euclidean space, and Ai(x)A_{i}(\vec{x}) must be a tangent vector of Γ\Gamma at x\vec{x}. The axis functions must fulfill the following condition in order for the resulting filter to locally resemble a nn-dimensional Gaussian filter:

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

In addition, from a practical perspective, A1(x),,Ak(x)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 (Γ=S2\Gamma=S^{2}, n=2n=2), there exists no A1(x),A2(x)A_{1}(\vec{x}),A_{2}(\vec{x}) that satisfies this condition on every xΓ\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 nn-spheres. Therefore, at least 3 axis functions are required to realize a spherical Gaussian blur using this technique.

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

A1(x)=σ(a1x(xa1))A2(x)=σ(a2x(xa2))A3(x)=σ(a3x(xa3))\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 ±ai\pm\vec{a_{i}} are the north and south poles of the sphere. If {a1,a2,a3}\left\{ \vec{a_{1}},\vec{a_{2}},\vec{a_{3}}\right\} is substituted with the standard basis, they can be written more neatly as:

A1(x)=σ(exxxx)A2(x)=σ(eyxyx)A3(x)=σ(ezxzx)\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{1,2,3}i\in\{1,2,3\}) and each cube face, there are two cases to handle:

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

  2. ±ai\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)(u,v), the corresponding point xS2\vec{x}\in S^{2} is found as:

x=11+u2+v2(uv1) \vec{x}=\frac{1}{\sqrt{1+u^{2}+v^{2}}}\begin{pmatrix}u\\ v\\ 1 \end{pmatrix}

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

Ai(x)=σ(u1+u2+v2v1+u2+v2111+u2+v2) 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=1z=1 we obtain:

ddtx+Ai(x)tez(x+Ai(x)t)t=0=(uσ1+u2+v2vσ1+u2+v20) \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 ±ai\pm\vec{a_{i}} is inside the face, assuming ai=ex\vec{a_{i}}=\vec{e_{x}}

Ai(x)=σ(1u21+u2+v2uv1+u2+v2u1+u2+v2) 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=1z=1 we obtain:

ddtx+Ai(x)tez(x+Ai(x)t)t=0=(σ1+u2+v200) \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.