Games202学习小结 | wyzwzz 的博客

这一周刚好七天重新写了一边Games202的前四个作业，把以前有漏洞或者没有涉及的Bouns部分也全部解决和实现，显示的效果还不错，代码量（包括Shader）估计在3500行左右，这还是OpenGL包装过以及使用一些小轮子的情况下…lambda函数嵌套用起来真优雅，不过编译报错也是看的人废了。

仓库链接

最终实现效果

ShadowMap
PRT
SSRT
KullaConty BRDF OFF
KullaConty BRDF ON

作业1 ShadowMap

ShadowMap的基本思路为：将相机放在光源处看向场景，记录深度信息，一般方向光对应正交投影，而聚光灯或者点光源对应透视投影，以及此时的View和Proj矩阵，记为LVP。有了深度信息后，正常渲染时，在片段着色器中，可以获取片段的位置变量pos，使用LVP对其进行变换。

$clip\_pos = LVP * vec4(pos,1.0)\\ ndc = clip\_pos / clop\_pos.w \\ uv = ndc.xy * 0.5 + 0.5 \\ z = ndc.z * 0.5 + 0.5$

注意OpenGL下NDC的xyz范围都是-1~1
如果uv都是0到1之间，说明在相机视角下的这个片段位置，也可以被灯光看见，可以根据uv从之前的深度信息图中查询到该片段在灯光视角下的深度值lz，将其与计算得到的z值做比较，如果z > lz，说明该片段被遮挡住了。

但是这样子会造成一个问题，因为记录深度信息的纹理分辨率是有限的，因此图中的一个像素，代表一个投影区域中心的深度值，任何在该区域内的其它位置，在第二步渲染查询深度值时，都将是一样的，这样子就会导致物体表面出现自遮挡现象，像是黑色的摩尔纹，放一下课件里的图。
self shadow

解决这个问题，一般思路是满足 $z - eps > lz$ 才被认为遮挡住，eps可以根据片段法向量和实现的夹角来缩放，
实际中，可以将原来的pos沿着normal移动eps的距离再进行变换。

解决自遮挡的问题后，还是有一定的问题，比如受限于深度图的分辨率，阴影会出现锯齿，一个简单而自然的想法就是对阴影进行过滤、模糊，即计算一个片段的阴影系数时，不再是0或者1，而是0~1之间，这就是PCF的思路。最简单的做法，就是取固定的滤波核，比如5x5，滤波核取得越大，阴影越模糊，反之阴影越尖锐。滤波范围理论上取圆盘区域可能更加符合实际，因此也可以划定一个半径，然后均匀采样深度信息图上圆盘内的其它深度值，与z作比较，求出阴影系数的平均值。作业还提供了圆盘泊松采样的方法，类似于一种由中心向外的螺旋线。

然而现实中的阴影其实是hard和soft兼具的，被遮挡的点里遮挡物越近，阴影越尖锐，反之则越模糊。放一张闫老师特别称赞的一张图。
percentage shadow

可以观察到，笔尖的阴影边界特别清晰，而远一点的阴影边界则自带模糊，我觉得这张图对形成原因解释地比较清晰。
模糊的阴影形成是因为部分光线被遮挡。因此PCSS的基本思路就是找出适当的滤波区域大小，再进行PCF过滤。
sun earth shadow

PCSS的原理图如下。
pcss
阴影的滤波区域大小与三个变量有关，灯光区域大小、遮挡物的深度和阴影处的深度。其中唯一不知道的就是遮挡物的深度，但还是要用平均深度替代，因此为了求这个平均深度，我们需要PCF求出平均遮挡深度，再根据这个平均遮挡深度求出滤波半径，再进行一次PCF，相当于两次的PCF，因此PCSS一般比PCF更消耗性能，但是获得的效果也会更好。另外，求平均遮挡深度时的PCF的滤波范围大小怎么确定呢？可以是固定的，不过更好的做法是根据光源的大小和遮挡处的深度。
search radius
将深度图放在灯光的近平面（假设不是方向光源，因为方向光源理论上不会产生软阴影），根据相似三角形关系可以求出滤波核在深度图上的对应大小。然后套用PCF求出遮挡的平均深度即可。实际应用时，光源的大小和近平面的大小需求调整到合适的范围。

