GAMES202笔记：实时阴影

Shadow Mapping

Shadow Mapping是一个Two-Pass的算法。
Pass1从光源出发进行“渲染”（此处并非普遍意义上的渲染，其Buffer中需要记录的只是深度），记录从光源向各方向发射光线打到场景物体的最小深度，由此得到Shadow Map
ShadowMappingPass1

Pass2从相机出发看向场景，针对场景中每一个被Camera Ray打到的点，将其连向相机得到该点到相机的深度。再与Pass1中得到的Shadow Map对应方向的深度进行比较，如果两个深度相等则从光源到该Shading Point之间无遮挡，如果Shading Point到Light的深度大于Shadow Map中对应方向的深度，则光源到该Shading Point之间有遮挡，该Shading Point即可视为在阴影中。
ShadowMappingPass2

PS：如GAMES101中对变换部分的叙述，在渲染时需要先将场景的Frustum变换为Cuboid，然后变换至归一化设备坐标(NDC)。介于近平面和远平面之间的点，进行视锥体挤压后，其坐标Z的值会变小，向远处移动。因此Shadow Mapping的Pass2中深度的比较需要统一坐标，都用变换后的或都用实际深度。
FrustumToCuboid

Shadow Mapping存在的问题

自遮挡现象(Self Occlusion)

Shadow Map中记录的场景深度，其中每一个“像素”深度相同，也就是说，即便某个场景物体表面是平滑的连续的Shadow Map记录的深度也并非是连续的，如下右图可见，其记录的深度是阶梯变化的。
那么，如下右图所示的Camera Ray打到的Shading Point的深度便大于Shadow Map中对应方向的深度，根据前文所述的原理，该Shading Point应被判定为在阴影中，这显然与事实不符。下左图中地面上纹路状的Artifacts便是这样产生的。
ShadowMappingIssues1

当光源在正上方垂直照向地面时，上述问题相对较小；相应地，当光源接近掠射角(Grazing Angle)时，上述问题最为严重。

解决办法
可以加入一个Bias，即比较上述两个深度时候可容忍的偏差阈值，当偏差小于阈值时就不算做被遮挡，如下图
这个Bias可以设置为根据角度变化，比如垂直角度Bias较小，Grazing Angle下Bias较大
ShadowMappingIssues2

但这随之造成了一个新的问题——Detached Shadow
由于深度偏差低于阈值的Shading Point被视为未遮挡，则场景中一些距离很近的物体或是相互接触的部分会因此被误判为没有被遮挡，如上左图中人物鞋子与地面接触的位置。

针对上述问题，有一种叫做Second-Depth Shadow Mapping的方法
其在Pass1生成Shadow Map的过程中不仅记录最小深度，也同时记录第二小的深度，根据这两个深度计算出对应方向的一个中间深度（平均深度），在后续Pass2中就与这个中间深度进行比较，深度小于该中间深度的都视为未遮挡，大于该中间深度的都视为被遮挡。（此时不再需要Bias）
SecondDepthShadowMapping

Second-Depth Shadow Mapping的想法很好，然而并不实用
首先，Second-Depth Shadow Mapping要想得以实现，必须保证所有物体的模型都是Watertight的，即某个物体有正面就必须有反面，这就要求即便是薄如纸张的物体也必须要有一个确切的厚度。因此在实际渲染中，很难保证所有物体的模型都是Watertight的。
再者，Second-Depth Shadow Mapping在实现的过程中需要记录最小深度和次小深度。时刻记录最小深度是很容易的，但是要同时记录次小深度就会涉及排序的问题，代码的实现中可能会多出许多诸如swap的操作。最终导致渲染的实时性难以保证。
PS：即便涉及排序，Second-Depth Shadow Mapping算法的复杂度仍然是O(n)，也就是说不需要再多的Pass。但是这依然可能严重拖慢RTR的速度，这也就是所谓的"RTR does not trust in Complexity"，即实时渲染不相信复杂度，只相信绝对速度。通俗地理解，可以想象一个O(2n)的算法，它同样是O(n)的复杂度，但是相比真正O(n)的算法它却差不多慢了一倍，这在RTR中是绝对无法接受的。

