6.4 Visibility的烘焙、存储与使用

在前几节中，我们烘焙出了用Volume Lightmap存储的Irradiance，不过这还不够，我们还需要Visibility信息。Visibility就是周围物体对某一点的遮挡信息，因为光照是从四面八方而来的，所以Visibility数据呈现球面分布，即球面Visibility信息。什么是球面Visibility信息？球面Visibility信息可以可视化出来，模型腋下一个点X均匀发射出很多射线（一共720条射线），如果某条射线R没有被任何三角面挡住，则X点在R方向是可见的，否则不可见。图6.8和图6.9分别从不同角度、不同距离画出了X点的可见射线，若某条射线R被挡住了，则标为红色，否则标为黑色。

本节要解决的问题就是如何在渲染时知道模型身上任何一点的球面Visibility信息。由图6.8和图6.9可以知道，任何一点的Visibility信息数据量是很大的，共发射了720条射线，每条射线即使用一个bool值来表达，也是90 byte=22.5个32 bit float的数据量。更严重的是，这些数据的使用，将它们编码传递给Shader并在Shader里解码，是很复杂的。本节要解决的核心问题就是离线计算、编码Visibility信息，在线解码任何一点的球面Visibility信息，同时保证编码占用的存储空间尽量少，解码耗费的计算量尽量小。Visibility补全了渲染的信息，迈向了更为正确的Global Illumination。看一下Global Illumination的Rendering Equation，公式中的V就是Visibility，本节要做的就是预计算V（Ω代表半球，ω是半球的任意单位向量，Li（ω）可以看作ω方向的入射光强度，Lo（ω）是算出来的o方向的出射光强度）。

6.4.1 Visibility数据的生成

Visibility本质上是一个球面数据，为了尽可能准确地获得数据，这里采用光线追踪来计算，大致步骤如下。

（1）对物体及周围场景进行Voxelization。

（2）在整个Mesh表面均匀生成采样点。

（3）对于每个采样点，发射呈球面分布的射线进行光线追踪，并记录下相交数据（是否相交及相交距离）。

（4）对相交数据进行处理并存储。

需要说明的是，步骤（1）在进行Voxelization的时候需要高精度的Voxel，因为后续的光线追踪需要比较精确的碰撞。作者采用的是0.1～1mm大小的Voxel。

步骤（3）在进行光线追踪的时候和前面计算Volume Lightmap没有本质区别，最大的区别就是Voxel精度高得多。另外，需要对射线初始位置进行一定的偏移，因为即使Voxel精度再高，也会有一定的误差，不如BVH来得精确。

步骤（4）是必须的，前面的步骤得到了大量的数据，需要降频存储才能高效利用。下面对步骤（4）进行详细的说明。

6.4.2 Visibility数据降频

球面函数例子如图6.10所示。从数学上来描述，数据降频就是要计算任何一点的一个球面函数：Visibility=f（θ，ϕ），这里采用球面坐标来描述函数。

宏观来说，Visibility其实是一个无限频率的数据，因为它可以不连续：在某个方向，R没有被挡住，但是该方向旁边的方向，R被挡住了。这种情况比比皆是。所以要对Visibility降频必然要损失效果，本节要做的就是保留大致的区域信息。实时渲染使用的一般是Spherical Harmonic、Spherical Gaussian、Besier Surface等。这里选用Spherical Harmonic，因为Spherical Gaussian并不是严格意义的球面基函数，并且拟合困难（为了减少拟合难度，目前的游戏基本都固定了几个方向来进行Nonlinear Fitting，解法不够优雅）；Besier Surface解码烦琐。虽然Spherical Harmonic（下文简称SH）也有问题，如会拟合出负值，会有Ringing现象（数值振荡），计算需要Global Support（每个点的值都需要计算所有项），不过它最大的优点在于它是一个严格的球面基函数，而且频率区分很好，投影简单。用SH来表达Visibility信息已经在论文[11]中描述过，SH只能表达低频的信息，这里选择Band-4的SH，因为Band-3的SH过于低频，丢失的信息太多。依然考虑腋下的位置点，首先用Ray Tracing按照球面均匀发射射线得到Visibility数据，然后把数据投影到Band-4的SH空间中，得到16个SH参数，接着把16个SH参数传给Shader，并在Shader中重构出低频的Visibility信息，如图6.11和图6.12所示，可以看出该点左右分别被手臂和躯干挡住，所以画出来的Visibility左右和上边都是黑的。

