Java计算菲涅耳积分 菲涅尔积分的结论

Lec5~7

• 1 环境光贴图下的着色计算（不考虑阴影）

• 1.1 Image-Based Lighting(IBL)
• 1.2 The Split Sum

• 2 环境光贴图下的着色计算（考虑阴影）

• 2.1 球谐函数（Spherical Harmonics）
• 2.2 预计算辐射传递 （Precomputed Randiance Transfer，PRT）

• 2.2.1 diffuse物体计算案例
• 2.2.2 glossy物体计算案例

1 环境光贴图下的着色计算（不考虑阴影）

环境光贴图： 一张图，图上记录着场景中任意方向的来自无限远的光照是多少。场景中任何位置的物体接收的光照都是一样的，因为应用同一个贴图。

场景中已经应用了一张环境光贴图，如何计算物体的着色？（不考虑可见性项）

1.1 Image-Based Lighting(IBL)

IBL:基于图像的光照计算，即光照信息是从图像中获取，在渲染方程中的体现为入射光Li。并且我们认为，着色点不存在遮挡问题即不考虑可见项。

一般方法：蒙特卡洛（不适用），数值解、需要大量采样。缺点：特别慢，实时渲染中不可用。只要涉及采样，都会特别慢。

如何解上述方程并且避免采样？？？

观察

• 如果BRDF是glossy —— 球面上积分域很小，brdf返回值最值差距比较大 —— 不太smooth，但也还能接受
• 如果BRDF是diffuse —— 积分域很大，为整个球面，brdf返回值为常数 —— smooth

马上考虑到把 乘积的积分 拆成 积分的乘积

---

1.2 The Split Sum

就是在一定条件下，把积分拆开计算，得到近似的结果

• 此方法在g(x)的积分域比较小 or 值比较smooth（最值相差不大）的情况下比较准确。注意：$Ω_G$为g(x)的积分限
• 因此我们可以直接把光照项Li拆出来，把BRDF项当做g(x)留在里面，因为很满足积分域较小（glossy）、smooth（diffuse）再次注意： 左边是Li在BRDF范围的积分，并且归一化

---

为了得到着色点的环境光计算结果，这个公式分为了两部分

Stage 1：左半部分计算

• 如果是glossy的BRDF，这个公式的意义就是在对环境光贴图上位于镜面反射方向，周围的所有Li进行求和，然后再求均值
• 如果是diffuse，就不用算镜面反射方向了，直接对整张环境光贴图做均值
• 因此，整个过程就是：采样BRDF范围的贴图上的光照强度，然后求均值

采样太慢了，不想采样

这里就引出一个思想：预计算(Pre filtering)环境光照

• 根据原始环境光贴图，预先计算生成一系列 应用不同大小的卷积核后的环境光照图像
• 查询时，如果BRDF的范围处于两个卷积核大小的层的中间 则进行三线性插值得出（类似于MIPMAP分层）
• 然后，就能很方便的 查询 上面公式的左半部分的值，不用采样，更不用计算
• 以glossy的brdf为例

• 本来应该从着色点向镜面反射方向的附近发射N根光线，去采样环境光贴图然后均值，作为着色点的颜色
• 现在只需要发射一根光线，就能查询到已经预计算好的均值

---

Stage 2：右半部分的计算

现在虽然已经有更好的解决方案，但是接下来的方法很经典 必须掌握

如何避免采样 —— 同stage1的思想，预计算

如果此处的BRDF采用微表面模型（Microfacet BRDF）来描述，如何用预计算来减少工作量呢？

• 预计算所有参数变量的所有可能的组合（roughness(NDF), color(Fresnel term)…）
• 但是这样的预计算参数空间太大，每一项基本都有1~2个变量，导致存储空间太大

Fresnel term —— 颜色

• 菲涅尔项可以用Schlick’s approximation来表示
• 两个变量：入射角度θ、基础反射率Ro（入射角度为0的）

NDF term —— 表面法线分布

• 有很多不同的分布函数，这里以beckmann distribution为例
• 两个变量：粗糙度α、半程向量与法线的夹角