Aliasing

如前文所述，生成Shadow Mapping的过程也可看作一种“渲染”，因此Shadow Mapping同样也有分辨率，当其分辨率不足时，最终渲染产生的阴影便会出现走样（锯齿）现象。
ShadowMappingAliasing

RTR中的重要数学原理

重要不等式

Schwarz不等式

[\int_{a}^bf(x)g(x)\mathrm{d}x]^2\leq\int_{a}^bf^2(x)\mathrm{d}x\cdot\int_{a}^bg^2(x)\mathrm{d}x

Minkowski不等式

\{\int_{a}^b[f(x)+g(x)]^2\mathrm{d}x\}^{1/2}\leq[\int_{a}^bf(x)^2\mathrm{d}x]^{1/2}+[\int_{a}^bg(x)^2\mathrm{d}x]^{1/2}

RTR中的近似

在RTR中，相比不等关系，我们更关系近似相等
RTR中有一个最为常用的近似关系

\int_{\Omega}f(x)g(x)\mathrm{d}x\approx\int_{\Omega}f(x)\mathrm{d}x\bigg/\int_{\Omega}\mathrm{d}x\cdot\int_{\Omega}g(x)\mathrm{d}x

上式在数学上严格来说绝对是错误的，但是在一定条件下可以使这一近似相等关系相对比较准确（或者说使其误差在可接受范围内），其条件如下（在RTR中，下列条件满足其一便视为可以使用上述近似）：

当 $g(x)$ 的积分域很小时；
当 $g(x)$ 比较smooth时（此处的smooth并非所谓的连续，而是要求其在积分域中变化比较小、比较低频）

RTR中对渲染方程的近似

渲染方程如下：

L_o(p,\omega_o)=\int_{\Omega^+}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)\cos\theta_iV(p,w_i)\mathrm{d}\omega_i

式中， $L_o$ 为出射的Radiance， $L_i$ 为入射的Irradiance， $f_r$ 为BRDF， $V$ 为Visibility

利用上一节中的近似关系，对渲染方程进行如下近似：

L_o(p,\omega_i)\approx\int_{\Omega^+}V(p,\omega_i)\mathrm{d}\omega_i\bigg/\int_{\Omega^+}\mathrm{d}\omega_i\cdot\int_{\Omega^+}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)\cos\theta_i\mathrm{d}\omega_i

上式相当于将Rendering Equation中的Visibility项拆了出来，而右侧除了Visibility项剩下的积分做的正是Shading(即不考虑Shadow)，渲染方程就变成了先做Shading然后再乘上Visibility项

要使得上述近似较为准确，需要满足如下条件之一：

积分域较小（常见于Light/Directional Lighting，此时相当于Visibility项不用积分，做硬阴影）
除Visibility项以外剩下的积分中被积函数比较Smooth（常见于Diffuse BRDF或者Constant Radiance Area Lighting）

当然，这本身就是近似，即便上述条件一个都不满足，强行使用也不是不行，比如Ambient Occlusion(环境光遮蔽)就会用到

Percentage Closer Soft Shadows

前文所述的Shadow Mapping使我们得以实现硬阴影，然而，现实中的光源通常为面光源，而因此生成的阴影通常为软阴影，如太阳光照向地球便是最为典型的例子
SoftShadows

在正式进入PCSS前，需要先引入一个常用的“工具”——Percentage Closer Filtering

Percentage Closer Filtering(PCF)

PCF在发明之初并非是为了生成软阴影，而是对阴影边缘进行抗锯齿，只是人们后来发现利用它的思想方法也可以也用来生成软阴影

PCF是对Shadow Mapping的Pass2中深度比较的结果进行Filter
特别注意：

PCF并非Filter Shadow Mapping生成的结果（即生成的阴影），正如在反走样中的道理，不能先得到一个走样的结果，再在这个结果上进行模糊、平均
PCF也不是Filter生成的Shadow Map，因为Shadow Map中的深度被用于在Pass2中与Shading Point连向光源所得的深度进行比较，这个比较结果无论如何都是非零即一的，因此对Shadow Map进行Filter是没有意义的