当然可以计算更多的例子出来，从实验来看，Visibility用Band-4的SH来编码解码基本能满足要求。但SH的问题还是比较明显的，那就是每个点需要存储16个float值，并且在Shader里需要计算Band-4的SH，计算开销比较大。接下来描述如何继续压缩所需的存储数据和减小还原的计算代价。SH方法能压缩到16个float值，Ambient Dice方法能压缩到12个float值，但是数据量依然太大，想要继续降低存储代价的话，就要考虑更为Lossy的压缩方式了。最好能把数据压缩到每个点4个float值。论文[12]在SH方法的基础上进行了进一步压缩，首先把Ray Tracing的Visibility数据投影到Band-4的SH，然后将其压缩为一个Cone，Cone内部的Visibility不为0，Cone外部的Visibility为0。为什么可以用Cone呢？因为需要计算的是Mesh表面任意点X的Visibility数据，X的背面都是0，如果画出来半球的SH Visibility会发现大部分半球的Visibility都很像一个Cone。用Cone表达Visibility的话，用4个float值就够了，因为Cone的方向用2个float值即可，Cone的Angle用1个float值，Cone的Scale也用1个float值。在SH投影中已经压缩到了Band-4 SH空间，剩下的就是要用一个Cone去拟合一个Band-4 SH。采用Least Square拟合算法，需要拟合的参数过多，既有Cone的中心轴，又有Cone的Angle和Scale，即使使用简单的Least Square算法，这个拟合过程也是复杂的Nonlinear Optimization，而且具有很多局部最优解和鞍点等，很难短时间找到全局最优解。幸好SH有个非常好的性质，就是投影到SH空间之后，SH的最优解方向已经算出来了，即OptimalDir=（-SH[3]，SH[2]，-SH[1]）（参见论文[13]）。如图6.13所示，绿色射线就是计算出来的SH的最优解方向，红色射线和蓝色射线是根据这个最优解方向重新生成的x轴和z轴，可以看出SH的最优解方向只需要SH的线性系数部分就够了。可以用OptimalDir作为Cone的中心轴，这样拟合的时候Cone的中心轴就可以当作常量了，只需要拟合Angle和Scale，大大降低了拟合难度。其实SH的最优解方向（Cone的中心轴方向）可以看作Bent Normal。

Cone的中心轴可以直接用SH的最优解方向，剩下需要拟合的数据就只有Angle和Scale了。这里不能直接用Cone的公式来和SH拟合，因为前者在时域空间，后者在频率空间，拟合是很困难的，需要将Cone投影到Band-4 SH空间来拟合。假设Cone的Angle是α，Scale是S，现在要在新的坐标系中拟合，这个坐标系的正z轴和Cone的中心轴重合（如果不做这个旋转，则Cone投影到SH空间的系数很难算），那么Cone投影到Band-4 SH空间可以写作下式，其中ci（α）是Cone投影后的各项SH值。需要注意的是，从这里开始，之后所使用的坐标轴是右手系的，因为SH的计算使用的是右手系坐标系。

式中，ω是方向向量。每一个ci（α）都可以通过Symbolic Integration工具（参见论文[14]）计算出来，这里在投影Cone的时候要假设Cone内部的值都为1（之后会说明为什么要假设内部值都是1，而不假设离Cone中心轴越远的值会线性下降到0）。投影之前的Cone的公式可以写作：

VCone(x)=S(0≤x≤α)

那么投影之后的系数计算为：

类似地，可以算出：

c1,3,4,5,7,8,9,10,11,13,14,15(α)=0

45°Cone投影前后如图6.14所示，左侧是45°的投影前的Cone，右侧是投影后还原的Cone，右侧下半球的白色条带就是SH固有的Ringing缺点，不过位于下半球，没有关系。

图6.15给出了10°、30°、45°、60°的Scale为1的Cone投影到SH空间重建出来的结果，所以用Cone拟合可以直观看作用不同Angle的Cone投影后的结果比对原始SH数据，从而选择一个球面所有方向的数据和SH数据最相像的Angle。

