前言

好久没写博客了，今天来共享一下最近研读的一篇论文吧，文章中[***.pdf]为参考文献可自行谷歌学术下载。因为本人还是个图形学菜鸟，所以有什么问题希望大家多多指正。

Stochastic Light Culling for VPLs on GGX Microsurfaces

词汇解析

VPL

VPL - virtual point light 虚拟点光源

以虚拟的点光源为间接光提供光照。可以在场景中分布多个虚拟点光源，对每个点光源的亮度都设置最大限度，这样就可以模拟在空间中传播的间接光。核心观点是，一般的光传输问题可以近似为计算许多虚拟点光源的直接照明的简单问题。

路径跟踪通过从相机开始跟踪光线来构建路径。在与场景的每个交点上，通过在交点表面上方的半球上随机采样一个出射方向(通常在每个交点上也估计直接光照)来继续路径。基于路径跟踪，双向方法跟踪摄像机子路径和光子路径，然后(确定地)将它们连接起来形成完整的路径。

Instant Radiosity (IR)是一种以特定方式构造这两组子路径的双向方法，如上图2所示。首先，生成并存储一组任意长度的光子路径[问题6]。为子路径上的每个顶点(不是路径上的每个点)，会完整地储存局部环境：位置、法线、入射方向、BRDF和当前“通量”(即发光光源的辐射通量乘以顶点的吞吐量除以该点的路径的概率密度)。其目的是，为每个顶点存储足够的数据以计算从该顶点向任何方向散射的出射照度。如果按照这种思维，我们可以放弃原始路径的概念，而将顶点建模为一种不同寻常的点光源。由于这些不是场景中的物理光源，我们称之为虚拟点光源(VPLs)。

为了完成IR算法，在第二阶段为每个像素构造相机子路径。因为光的子路径是任意长的，所以只考虑长度为1的摄像机路径就足够了。然后，像所有的双向算法一样，IR将这些摄像机子路径的顶点连接到光路径的顶点，形成完整的路径。但是，由于VPL的生成是预处理过程，所以这个连接步骤与VPL对第一个摄像机命中点的直接照明相当。

一个基本的VPL照明方法总结如下：

第一阶段：生成VPL

在场景中随机选择一个主光源，随机采样一个位置 x \boldsymbol{\mathbf{x}} x和方向 ω \boldsymbol{\omega} ω(如果不处理直接照明，在此位置创建一个VPL)。

追踪射线 x + t ω \boldsymbol{\mathbf{x}} + t\boldsymbol{\omega} x+tω。如果它与一个表面相交，那么在这个相交位置创建一个VPL

使用俄罗斯轮盘随机决定是否终止路径。如果继续，则采样输出方向，基于BRDF和方向更新路径吞吐量，并继续跟踪。

第二阶段：使用VPL渲染

为一个表面点着色，只需遍历所有VPL，测试其照明是否被阻挡，如果没有，计算各自的贡献。

[Scalable Realistic Rendering with Many-light Methods.pdf]

GGX

GGX是一种更精确的微表面分布函数，是一种描述微表面法线方向分布的函数，可对于表面粗糙的半透明物体和表面粗糙的不透明物体均可适用。

https://www.jianshu.com/p/be4f025aeb3c

Caustics

Caustics(光焦散)是由曲面或物体反射或折射的光的聚合，或该光的聚合在另一表面上的投影。光焦散是每条光与之相切的曲线或表面，将光的聚合线的边界定义为聚光曲线。因此，在相邻图像中，焦散可以是光斑或它们的明亮边缘。这些形状通常具有尖锐的奇点。光焦散的常见情形就是水波纹。

https://en.wikipedia.org/wiki/Caustic_(optics)

Helmholtz Reciprocity

Helmholtz互异性:入射角和出射角互换，函数值保持不变。

https://blog.csdn.net/weixin_41036461/article/details/79846176

Lightcuts

Lightcuts(光剪切)是第一种实用的、可扩展的多光源方法。在需要计算光照的每个接收点，Lightcuts会根据分析的每个群集的误差范围和感知指标生成自定义的切割。首先，基于空间和方向相似性，将VPL先组织成二叉树。要为接收点选择一个片段，我们从一个简单的粗略聚类开始，例如将所有灯光放在一个聚类中。这对应于仅由光树的根节点组成的剪切。然后，我们迭代地选择当前剪切中误差范围最大的聚类，并对其进行精炼以替换其在光树中的孩子。重复此过程，直到剪切中所有聚类的误差范围低于基于感知的阈值（通常设置为总数的2％）为止。剪切的大小通常仅取决于VPL的数量，与完全评估相比，随着VPL数量的增加，剪切速度会大大提高。使用分析误差界限还可以确保始终找到并评估最重要的VPL，从而使估算更加可靠。 Davidovic描述了渐进的、对GPU友好的Lightcuts变体。

然而，使用误差界限也有一些缺点。误差边界忽略了遮挡，并且通常过于保守，导致比最优剪切更大。此外，为新的材质模型开发良好的误差界限是具有挑战性和困难的。Lightcuts也没有利用附近接收点的剪切通常非常相似这一事实，而是为每个点从头生成一个剪切。