NDF中$θ_h$可以用入射角θ通过一定的方式来表示，从而降低至三维，对预计算来说，三维依然太高，如何能将参数的纬度降低至2维？

• 把积分中的BRDF项$f_r$除以F再乘以F，乘以F的这一部分菲涅尔项用公式表达$R_0+(1-R_0)(1-cosθ)^5$，整理后得到下面这个式子。注意$\Large \frac{fr}{F}$可以看做约分，fr只剩下NDF项和几何项(忽略)
• 因此stage 2又被拆成2个积分，这个操作消除了对于基础反射率的依赖，即把Ro给拿到积分外面去了
• 此时对于这两个积分，参数空间仅剩2维：入射角$\color{red}θi$、NDF中的粗糙度$\color{red}α$
• 预计算拆分后的这两个积分的计算结果的所有组合，由于其只包含2个变量（入射角、吐槽度），所以只需要两张texture(分别对应两个积分的预计算结果)，甚至省一点，一张texture的两个通道存放预计算结果。渲染时根本不用计算，直接查询

以上便是借助Split Sum方法进行预计算，从而完成IBL下物体的着色计算

UE4中的Split Sum应用效果

• 通过把渲染方程拆分成两个积分（stage1、stage2），再把stage2拆分成两个新的积分
• 通过预计算的技术即可 避免采样，效果非常好，很接近离线渲染的效果了

---

2 环境光贴图下的着色计算（考虑阴影）

不考虑阴影：使用Split Sum方法进行IBL下的着色计算，确实能得到很好的效果，并且与计算好后，渲染时能够极快的查询到结果

考虑阴影，则计算难度就不同了，几种不同的视角看待环境光贴图中的光源

• Many-light
  把环境光贴图上比较白的部分用一个点光源代替，而多少个光源就产生多少shadow map ，代价高不好用。工业界通常处理办法就是给环境贴图上最亮的那几个地方放光源，基本够用
• 采样
  解渲染方程，以一定的频率采样位于着色点的各个方向的光照，最难的是每个着色点每次采样的Visibility项也是一个变量，V项不能简单的对环境光贴图进行采样，因为中间可能有各种阻挡问题

环境光贴图下的阴影计算相关研究：imperfect shadow maps 、light cuts（离线）、实时光线追踪、PRT

这里主要学习PRT，而学习PRT之前必须了解什么是 球谐函数

2.1 球谐函数（Spherical Harmonics）

傅里叶变换、卷积、滤波等知识回顾看这里

前置基本概念

• 把任何两个函数乘积的积分看做filtering操作
  $\large \int_{Ω}{f(x)g(x)}dx$

• 低频 —— 代表函数比较smooth，函数的最值相差小，变化缓慢
• 积分运算的结果的频率取决于被积函数的 频率最低的那个函数，相当于给另一个函数应用低通滤波器

• 基函数（Basis Function）： 一系列用来表示其他任意函数的函数$\mathbf{B_i(x)}$。比如傅里叶展开、泰勒多项式展开
  $\large f(x) = \sum_{i}c_i·B_i(x)$

---

球谐函数SH：一系列定义在球面上的2D基函数$\mathbf{B_i(\omega)}$，可以理解为关于方向的函数，三维空间的方向可以用两个角度描述（θ、φ）。

球谐函数的可视化

• 蓝色为正，黄色为负，离中心越远绝对值越大，值变化越快，频率越高
• $l$为SH函数的阶，同阶的基函数频率相同；阶数越高，函数频率越高，函数个数越多
• 不同阶的基函数个数 = $2l+1$
• 每阶的基函数编号m从$-l$ 到 $l$
• 前n阶基函数个数：$n^2$

对于任何一个原始的函数，都能用一堆基函数和系数来描述

• 比如用4阶，用16个基函数存储，则只能还原出原函数的较为低频的信息，高频信息丢失，相当于应用了一个低通滤波
• 基函数用得越多，越能记录原函数的高频的信息