上述函数的导数很容易求出，这里就省略了，表示为。为了拟合，之前计算出来的Visibility的SH值也要转到Cone所在的坐标系，该坐标系的正z轴即（-SH[3]，SH[2]，-SH[1]）。这就需要用到SH的Rotation矩阵，SH的Rotation矩阵虽然是稀疏矩阵，不过并不好计算，一般都是通过低Band的SH Rotation矩阵递归计算出来的。工程上可以通过之前计算出来的SH参数重建Visibility，然后重新投影到新的坐标系来绕开计算Rotation矩阵。计算出新坐标系下的SH参数之后，SH空间的Visibility可以表示为：

接下来要用Least Square来最小化VCone（S，α，ω）和VSH（ω）之间的差值。数学上就是最小化E：

对α和S分别求导：

由SH的Orthonomality：

可以推出：

把它们都设为0，可以求得极值点条件等式为：

只要求得了α，由式（2）就可以马上算出S的最优解。式（1）看起来很复杂，其实只是一个一维的Nonlinear Equation，最简单的解法是Bisection算法（参见论文[15]），先找到一个区间[x1，x2]，使得式（1）在x1和x2处取不同的正负号，再二分缩小这个区间，保证式（1）在区间两端点处一直取不同的正负号，因为式（1）是连续的，因此必然跨过零点，这样就可以得到任意精度的解了。图6.16所示为SH拟合到Cone的例子。

这里说下为什么投影前Cone采用的是VCone（x）=S（0≤x≤α）而不是，后者直观来说就是Cone中心轴位置Visibility为1，随着角度变大，Visibility线性下降，当和中心轴角度为α时降为0。我们尝试重新投影到SH空间，如c2（α）：

对比之前的c2（α），可以看出复杂很多，更高频的项会更复杂，拟合起来更困难。在拟合的时候最重要的是Cone的形状，即Angle的计算，所以在拟合的时候假设Cone内部都是S×1是可行的，算出来后，Shader在用Cone数据还原Visibility的时候可以进行衰减（不过要注意能量守恒）。

6.4.3 Visibility数据的存储

有了任何一点的Visibility的计算算法，剩下的就是如何存储的问题。一共有3种存储选择，如下。

● 为每个顶点计算Visibility，并存储在顶点数据里。

● 在Mesh表面计算像素级数据，并存储在贴图里。

● 按照Voxel计算，存储在Volume Texture里。

为了尽量减少带宽压力，我们选择存储在顶点数据里，这里用Cone编码进行实验。本节尝试了直接针对每一个顶点计算Visibility数据，并存放在顶点数据里，如图6.17所示，画出了插值出的每个像素的Cone.Scale×Cone.Angle/π，越黑的地方，Cone的Scale和Angle越小。结果发现，数据不连续的地方很多，表现出来就是画面比较脏。这是一个比较经典的问题，Mesh的面数越低，问题越明显，本质上这是因为Vertex的数据是对像素级数据的采样，低于了Nyquist极限。从工程上来说，这是因为三角面内的Visibility数据都是通过Barycentric Coordiante插值3个顶点的Visibility数据得到的，而顶点的Visibility的计算是局部的，并没有考虑到自己的数据会被拿来插值，所以出现错误数据是不可避免的。

论文[16]描述了解决方案，它研究了如何将像素级任意Signal存入Vertex数据使得三角面插值还原的数据与原始像素Signal数据之间的差距最小。这是个最优化问题，论文采用的是Least Square算法。假设Mesh的顶点集合是vi，…，vN，硬件插值用数学表达出来就是下面的公式，其中vivjvk是Mesh上的任意三角形；xi，xj，xk是对应顶点vi，vj，vk上的值，可以是任意维度的向量；Bi（p），Bj（p），Bk（p）是插值函数，叫作Linear Hat Function，如图6.18和图6.19所示。其实很简单，例如，对于Bi（p）来说，它的值只有在三角形vivjvk内部才非0，并且Bi（vi）=1，Bi（vj）=0，Bi（vk）=0，Bi（p）在三角形内部是线性下降的。p是三角形vivjvk内部任意一点。注意，对于一个有N个顶点的Mesh来说，一共存在N个对应的Bi函数。所以下面的公式中∑的上限是N。对三角形vivjvk来说，除Bi（p），Bj（p），Bk（p）外，其他函数的值都是0。为什么采用这种表达式呢？因为这种全局表达的式子可以利用矩阵运算的能力规约为矩阵问题。

