基于预计算的实时全局光照(Real-time Global Illumination Based On Precomputation)

2021年9月30日 5点热度 0条评论 来源: KillerAery

目录

基于图像的光照(Image Based Lighting,IBL)

基于图像的光照(IBL),简单的说就是一类通过 环境贴图(Environment Map) 来保存某个物体的环境光信息,从而实现该物体的基于物理的物体渲染方法。

  • IBL 可用于 diffuse/glossy/specular 物体(基本上包含大部分物体了)渲染:这取决于环境贴图,环境贴图分辨率越大,那么所能表示高频信息就越多,从而越适合 specular 物体的渲染;而分辨率越低,所表示的信息更多是低频信息,则 diffuse 物体的渲染也足以满足,而且还能节省一定存储空间。

  • IBL 可以实现成动态环境光照:实时渲染出动态的环境贴图。有一定开销,一般只用于少量的specular物体。

  • IBL 可以实现成静态环境光照:预渲染环境贴图。需要环境光是静态的。

用 IBL 实现静态环境光照就有点类似 Light Map的做法了,区别在于 IBL 存的是静态环境光信息,Light Map 存的是接受静态环境光信息后的着色结果。

为了表示来自四面八方的环境光信息,IBL 使用 Spherical Map 或者 Cube Map 的方式来存储:


The Split Sum Approximation

IBL 中最常见的算法便是基于 The Split Sum Approximation 的算法。

我们知道,一般的渲染方程如下:

\(L_{r}\left(x, \omega_{r}\right)=\int_{\Omega^{+}} L_{i}\left(x, \omega_{i}\right) f_{r}\left(x, \omega_{i}, \omega_{r}\right)\left(n \cdot \omega_{i}\right) \mathrm{d} \omega_{i}\)

在实际渲染的时候,我们当然可以使用蒙特卡洛方法实现该渲染方程,然而这样的开销是巨大的(每个shading point都要做多重采样,而且结果很容易是noisy的)。

为了进一步加快式子,这里有一个很经典的近似公式:

\(\int_{\Omega^+} f(x) g(x) \mathrm{d} x \approx \frac{\int_{\Omega_{G}} f(x) \mathrm{d} x}{\int_{\Omega_{G}} \mathrm{~d} x} \cdot \int_{\Omega^+} g(x) \mathrm{d} x\)

积分域 \(\Omega_G\) :在原本 \(\Omega^+\) 的积分域范围内剔除掉 \(f(x) = 0\) 的地方而剩余的范围。

想让这个公式近似效果比较精确,那么需要满足以下一种或两种条件:

  • 积分域 \(\Omega_G\) 比较小
  • \(g(x)\) 比较光滑,即变化不是很大

而我们的观察是:

  • 如果 BRDF 是 glossy/specular 的,那么它的 lobe 往往是花瓣状,即只有很小的积分域才能接受环境光。
  • 如果 BRDF 是 diffuse 的,那么它的 lobe 往往是均匀的半球状,即无论哪个方向的环境光打进来, \(f_r\) 函数的输出几乎没多少变化(甚至是个常数)。

于是基于以上理论,Split Sum 方法对渲染方程改造成这样的近似公式来获得渲染的加速:

\(L_{r}\left(x, \omega_{r}\right) \approx \frac{\int_{\Omega_{f_{r}}} L_{i}\left(x, \omega_{i}\right) \mathrm{d} \omega_{i}}{\int_{\Omega_{f_{r}}} \mathrm{~d} \omega_{i}} \cdot \int_{\Omega^{+}} f_{r}\left(x, \omega_{i}, \omega_{r}\right) \cos \theta_{i} \mathrm{~d} \omega_{i}\)

对渲染方程拆分成两个部分(环境光积分、BRDF积分)后就可以通过预计算的方式(后面两节会介绍如何预计算)分别减少这些积分的运行时开销,总结这种方法的好处是:

  • Split Sum 方法和原始蒙特卡洛方法的图像效果几乎一模一样
  • 由于不用对环境贴图进行多重采样,性能开销大大减低了

过滤环境贴图

环境光积分:\(\frac{\int_{\Omega_{f_{r}}} L_{i}\left(x, \omega_{i}\right) \mathrm{d} \omega_{i}}{\int_{\Omega_{f_{r}}} \mathrm{~d} \omega_{i}}\)

因为我们已经拥有一张环境贴图(无论是实时的还是预渲染的)来存储环境光信息了,为了计算环境光部分的积分,需要在 \(\Omega_{fr}\) 范围内做多次光线采样。但是,可以有一个几乎等价但避免运行时多重采样的方式:

预先对纹理进行滤波操作(模糊),我们只需要对滤波后的环境贴图采样一次光线方向就能得到积分。

滤波后的环境贴图实际上称为 Irradiance Environment Map :原来环境贴图中的单位是 radiance \(L(\mathrm{x}, \omega)\),贴图积分后的单位则变成了 Irradiance \(E(x)\)

当然,BRDF Lobe 的形状越尖锐,即环境光积分范围越小,就需要使用模糊程度更低的环境贴图;反之 BRDF Lobe 的形状越粗壮,即环境光积分范围越大,就需要使用模糊程度更高的环境贴图。