作业2 PRT

这部分强数学相关…很无奈自己没有学过离散数学和信号之类的知识，靠着微弱的微积分基础，勉强看懂了吧orz…
这一切都要从傅里叶展开说起，任意一个y(x)函数都可以由一系列的正弦和余弦函数组成。函数的级数展开也是类似的，不过傅里叶的展开有其独特的性质，比如基函数之间是相互正交的，不同基之间的乘积为0，这一点十分重要，prt里就是用到这一点性质。公式可以写作如下：

$f(x) = \sum_{i}c_i \cdot B_i(x)$

基函数除了正弦和余弦函数其实可以有其它的选择，而且这个也可以拓展到二维：

$f(\omega) = \sum_{i}c_i \cdot B_i(\omega)$

比如这里使用到的球弦基函数可视化出来是这样子的：

将渲染方程分为两部分，光照和光照传输：
render equation
对于环境光，将其看作一个二维函数 $L_i=f(\theta_i,\phi_i)$ ，因此可以投影到基函数中：

$L(\omega_i) \approx \sum^N_{j = 0} l_j \cdot B_j(\omega_i)\\ l_j = L(\omega_i) \cdot B_j(\omega_j)$

将上面的结果带入原来的渲染方程可以化简得到

$L(o) \approx \rho \sum^N_{j=0} l_j\int_{\Omega}B_j(i)\cdot V(i)\cdot (\omega_n \cdot \omega_i) d\omega_i \\ \approx \rho \sum^N_{j=0}l_i \cdot T_i$

其中 $T_i$ 是light transport投影到基函数的系数。在这里把brdf视为constant变量，所以可以把 $\rho$ 提到积分外面，
这样子预计算的部分就和具体的 $\rho$ 无关了，可以任意调节。不过如果有albedo map贴图的话，就需要在把light transport部分投影的时候也把 $\rho$ 带入，但是最后预计算的结果是基于这个贴图的，渲染时无法改变反射率。
另外light transport也可以考虑多次bounce之后的入射光，也可以考虑阴影的遮挡，和path tracing的流程是差不多的。简单来说，就是把light和light transport投影到球弦基函数上，算出投影系数，然后再实际渲染时相乘一下即可。而投影到某一坐标系，就是与坐标系的所有基函数做点乘。
至于环境光的旋转部分，主要用到了球弦基函数旋转后仍然可以用原来的基函数线性表示，那么旋转环境光，等价于在旋转之后的基函数上重新投影，而旋转后的基函数可以由原来的基函数线性表示，那么就可以求出两者之间的变换关系，具体的操作方法这里不多介绍了，因为我也是看了一个大佬的博客明白的，链接在这。

作业3 SSRT

屏幕空间反射（光线追踪），建立在延迟渲染的基础上，也就是说进行光线追踪的信息来源于GBuffer，并不是整个场景的信息，所以它并不是那么准确，比如对于一个处于阴影处的物体，它的光照只能来源于间接光，但是如果屏幕中没有其它被直接光照射的物体，那么它就没有光照来源，只能为黑色，或者靠着ambient环境光证明自己的存在。但是SSR当然尤其好处，就是可以做到光栅化无法解决的反射求交问题。在光栅化流水线中，如果想要对一根射线与场景三角形求交是十分困难的，如果是OpenGL只能是自己动手实现一下BVH？见过大佬实现过orz…而对于SSR，可以很简单地使用ray marching求出反射光线的交点，当然ray marching这一方法本身是比较消耗性能的，不过也有优化加速的方法，因此ssr在特定场景（比如限制相机视角）效果还是十分好的。ssr的原理和流程图如下，其实还是十分简单的。
ssr