环境光贴图就是一个二维的函数，光照值为立体角的函数，因此环境光贴图可以用一系列基函数表示！可以把这个函数投影到任意阶数的HS的基函数上，用前n阶的球谐函数来恢复这个二维函数的低频信息。越高阶占用的存储空间越多，能恢复出来的频率越高

投影（Projection）：已知任何一个2D函数$f(\omega)$都能用一系列基函数的和来表示，而计算每个基函数的系数$\large{c_i}$的过程就是投影。求解公式为
$\large \color{red} c_i=\int_{Ω}{f(ω)B_i(ω)dω}$

在明白球谐函数后，得到系数的过程叫投影，思考一下空间中的一个点如何表示的？

• 空间中有一组相互垂直的xyz轴，任意一个向量的坐标表示是该向量在三个坐标基上的 投影
• 投影是 点乘，，积分本质是连续求和，如果用离散的思维来考虑，每个$ω$对应一个$f(ω)$和$B_i(ω)$，遍历所有ω，得到的其实是两个向量，而两个向量的乘积就是 点乘；因此可以很容易的看出上面的基函数系数计算本质就是 点乘再连续求和
• 再看一个正交坐标系中，xyz坐标基相互正交，任意基投影到另一个基上值都是0。其实SH表示的基函数投影到另一个基函数上也是0，而基函数自己对自己的投影（点乘）结果为1

---

应用球谐函数：不考虑阴影的情况下，环境光照下diffuse物体的着色计算

Prefiltering + single query == no filtering + multiple queries

• 左图：不对环境光贴图做Prefiltering，直接按照反射方向去贴图上做1次查询，结果为镜面小球
• 如果先对环境光贴图Prefiltering，再沿着某个方向进行1次查询，结果如右图所示，像一个diffuse小球，等价于不Prefiltering，往任何一个着色点的法线方向的球面多次查询，取平均

再看Render Equation。被积函数有两个，逐点相乘后积分起来，这就是 点乘，$L_i$是环境光贴图提供，而环境光贴图相当于一个 二维的球面函数，初步考虑能用SH函数来表示，问题是存多少阶合适呢？

• 环境光$L_i$项可以有低频、高频信息，BRDF项是定义在整个半球上的很光滑的函数（diffuse的brdf是个常量，处处相等当然光滑），此时它相当于一个 低通滤波器。所以Li 不管高频有多少，只要乘上漫反射BRDF项后，根本就没有高频信息留下来！
• 有实验表明，用球谐函数来描述漫反射BRDF只需要3阶就完全足够。因此根本没必要用高阶SH来存储环境光贴图

因此，利用这个信息，既然任何环境光贴图在漫反射物体上的着色结果都只有低频信息，那 直接用低阶SH来表示环境光贴图不就好了

• 一阶SH存储环境光的效果（左图为正常计算环境光映射的图像，右图为用一个基函数表示的SH计算出来的图像）
• 取前2阶SH基函数来描述环境光
• 取前3阶SH基函数，几乎一模一样，误差为1%
  
  结论：任何光照，只要是照亮漫反射物体，用3阶SH就足够

SH的美丽属性（下面diffuse案例结尾中还有补充）

• 正交性 任何两个基函数互相垂直，自己投影自己为1，不同的基函数相互投影为0



• 投影计算简单，任意一个函数可以投影到任何基函数上，做product integral即可
• SH预计算好的光照可以任意旋转，重新计算基函数的系数代价很低。如果旋转光照，等于旋转SH的所有基函数，而旋转任何一个基函数，都可以用同阶的基函数线性表示，所以只要查表就能得到旋转后的光照所对应的SH基函数的线性组合
• 少量几个基函数，就能很好的还原球形函数的低频信息。下图N代表基函数个数（高频信息想要很好地还原需要很大代价）



• 只要确定阶数，则基函数就确定了，不确定的只有系数，可通过投影得到。（个人猜测）
• 球谐函数跟平时常用的XYZ坐标系一样，每个基函数都是互相垂直的坐标基，因此SH基函数也构成了一个SH space（个人猜测）
• 足够多的基函数能表示任何函数
• 保留某个频率以下的内容，只需要取前面的一部分基函数就能重建原函数，得到原函数的低频信息
• 如果两个函数都是SH的基函数表示，那么他们其实就是点乘