PCF实际上进行的工作是，对于一个Shading Point，将其投影会光源处后，对其周围的 $n\times n$ 个像素进行深度比较（每一个比较的结果都是非零即一的），然后最这些结果进行平均（也可加权平均）作为该Shading Point深度比较的结果。
PCF

Percentage Closer Shadow Mapping(PCSS)

根据PCF的原理，可以自然而然地想到，当进行PCF的Kernel(不确定是否可以称为Kernel，但由于很像CV里的卷积核的概念，故暂且这么叫)为 $1\times1$ ，就相当于没有进行PCF而做传统的Shadow Mapping，此时生成的是硬阴影；而当Kernel很大时，阴影的锯齿变少了，但是阴影也变模糊了，就形成了类似软阴影的效果。
也就是说，Kernel越小，阴影越“硬”；Kernel越大，阴影越“软”

利用这种特性，可以对需要产生硬阴影的位置使用较小的Kernel，对需要产生软阴影的位置使用较大的Kernel。这便是PCSS的基本思想

那么，何处应当形成硬阴影，何处应当形成软阴影呢
考察下面的例子，可以得到一个基本观察
Pen
笔尖处的阴影更“硬”，笔杆处的阴影更“软”
推而广之，也就是说当阴影的接收物与投射物之间的距离较近时，生成较硬的阴影；当阴影的接收物与投射物之间的距离较远时，生成较软的阴影；

结合上述的几个观察可知，所使用的Kernel Size与Blocker Distance(即上述的阴影接收物与投射物之间的距离)相关，更确切地说，是与Blocker相对的、平均的投影深度有关（这是因为阴影具体的软硬程度还与光源与Blocker以及光源大小等因素相关）

如下图所示，可以有这样的关系：

\omega_{Penumbra} = (d_{Receiver}-d_{Blocker})\cdot\omega_{Light}/d_{Blocker}

DepthPenumbraMathematicalTranslation

但是，上图所示的是一种比较取巧的情况，即Blocker是一个平面（没有厚度）且完全与Shading Point所在平面平行。现实中的Blocker通常有一定厚度且其整体不一定与Shading Point所在平面平行。
所以通常我们说Blocker Depth是一个Average Blocker Depth
准确地说，对于一个Shading Point，考察其投影至Shadow Map后周围一小块范围内记录的有遮挡的深度的平均值，该值即为Average Blocker Depth
注意：上面特别标出了“有遮挡”，即该点在Shadow Map对应像素的邻域中，记录的深度小于Shading Point的深度（即被遮挡）的值才被纳入平均值计算，因为没有遮挡的点本身就没有Blocker和Blocker Depth，将其也纳入计算就不准确了

根据以上思路，将PCSS的步骤总结如下：

Blocker Search
对于Shading Point对应的Shadow Map上的像素，在一小块范围内计算记录的Average Blocker Depth
Penumbra Estimation
根据第1步得到的Blocker Depth计算伴影大小（即PCF的Kernel的大小）
进行PCF

到此，我们还剩最后一个问题，那就是应该在多大的范围内进行Blocker Search
这里可以定义一个固定大小的Blocker Search Region，但是有更好的方法：
首先说明，Shadow Map的生成就如我们进行渲染，也需要定义一个Viewport。由此，从光源（对于面光源则是其中心点）看向的场景也会在被远近两个平面切割出来的Frustum内，而Shadow Map就成像于Viewport(即近平面)上。
将Shading Point连向面光源，也会形成一个锥体，该锥体被Shadow Map所在平面截断，在Shadow Map上截出一块区域，该区域就可作为我们所需的Blocker Search Region
BlockerSearchRegion

根据上述原理可见，PCSS方法的开销是十分大的

A Deeper Look

PCF可以表示成如下的卷积形式：

[\omega*f](p)=\sum_{q\in N(p)}\omega(p,q)f(q)

其中， $\omega$ 为Region内各点的权重

于是，在PCSS中，Visibility项就可以表示如下：

V(x)=\sum_{q\in N(p)}\omega(p,q)\cdot\chi^+[D_{SM}(q)-D_{scene(x)}]