有了上面的理解，现在要解决的问题用数学表达出来就是：给定一个定义于Mesh表面S的函数f（p），其中p是S中的任意一点，x=（x1，…，xN）T使得根据g（p）插值出来的所有值的和f（p）差距最小。怎么定义差距呢？这里使用Linear Least Square算法中的平方差距离，也就是要最小化平方差距离。

需要注意的是，f（p）是已知的，如对于烘焙颜色贴图到顶点上来说，f（p）就是p点的颜色值；对于烘焙AO项来说，f（p）就是离线光线追踪算出来的p点的AO值。

代入g（p）得：

可以看出，式（1）中第一项是关于x的二次项，第二项是关于x的一次项，第三项是常数，用矩阵形式写出来就是：

其中：

式（2）看起来很复杂，但其实就是一个二次函数，可以证明矩阵A是一个稀疏对称正定矩阵，所以这个二次函数是一个单调递增函数，并且当它的Gradient等于0时取得最小值，也就是当所有都取0时得到最优解，所有的式子合并在一起即：

Ax= b

因为对于vi，…，vN来说，Bi（p），Bj（p），Bk（p）都是线性的，而且当l！=i，j，k，Bl（p）都是0，所以可以对每个三角形分开求积分再求和。其中，Ai，j和f（p）无关，可以手算或用积分工具算出，可以求得Ai，i就是包含vi的所有三角形面积和的1/6，如图6.20中围绕vi的一圈三角形。Ai，j（i！=j）就是包含边vivj的三角形面积和的1/12，如图6.20中的三角形1和三角形2。

对于bi和c来说，因为都包含f（p），所以需要用蒙特卡罗方法对每个三角形进行f（p）采样，进而求得积分：

式中，iμ是包含vi的所有三角形面积和；Xi是这些三角形内的采样点集；|Xi|是采样点数量。c可以用同样的方法算出。

众所周知，用Least Square来进行Regression会造成Overfitting，需要加入Regularization项。形式就是：

(A+αR)x=b

R矩阵就是Regularization项。R矩阵的选取不是那么容易的，不能影响到算出来的x的绝对值大小，如果用最简单的矩阵，即A自己，这样会造成x的绝对值变小，考虑极端情况，当参数α很大时，x趋向于0，明显不对。想一下我们的需求，一是要让插值出来的值尽量逼近原始值，二是要比较“平滑”。怎么表达“平滑”这个需求呢？考虑两个相邻的三角形，它们公用两个顶点，当硬件插值时，这两个三角形有两个值是公用的，为了使得这两个三角形在交接的地方不出现明显的硬边和不连续的现象，要求这两个三角形的g（p）函数的变化率尽量一致，也就是g（p）的Gradient，即∇g（p）要尽量一致。之前的Error Function就应该是：

式中，t，u表示两个相邻的三角形。第二项表明相邻的三角形的Gradient差值越大，Error越大，在逼近原始值的同时要尽量减小全局Gradient差值大小。

现在的问题就是怎么求∇g（p），考察某个三角形t，先展开vi，vj，vk：

∇g( p)=x∇iBi(p)+x∇jBj(p)+xk∇Bk(p)

因为Bi，Bk，Bk都是线性函数，所以∇Bi，∇Bj，∇Bk都是常量Vector。现考察Bi，如图6.21所示，通过形状来看能猜到∇Bi应该和vjvk垂直并指向三角形内部，因为这是数值上升最快的方向。

严格推导如下。如图6.22所示，因为Bi（p）是相对于vi的重心值，所以：

式中，at是三角形t的面积。这个公式中的变量只有hi，所以Bi（p）是hi的函数：

式中，ui是∇Bi（p）的方向。所以：

为了将积分表示成矩阵形式，这里把∇g（p）写成矩阵形式，由上式可知∇g（p）是xi，xj，xk的线性函数，所以可以写成：

其中，第一个矩阵的大小是3×N，并且是个非常稀疏的矩阵，只有和xi，xj，xk相关的值才非0。简化为：

∇g ( p)= Fx

现在可以计算Regularization了：

所以Regularization矩阵就是：

R=(Ft -Fu)T(Ft -Fu)

实现时并不需要真正计算每个很大的稀疏矩阵F，因为每个F只和两个相邻的三角形有关，可以只计算一个很小的4×4矩阵，然后填入矩阵R。