[Scalable Realistic Rendering with Many-light Methods.pdf]

[Lightcuts A Scalable Approach to Illumination.pdf]

阴影遮罩函数

masking-shadowing function即阴影遮罩函数(也称为可见性函数)，是描述光线从入射到出射的过程中，有多少比例被微表面自身的凹凸不平遮挡主了。这部分和法线分布函数很像，依据的是表面的粗糙程度。为了准确描述，首先需要这样的一个函数：给一个方向，获得微表面在这个方向上被自遮挡的面积的比例。

[Geometrical shadowing of a random rough surface.pdf]

[Understanding the masking-shadowing function in microfacet-based BRDFs.pdf]

主要目标

由于现存的光剔除技术需要限制光源的影响范围，所以会造成场景明显偏暗。

本文的方法为利用GGX微表面生成光滑的单反弹光焦散的实时绘制方法。该方法是基于VPL的随机光剔除，这是一种无偏差的剔除方法，可以随机确定每个VPL的光影响范围，并且是对每个VPL的包围椭球进行分块剔除。包围椭球的计算使用的是一种推导出的简单分析方法，可以为GGX微表面上的VPL生成更紧密的边界椭球。对于成千上万的VPL时，该方法结合交错采样比基于保守光栅化的聚类着色的速度更快(聚类着色法是最新的支持剔除边界椭球的技术)。

本文主要有如下贡献：

将随机光剔除方法推广到光滑VPL的方向相关影响范围。

从GGX微平面BRDF反射的VPL中解析地推导出光范围的包围椭球体。

本文提出了一种有效的包围椭球体的分块剔除算法。

本文的绘制管线如下图2所示。上诉贡献在绘制管线中的橙色部分实现了。

算法细节

3 光滑虚拟点光源的随机光剔除

3.1 随机光剔除

随机光剔除算法是通过俄罗斯轮盘确定不重要[问题5]光源的影响范围来随机地删除不重要的光源。对于每一个光源，俄罗斯轮盘通过其在着色点的概率来决定是否剔除这个光源，这个概率与辐射率成正比。对于每一个点光源，其在某个着色点的辐射率如下：

L ( ω o , l ) = I ( ω o ) l 2 L\left(\boldsymbol{\omega}_{o}, l\right)=\frac{I\left(\boldsymbol{\omega}_{o}\right)}{l^{2}} L(ωo,l)=l2I(ωo)

其中， I ( ω o ) I(\boldsymbol{\omega}_{o}) I(ωo)表示辐射强度， ω O ∈ S 2 \boldsymbol{\omega}_{O}\in\mathcal{S}^{2} ωO∈S2表示从点光源到着色点的方向， l ∈ [ 0 , ∞ ) l\in[0, \infty) l∈[0,∞)表示从点光源到着色点的距离。

在这篇文章中，使用如下的概率函数：

p ( ω o , l ) = min ⁡ ( I ( ω o ) δ l 2 , 1 ) p\left(\boldsymbol{\omega}_{o}, l\right)=\min \left(\frac{I\left(\boldsymbol{\omega}_{o}\right)}{\delta l^{2}}, 1\right) p(ωo,l)=min(δl2I(ωo),1)

其中 δ ∈ ( 0 , ∞ ) \delta \in(0, \infty) δ∈(0,∞)是一个用户自定义的参数(这篇文章用的 δ = 0.001 \delta = 0.001 δ=0.001)。

如果某个点没有被剔除的话(概率超过某个阈值)，就会将它的辐射率除以概率值[问题1]，如下所示：