---

2.2 预计算辐射传递 （Precomputed Randiance Transfer，PRT）

实时渲染中，环境光照下考虑可见性(阴影)的Render Equation可简写成下面这样

• 以右下角框出来的着色点为例，它的环境光照、可见性、BRDF分别可以球面函数（6张贴图）表示。（由于相机位置固定（o为常数），BRDF就只剩两个参数了（i入射角拆为天顶角，方位角），所以也能看做是球面函数）
• 就一张cube map而言，假如一个面分辨率64x64，六个面就是6x64x64，三张cube map需要3x6x64x64，这还只是一个着色点，所以预计算的代价相当之大。可见，没有Spherical Hamonics的情况下，阴影的预计算环境光着色计算是多么困难的事情

SIGGRAPH 2002中的这篇论文提出了PRT技术

PRT基本思想：假设整个场景 只有光照可以发生变化（光照可以旋转、可以更换），该着色点的其他部分如BRDF、可见性、相机位置都不变。总之被积函数分为两部分，lighting部分以及与lighting无关的light transport部分，并且这两部分都能用球谐函数表示。

Light Transport部分就像其命名一样，描述光如何从入射$L(i)$变化到出射$L(o)$

• Precomputation阶段

• lighting：用SH来表示 $\mathbf{ L(i)≈\sum l_iB_i(i)}$ 。旋转、替换环境光贴图后，需要重新预计算SH组合
• light transport：这一部分对于每个着色点来说，可以看做是着色点本身 不随光源变化的 性质

• $V(i)$是关于$i$的球面函数，因为着色点是否被某个物体阻挡，只与着色点自身所在位置有关
• Brdf项，也是着色点本身的材质属性，并且相机位置还是固定的，是2D球面函数
• 后面的cos项也是球面函数

• Light Transport整个部分，合并在一起用一套SH来表示（只与入射光的两个角度有关）

• Runtime阶段：因为两部分都是预计算的，直接查询后做乘法

• 如果diffuse材质，做 向量点乘
  积分就是求和，被积函数的左右两部分对于每个不同入射角度i都把查询到的值相乘，最终把每次相乘结果加起来，线性代数表达方式：$\mathbf {(l_1,l_2,...,l_n)^T ·(lt_1,lt_2,...,lt_n)^T}$。$l_1和lt_1$相当于入射角度为1时查询到的预计算的lighting 和 light transport的值
• 如果glossy材质，做 矩阵-向量乘法
  没啥区别，就是Lighting部分每个入射角度i查出来是一个矩阵

2.2.1 diffuse物体计算案例

• 为了与SH基函数下标相区分，文字描述中我把积分变量$i$用$ω_i$表示
• 如下图所示，BRDF项$ρ$是常数，可以提前。
• lighting部分可以用SH表示，然后可以把系数部分的求和提出来，因为系数$l_i$与积分变量$ω_i$无关

实时渲染中积分符号与求和符号可以认为是等同的（连续和离散的区别嘛）

• 然后再看被积函数剩下基函数部分$B_i(ω_i)$以及light transport部分，而light transport可以看做一个函数球面函数$f(ω_i)$，即积分变成了$\displaystyle\int_ΩB_i(ω_i)·f(ω_i) dω_i$，可以理解为light transport函数在基函数$B_i(ω_i)$上的 投影，投影就是在算基函数$B_i$的系数（这里注意，外面的$l_i$是Lighting部分用环境光函数投影出来的系数，假设是前3阶，则应该是9个系数求和）。积分内预计算结果就是用不同于$l_i$的另一组 系数
• 两部分预计算做好后，着色计算只需要两套系数之间做 向量点乘

下一节课中有提到另一个思路：

• 把lighting和transport部分都用SH替代

• 然后渲染方程可以变成这样两部分