ssr建立在延迟渲染的基础上，那么就需要GBuffer，如上图，需要pos、normal、albedo、view depth，前三者是vec3类型，最后一个float类型，如果简单的话，我们需要四张二维纹理，其中pos应该需要rgb32f，normal和albedo可以使用rgb8，而view depth要使用r32f。综合下来，虽然内存使用是不大的，但是每次从GBuffer中读取这些信息，是有代价的，现在需要从四张纹理中读取，还是可以进行压缩的。其中normal和albedo可以进行压缩，从vec3变为vec2，因此所有的GBuffer最终只需要两张rgba32f即可。

法向量的压缩算法可以参考这个

vec2 octWrap(vec2 v)
{
    return (1.0 - abs(v.yx)) * (all(greaterThanEqual(v.xy, vec2(0.0))) ? vec2(1.0) : vec2(-1.0));
}

vec2 encode(vec3 n)
{
    n /= (abs(n.x) + abs(n.y) + abs(n.z));
    n.xy = n.z >= 0.0 ? n.xy : octWrap(n.xy);
    return n.xy;
}
vec3 decode(vec2 f)
{
    vec3 n = vec3(f.x, f.y, 1.0 - abs(f.x) - abs(f.y));
    float t = clamp(-n.z, 0, 1);
    n.xy += all(greaterThanEqual(n.xy, vec2(0.0))) ? vec2(-t) : vec2(t);
    return normalize(n);
}

对于颜色，直接将r和b通道压缩为float16当作一个float32进行存储即可，代码如下。

//encode
vec3 color  = albedo / 3.14159265;
float color1 = uintBitsToFloat(packHalf2x16(color.rb));
float color2 = color.g;
//decode
vec2 color1 = unpackHalf2x16(floatBitsToUint(p1.x));
vec3 albedo = vec3(color1.x, p1.g, color1.y);

直接光的光照计算直接采用公式，计算完后存储到一张纹理中，作为计算间接光时的输入。

$L_{direct} = BRDF * L_{light} * Vis * \cos\theta$

间接光的计算采用蒙特卡洛采用，使用z-weight半球采样，对应diffuse材质，公式为：

$L_{indirect} = \frac{1}{N}\sum^{N}_{i = 1}\frac{BRDF * L_i * (\omega_{i}\cdot \omega_{n})}{pdf(\omega_{i})}\\ =\frac{\pi}{N}\sum^{N}_{i=1}BRDF * L_{i direct}$

这里去掉了Vis项，因为通过ray marching找到的都是可见的，否则直接跳过这次采样。
对于线性的ray marching，也就是说光线在view space中按照固定步长前进，每次都要判断当前的位置是不是小于通过纹理查得的z值，如果小于则说明光线与该片段相交，即该被直接光照射的片段，可以被计算间接光这一片段看见，作为间接光的入射光源。这里其实会存在一个问题，那就是当判断当前光线所处的片段比z值小时，还有一个情况就是光线当前位置处于直接光照射的阴影处，这样子就会造成虚假的间接光入射源，特别是断层附近，低处会突然比较亮，如下图。
depth threshold off
depth threshold on
右图是进行了DepthThreshold判断后的结果，即光线位置的深度不可以小于z值太多，因此光线每一步前进的步长Step和DepthThreshold这两个变量的值是相互关系的。