我们可以使用 MIPMAP技术 来生成不同Level的环境贴图,通过三线性插值(Trilinear Interportion)的方式来得出任何模糊程度且任何2D位置的环境光滤波结果。

预计算BRDF积分

回顾 Microfacet BRDF 的组成:

  • Fresnel项(举例 Schlick's approximation)

\(F=F_0 +(1-F_0)(1-(h \cdot \omega_{r}))^{5} = F_0 +(1-F_0)(1-(\cos\theta_{vh}))^{5}\)

  • NDF项(举例 Beckmann NDF)

\(D(h)=\frac{1}{\pi \alpha^{2}({n} \cdot {h})^{4}} \cdot \exp {\left(\frac{({n} \cdot {h})^{2}-1}{a^{2}({n} \cdot {h})^{2}}\right)} = \frac{1}{\pi \alpha^{2}\cos^4 \theta_h } \cdot \exp {\left(\frac{\cos^2 \theta_h -1}{a^{2}\cos^2 \theta_h }\right)}\)

可以预想到,该BRDF 的积分结果依赖于三个参数:

  • \(F_0\)(Fresnel项系数)
  • \(\alpha\) (粗糙度 roughness)
  • \(\theta_v\) (反射方向与法线的夹角,实际上决定了 \(\theta_{vh}\)\(\theta_{h}\)

我们可以把 \(F_0\) 项拆出来,让BRDF 积分拆成两个积分,但是这些积分都减少了一个依赖的参数:

\(\begin{align} \int_{\Omega^{+}} f_{r}\left(p, \omega_{i}, \omega_{o}\right) \cos \theta_{i} \mathrm{~d} \omega_{i}\approx & F_{0} \int_{\Omega^{+}} \frac{f_{r}}{F}\left(1-\left(1-\cos \theta_{vh}\right)^{5}\right) \cos \theta_{i}\mathrm{~d} \omega_{i} \\ & + \int_{\Omega^{+}} \frac{f_{r}}{F}\left(1-\cos \theta_{vh}\right)^{5} \cos \theta_{i} \mathrm{~d}\omega_{i} \end{align}\)

这样,对于相同的材质(相同的BRDF),我们就可以针对剩余两个依赖的参数 \(\alpha\)\(\theta_v\) 建立一张 二维查询表

预计算辐射度传输(Precomputed Radiance Transfer,PRT)

预计算辐射度传输(Precomputed Radiance Transfer,PRT) 是一类通过预计算辐射度传输的基于物理的物体渲染方法。所谓辐射度传输,可以理解成物体自身的阴影/AO(AmbientOcclusion)和表面互反射等有关于光路传输的信息。因此 PRT 往往适用于动态光照下的静态物体/静态材质。

  • PRT 可以实现动态环境光照:仅预计算 transfer 部分
  • PRT 可以实现静态环境光照:不仅预计算 transfer,还顺便预计算环境光部分(可以称为 radiance);这样就能以更低性能开销实现运行时基于物理的渲染,只是光照动态性会有所限制(甚至得完全静态光照)。

顺便一提,PRT 方法在2002年 Siggraph 会议被 Peter-Pike Sloan 首先提出,并从此掀起了一波 PRT 方法的研究浪潮。当时 Peter-Pike Sloan 论文所实现的 PRT 算法便是基于球谐(SH)的,实际上 PRT 不仅可以通过 SH 去表示低频环境光场景,还可以有其它方法(如Wavelet)。不过本文仍将主要介绍球谐光照(SH Lighting)方法。

球谐(Spherical Harmonics,SH)

对于任意函数 \(f(x)\) ,不管是连续的还是不连续的,我们可以展开成一系列基函数(每项都带某个系数)的线性组合:

\(f(x)=\sum_{i} c_{i} \cdot B_{i}(x)\)

例如,多项式展开可以看成是一系列基函数(多项式)的线性组合:\(f(x)=c_{0} +c_1\cdot x^1 + c_2 \cdot x^2 + ...\)

再例如,傅立叶变换也可以将 \(f(x)\) 表示成另一系列基函数(各种频率的正弦谐波)的线性组合:

球谐(Spherical Harmonics,SH) 便是定义在球面上的一系列2D基函数,它与2D傅里叶序列有点相似,但非常适合球面函数 \(f(\omega)\)(即参数为单位球面向量)。

SH 是分阶数的:在第0阶(l=0)有1个基函数(m=0);第1阶(l=1)有3个基函数(m=-1,0,1)...第n阶(l=n)有2n+1个基函数(m=-l,...,0,...,l)。实际上,阶数越低的基函数所代表的信息就越是低频。如果要完全复原一个任意函数,我们需要无穷阶的 SH;而如果我们只要复原一个任意函数的近似(换句话说只重建出该函数的低频信息),那么我们完全可以只需要前几阶的 SH(包含 l=0,l=1,...,l=n 每一阶的所有基函数)。

SH 的实数形式的基函数如下:

SH 基函数的形式比较复杂,只建议看一眼就跳过。

\(y_{l}^{m}(\theta, \varphi)=\left\{\begin{array}{rl}\sqrt{2} \operatorname{Re}\left(Y_{l}^{m}\right) & m>0 \\ \sqrt{2} \operatorname{Im}\left(Y_{l}^{m}\right) & m<0 \\ Y_{l}^{0} & m=0\end{array}=\left\{\begin{array}{cc}\sqrt{2} K_{l}^{m} \cos m \varphi P_{l}^{m}(\cos \theta) & m>0 \\ \sqrt{2} K_{l}^{m} \sin |m| \varphi P_{l}^{|m|}(\cos \theta) & m<0 \\ K_{l}^{0} P_{l}^{0}(\cos \theta) & m=0\end{array}\right.\right.\)

其中,

\(K_{l}^{m}=\sqrt{\frac{(2 l+1)(l-|m|) !}{4 \pi(l+|m|) !}}\)

\(p^m_l\) 为勒让德多项式(Legendre polynomials):

  • $P_0^0= 1 $
  • $P_{m}^{m} =(1-2 m) P_{m-1}^{m-1} $
  • $P_{m+1}^{m} =(2 m+1) z P_{m}^{m} $
  • $ P_{l}^{m} =\frac{(2 l-1) z P_{l-1}^{m}-(l+m-1) P_{l-2}^{m}}{l-m} $

z为球面向量对应球面坐标的z值,该多项式不依赖x值y值。

以下是前3阶的SH实数形式的基函数:

\(l=0\)

  • \(Y_{00}=s=Y_{0}^{0}=\frac{1}{2} \sqrt{\frac{1}{\pi}}\)

\(l=1\)

  • \(Y_{1,-1}=p_{y}=i \sqrt{\frac{1}{2}}\left(Y_{1}^{-1}+Y_{1}^{1}\right)=\sqrt{\frac{3}{4 \pi}} \cdot \frac{y}{r}\)
  • \(Y_{1,0}=p_{z}=Y_{1}^{0}=\sqrt{\frac{3}{4 \pi}} \cdot \frac{z}{r}\)
  • \(Y_{1,1}=p_{x}=\sqrt{\frac{1}{2}}\left(Y_{1}^{-1}-Y_{1}^{1}\right)=\sqrt{\frac{3}{4 \pi}} \cdot \frac{x}{r}\)

\(l=2\)

  • \(Y_{2,-2}=d_{x y}=i \sqrt{\frac{1}{2}}\left(Y_{2}^{-2}-Y_{2}^{2}\right)=\frac{1}{2} \sqrt{\frac{15}{\pi}} \cdot \frac{x y}{r^{2}}\)
  • \(Y_{2,-1}=d_{y z}=i \sqrt{\frac{1}{2}}\left(Y_{2}^{-1}+Y_{2}^{1}\right)=\frac{1}{2} \sqrt{\frac{15}{\pi}} \cdot \frac{y z}{r^{2}}\)
  • \(Y_{2,0}=d_{z^{2}}=Y_{2}^{0}=\frac{1}{4} \sqrt{\frac{5}{\pi}} \cdot \frac{-x^{2}-y^{2}+2 z^{2}}{r^{2}}\)
  • \(Y_{2,1}=d_{s z}=\sqrt{\frac{1}{2}}\left(Y_{2}^{-1}-Y_{2}^{1}\right)=\frac{1}{2} \sqrt{\frac{15}{\pi}} \cdot \frac{z x}{r^{2}}\)
  • \(Y_{2,2}=d_{x^{2}-y^{2}}=\sqrt{\frac{1}{2}}\left(Y_{2}^{-2}+Y_{2}^{2}\right)=\frac{1}{4} \sqrt{\frac{15}{\pi}} \cdot \frac{x^{2}-y^{2}}{r^{2}}\)

SH 与一般的基函数相比,具有以下性质:

  • 基函数之间具有正交性(orthonormal)

\(\int_{\Omega} B_{i}(\omega) \cdot B_{j}(\omega) \mathrm{d} \omega=\mathbf{1} \quad(i=j)\)

\(\int_{\Omega} B_{i}(\omega) \cdot B_{j}(\omega) \mathrm{d} \omega=\mathbf{0} \quad(i\neq j)\)

  • 通过 投影(Projection) 可以很方便得到 SH 系数(SH coefficients)

\(f_i = \int_{\Omega}f(\omega)\cdot B_i(\omega)d\omega\)

  • 通过系数向量组与基函数组的点积(前n阶的基函数/系数共有 n*n 个)可以很方便重建球面函数

\(f(\omega)\approx \sum^{n^2}_{i=1} f_{i} B_{i}(\mathbf{\omega})\)

  • product projection:\(c(\omega)=a(\omega)b(\omega)\) ,已知 \(a\) 的形式而 \(b\) 未知,有

\(\vec{c} = M\cdot \vec b\)

其中,\(\vec{c}\)\(c(\omega)\) 的SH系数向量组,\(\vec{b}\)\(b(\omega)\) 的SH系数向量组,\(M\)\(n^2\) X \(n^2\) 的矩阵,其元素为

\(M_{ij}=\int_{\Omega}a(\omega)B_i(\omega)B_j(\omega)d\omega\)

这样我们就可以预计算好 \(M\) ,等到知道 \(\vec{b}\) 后,乘起来就能得到 \(\vec{c}\)

推导:

\(令 M_{i}(\omega)=a(\omega)B_i(\omega)\)

\(\begin{align} c_i &= \int_{\Omega} a(\omega)b(\omega)B_i(\omega) d\omega \\ &= \int_{\Omega} M_i(\omega)b(\omega)d\omega \\ &= \int_{\Omega} \sum_{j=1}^{n^2} (M_{ij} B_j(\omega)) b(\omega) d \omega \\ &= \sum_{i=1}^{n^2} M_{ij} \int_{\Omega} B_j(\omega)b(\omega) d \omega \\ &= \sum_{j=1}^{n^2} M_{ij}\cdot b_j \end{align}\)

  • 支持插值,对 SH 系数的插值相当于对重建的函数值的插值。

  • 旋转不变性(rotational invariance),对函数 \(f\) 的旋转 \(R_{SH}\) 等价于对 \(f(\omega)\) 的自变量的旋转 \(R_{3D}\)

\(R_{SH}(f(\omega))=f(R_{3D}(\omega))\)

球谐光照(Spherical Harmonic Lighting)

纹理本质上也是一个函数(信号),输入二维坐标,输出对应纹素的RGBA值。传统的环境光信息往往使用环境贴图(Environment Map)表示,这往往需要一个庞大的二维数组存储各个纹素的值,十分耗费空间,而且采样和纹理I/O也有一定开销。

球谐光照(SH Lighting) 的思路:将渲染方程分为两个球面函数(即 lighting function 和 transfer function),这些球面函数将使用 SH 方法来表示:

  • 对于环境光(lighting)部分,只需预先存储若干个 SH 系数而不必存储一整张环境贴图。
  • 对于传输函数(transfer function) 部分,只需预先存储若干个 SH 系数(diffuse情况下)或者矩阵(glossy情况下),预计算好传输率可以在实时渲染中以极低代价实现自阴影、互反射的效果。

不过 SH Lighting 一般采用3阶SH,因此 SH 所能表示的 lighting 和 transfer 信息是低频的,只适用于 diffuse 和 glossy 的物体而不适用于 specular 的物体。

球谐光照(Spherical Harmonic Lighting)效果:

图分别为:diffuse物体的unshadowed效果、diffuse物体的shadowed&interreflection效果、glossy物体的unshadowed效果、glossy物体的shadowed&interreflection效果

Diffuse 物体的球谐光照

物体的渲染方程:

\(L(r)=\int_{\Omega^+} L(\mathbf{\omega})\cdot \rho \cdot V(\mathbf{\omega}) \max (0, n\cdot \mathbf{\omega}) \mathrm{d} \mathbf{\omega}\)

其中,\(\rho\) 为 BRDF 项;\(V\) 为 Visibility,表示不被遮挡的程度,往往表现为自遮挡产生的阴影现象。

  • 由于物体是 diffuse 的,因此它的 BRDF 将是一个常数 \(\rho\) (无论从哪个方向观察都得到相同的BRDF值)
  • 对于 lighting 部分 \(L(\mathbf{\omega})\),原本需要对 Environment Map 进行查询,而现在可以换成使用 SH 函数去表示:\(L(\mathbf{\omega}) \approx \sum l_{i} B_{i}(\mathbf{\omega})\)
  • 对于 transfer function 部分 \(T(\omega)=V(\mathbf{\omega}) \max (0, n\cdot \mathbf{\omega})\) , 也可以换成使用 SH 函数去表示:\(T(\omega) \approx \sum T_j B_j(\omega)\)

代入渲染方程整理后得:

\(\begin{align} L(r) &\approx \frac{\rho}{\pi} \int_{\Omega^+} \sum_{i=1}^{n^2}l_{i}B_i(\omega) \sum_{j=1}^{n^2}T_j\mathrm{B}_{j}(\mathbf{\omega}) \mathrm{d} \mathbf{\omega} \\ &= \frac{\rho}{\pi}\sum_{i=1}^{n^2} \sum_{j=1}^{n^2} l_{i}T_j\int_{\Omega^+} B_i(\omega) \mathrm{B}_{j}(\mathbf{\omega})d\omega \\ &= \frac{\rho}{\pi}\sum_{i=1}^{n^2} l_{i}T_i\int_{\Omega^+} B_i(\omega) \mathrm{B}_{i}(\mathbf{\omega})d\omega \\ &= \frac{\rho}{\pi}\sum_{i=1}^{n^2} l_{i}T_i \end{align}\)

这样,我们只需要预计算出 \(l_i\)\(T_i\),就能让运行时的diffuse物体渲染速度大大提升。

预计算过程,对物体模型的每个顶点:

  1. 预计算 lighting 的 SH 系数向量组 \(\vec{l}\),其中一个元素为:\(l_{i}=\int_{\Omega^+} L(\mathbf{\omega}) B_{i}(\mathbf{\omega}) \mathrm{d} \mathbf{\omega}\)

  2. 预计算 light transfer 的 SH 系数向量组 \(\vec{T}\),其中一个元素为:\(T_i=\int_{\Omega^{+}}V(\mathbf{\omega}) \max (0, n\cdot \mathbf{\omega}) \mathrm{B}_{i}(\mathbf{\omega}) \mathrm{d} \mathbf{\omega}\)

此时,我们可以理解成系数 \(l_i\) 代表了环境光照(lighting)的信息,而系数 \(T_i\) 代表了光路传输(light transfer)的信息

运行时渲染过程:

  1. 在 vertex shding 阶段计算顶点的 SH 颜色 \(L(r) \approx \frac{\rho}{\pi} \sum_{i=1}^{n^2} l_{i} T_{i} = \frac{\rho}{\pi}( \vec{l}\cdot \vec{T})\)

  2. 在 pixel/fragment shading 阶段得到插值后的 SH 颜色即为该像素的颜色

Glossy 物体的球谐光照

物体的渲染方程:

\(L(r)=\int_{\Omega^+} L(\mathbf{\omega})\cdot \rho(r,\mathbf{\omega}) \cdot V(\mathbf{\omega}) \cdot \max (0, n\cdot \mathbf{\omega}) \mathrm{d} \mathbf{\omega}\)

对于 glossy 物体,不同视角观察物体表面同一点会有不同的光照。因此 glossy 的 BRDF 将是四维的函数(参数不仅包含\(\mathbf{\omega}\),还包含\(r\)),这次将 BRDF 算入 transfer,则 transfer function 将额外增加一个二维参数 \(r\)

  • 对于环境光(lighting)即 \(L(\mathbf{\omega})\),原本需要对 Environment Map 进行查询,而现在可以换成使用 SH 函数去表示:\(L(\mathbf{\omega}) \approx \sum l_{i} B_{i}(\mathbf{\omega})\)
  • 对于传输函数(transfer function) \(T(r,\omega)=\rho(r,\mathbf{\omega}) V(\mathbf{\omega}) \max (0, n\cdot \mathbf{\omega})\) ,换成使用 SH 函数去表示:\(T(r,\omega) \approx \sum_{i=1}^{n^2} T_{i}(r)B_i(\omega)\approx \sum_{i=1}^{n^2} \sum_{j=1}^{n^2} T_{ij} B_j(r)B_i(\omega)\)

由于额外多了一个二维的参数,transfer 系数将是一个矩阵而非之前 diffuse 情况下的系数向量。

代入渲染方程整理后得:

\(\begin{align} L(r) & \approx \int_{\Omega^+} \sum_{i=1}^{n^2} l_{i} B_{i}\sum_{j=1}^{n^2}\sum_{k=1}^{n^2} T_{jk} B_j(\omega)B_k(r)\mathrm{d} \mathbf{\omega} \\ &= \sum_{i=1}^{n^2}\sum_{j=1}^{n^2}\sum_{k=1}^{n^2} l_iT_{jk }B_k(r) \int_{\Omega^+} B_{i}(\omega)B_j(\omega)\mathrm{d} \mathbf{\omega} \\ &= \sum_{i=1}^{n^2}\sum_{k=1}^{n^2}l_iT_{ik }B_k(r) \end{align}\)

预计算过程,对物体模型的每个顶点:

  1. 预计算 lighting 的 SH 系数向量组 \(\vec{l}\),其中一个元素为:\(l_{i}=\int_{\Omega^+} L(\mathbf{\omega}) \cdot B_{i}(\mathbf{\omega}) \mathrm{d} \mathbf{\omega}\)

  2. 预计算 light transfer 矩阵 \(T\),其中一个元素为:\(T_{ij}=\int_{\Omega^{+}}\mathrm{B}_{j}(\mathbf{\omega}) T_i(\omega) \mathrm{d} \mathbf{\omega}\)

运行时渲染过程:

  1. 在 vertex shding 阶段计算顶点的 SH 颜色 \(L(r) \approx \sum_{i=1}^{n^2}\sum_{j=1}^{n^2} l_{i} T_{ij} B_j(r) = T\vec{l}\cdot \vec{B(r)}\)
  2. 在 pixel/fragment shading 阶段得到插值后的 SH 颜色即为该像素的颜色

原论文实际上还有个 BRDF 卷积步骤,具体还没细看。

球谐光照预计算 Transfer

PRT 的 SH Lighting 预计算是针对单个模型物体的预计算,即输入一个 模型物体 + Environment Map,然后输出模型中每个顶点的 Lighting SH 系数和 Transfer SH 系数或矩阵。这样也意味着这种 GI 方式的 Transfer 只考虑了模型自身的几何遮蔽关系,并不会考虑其它物体。

TODO

球谐的旋转(SH Rotation)

PRT 的一个问题是如果 lighting 部分是预计算的,那就只适用于静态环境光下的静态物体渲染;环境光或者物体只要有变化,PRT 就不得不进行重新预计算;但得益于 SH 的旋转不变性,我们至少可以让 SH Lighting 适用于动态旋转的情形而不必重新预计算。

环境光旋转往往用来实现场景的昼夜轮转效果,物体旋转就比较常见了。顺便注:环境光旋转和物体旋转在 PRT 渲染中是等价的(只是说看相对于哪个东西来看待旋转而已)

现在给定旋转 \(R\) ,然后有SH投影函数 \(P\)(输入一个球面向量,输出第 \(l\) 层band的 SH 系数向量组)。假设我们有 \(2l+1\) 个任意的球面向量,我们需要想办法求出能等价于旋转 \(n_i\) 的矩阵 \(M\) 来旋转 SH投影函数 \(P\)

\(M P\left(n_{i}\right)=P\left(R\left(n_{i}\right)\right), i \in[-l, l]\)

整理得:

\(M\left[P\left(n_{-l}\right), \ldots, P\left(n_{l}\right)\right]=\left[P\left(R\left(n_{-l}\right)\right), \ldots, P\left(R\left(n_{l}\right)\right)\right]\)

\(A = [P_{(n−l)}, ..., P (n_{l})]\) ,如果矩阵 \(A\) 是可逆的,则:

\(M=\left[P\left(R\left(n_{-l}\right)\right), \ldots, P\left(R\left(n_{l}\right)\right)\right] A^{-1}\)

这样,无论 \(R\) 怎么变化,我们都可以即时根据 \(R\) 求出 \(M\) ,然后把最新的 \(M\) 应用到所有预计算好的 \(P(\omega)\) 上(即预计算好的SH系数上),而不必重新预计算 \(P(R(\omega))\)

SH 的快速旋转实现:

对于第 \(l\) 层 band 需要做如下处理:

  1. 选取共 \(2l + 1\) 个 单位向量 \(n_i\),这些向量 \(n_i\) 投影在该层 band 上的基函数就能得到系数 \(P (n_i)\)。共 \(2l + 1\)\(2l+1\) 维向量 \(P (n_i)\) 构成了矩阵 \(A\),并求出 \(A^{−1}\)

如何选取单位向量:要保证投影后构成的 \(A\) 矩阵可逆。

  1. 给定旋转 \(R\) ,对所有 \(n_i\) 依次做旋转 \(R(n_i)\) ,这些旋转后的单位向量同样投影在该层 band 上的基函数就能得到系数 \(P(R(n_i))\)。 共 \(2l + 1\)\(2l + 1\) 维向量 \(P(R(n_i))\) 构成了矩阵 \(S\)

  2. 求出该层 band 上球谐函数的旋转矩阵 \(M = SA^{−1}\)

  3. \(M\) 乘以该层 band 上的 SH 系数向量组就可以得到旋转后的 SH 系数向量组。

\(l=0\) 时只有1个系数,故不需要处理;\(l=1\) 时需要处理3个系数,其中的 \(A、S、M\) 将会是 3X3矩阵;\(l=2\) 时需要处理5个系数,其中的 \(A、S、M\) 将会是 5X5矩阵......

最后,将每一层 band 的结果重新拼接起来即可得到完整的旋转后的 SH 系数结果。

// 伪代码:三阶SH旋转
SHCoeffientsAfterRotation(Array coeffients,Rotation rotation){
    Array res;
    // 处理 l = 0
    res[0] = coffients[0];
    // 处理 l = 1
    // step 1
    n1 = [1,0,0],n2 = [0,1,0],n3 = [0,0,1];
    A = [ // SHProject(n,l,m):向量n投影到SH中l层第m个基函数,输出对应系数
    [SHProject(n1,1,-1),SHProject(n2,1,-1),SHProject(n3,1,-1)],
    [SHProject(n1,1,0),SHProject(n2,1,0),SHProject(n3,1,0)],
    [SHProject(n1,1,1),SHProject(n2,1,1),SHProject(n3,1,1)]
    ];
    A_inv = InvMatrix(A);
    // step 2
    rn1 = rotation*n1,rn2 = rotation*n2,rn3 = rotation*n3;
    S = [
        [SHProject(rn1,1,-1),SHProject(rn2,1,-1),SHProject(rn3,1,-1)],
        [SHProject(rn1,1,0),SHProject(rn2,1,0),SHProject(rn3,1,0)],
        [SHProject(rn1,1,1),SHProject(rn2,1,1),SHProject(rn3,1,1)]
    ];
    // step 3
    M = S*A_inv;
    // step 4
    coeff_l1 = M*[coeffients[1],coeffients[2],coeffients[3]];
    res[1] = coeff_l1[0];
    res[2] = coeff_l1[1];
    res[3] = coeff_l1[2];
    
    // 处理 l = 2    
    // step 1
	k = 1/sqrt(2);
    n1 = [1,0,0],n2 = [0,0,1],n3 = [k,k,0],n4 = [k,0,k],n5 = [0,k,k];
    A = [
[SHProject(n1,2,-2),SHProject(n2,2,-2),SHProject(n3,2,-2),SHProject(n4,2,-2),SHProject(n5,2,-2)],
[SHProject(n1,2,-1),SHProject(n2,2,-1),SHProject(n3,2,-1),SHProject(n4,2,-1),SHProject(n5,2,-1)],
[SHProject(n1,2,0),SHProject(n2,2,0),SHProject(n3,2,0),SHProject(n4,2,0),SHProject(n5,2,0)],
[SHProject(n1,2,1),SHProject(n2,2,1),SHProject(n3,2,1),SHProject(n4,2,1),SHProject(n5,2,1)],
[SHProject(n1,2,2),SHProject(n2,2,2),SHProject(n3,2,2),SHProject(n4,2,2),SHProject(n5,2,2)]
    ];
    A_inv = InvMatrix(A);
    // step 2
    rn1 = rotation*n1,rn2 = rotation*n2,rn3 = rotation*n3,rn4 = rotation*n4,rn5 = rotation*n5;
    S = [
[SHProject(rn1,2,-2),SHProject(rn2,2,-2),SHProject(rn3,2,-2),SHProject(rn4,2,-2),SHProject(rn5,2,-2)],
[SHProject(rn1,2,-1),SHProject(rn2,2,-1),SHProject(rn3,2,-1),SHProject(rn4,2,-1),SHProject(rn5,2,-1)],
[SHProject(rn1,2,0),SHProject(rn2,2,0),SHProject(rn3,2,0),SHProject(rn4,2,0),SHProject(rn5,2,0)],
[SHProject(rn1,2,1),SHProject(rn2,2,1),SHProject(rn3,2,1),SHProject(rn4,2,1),SHProject(rn5,2,1)],
[SHProject(rn1,2,2),SHProject(rn2,2,2),SHProject(rn3,2,2),SHProject(rn4,2,2),SHProject(rn5,2,2)]
    ];
    // step 3
    M = S*A_inv;
    // step 4
    coeff_l2 = M*[coeffients[4],coeffients[5],coeffients[6],coeffients[7],coeffients[8]];
    res[4] = coeff_l2[0];
    res[5] = coeff_l2[1];
    res[6] = coeff_l2[2];
    res[7] = coeff_l2[3];
    res[8] = coeff_l2[4];
    
    return res;
}

这样我们对 lighting 的 SH 系数做 Rotation,就可以实现支持环境光旋转或物体旋转的实时 SH Lighting。

但是注意,不应当对 transfer 的 SH 系数做 Rotation,因为一个物体的 transfer function 只考虑自身的几何遮蔽(自身与自身的遮蔽关系),而非像物体与环境光那样存在相对关系。

SH Rotation 效果图:

光照探针(Light Probe)

光照探针(Light Probe) 是一类场景GI方案。它在场景里设置若干个Probe点(可以理解成向四面八方探测光照的点),为每个 Probe 预计算出环境光信息后,在运行时通过物体周围的 Probe 信息来插值来得到物体此时受到的环境光。

实际上,PRT也算是一种基于Probe的GI方案:PRT物体上的每个顶点都是一个Probe,它所探测的环境光照信息将预计算为 SH 系数,然后 shading point 可以根据三角形的顶点对 SH 系数做重心插值从而得到该点的 SH 系数(即代表了环境光)。只是 PRT 的每个三角形顶点都是一个Probe,这些Probe综合起来完成了单个物体的 GI 效果;而一般的 Light Probe 方案是在场景中每隔一段空间放置一个Probe,这些Probe综合起来提供了场景中所有物体的 GI 效果。

不过 Light Probe 经常会遇到错误的 Bleeding 问题:受到几何上不应该存在光照关系的 Probe 影响,常见于墙壁遮挡的内外侧。

例如,下面的一个 Probe 生成在室内:

结果由于 Light Probe 方法没有考虑遮挡,直接且错误地进行了对相邻的 Probe 插值,从而室外本该明亮的一侧变暗,室内本该黑暗的一侧发生漏光:

因此一个好的 Probe 放置方案可以尽量减少这种 Bleeding 问题。

Unity 四面体镶嵌(Tetrahedral Tessellations)

Unity 允许在场景中放置自定义位置的 Probe ,而且这些Probe将相连成一个个四面体。

为了增加这种GI方案的真实感和避免过多的Probe带来的存储开销,还应当把 Probe 在光照发生明显变化的地方(如明暗交接处)放密集些,而在光照不怎么变化的地方可以稀疏地放置。

预计算过程:

  1. 每个 Probe 从自身往四面八方看到的环境光信息烘焙成3阶SH系数(SH Lighting方法)并记录之

Unity 曾经提供了一种实时GI方案,只不过现在被废弃了:它的 Probe 是运行时采样并实时计算出最新的2阶 SH 系数。

  1. 将这些 Probe 进行三角化(将空间切分成四面体),切分的算法是使用Delaunay Triangulation 完美三角剖分。

运行时渲染过程:

  1. 在渲染物体的 shading point 时,通过四面体插值的方法取周围4个 Probe 的结果的加权平均作为该点受到的 lighting

四面体插值:假设 \(a,b,c,d\) 分别为四个 Probe 的权重,则

\(\left[\begin{array}{l}a \\ b \\ c\end{array}\right]=\left[\begin{array}{lll}\vec{P}_{0}-\vec{P}_{3} & \vec{P}_{1}-\vec{P}_{3} & \overrightarrow{P_{2}}-\vec{P}_{3}\end{array}\right]^{-1}\left[\vec{P}-\overrightarrow{P_{3}}\right]\)

\(d=1-a-b-c\)

渲染效果:

Unity Tetrahedral Tessellations 的弊端:

  • 运行时,查找所处四面体的CPU运算量较高:四面体的分布并不是均匀的,GPU很难做到查找,所以Unity需要通过CPU把每个使用GI的物体所处的四面体找到并将查找结果传递给GPU,这就带来了性能瓶颈。

UE4 间接光缓存(ILC) & 体积光照贴图(VLM)

间接光缓存(Indirect Light Cache):与 Unity Light Probe 方案基本一样,Probe 也是采用烘焙 SH 系数和插值的方案,只是放置 Probe 的方式有所不同。ILC 基本上是在静态物体表面法线向上均匀放置 Probe ,这是基于假设受环境光影响最大的地方都是在靠近静态物体的地方。

当然实际上 ILC 还包含两种 Probe 放置方式:通过手动放置 Lightmass Importance Volume(重要光照范围) 来限制烘焙光照采样范围,或 Lightmass Character Indirect Detail Volume,来增加一段均匀放置 Probe 的区域。

ILC 通过 CPU 寻找物体周围的 Probe:所有的 Probe 的数据将保存到八叉树中以方便物体找到周围的 Probe。

体积光照贴图(Volumetric Light Map,VLM):还是采用烘焙 SH 系数和插值的方案,也还是放置 Probe 的方式有所不同。VLM 使用网格(Grid)采样点保存 Probe 数据,在静态物体表面附近,会对网格进一步细分。

VLM 将通过 GPU 来寻找物体周围的 Probe :所有 Probe 的数据将烘焙至贴图中,不同细分层度的网格将使用 level 不同的贴图,这样可以方便地在 GPU 中进行逐像素的三线性插值。

VLM 与 ILC 比较:

  • VLM 比 ILC 的 Probe 要多更多,从而 VLM 的 GI 效果更佳,但要存储的数据更多
  • VLM 的存储结构决定了可以通过 GPU 算法来寻找 Probe,这比传统的 CPU 算法性能更加客观(即是 VLM Probe 数量要多得多)
  • VLM 更适合将 Probe 应用到体积雾效果中,这是因为 ILC 往往不在几乎没什么物体的空间放置 Probe

UE 4.18 以后使用 VLM 代替了 ILC 方法作为 UE4 默认的 Light Probe 方案。

UE4 ILC & VLM 效果图:

光照贴图(Light Map)

有关 Light Map 的初步介绍可见 实时渲染基础(4)纹理-光照贴图

Light Map 结合动态光照

在 UE4 和 Unity 中,静态物体的全局光照效果都采用了烘焙 Light Map 的方法。

Light Map 方法实际上就是把环境光被分成了静态部分和动态部分,而 Light Map 记录的正是受静态部分环境光后着色结果,因此参与烘焙的物体(位置、形状、材质)和光照都应当是静态的,这样的局限性很大。

为了让静态物体也能受到动态部分环境光的影响,将 Light Map 中对应位置的着色信息解析出来,并与动态光照部分计算出来的着色进行混合。

这里对Unity、UE4所实现的 Light Map 做一个更深入分析:

  • Light Map 一般记录的是受静态部分环境光后的着色结果,着色方法不限于Ray Tracing、辐射度、Shadow、环境光遮蔽(Ambient Occlusion,AO)等算法。

比较有意思的是,对于颜色存储的方案,Unity采用了直观的线性颜色存储方式,而UE4则采用了非线性映射的颜色存储方式:暗色区域的颜色可以有更多细节(更高精度)而明色区域的颜色则更少细节(更低精度)。这是因为在实际画面中,我们更容易注意到暗色区域的细节而不容易注意到明亮区域的细节。

  • Light Map 也可以直接混合写入正常纹理,这样物体着色只需要采样一个纹理,但会极度影响纹理的重用(有别的物体使用了相同的纹理,但是光照结果却是截然不同的;有些着色算法依赖原生纹理而非混合纹理的数据)

  • Light Map 对于带 Normal Map 的物体,可能无法得到正确的结果(Normal Map没有起作用)。这是因为光照烘焙系统只使用了物体的顶点数据(而没有使用 Normal Map)。而 Unity 和 UE4 的 Light Map 在启动了高质量选项后会多额外一张纹理(称为 Directional Light Map )来存储入射光的主要方向(世界坐标空间):在对某个 shading point 着色时,通过解析出来的入射光主要方向与 shading point 的 normal(Normal Map 解压出来的,且需变换成世界坐标空间下的)进行点乘,便可以得到对应 Light Map 颜色的贡献权重。

Lighting Scenarios 实现昼夜天气系统

Lighting Scenarios :一个物体使用多套 Light Map,常用于模拟天气系统以及一天中不同时段的光照。实际上,就是烘焙几个关键时间点的 Light Map, 根据时间进行线性插值。

UE4 的 Lighting Scenarios 效果:

总结

GI 总体上分为两派:一派是通过一套统一标准的渲染流程把任何物体的 GI 计算出来,往往计算量极大但效果更加物理更加真实,典型的例子便是离线光线追踪;另一派是把 GI 的效果看成一个个部分来组成,这样我们可以选择其中一些GI效果组合使用,适应不同的性能(往往计算量相对低些,尤其是实时渲染)的同时也能带来能接受的 GI 效果(虽然往往不是严谨的物理正确),典型的例子就是这篇博客所介绍的方法和另一些完全动态的 GI 方法(如 SSAO)。

个人归纳:

  • IBL 是半动态的 GI 方法(环境光可动可静、动态物体、材质仅粗糙度参数可变),常用于 Specular 物体的渲染
  • PRT 是半动态的 GI 方法(环境光可动可静、即使静态环境光还可支持旋转、物体不可形变、静态材质), 常用于 Diffuse/Glossy 物体的渲染
  • Light Probe 是半动态的 GI 方法(静态环境光、动态物体), 实时给场景中任何物体提供合适的静态环境光信息
  • Light Map 是全静态的 GI 方法(静态环境光、静态物体),预先给场景中静态物体提供受环境光静态部分影响后的着色结果

参考

    原文作者:KillerAery
    原文地址: https://www.cnblogs.com/KillerAery/p/15335369.html
    本文转自网络文章,转载此文章仅为分享知识,如有侵权,请联系管理员进行删除。