• 积分外：两套系数的双重求和，乍一看是一个矩阵的所有元素的和（每个元素都是$c_pc_q$的乘积）
• 积分内：俩基函数的互乘。p==q时，结果为1；p≠q时，结果为0。
  这里也间接印证了我的前面的在SH美丽的属性部分的猜测吧？两部分的原函数不同，但是都转换成SH表示的话，下标相同的基函数是同一个基函数，所以点乘为0。百度了一下确实存在基函数表这种东西

• 因此：实际上仅当p==q时，结果不为0，也就是说除了主对角线，其他部分都是0，因此依然是两个向量/两套系数做 点乘

PRT Diffuse总结（SH函数性质的部分补充）：

• 任何球面函数，用SH来表示后，其实就是一个 系数向量



• 系数向量计算方式为投影，即把原函数投影到SH空间中，即product integral
  $l_i=\int_ΩL(\omega)·B_i(\omega)d\omega$
• 对原函数的还原只需要用每个系数与对应的基函数相乘后相加即可，也就是个点乘，一个系数向量，一个基函数向量
  $\displaystyle L(w)\approx \sum l_iB_i(w)$

效果

• 环境光照是直接光照，即光源打到物体进入摄像机
• 可见性项预计算是每个着色点往球面上发射光线 判断可见性 算出来的，因为是渲染前的预计算，所以不考虑花费多少时间

diffuse部分结束

---

2.2.2 glossy物体计算案例

与diffuse不同的地方在于，glossy的brdf不是一个常数，是一个4D函数(io分别两维)，它的o与视点是有关的，不同视角/入射角度的反射率都不同

• 把$L(i)$用SH表示很简单，跟diffuse案例中是一样的
• 而Light transport部分就不一样了，因为diffuse的brdf是常数，不用考虑，而glossy部分的brdf既是i的函数，也是 o 的函数，这就很麻烦了。
• 如果仅仅是i的函数也还好，因为transport部分依然为2D球面函数，可以用球谐函数表示。关键还是o的函数，即每更换一个相机位置 o，就要重新投影得到一组系数，因此最终着色点的计算是向量·矩阵。
• 注意下面公式省略了中间一些过程，即$L(i)$用SH表示，系数求和提到积分外，积分内的基函数用于transport部分的投影计算，最终得到下面的$\displaystyle L(o)\approx\sum_il_iT_i(o)$，必须看懂这个看起来简单的式子（$\sum T_i(o)$ 相当于二维矩阵，$T_{ij}$，即不同相机角度o对应不同的$T_i$向量，L(o)最后得到的是一个向量，确定o才能确定具体着色结果）



• $T_i(o)$又是一个球面函数，又能投影到SH 空间中去，但是这个过程略有不同，
  因为此时的球面函数$T_i(o)$是位于$\sum_i$的作用域下的
  $\displaystyle L(o) \approx\sum_il_i(\sum_jt_j·B_j(o))=\sum(\sum l_it_{ij})B_j(o))$

最终着色点结果就是计算 向量·矩阵

代价：任何一个着色点 or 顶点 都要存一个矩阵（i、o的组合）效果图

---

总结

如果选择SH基函数为16个，则着色计算时

• Diffuse Rendering：vector(16)·vector(16)
• Glossy Rendering：vector · matrix(16x16)

Light Transport部分预计算过程的理解方式

• 整个Transport部分看做一个球面函数，向第i个基函数上投影，结果记作$T_i$
• 也能把每个基函数都看做一种光照，brdf是常数不用管，而预计算的工作就是把N种基函数光照下的物体每个物体的着色结果都计算出来
• 与计算好后，Run-time 渲染计算，只需要做向量点乘，这在shader中是很好写的
  
  不同BRDF的PRT渲染结果

接下来的学术界研究内容

• 新的基函数
• 更多的预计算内容，比如把Light transport部分拆开，然后点乘从两个变成多个
• 动态场景
• 动态材质
• 其他效果，如：半透明材质、毛发

---

还有很多基函数种类，其中具体讲了Wavelet，效果很好，不仅能还原低频信息，还能还原高频，但是不能旋转光源，具体内容这里就不写了。