看一个Diffuse贴图烘焙到顶点上的例子，如图6.23所示，左侧是原始效果，中间是顶点Point Sample Diffuse贴图的结果，右侧是使用α=0.01拟合出来的结果。和Point Sample的结果对比，可以看出，肩甲和脸部的错误已经被修正了，同时大大减弱了不连续的Mach Banding现象。

6.4.4 Visibility数据的使用

本节通过3个例子介绍Visibility在渲染中起到的一些作用。由6.4.3节可知，Visibility数据是存放在每个顶点里面的，形式为（Cone Angle，Cone Axis，Cone Scale），分别表示锥的开口角度、中心轴方向及缩放值，中心轴方向就是该点未被遮挡的方向。虽然锥是比较粗糙的遮挡数据，但是因为游戏计算时的性能问题，不可能存放更加精确的遮挡信息（如Spherical Harmonic编码的遮挡信息）。每个Pixel拿到的是插值过后的Visibility数据，从而能判断出大致的可见性信息，知道光照从何处而来。

6.4.4.1 Ambient Occlusion计算

这里采用Uniform-Weighted Ambient Occlusion的定义：

Cone内部的Visibility为1，代表没有遮挡；Cone外部的Visibility均为0。另外，假设Cone的Angle为a，Cone的scale为s，代入上面的积分可得：

图6.24所示为《王者荣耀》红方基地仅显示AO项的示意图。

6.4.4.2 Specular Occlusion计算

在PBR中，IBL存储了环境高光信息，Specular Occlusion（SO）就是对IBL的遮挡信息。那为什么不直接使用Ambient Occlusion（AO）作为SO呢？因为AO只是一个降维的数值，没有方向性，它的计算依赖于假设——物体是Diffuse的，而高光是有方向的，光源或视角变化，遮挡信息就会变化，相应地，SO就会变化，而AO不会变化。

来看怎么计算Ground Truth SO。按照论文[17]中的描述，fr是BRDF；V（wi）表示这条入射光线是否被遮挡，被遮挡为0，否则为1；Ω是整个球面。

考察这个积分可知，当Visibility用Cone来表达，并且BRDF用Phong或Microfacet BRDF来表达时，这个积分其实可以预计算成一张LUT表，输入是View Dir的方向、Cone Angle、Cone Axis和Reflection Dir的夹角。

这就是该论文采用的方式，LUT表最后可以存成一张3D纹理图，Runtime时通过View Dir及Cone Angle和Reflection Dir的夹角，就可以从这张图中取出预计算出的遮挡数据了。

从上面的GTSO描述可以看出，虽然SO计算得比较精确，但是Runtime时需要采样一张3D纹理图，这对手机端来说是不小的开销。那怎么省下这张图呢？必须要抛弃一部分精度。

如图6.25所示，仔细看BRDF的形状，其实和一个Cone的形状很类似，如果BRDF可以用Cone来表示，那么就可以在Runtime时求Visibility的Cone和BRDF的Cone相交的Solid Angle，再除以BRDF的Solid Angle，就能算出SO了。

怎么求BRDF拟合出来的Cone Angle呢？这个Angle大小其实和该点的粗糙度有关，粗糙度越大，Angle越大。对于GGX，论文[18]给出了一个简单的拟合公式。

那怎么求Cone和Cone相交的Solid Angle呢？精确计算Cone与Cone相交的Solid Angle用论文[19]中的算法虽然办得到，但是计算量过大，使用了大量三角函数，为了进一步减小开销，可以舍弃一部分的精度，求一个拟合的Cone和Cone相交的Solid Angle的算法。这里主要采用了论文[20]里面的算法，虽然是拟合解，但是拟合得很好，损失是很小的，最大的Error来自前面用Cone拟合BRDF而不来自这里，具体的相交代码在原始论文里比较详细，读者可以参考。

图6.26对有无Specular Occlusion进行了对比，左侧没有Specular Occlusion，右侧有Specular Occulusion。为了对比更清晰，这里把整个人都设置为了“金属”。可以看出，加入Specular Occlusion后，解决了大部分的漏光问题，模型更为真实了。

6.4.4.3 Irradiance计算

前面介绍了采样Volume Lightmap的方法，即SampleVolumeLightMap3D函数，之前使用Normal去采样，有了Visibility Cone Axis之后，可以使用Axis去采样，这样得出的Irradiance更为准确，因为考虑了可视性信息。