线性ray marching还是很简单的，在间接光着色时，根据片段的当前位置，以及半球采样的方向发出一条光线，每一步前进固定大小的距离，每次都判断当前位置对应的view depth是不是比纹理记录的要小，以及如果当前光线已经跑出view空间范围了，如果满足所有要求，则返回true和直接光信息，累加起来就好。至于效率，当然是不太好，因为每次前进的距离都是固定的，在找到满足条件的位置前，有很多不必要的循环计算和查表，因此可以使用mipmap进行加速，建立view depth的mipmap，使用最小值策略。在光线前进时，如果当前位置深度大于某一lod下对应的深度，那么可以直接跳过这一lod这一像素对应的范围。这里就有一个重要的点，在view space里的以均匀步长前进是一种浪费，因为最后的深度纹理还要经过project变换才能得到，也就是说，view space里某个方向上两个点之间的连线经过project后，它的长度不是线性变化的，在投射投影下。另外它所代表的像素范围可能很小，即我们在view space里线性前进的步长设置应该根据view depth图像的分辨率来调整。不过由于投射摄影下，即便是同一方向上固定距离的两个点，随着它们与相机距离的变化，它们经过投影变化后在图像上代表的范围也会发生变化。比如根据近大远小，在远处前进距离x，对应查询的纹理范围更小，假设占据1个像素，而在近处前进距离x，对应查询的纹理范围就更大，可能为6个像素，因此为了有效利用，显而易见应该在远处以更大的步长前进，这个比例理论上可以找到的，根据z值和view space下方向的方向，不过很麻烦，简单的做法是在ncd坐标系下进行。在ndc坐标系下，计算view space下光线在ndc坐标系中的方向，以在ndc坐标系下固定的步长前进，对应到view space就是一种自适应采样步长的方法，这里都没有涉及到mipmap。

ray mipmap test

层次光线投射的原理是，如果光线没有与较大lod下的node相交，那么也就肯定不会与其子节点相交。伪代码为如下。

level = 0;
level_advance_dist[max_level];
hit = 4;
while(level > -1)
    step through current cell
    update level advance distance
    if(above Z plane ) 
        ++level;
        increase step
        continue;
    if(level == 0)
        test if satisfy depth threshold
    level--
    decrease step
    update level advance distance

这里的hit是用来提高效率的。如果光线与一个非叶节点相交，显然易见并不是说它一定会与某个子节点相交，因此我们需要递归测试其子节点，比如该节点有两个子节点，node0和node1，如果不适用hit，当测试node0不相交时，会马上又进行level++，即便它原本可能会与node1相交，因此hit就是保证当光线与某个level的node相交后，可以保证测试完其level-1的子节点。

感觉自己的算法还是有那么些问题的，算法没有那么完整、准确或者鲁棒，不过最终似乎对于这个场景work了，虽然是在漫反射的材质上，还没有测试过镜面反射的材质，原方法叫做 stochastic screen-space reflections，还结合了tile based的方法，性能和效果都很好。

此外对于间接光的采样，采样数可以设置地少一点，比如8个，这样子每一帧的耗时就会少很多，但是可以通过时间上的累加最终达到收敛的效果。

小结

SSR的优点

对于glossy和specular的材质可以效果好并且速度快

SSR的缺点

对于diffuse的材质，需要较多的采样点才能收敛到比较好的效果
丢失了屏幕外的信息
只能处理二次反射

作业4 KullaConty BRDF

这个首先建立在微表面模型的基础上，即BRDF = D * F * G，由法线分布、菲涅尔项和几何遮挡项组成，如果不太了解，可以先看看这篇文章。

在IBL中，把BRDF分成了specular和diffuse项，specular项使用严格的DFG计算，而diffuse项则是根据所谓的能量守恒原理，即根据菲涅尔公式求出F后，kSpecular = F，那么 kDiffuse = 1 - F 。这样子的做法似乎看起来是正确的，显示的效果也不错，至少在IBL的结果中是很nice的。然而… 闫老师指名批评了一下，这个做法是一种undesirable hack，完全的错误的。
consider diffuse in microface model is totally wrong
Diffuse项是什么呢？之前认为是光线投射到物体内部再次弹射出来的部分，为了简单考虑认为是均匀分布的。感觉这个应该是有那么点物理可靠的，对于电介质而言，光线进入内部后还是有很大概率从表面射出的，而金属会把投射进来的光线能量几乎全部吸收，因此对于电介质而言，考虑简单的diffuse应该是挺可靠的。但是对于金属而言，应该是没有diffuse项的，那么随着金属表面的粗糙度增加，金属表面的能量损失就会越严重，表现出来越暗，在这个时候加一个diffuse项那就是物理上不正确的，感觉闫老师指的应该是这种情况，这个时候表面粗糙度较大，如果只考虑单次的bounce，就会有很多的能量损失，看上去很暗，但实际生活中看到的金属粗糙度虽然大，也不会很暗，因为还有多次bounce的情况，即光线被物体表面第一次反射后，还会与其它微表面继续相交，经过多次这个过程后再离开物体表面最后被观察到。
energy loss