其中， $\chi^+$ 是符号函数， $D_{SM}(q)-D_{scene(x)}$ 是进行深度比较
从上式也可以清晰地看出，PCF并非Filter Shadow Map，也不是在图像空间Filter生成的阴影结果，而是Filter深度比较后的结果，再以此生成阴影

\begin{aligned} V(x)&\neq\chi^+\{[\omega*D_{SM}](q)-D_{scene}(x) \} \\ V(x)&\neq\sum_{q\in N(p)}\omega(p,q)V(q) \end{aligned}

补充说明

PCSS的第一步和第三步都需要在一定大小的Region里进行Filter，因此开销比较大。需要生成越软的阴影，所需的Region越大，开销也越大。

针对这个问题，可以对Region里的像素进行稀疏采样，采样某些点进行加权平均，而不是依次加权其中所有像素求平均。但是这种采样会使得结果相对较Noisy，因此工业界通常将上述的采样与最终在图像空间对阴影结果降噪相结合，以得到较好的结果

Variance Soft Shadow Mapping(VSSM)

如前文所述，PCSS主要慢就慢在1、3两步中在像素的邻域进行filter，而VSSM就是针对性地解决这个问题

VSSM的基本思想

首先考虑这样一个情形：我所在的班级进行了一场考试，我想要知道本次考试我的成绩在班里排百分之几，那么我便需要知道有多少同学考得比我好（或差），最简单直接的办法就是与班里其他每位同学一一比较
这也正是PCSS中计算Average Blocker Depth所做的事情

在这种情形下，如果我能得到一个班级同学考试成绩的统计直方图，那么我便不再需要与同学一一比较成绩，直接根据自己的成绩在直方图中一对照就可以知道我排在百分之几
再进一步，如果对精度的要求不那么高，我甚至可以将班级同学的成绩看作是正态分布
NormalDistribution

套用至PCSS中，就是将Region内所有像素的深度值当成是成绩，以此拟合一个正态分布，来得到当前Shading Point对应的Shadow Map上像素的深度值在其中排在百分之几

VSSM的具体实现方法

为了拟合出正态分布，我们需要确定其均值与方差

Mean

对于均值，就是在Region中求出所有像素的深度的平均，这种对于一块方形区域，欲快速得到其所有像素的平均，这便是Mipmap适合做的事情
但是正如GAMES101中所述，Mipmap并不太准确，因为它需要在不同层级间插值，并且Mipmap只能用于正方形区域，而不适合用于长方形区域

而对于一个2D矩形区域，有一种用于求其均值的、更为精准的数据结构，称为Summed Area Tables(SAT)

Variance

对于2D矩形区域，如何快速求其方差呢？
这就需要用到概率论当中非常经典的一个公式：

Var(X)=E(X^2)-E^2(X)

其中， $E(X)$ 项在前一步计算均值时就已经得到。而对于 $E^2(X)$ 项，可以在生成Shadow Map时同时生成一张记录 $Depth^2$ 的Shadow Map
注：这里看似用了两倍的空间以换取时间，但是实际上在OpenGL中，Shadow Map也是一张图像，由各个像素组成，而像素值就像在RGBA空间一样用四个通道表示，也就是说我们先前生成的Shadow Map记录深度值只用了其中的一个通道，对于记录 $Depth^2$ 的Shadow Map，我们完全可以再使用一个通道来记录其值，所以本质上只是用了一张图的两个通道来分别表示两种Shadow Map，并不用额外增加空间。所以额外生成记录 $Depth^2$ 的Shadow Map只需要顺便算一下各个深度值的平方这一点点额外开销

PDF->CDF

经过前述步骤，我们可以得到一个正态分布的概率分布函数(Probability Distribution Function, PDF)，而需要知道当前像素在该Region内排在百分之几，就是计算其积分，即累积分布函数(Cumulative Distribution Function, CDF)

但是实际使用时，我们并不可能真的一一计算其积分，那样将非常耗时
对于正态分布有一个重要工具——误差函数/高斯误差函数(Error Function/Gauss Error Function, erf)
正态分布的累积分布函数通常可以表示为：