L ( ω o , l ) ≈ { L ( ω 0 , l ) p ( ω o , l ) ( p ( ω o , l ) > ξ ) 0 ( otherwise ) L\left(\boldsymbol{\omega}_{o}, l\right) \approx\left\{\begin{array}{ll} \frac{L\left(\boldsymbol{\omega}_{0}, l\right)}{p\left(\boldsymbol{\omega}_{o}, l\right)} & \left(p\left(\boldsymbol{\omega}_{o}, l\right)>\xi\right) \\ 0 & (\text { otherwise }) \end{array}\right. L(ωo,l)≈{p(ωo,l)L(ω0,l)0(p(ωo,l)>ξ)( otherwise )

其中， ξ ∈ [ 0 , 1 ) \xi \in[0,1) ξ∈[0,1)是一个统一的随机数。

随机光剔除算法中，对每一个光源都使用单个随机数 ξ \xi ξ；并且对于一个给定的光源，所有的着色点都会使用相同的 ξ \xi ξ。这样做的好处就是可以限制每个光源的影响范围，并且能够以无偏差的方式时利用现有的剔除方法。

3.2 虚拟点光源的BRDF相关影响范围

Tokuyoshi为了使用包围球和现有的块剔除方法，但却忽略了概率函数 p ( ω o , l ) p\left(\boldsymbol{\omega}_{o}, l\right) p(ωo,l)的方向性。这篇文章考虑了光滑虚拟点光源的方向性，如下图4所示：

所以，某个VPL的辐射强度计算如下：

I ( ω o ) = Φ f ( ω i , ω o ) max ⁡ ( ω o ⋅ n , 0 ) I\left(\boldsymbol{\omega}_{o}\right)=\Phi f\left(\boldsymbol{\omega}_{i}, \boldsymbol{\omega}_{o}\right) \max \left(\boldsymbol{\omega}_{o} \cdot \mathbf{n}, 0\right) I(ωo)=Φf(ωi,ωo)max(ωo⋅n,0)

其中 Φ \Phi Φ表示光子到达VPL处的辐射通量， f ( ω i , ω o ) f\left(\boldsymbol{\omega}_{i}, \boldsymbol{\omega}_{o}\right) f(ωi,ωo)是BRDF， ω i ∈ S 2 \boldsymbol{\omega}_{i} \in \mathcal{S}^{2} ωi∈S2是光子的入射方向， n ∈ S 2 \mathbf{n} \in \mathcal{S}^{2} n∈S2是VPL处的几何法线。因此，上图4所示的VPL的影响范围如下：

l max ⁡ ( ω o ) = p − 1 ( ξ ) = Φ f ( ω i , ω o ) max ⁡ ( ω o ⋅ n , 0 ) δ ξ l_{\max }\left(\boldsymbol{\omega}_{o}\right)=p^{-1}(\xi)=\sqrt{\frac{\Phi f\left(\boldsymbol{\omega}_{i}, \boldsymbol{\omega}_{o}\right) \max \left(\boldsymbol{\omega}_{o} \cdot \mathbf{n}, 0\right)}{\delta \xi}} lmax(ωo)=p−1(ξ)=δξΦf(ωi,ωo)max(ωo⋅n,0)

其中 l max ⁡ ( ω o ) l_{\max }\left(\boldsymbol{\omega}_{o}\right) lmax(ωo)是BRDF( f ( ω i , ω o ) f\left(\boldsymbol{\omega}_{i}, \boldsymbol{\omega}_{o}\right) f(ωi,ωo))所反射的辐射率的等值面。为了绑定这个等值面，我们使用了与Dachsbacher溅射算法中使用过类似的椭球面。但与之不同的是，我们为GGX微平面的BRDF引入了包围椭球。

公式推导如下:

(字丑，见谅)

4 GGX反射的包围椭球

4.1 GGX微平面的BRDF

建模微平面的BRDF可以用来表示粗糙平面光的反射，BRDF计算如下：

f ( ω i , ω o ) = F ( ω i ⋅ ω h ) G 2 ( ω i , ω o ) D ( ω h ⋅ n ) 4 ∣ ω i ⋅ n ∥ ω o ⋅ n ∣ f\left(\boldsymbol{\omega}_{i}, \boldsymbol{\omega}_{o}\right)=\frac{F\left(\boldsymbol{\omega}_{i} \cdot \boldsymbol{\omega}_{h}\right) G_{2}\left(\boldsymbol{\omega}_{i}, \boldsymbol{\omega}_{o}\right) D\left(\boldsymbol{\omega}_{h} \cdot \mathbf{n}\right)}{4\left|\boldsymbol{\omega}_{i} \cdot \mathbf{n} \| \boldsymbol{\omega}_{o} \cdot \mathbf{n}\right|} f(ωi,ωo)=4∣ωi⋅n∥ωo⋅n∣F(ωi⋅ωh)G2(ωi,ωo)D(ωh⋅n)

其中， ω h = ω i + ω o ∥ ω i + ω o ∥ \boldsymbol{\omega}_{h}=\frac{\boldsymbol{\omega}_{i}+\boldsymbol{\omega}_{o}}{\left\|\boldsymbol{\omega}_{i}+\boldsymbol{\omega}_{o}\right\|} ωh=∥ωi+ωo∥ωi+ωo表示半程向量， F ( ω i ⋅ ω h ) ∈ [ 0 , 1 ] F\left(\boldsymbol{\omega}_{i} \cdot \boldsymbol{\omega}_{h}\right) \in[0,1] F(ωi⋅ωh)∈[0,1]是菲涅尔因子。

D ( ω h ⋅ n ) D\left(\boldsymbol{\omega}_{h} \cdot \mathbf{n}\right) D(ωh⋅n)是NDF，表示微平面法线的分布。在GGX中，NDF是一个钟形的函数，定义如下：

D ( cos ⁡ θ m ) = α 2 χ + ( cos ⁡ θ m ) π ( α 2 cos ⁡ 2 θ m + sin ⁡ 2 θ m ) 2 D\left(\cos \theta_{m}\right)=\frac{\alpha^{2} \chi^{+}\left(\cos \theta_{m}\right)}{\pi\left(\alpha^{2} \cos ^{2} \theta_{m}+\sin ^{2} \theta_{m}\right)^{2}} D(cosθm)=π(α2cos2θm+sin2θm)2α2χ+(cosθm)

其中， α \alpha α是粗糙度参数， χ + ( cos ⁡ θ m ) \chi^{+}(\cos \theta_{m}) χ+(cosθm)是海维赛德函数：如果 cos ⁡ θ m > 0 \cos \theta_{m} > 0 cosθm>0，则表示1；如果 cos ⁡ θ m ≤ 0 \cos \theta_{m} \leq 0 cosθm≤0，则表示0。本文假定 α ∈ ( 0 , 1 ] \alpha \in(0,1] α∈(0,1]，这在计算机图形学的数据压缩中很常用，并且很容易被艺术家控制。如果 α ∈ ( 0 , 1 ] \alpha \in(0,1] α∈(0,1]，那么GGX的NDF在 cos ⁡ θ m ∈ [ 0 , π / 2 ] \cos \theta_m \in[0,\pi/2] cosθm∈[0,π/2]时单调递减。

G 2 ( ω i , ω o ) ∈ [ 0 , 1 ] G_{2}\left(\boldsymbol{\omega}_{i}, \boldsymbol{\omega}_{o}\right) \in[0,1] G2(ωi,ωo)∈[0,1]是阴影遮蔽函数。在本文中，我们使用Smith微平面模型，因为这个遮蔽函数是可分离的[问题2]：

G 1 ( ω i , ω h ) = χ + ( ω i ⋅ ω h ) G 1 dist ( ω i ) G_{1}\left(\boldsymbol{\omega}_{i}, \boldsymbol{\omega}_{h}\right)=\chi^{+}\left(\boldsymbol{\omega}_{i} \cdot \boldsymbol{\omega}_{h}\right) G_{1}^{\text {dist }}\left(\boldsymbol{\omega}_{i}\right) G1(ωi,ωh)=χ+(ωi⋅ωh)G1dist (ωi)

其中, χ + ( ω i ⋅ ω h ) \chi^{+}(\boldsymbol{\omega}_{i} \cdot \boldsymbol{\omega}_{h}) χ+(ωi⋅ωh)表示的是正面微平面的二元可见性，并且 G 1 d i s t ( ω i ) G_{1}^{dist}(\boldsymbol{\omega}_i) G1dist(ωi)与 ω h \boldsymbol{\omega}_h ωh无关。根据可见法线分布的约束， G 1 d i s t ( ω i ) G_{1}^{dist}(\boldsymbol{\omega}_i) G1dist(ωi)定义如下：

G 1 dist ( ω i ) = ∣ ω i ⋅ n ∣ ∫ S 2 D ( ω ⋅ n ) max ⁡ ( ω i ⋅ ω , 0 ) d ω = 2 ∣ ω i ⋅ n ∣ ∣ ω i ⋅ n ∣ + ( 1 − α 2 ) ( ω i ⋅ n ) 2 + α 2 \begin{aligned} G_{1}^{\text {dist }}\left(\boldsymbol{\omega}_{i}\right) &=\frac{\left|\boldsymbol{\omega}_{i} \cdot \mathbf{n}\right|}{\int_{\mathcal{S}^{2}} D(\boldsymbol{\omega} \cdot \mathbf{n}) \max \left(\boldsymbol{\omega}_{i} \cdot \boldsymbol{\omega}, 0\right) \mathrm{d} \boldsymbol{\omega}} \\ &=\frac{2\left|\boldsymbol{\omega}_{i} \cdot \mathbf{n}\right|}{\left|\boldsymbol{\omega}_{i} \cdot \mathbf{n}\right|+\sqrt{\left(1-\alpha^{2}\right)\left(\boldsymbol{\omega}_{i} \cdot \mathbf{n}\right)^{2}+\alpha^{2}}} \end{aligned} G1dist (ωi)=∫S2D(ω⋅n)max(ωi⋅ω,0)dω∣ωi⋅n∣=∣ωi⋅n∣+(1−α2)(ωi⋅n)2+α2 2∣ωi⋅n∣

Smith阴影遮蔽函数有几种形式(比如高度相关的形式)，但任何形式都满足 G 2 ( ω i , ω o ) ≤ G 1 dist ( ω i ) G_{2}\left(\boldsymbol{\omega}_{i}, \boldsymbol{\omega}_{o}\right) \leq G_{1}^{\text {dist }}\left(\boldsymbol{\omega}_{i}\right) G2(ωi,ωo)≤G1dist (ωi)

4.2 包围椭球

因为微平面的BRDF是非常复杂的，所以我们考虑使用一个简化的函数来表示反射波瓣(简化函数大于等于反射波瓣)，定义如下：

f ( ω i , ω o ) max ⁡ ( ω o ⋅ n , 0 ) ≤ F max ⁡ ( ω i ) G 1 d i s t ( ω i ) D ( ω h ⋅ n ) 4 ∣ ω i ⋅ n ∣ f\left(\boldsymbol{\omega}_{i}, \boldsymbol{\omega}_{o}\right) \max \left(\boldsymbol{\omega}_{o} \cdot \mathbf{n}, 0\right) \leq \frac{F_{\max }\left(\boldsymbol{\omega}_{i}\right) G_{1}^{\mathrm{dist}}\left(\boldsymbol{\omega}_{i}\right) D\left(\boldsymbol{\omega}_{h} \cdot \mathbf{n}\right)}{4\left|\boldsymbol{\omega}_{i} \cdot \mathbf{n}\right|} f(ωi,ωo)max(ωo⋅n,0)≤4∣ωi⋅n∣Fmax(ωi)G1dist(ωi)D(ωh⋅n)

其中， F max ⁡ ( ω i ) F_{\max}(\boldsymbol{\omega}_i) Fmax(ωi)是对于 ω i \boldsymbol{\omega}_i ωi来说最大的菲涅尔因子。在不等式右边，与 ω o \boldsymbol{\omega}_o ωo有关的就只有NDF的 D ( ω h ⋅ n ) D(\boldsymbol{\omega}_h \cdot \boldsymbol{n}) D(ωh⋅n)。当出射方向 ω o \boldsymbol{\omega}_o ωo位于以理想镜面反射方向为中心的球形圆上时， ω h \boldsymbol{\omega}_h ωh是一个位于半短轴为 θ 2 \frac{\theta}{2} 2θ的球形椭圆上(如下图5.a和5.b)或半长轴为 θ 2 \frac{\theta}{2} 2θ的球形双曲线上(如下图5.c)，其中理想镜面反射方向为 ω u = 2 ( ω i ⋅ n ) n − ω i \boldsymbol{\omega}_{u}=2\left(\boldsymbol{\omega}_{i} \cdot \mathbf{n}\right) \mathbf{n}-\boldsymbol{\omega}_{i} ωu=2(ωi⋅n)n−ωi，球形圆的半径为 θ = arccos ⁡ ( ω o ⋅ ω u ) \theta=\arccos \left(\boldsymbol{\omega}_{o} \cdot \boldsymbol{\omega}_{u}\right) θ=arccos(ωo⋅ωu)， θ \theta θ表示出射方向 ω o \boldsymbol{\omega}_o ωo与理想镜面反射方向 ω u \boldsymbol{\omega}_u ωu之间的夹角。

因此，我们可以获得一个在 ω h \boldsymbol{\omega}_h ωh和 n \boldsymbol{n} n之间的角度下界[问题3]：

arccos ⁡ ( ω h ⋅ n ) ≥ θ 2 \arccos \left(\boldsymbol{\omega}_{h} \cdot \mathbf{n}\right) \geq \frac{\theta}{2} arccos(ωh⋅n)≥2θ

当出射方向 ω o \boldsymbol{\omega}_o ωo位于经过入射方向 ω i \boldsymbol{\omega}_i ωi和法线 n \boldsymbol{n} n的大圆上时，不等式左右相等。对于 α ∈ ( 0 , 1 ] \alpha \in(0,1] α∈(0,1]，因为GGX的NDF是单调递减的，根据上式就有：

D ( ω h ⋅ n ) ≤ D ( cos ⁡ θ 2 ) D\left(\boldsymbol{\omega}_{h} \cdot \mathbf{n}\right) \leq D\left(\cos \frac{\theta}{2}\right) D(ωh⋅n)≤D(cos2θ)

因此，基于GGX的光滑反射的等值面就被如下的面 s ( ω o ) s(\boldsymbol{\omega}_o) s(ωo)所包围：

l max ⁡ ( ω o ) ≤ s ( ω o ) = Φ F max ⁡ ( ω i ) G 1 dist ( ω i ) D ( cos ⁡ θ 2 ) 4 δ ξ ∣ ω i ⋅ n ∣ = r π D ( cos ⁡ θ 2 ) \begin{aligned} l_{\max }\left(\boldsymbol{\omega}_{o}\right) \leq s\left(\boldsymbol{\omega}_{o}\right) &=\sqrt{\frac{\Phi F_{\max }\left(\boldsymbol{\omega}_{i}\right) G_{1}^{\text {dist }}\left(\boldsymbol{\omega}_{i}\right) D\left(\cos \frac{\theta}{2}\right)}{4 \delta \xi\left|\boldsymbol{\omega}_{i} \cdot \mathbf{n}\right|}} \\ &=r \sqrt{\pi D\left(\cos \frac{\theta}{2}\right)} \end{aligned} lmax(ωo)≤s(ωo)=4δξ∣ωi⋅n∣ΦFmax(ωi)G1dist (ωi)D(cos2θ) =rπD(cos2θ)

其中， r = Φ F max ⁡ ( ω i ) G 1 dist ( ω i ) 4 π δ ξ ∣ ω i ⋅ n ∣ r=\sqrt{\frac{\Phi F_{\max }\left(\boldsymbol{\omega}_{i}\right) G_{1}^{\text {dist }}\left(\boldsymbol{\omega}_{i}\right)}{4 \pi \delta \xi\left|\boldsymbol{\omega}_{i} \cdot \mathbf{n}\right|}} r=4πδξ∣ωi⋅n∣ΦFmax(ωi)G1dist (ωi) 。意外的是，这个平面 s ( ω o ) s(\boldsymbol{\omega}_o) s(ωo)是一个椭球体，如下图6所示(推导见附录A)。

这个椭球体的半轴是( α \alpha α表示粗糙度参数)：

( r u , r v , r w ) = ( 1 + α 2 2 α r , r , r ) \left(r_{u}, r_{v}, r_{w}\right)=\left(\frac{1+\alpha^{2}}{2 \alpha} r, r, r\right) (ru,rv,rw)=(2α1+α2r,r,r)

这个椭球体的旋转矩阵是：

R = [ ω u ω w × ω u ω w ] \mathbf{R}=\left[\begin{array}{lll} \boldsymbol{\omega}_{u} & \boldsymbol{\omega}_{w} \times \boldsymbol{\omega}_{u} & \boldsymbol{\omega}_{w} \end{array}\right] R=[ωuωw×ωuωw]

其中， ω w ∈ S 2 \boldsymbol{\omega}_w \in\mathcal{S}^{2} ωw∈S2是一个单位向量，且与 ω u \boldsymbol{\omega}_u ωu和 n \boldsymbol{n} n正交(也就是， ω w = ω u × n ∥ ω u × n ∥ \boldsymbol{\omega}_{w}=\frac{\boldsymbol{\omega}_{u} \times \boldsymbol{n}}{\|\boldsymbol{\omega}_{u} \times \boldsymbol{n}\|} ωw=∥ωu×n∥ωu×n，当 ω u ≠ n \boldsymbol{\omega}_u \neq \boldsymbol{n} ωu=n)。

这个椭球体的中心是：

c = q + 1 − α 2 2 α r ω u \mathbf{c}=\mathbf{q}+\frac{1-\alpha^{2}}{2 \alpha} r \omega_{u} c=q+2α1−α2rωu

其中， q ∈ R 3 \mathbf{q} \in \mathbb{R}^{3} q∈R3是VPL的位置。

5 使用包围椭球的分块剔除

尽管基于光栅化的剔除方法支持包围椭球，但本文扩展基于计算的分块剔除方法来使用包围椭球。这种基于计算的剔除方法比基于光栅化的剔除方法更高效、代价更低，因为对于成千上万的光源来说，随机光剔除是与交错采样结合使用的。这种扩展在数学上是微不足道的，但是它的计算代价却高于包围球。因此，本文引入了一个优化的方法。

5.1 椭球-截椎体的相交测试

在视图空间中，对于每一个分块，分块剔除都会执行包围体和截椎体的粗略的相交测试(如下图7.a)。

这个视图空间的椭球体使用一个旋转矩阵 R ^ = V R \hat{\mathbf{R}}=\mathbf{V R} R^=VR和一个中心位置 c ^ = V c + o \hat{\mathbf{c}}=\mathbf{V} \mathbf{c}+\boldsymbol{o} c^=Vc+o来表示，其中，3x3的矩阵 V \mathbf{V} V表示从世界空间到视图空间的旋转， o ∈ R 3 \boldsymbol{o} \in \mathbb{R}^{3} o∈R3表示从世界空间到视图空间的平移。该椭球体与截锥体的相交测试可等价地表示为球面与截锥体的相交测试，其方法是利用以下变换矩阵对空间进行拉伸：

S = R ^ [ 1 r u 0 0 0 1 r v 0 0 0 1 r w ] R ^ T \mathbf{S}=\mathbf{\hat{R}}\left[\begin{array}{ccc} \frac{1}{r_{u}} & 0 & 0 \\ 0 & \frac{1}{r_{v}} & 0 \\ 0 & 0 & \frac{1}{r_{w}} \end{array}\right] \mathbf{\hat{R}}^{T} S=R^⎣⎡ru1000rv1000rw1⎦⎤R^T

使用这个变换矩阵，每一个包围椭球体都会被拉伸到单位球(如图7.b)。因此，通过变换每一个光源的测试空间，我们能够对包围椭球做一些修改就可以重用现有的基于包围球的剔除方法。本文采用了Modified HalfZ剔除法，这个方法使用两个深度聚类和一个截椎体的轴对齐包围盒(AABB)来对包围球-截椎体进行粗糙的相交测试。

5.1.1 Modified HalfZ剔除法

图1.4显示了一个分块包围体的2D表示。如图1.4(a)所示，背景物体前面的前景物体会导致分块包围体的深度范围变大。光源可以位于前景和背景之间的空白空间，这样可以通过相交测试，但实际上不影响任何分块中的任何像素。也就是说，深度不连续会导致假阳性交点的增加。

图1.4(b)显示了称为Half Z方法，它是可以更好地处理深度不连续的策略。它只是在中点将深度范围划分为两个，并根据两个深度范围进行筛选：一个是从最小Z到半Z，另一个是从半Z到最大Z。并且逐块为每个深度范围维护一个单独的列表，该方法只需要两个额外的平面测试。

图1.4©显示了第二种策略，称为Modified HalfZ方法。它执行附加的原子操作以在最小Z和半Z之间找到第二个最大值(最大值Z2)，在最小Z和半Z之间找到第二个最小值(最小Z2)。相比Half Z方法，这样可以形成更紧密的包围体，但是，由于需要额外的原子操作，计算附加的最小值和最大值比简单地计算半Z要昂贵得多。

5.2 通过旋转测试空间来加速

对于代码优化，大多数现有的剔除方法实现都假定了截椎体的深度平面是与z轴垂直的。然而，我们的拉伸变换却打破了这个假设。这就增加了相交测试的AABB计算的代码复杂度。另外，由于拉伸后的截椎体与AABB的不匹配(如图7.b)，就会捕获很多的误报。为了缓和这些问题，本文使用如下的旋转矩阵对测试空间进行额外的旋转：

B = { [ b x ∥ b x ∥ b z × b x ∥ b z × b x ∥ b z ∥ b z ∥ ] T ( ∥ b x ∥ > ∥ b y ∥ ) [ b y × b z ∥ b y × b z ∥ b y ∥ b y ∥ . b z ∥ b z ∥ ] T (otherwise) \mathbf{B}=\left\{\begin{array}{ll} {\left[\frac{\mathbf{b}_{x}}{\left\|\mathbf{b}_{x}\right\|} \quad \frac{\mathbf{b}_{z} \times \mathbf{b}_{x}}{\left\|\mathbf{b}_{z} \times \mathbf{b}_{x}\right\|}\right.} & \left.\frac{\mathbf{b}_{z}}{\left\|\mathbf{b}_{z}\right\|}\right]^{\mathrm{T}} \quad\left(\left\|\mathbf{b}_{x}\right\|>\left\|\mathbf{b}_{y}\right\|\right) \\ {\left[\frac{\mathbf{b}_{y} \times \mathbf{b}_{z}}{\left\|\mathbf{b}_{y} \times \mathbf{b}_{z}\right\|} \quad \frac{\mathbf{b}_{y}}{\left\|\mathbf{b}_{y}\right\|}\right.} & .\left.\frac{\mathbf{b}_{z}}{\left\|\mathbf{b}_{z}\right\|}\right]^{\mathrm{T}} \quad \text { (otherwise) } \end{array}\right. B=⎩⎪⎨⎪⎧[∥bx∥bx∥bz×bx∥bz×bx[∥by×bz∥by×bz∥by∥by∥bz∥bz]T(∥bx∥>∥by∥).∥bz∥bz]T (otherwise)

其中， b x = S [ 1 0 0 ] T \mathbf{b}_{x}=\mathbf{S}\left[\begin{array}{lll}1 & 0 & 0\end{array}\right]^{\mathrm{T}} bx=S[100]T， b y = S [ 0 1 0 ] T \mathbf{b}_{y}=\mathbf{S}\left[\begin{array}{lll}0 & 1 & 0\end{array}\right]^{\mathrm{T}} by=S[010]T，并且 b z = b x × b y \mathbf{b}_{z}=\mathbf{b}_{x} \times \mathbf{b}_{y} bz=bx×by。旋转矩阵 B \mathbf{B} B用于使变换后的截椎体的深度平面与z轴保持垂直(如图7.c)。这就能让我们简化AABB的计算代码，并且计算一个对剪切后的截锥体包裹更紧密的AABB。在本文中，每一个光源的矩阵 B S \mathbf{BS} BS和位置 B S c ^ \mathbf{BS\hat{c}} BSc^被计算后都会存储到内存当中。然后，在光剔除阶段的椭球体-截锥体相交测试时，都会加载这些 B S \mathbf{BS} BS和 B S c ^ \mathbf{BS\hat{c}} BSc^。

实现细节

交错采样

除了随机光剔除，本文使用了对GPU友好的交错采样技术来减少每个像素中VPL的数量(如下图8所示)。首先，该技术将规则采样模式中的像素去交错到屏幕的子区域中。然后，对于每一个子区域，使用不同的VPL的子集来进行着色。比如，对65536个VPL使用8x8的交错采样，每一个子区域的VPL数量就会减少到1024个。在本文中，我们的随机光剔除在每个子区域中都会执行。所以，这种交错采样技术不仅减少了着色的时间，还减少了剔除的时间。我们基于分块延迟着色的方法在单通道下实现了这个算法交错样本模式的非交错延迟着色。尽管方差在渲染的图像中可见为噪声，但是在后期处理中使用交叉双边滤波器可以消除此噪声。

漫反射VPL

尽管包围椭球也可以用于漫反射VPL，这个形状(由Dachsbacher导出的)几乎是一个球体。因此，对于漫反射VPL，我们使用的是包围球，其半径等于椭球的最长半轴。在本文中，包围球的中心点位置和半径分别由 q + ( 1 3 ) 3 4 Φ k π δ ξ n \mathbf{q}+\left(\frac{1}{3}\right)^{\frac{3}{4}} \sqrt{\frac{\Phi k}{\pi \delta \xi}} \mathbf{n} q+(31)43πδξΦk n和 ( 4 27 ) 1 4 Φ k π δ ξ \left(\frac{4}{27}\right)^{\frac{1}{4}} \sqrt{\frac{\Phi k}{\pi \delta \xi}} (274)41πδξΦk 随机给定，其中，k是漫反射系数。

不足和改进

高度镜面的表面

尽管我们的方法通过每一个VPL的贡献来剔除VPL，但渲染的质量会受到剔除前VPL密度的限制，类似于lightcuts。因此，对于高频的BRDF有必要生成更多的VPL(如下图11所示)。对于这种大量VPL的情况来说，着色之前的剔除阶段将会是瓶颈。另外，后处理的去噪会模糊细节处的光焦散。我们的方法不能渲染由完美镜面反射的光焦散(即 α = 0 \alpha = 0 α=0时)。

光滑对光滑的相互反射

尽管本文提高了光滑VPL的采样的效率，但由于采样概率的问题，仍然忽略了着色点处的BRDF。如果着色点处是高度镜面的BRDF(如下图12所示)，这个限制会造成显著的差异(即噪声)。在实际中，对于漫反射对光滑和光滑对光滑的相互反射，应该使用有偏差的近似方法(如屏幕空间的反射，screen-space reflection)，而不是我们的方法。考虑到光滑对光滑的相互反射，未来的研究可能会研究一种无偏差的剔除方法。

各向异性反射瓣

对于由入射方向 ω i \boldsymbol{\omega}_i ωi和法线 n \boldsymbol{n} n定义的平面来说，我们的包围椭球体可以很好地拟合。另一方面，对于方向 ω w \boldsymbol{\omega}_w ωw(与这个平面正交)来说，包围椭球会非常松弛。这是因为微平面BRDF模型的反射瓣对于掠射入射方向(a grazing incoming direction[问题4])是各向异性的。因此，可能要减少半轴 r w r_w rw的空间。此外，我们的包围椭球体没有考虑NDF的各向异性。尽管对于一个各向异性的NDF，包围椭球体可以通过最大的粗糙度来计算，但是这样会形成一个松弛的包围体。根据这种各向异性的反射来缩小包围椭球体也是我们未来的工作。

问题

问题1

如果用辐射率除以概率值的话，不就增强了光源的辐射率了吗？(概率值小于等于1)。这样是否会违背能量守恒定理？且分别将 L ( ω o , l ) L\left(\boldsymbol{\omega}_{o}, l\right) L(ωo,l)和 p ( ω o , l ) p\left(\boldsymbol{\omega}_{o}, l\right) p(ωo,l)带入后计算得出的 L ( ω o , l ) ≈ δ L(\boldsymbol{\omega}_{o}, l) \approx\delta L(ωo,l)≈δ！

a: 1、除以概率值是为了补偿被轮盘截断的光强；2、见上面公式推导

问题2

这个可分离是怎么分离的： G 1 G_1 G1和 G 2 G_2 G2有什么关系，或者怎么将 G 1 G_1 G1用于GGX的BRDF？

a：详见阴影遮罩函数解析。

问题3

为什么是下界呢？

a：因为 ω o \boldsymbol{\omega}_o ωo和 ω u \boldsymbol{\omega}_u ωu之间的夹角为 θ \theta θ，所以当 ω i \boldsymbol{\omega}_i ωi、 ω o \boldsymbol{\omega}_o ωo和 ω u \boldsymbol{\omega}_u ωu三者共面的时候， ω h \boldsymbol{\omega}_h ωh和 n n n之间的夹角最小(为 θ 2 \frac{\theta}{2} 2θ)；当 ω o \boldsymbol{\omega}_o ωo与 ω i \boldsymbol{\omega}_i ωi、 ω u \boldsymbol{\omega}_u ωu不共面时， ω h \boldsymbol{\omega}_h ωh和 n n n之间的夹角就大于三者共面时的 θ 2 \frac{\theta}{2} 2θ，所以 θ 2 \frac{\theta}{2} 2θ就是下界。

思维错误原因：以前认为 θ \theta θ是一个固定值，所以 ω h \boldsymbol{\omega}_h ωh和 n n n之间的夹角最小值是能够为0的，因此 ω h \boldsymbol{\omega}_h ωh和 n n n之间的夹角肯定是可以小于这个固定值的。但是事实上 θ \theta θ是由 ω o \boldsymbol{\omega}_o ωo和 ω u \boldsymbol{\omega}_u ωu确定的，其是一个可变量。但一旦确定了 ω o \boldsymbol{\omega}_o ωo和 ω u \boldsymbol{\omega}_u ωu的方向， θ \theta θ也就确定了，这个时候 ω h \boldsymbol{\omega}_h ωh的方向也就确定了，但是由于 ω o \boldsymbol{\omega}_o ωo可以在球形圆上的任意位置，所以 ω h \boldsymbol{\omega}_h ωh和 n n n之间的夹角大小是处在一个范围内的，这个下界就是这个范围的最小值。

问题4

a grazing incoming direction是啥？

a: 掠射入射角。光从一种介质向另一种介质传播，入射角接近于90度时称之为掠射。注意：一定要从光疏介质（折射率小）向光密介质（折射率大），入射角一定要极其接近于90°。入射角为90°（事实上略小于90°，但在计算时完全可以按90°算）

问题5

这个重要性是怎么判断的？

a: 概率是取1与辐射率采样的较小值，如果某个VPL对该着色点来说重要的话，概率值就为1；如果不重要，就小于1。

问题6

这些信息怎么储存，2D纹理或者立方体纹理都不足够吧？使用3D纹理储存，还是直接保存在OpenGL缓冲中然后直接传给绘制？

a: 应该不是是用纹理存储，类似普通数组的存储方式