有一个测试叫白炉测试，用于测试能量的损失，对于BRDF = DFG而言，需要假设F项始终为1，以及环境光从每个方向发射的辐射度均为1，那么根据渲染方程的定义，从每一个角度观察物体应该都是白色的。

可以看到，如果只考虑单次的boundce，随着粗糙度的增加，物体反射损失的能量越多。而KullaConty-BRDF的基本思想就是，将这一部分损失的能量以BRDF的形式弥补。假设固定观察方向o，那么看到的能量应该是

$E(\mu_o) = \int^{2\pi}_{0}\int^1_{0}f(\mu_o,\mu_i,\phi)\mu_i d\mu_i d\phi \\ \mu = \sin\theta$

因此损失的能量是 $1 - E(\mu_o)$ ，考虑BRDF的对称性以及不同入射角和出射角对应的损失能量不同，那么弥补能量的BRDF的形式应当是
$c(1-E(\mu_i))(1-E(\mu_o))$

因此所有的推导结果，直接放图了，不手敲公式了。
推导过程

因此我们只需要预计算 $E(\mu_o)$ 和 $E_{avg}$ 两张表就好，在实际渲染中根据上面的公式就可以计算出 $f_{ms}$ 。
在这里解释一下 $E_{avg}$

$E_{avg} = \frac{\int^1_0E(\mu)\mu d\mu}{\int^1_0 \mu d\mu} = 2\int^1_0E(\mu)\mu d\mu$

其实是求了 $E(\mu)$ 的平均期望，也就是说，任意出射方向的平均能量。
对于 $E(\mu_o)$ ，其与roughness和 $\mu$ 相关，因此是一张二维的查找表，而对于 $E_{avg}$ ，则只和roughness有关，因此是一张一维的查找表。

先前都是假设F=1，但实际上BRDF是有颜色，即会有能量的吸收和损失。这里定义了一个平均的菲涅尔项。

$F_{avg} = \frac{\int^1_0F(\mu)\mu d\mu}{\int^1_0\mu d\mu} = 2\int^1_0F(\mu)\mu d\mu$

这里提到的E其实并不是真正的能量，因为积分的时候并没有带入L项，其实是一种比例。
我们以及得到了fms，代表不考虑颜色的情况下，brdf应该要加回的能量。那么在考虑了颜色吸收后，
假设fms是反射的光线能量，其值为单位1，那么这个光线的能量可以直接打到人眼的比例为 $E_{avg}E_{avg}$ ，即需要同时考虑能量本身的损失和因为颜色的损失，
那么再经过一次反射后，能量比例为：

$F_{avg}(1-E_{avg})\cdot(F_{avg}E_{avg})$

再经过k次后：

$F^k_{avg}(1-E_{avg})^k\cdot F_{avg}E_{avg}$

然后进行等比数列求和，最终因为颜色损失后剩余的能量占比为：

$\frac{F_{avg}E_{avg}}{1 - F_{avg}(1-E_{avg})}$

$F_{avg}$ 有快速计算的公式，对于电介质来说：

$F_{avg}(\eta) \approx \frac{\eta - 1}{4.08567 + 1.00071\eta}, 1 < \eta < 400\\ F_{avg}(\eta) \approx 0.997118 + 0.1014\eta - 0.965241\eta^2 - 0.130607\eta^3, 0 < \eta < 1$

对于导体来说：

$F_{avg}(r,g) \approx 0.087237 + 0.0230685g - 0.0864902g^2 + 0.0774594g^3\\ +0.782654r - 0.136432r^2 + 0.278708r^3\\ +0.19744gr + 0.0360605g^2r - 0.2586gr^2\\ r = reflectance, g= edgetint$

关于KullaContyBRDF的更多内容，可以参考这里

那么IBL也可以校正一下，写了一个小demo，感兴趣的可以看看。