\Phi(x)=1/2(1+erf((x-\mu)/\sigma\sqrt{2}))

而其中的erf在工程上通常可以打表，用查表的方式快速得到，如下（数值解而非解析解）：
erf

Chebychev’s inequality

然而，VSSM比前文所述的更为粗暴，它使用了切比雪夫不等式：

P(x>t)\leq\sigma^2/[\sigma^2+(t-\mu)^2]

这样一来，甚至都不需要假设其为正态分布，只要得到其均值和方差即可得到当前像素值排在前百分之几，即下图的阴影区域面积
Chebychev
在RTR中，我们更关注相等而非不等关系，这里VSSM又做了一个大胆的假设，就是将切比雪夫不等式告诉我们的上界直接当成其估计值，即将不等关系直接改为近似相等关系

注：切比雪夫不等式在使用时需满足一个前提，就是 $t>\mu$ ，但是即便只有半边是好用的，它还是被十分普遍地使用了，管它呢，搞图形学就得敢做赌狗

VSSM性能分析

Mean of depth in range: O(1)
Mean of depth square in range: O(1)
Chebychev: O(1)
No samples/loops needed

Blocker Search优化

前述的内容只是解决了Shadow Mapping的Pass3中深度比较太慢的问题，Pass1中计算遮挡物平均深度的速度也需要进行优化

在如下的Region中，我们可以得到像素的平均深度( $Z_{Avg}$ )、遮挡物的平均深度( $\textcolor{RoyalBlue}{Z_{OCC}}$ )以及非遮挡物的平均深度( $\textcolor{red}{Z_{UNOCC}}$ )
VSSMblockerSearch

假设Region中像素数量为 $N$ ，非遮挡物像素数量为 $N_1$ ，遮挡物像素数量为 $N_2$ ，则可以得到如下关系：

\textcolor{red}{Z_{UNOCC}}\cdot N_1/N + \textcolor{RoyalBlue}{Z_{OCC}}\cdot N_2/N=Z_{Avg}

此时， $N_1/N$ 的值表示大于当前Shading Point深度t的像素的比例（亦即班里成绩高于我的同学的比例）
于是又可以利用切比雪夫不等式近似求得，即 $N_1/N=P(x>t)$ ，而 $N_2/N=1-P(x>t)$

$Z_{Avg}$ 如前文所述可以用MipMap之类的方法快速得到，此时要求得 $\textcolor{RoyalBlue}{Z_{OCC}}$ ，我们还需要知道 $\textcolor{red}{Z_{UNOCC}}$ ，此处再次进行一个大胆的假设，就是认为非遮挡物的平均深度等于Shading Point的深度，即 $\textcolor{red}{Z_{UNOCC}}=t$

注：RTR就是这样，很多方法其实非常Hack。但是记住，“图形学第一定律(The First Law of Computer Graphics)告诉我们：If it looks right, it is right. 我们做RTR更多追求的是视觉上看起来正确，而不需要做到数学和物理上严格的正确。所以，方法很Hack，假设很大胆，都没有关系，有用就完事了

Summed Area Tables for Range Query

关于MipMap的原理和使用，以及其在范围查询时候的缺点参见GAMES101，此处不再赘述

SAT最核心的就是使用了前缀和
一维情况如下所示
SAT1D

二维情况下如下所示，其当前像素记录的前缀和就是左上角矩形中所有元素的和
SAT2D

SAT其实就是增加了O(n)的空间并使用O(n)的时间来记录一个前缀和数组用以进行范围查询
其中在二维情况下，计算前缀和需要计算行与列各一次的一维前缀和，其中每一行之间、每一列之间的计算不会相互影响，所以可以并行

Moment Shadow Mapping

VSSM的问题

首先，我们来观察一下VSSM存在什么问题
根据前文所述，VSSM为了提高计算速度而做了很多大胆的假设，但是当实际情况与其假设相差较多时，便会出现问题
如下左图所示，遮挡情况较复杂时，假设遮挡物深度的分布为正态分布是没问题的；但是如下右图所示，遮挡情况较简单时，假设其为正态分布显然会相差较多
VSSMproblem

估计的分布与实际分布相差较多时，可能会因为估得的深度大于Shading Point深度的像素比实际少而导致阴影变暗，也可能会因为估得的深度大于Shading Point深度的像素比实际多而导致阴影过亮(Light Leaking/Light Bleeding)
LightLeakingDis

如下图所述，车身中间是镂空的，当光源在上方时，挡住车底的Shading Points的遮挡物深度分布显然会与正态分布相差甚远，因此车底可以看见Light Leaking现象
LightLeaking

除了漏光问题，如前所述，切比雪夫不等式只有右半边比较准确，这也是VSSM中的一个问题

Moment Shadow Mapping的原理

Moment Shadow Mapping就是为了解决VSSM中估计深度分布不准确的问题，以求得到一个相对更接近实际的分布
这里的Moment是矩的意思，Moment Shadow Mapping就是使用更高阶的矩来描述分布，以获得更好的估计（VSSM实际上就是使用了一阶矩和二阶矩）

这里有一个结论：使用前 $m$ 阶矩来表示一个用一系列阶跃函数堆叠起来的函数，最多可以表示 $m/2$ 个台阶

通俗来说，就是用前 $m$ 阶矩，将分布表示成展开式形式
通常，使用前4阶矩就可以得到相对较好的结果

Distance Field Soft Shadows

Distance Functions(Signed Distance Functions)

关于距离函数，在GAMES101中已经有介绍，这里仅做简单复习
距离函数定义了某一点到物体表面的最小距离
DistanceFieldA

距离函数可以用于Blend运动的边界，可以得到较好的几何过渡
SDFblend1
SDFblend2

Distance Filed的应用

Ray Marching

假设已经得到了场景的SDF，则对于任意一点，其SDF的值告诉我们一个“安全距离”，即从该点向任意方向Trace一个小于等于其“安全距离”的步长，都不会打到其他物体
RayMarching
注：SDF可以用于运动的刚体，但对于形变物体就不太适用

我们还可以利用SDF计算出Shading Point的Safe Angle，如下图所示，从Shading Point向某一方向看去，通过其路径上的点的SDF值，我们可以计算出从该路径向任意方向偏转，在多大的角度以内都不会打到遮挡物
SDFpercentageofocclusion

而对于每个Shading Point，其Safe Angle越小，意味着其Visibility越小，该Shading Point就越黑
SDF->Visibility并不是一个准确的方法，但是符合我们对Visibility的观察，因而可以获得较好的效果

Safe Angle的计算如下图所示：
CalculateSafeAngle

从Shading Point向某一方向出发，Trace与其SDF值相等的长度得到下一个点（图中p1），再根据p1的SDF继续Trace到p2，以此类推。
对于这条Path上得到的每一个点，都可以以该点为圆心，SDF为半径画出一个球面，从Shading Point向每一个球面作切线，每一条切线都与Path形成一个夹角 $\theta$ ，其中最小的 $\theta$ 即为Safe Angle

此处， $\theta=\arcsin[SDF(p)\big/\lVert p-o \rVert]$
虽然只需要计算一个 $\arcsin$ ，但是反三角函数通常计算量较大，因此用下式来近似表示Visibility

\min\{k\cdot SDF(p)\big/\lVert p-o \rVert,0\}

上式的逻辑就是， $SDF(p)\big/\lVert p-o \rVert$ 的值的大小已经可以说明Visibility的相对大小了，该方法本身就不是完全精确的，完全准确计算出 $\theta$ 其实也没必要。除此之外，上式中的 $k$ 还有一个重要作用——控制阴影的软硬程度，k越大，Visibility从0到1的过渡带范围越小，也就是说生成阴影的伴影会更早截止，该阴影就越硬

将SDF可视化后为下图所示效果：
SDFvisualization

Pros and Cons of Distance Field Soft Shadow
Distance Field Soft Shadow相对比较快速（此处忽略渲染前生成SDF的时间），且质量较高
但是Distance Field Soft Shadow需要大量的预计算，且需要较大的存储空间