光线追踪及相关数学原理与实践(二)

标签:
游戏渲染光线追踪路径追踪 |
分类: 游戏 |
接上文,若没看过请先阅读之。
二.路径追踪(Path Tracing)
首先要说明的是路径追踪并不是一门全新的技术,可以认为是光线追踪的一种。实际上路径追踪和光线追踪的概念是有点模糊的,路径追踪只是在DRT的基础上的做了些扩展,现如今还经常把两者混为一谈。路径追踪的特点是光线每次反弹后只产生一条光线,所以更加关心光线反弹路径的选择质量,显然这就需要从概率论中寻找相关的知识帮助。此特点也解决了上文提到的光线数量的指数爆炸问题。但是路径追踪诞生的初衷可并不是为了解决这个问题,而是跟当时的图形学历史有关。
1986年Kajiya提出了著名的渲染方程,这在基于物理的渲染(Physicallly-Based Rendering,简称PBR)的发展史上是一个具有里程碑意义的事件,人们普遍认为想要渲染的真实就得符合渲染方程,一直影响至今。但是此方程无解,或者说人类无法用已有的知识直接求解,因为渲染方程本质是个积分,求一个复杂被积函数的积分本来就是一大数学难题,更何况它还是个递归形式的积分。所以同年,Kajiya自己又提出了路径追踪的理论,它的目的甚至以后所有光线追踪算法的主要目的都是为了更好的得到渲染方程的近似解。
那么具体该怎么对于一个不能直接求解的积分进行估值呢?这就不得不先谈下蒙特卡罗了。
1.蒙特卡罗(Monte Carlo)
蒙特卡罗(Monte Carlo,简称MC)本质上是一种数学思想,就像转化、归纳等数学思想一样,或者说是一类方法的统称,并不是指某个具体的算法。而且本身很简单,就是通过一些随机数来解决问题。蒙特卡罗起源于1777年法国数学家在一个正方形和其内接圆组成的形状里用投针法求圆周率π的值。再举个类似的例子,求《蒙娜丽莎》画中人脸的面积——抓把小米随机散在画上,数一数位于脸部区域的小米数,求出占总小米数的比例再乘以画的总面积,这就算运用了蒙特卡罗。它的神奇之处就在于用愚者都能理解的方法解决了智者都无可奈何的难题。
设被积函数f(x)在长度为L的积分区间[a,b]上的积分可表示为以下形式:

例一. 求以下积分的值:

此例中L=2,被积函数f(x)=e^(-x³)。这是个超越积分,因为e^(-x³)的原函数并不是个初等函数,无法用常规的分部积分法直接求解,但我们可以对它就行估值。
(1)思路一:利用定积分的定义
根据一维积分的几何意义,积分值就是函数曲线在积分区间内所围区域的面积S,因此就转化成了求面积近似值的问题。可将所围区域分割成N个等宽的小矩形:

图2-1 将 所围面积分割成10个等宽的小矩形面积
则每个小矩形宽为L/N,第i个小矩形高为f(xi),求所有小矩形的面积和:
当N趋于无穷大时,就是积分值。这种方法当然可以,依据的是定积分的定义。
(2)思路二:利用平均值
依然是数形结合思想,不过是想用一个固定值代替积分区间内的所有的函数值,这样所围区域就变成了一个矩形:

图2-2 y=f(x)拔
如果要保证该矩形的面积和f(x)函数曲线围成的面积相等,那么这个固定值显然就是所有f(x)的平均值,记为f(x)拔,再根据矩形面积公式:
其结果与式2-2相同,即一个积分的近似值可以表示成以下形式:
我们运用蒙特卡罗思想用计算机模拟这个过程:在[0,2]内生成N个随机数,每个随机数带入e^(-xi³) 求值相加后的结果再乘以(2/N)。样本数N越大估值越准确,相应的计算开销也会增加。
然而,从例一的函数图像可以看出,在自变量接近2时的函数值都普遍很小,而式2-4的运算本质上是对一堆函数值求和,如果某个随机数经过一些计算后带入求和的是一个很接近0的数,那可以认为这次计算是浪费的。我们希望选取的自变量的函数值要尽可能大一些,这样对最后求和的贡献更大。由于在模拟时用计算机直接生成的随机数是服从于均匀分布的,所以在本例中使用这种分布的采样点合理性欠佳。
PS.数学中的随机变量和计算机中的随机数
在数学中,随机变量一般用大写字母表示,比如X,可以是X~U(a,b),也可以是X~E(λ)、X~P(λ)、X~N(μ,σ2)、X~La(α,β),甚至可以是X服从于某个自定义的分布,这样的X在数学上都叫做随机变量。
在计算机中,高级编程语言都自带API可以直接生成[0,1]内的随机数,这个过程用数学语言描述就是生成了随机变量X且 X~U(0,1),若要生成Y~U(a,b)只需对X做线性变换即可。所以若无特别说明,直接利用计算机生成的随机数一定都是指均匀分布的随机变量。而想用计算机生成其它分布的随机变量则需要人工进行算法设计,后文会讲。
现在在这谈论的采样点选择的问题对应到路径追踪里实际上就是反弹光线的方向选择问题,而每条光线将来还要继续不断反弹做递归计算,所以这个选择如果做的不合理的话对最终渲染质量和时间开销的打击都是很大的。而解决方案下面就会讲到。
下图是我用蒙特卡罗渲染的效果。这里我考虑了背景因素,可以看到玻璃兔子下边折射出了蓝天白云的景象而兔子上边则是地上植被的景象,并且能看到铝合金底座对黄金材质的人偶(Substance的吉祥物MeetMat)的反射现象,以及黄金这种带有一定粗糙度的材质对周围环境的模糊反射。DRT那种从材质表面发射多条射线求平均的思想,当射线在一个结构复杂的物体表面来回反弹或者在一个复杂场景中的多个物体表面之间来回反弹时就极易导致指数爆炸,而蒙特卡罗就可以解决该问题。至于WRT中提到的射线一变三的问题,还要再配合拒绝采样(Reject
Sampling,简称RS)来解决。

图2-3 MC的实践
下图是我回到室内场景用蒙特卡罗渲染的效果。由于加入了朗伯体在法线半球内随机反弹光线的机制,首先导致顶上也有光线到达,所以终于看到了天花板;然后初步实现了全局光照的效果,人偶不再有死黑的背光面,观察发现人偶左侧明显受到了左侧绿墙的影响,而右侧也较弱的反射了右侧红墙的颜色,物体在地面上的投影也有了墙壁反射的颜色。
图2-4 MC的实践
不过问题也是显而易见的。首先最大的问题就是图像噪点现象太突出;再者就是这个光照效果也不是太正确,从素描的角度看的话人偶头顶的高光不明显,下颌和腹部接受白地面的反光也不明显,还有弧形脸上的明暗交界线也不明显,这些不明显加起来其实就是画素描时常说的一个问题——“灰”。在以前我们学素描的过程中相信大家都经历过“脏、花、灰”之类的通病,对这种画面问题应该会有比较敏锐的直觉。
2.重要性采样(Importance Sampling)
上面两种积分估值的思路都是转换成了面积问题利用几何来求解,下面我们试图从代数上寻找答案。
思路三:利用期望和大数定律
(1)期望
根据连续型变量的期望的定义:若随机变量X在[a,b]内服从于PDF为p(x)的分布,则X的期望为:
p(x)表示概率密度函数(Probability Density Function,简称PDF),PDF总是会跟积分联系在一起,PDF的原函数就是累积分布函数(Cumulative Distribution Function,简称CDF),记为F(x)。
PS:有些人喜欢把CDF称为概率分布函数或简称为分布函数,我个人是很不推荐的,一来呢其英文缩写容易跟PDF搞混,二来呢我觉得累积两个字能提醒你这个函数在定义内是单调递增的,进而记住是P{x≤a},即含义是随机变量X小于等于某个值x的概率,比如假设人的身高X的CDF在180的值就是全世界矮于等于180的人数除以全世界总人数。与CDF相对的P{x≥a}叫做生存函数(Survival Function),也是个概率论的知识点,常用于生存分析。
设Y=f(X),则有定理:
这个定理的意义在于我们可以直接求出E(Y),而避开了按照期望的定义需要先求出Y的PDF这一复杂计算,然而求期望依旧是一个复杂的积分过程,因为f(X)还在被积函数里所以我们还是解不出来。
PS:严谨一点的话,在f(x)无法求出原函数的情况下,f(x)p(x)是有可能直接求出原函数的,比如:
此积分无法求解(不过这个积分比较特殊可以用极坐标法解出,假设你不知道)。取p(x)=x,可以自行验证这是一个合法的PDF,此时:

现在一眼就看出答案了——用换元法设u=x²,则有

这个得p(x)很碰巧时才能出现这种情况。由于渲染中我们的积分对象是渲染方程,其中含有递归项,所以是一定无法直接求解的,因此我们一直关注的是寻找其它方法来估计积分值。从数学上来讲,确切的说渲染方程是一个Fredholm积分方程,而不是像式2-1所示形式的Riemann积分。
(2)大数定律
根据大数定律(Law of Large Numbers),对于独立同分布(Independently Identically Distribution,简称IDD)的随机变量X进行N次大量采样得到平均值X拔即为X的期望E(X)的一个估计值,即:
大数定理告诉了我们,求一个随机变量的期望的近似值,不管这个随机变量服从于什么样的分布,就是对其进行大量采样的结果求平均值。所以对于E(Y)当然也成立:

结合式2-6,可得:
我们得到了一个估计积分的新方法,但这个公式左边的积分项和式2-1形式上还有点区别需要变换一下,不妨设g(x)=f(x)p(x) ,于是有:
所以得到:
结论:求一个积分的近似值,可以生成大量的服从于PDF为p(x)的分布的随机数代入右侧表达式求值。
PS:记法上,到底写成f(x)还是g(x)是无所谓的,所以式2-11也可以写成:

若X~U(a,b)时,则有p(x)=1/L,带入上式后就可得到式2-4,进一步验证了上一节一直在使用服从于均匀分布的随机变量来处理问题。
运用蒙特卡罗思想用计算机模拟这个过程:
第一步,构造一个PDF;
第二步,生成服从于该PDF的随机数;
第三步,将随机数带入被积函数和PDF分别求值后相除;
第四步,大量重复第二步和第三步,并对其结果求和后再求平均值。
首先第一步该构造一个怎样的PDF呢?
(1)必须合法。比如要满足在整个积分区间上对PDF的积分一定是1、PDF的值域下限必须大于等于0等;
(2)能使采样点选取的更合理。蒙特卡罗存在的问题是统一进行了均匀采样,改进的思想就是要针对具体的被积函数在对积分贡献更高的区域放置更多的采样点,已体现出这部分区域的重要性,这就是重要性采样(Importance Sampling,简称IS)的核心思想。具体方法就是选取一个与被积函数曲线形状相似的PDF,进而能按照此PDF进行采样。
(3)要简单。到底怎样的PDF就算简单了呢?至少要能方便的求出它的原函数CDF以及CDF的反函数,原因跟下一步有关。
然后第二步牵扯到一个问题:如何用计算机生成服从于某个PDF的随机数呢?
设函数y=F(x),根据CDF的定义:


图2-5 Y=F(X)
PS:函数y=F(x)表示的是两个变量之间的映射关系,具体应用在统计学里也可以理解为当随机变量X取数值x时,此时对应的随机变量Y的取值为数值y。
由于CDF一定是单调递增函数,若再保证CDF严格递增,则CDF的反函数一定存在,记为x=F^-1(y),所以X=F^-1(Y) ,代入上式得:

图2-6 X=F^-1(Y)
又因为反函数的单调性与原函数一致,所以函数x=F^-1(y)也单调递增。又因为当X为x时对应的值为F(x),所以当值域X≤x对应Y的取值范围为[0,F(x)],即X=F^-1(Y)≤x的概率等于Y≤F(x)=y的概率,即:

而根据CDF的定义:
所以由1234得:
结论:设随机变量X服从于CDF为F(x)分布,对于把X带入F(x)的反函数后形成的新的随机变量Y与X同分布,即随机变量Y也服从于CDF为F(x)的分布。
由于CDF的值域一定是[0,1],所以我们可以用计算机直接生成X~U(0,1),代入F^-1(x)后得到的随机变量Y=F^-1(X)即为所求。这种生成服从于某个分布的随机变量的方法也称为Smirnov Transform。
最后再进行第三步和第四步常规的代数运算得到积分的估值。
例二.在[0,1]内构造一条浴盆曲线的积分:

这个积分我们能直接求出为1.1111111……。通过用计算机模拟,分别用蒙特卡罗和PDF为y=2x的重要性采样,每次试验选取1w个采样点,重复100次的试验数据如下:
|
蒙特卡罗 |
重要性采样 |
样本期望 |
0.111262202 |
0.111190937 |
样本标准差 |
0.002151918 |
0.001420695 |
表2-1 两种方法对例二采样的统计数据对比
可见,重要性采样对例二中积分的估值质量更高一些。
重要性采样是基于蒙特卡罗做的改进,可以认为重要性采样是蒙特卡罗的一种。重要性采样给每次采样结果一个权重,把简单的算数平均值变成了加权平均值。根据式2-11,如果按照分割矩形的方法直观来看的话,PDF越大对应矩形宽度会越窄,矩形排列更密集,在有限的矩形数量时对曲线所围面积估计越准,如下图所示:

图2-7 重要性采样对应的矩形分割示意
事实上,重要性采样是蒙特卡罗进行方差缩减的方法中的其中一种,理想情况下,如果选取的p(x)恰巧其实是f(x)与其期望之比,那么方差为0。
下图是我用IS渲染的效果。朗伯体的反射光线在法线半球内不再进行随机采样。现在人偶的体积感较之图2-4有了明显提升,我们更能感受到在物体结构的影响下表面与光的互动关系。噪点问题也有了些许改善但仍然严重。

图2-8 IS的实践
这一篇文章先写到这,下一篇文章会结合几个重要性采样的例子来具体求解。到目前为止我们仍然得到不了一张能让人满意的高质量图像,不过不要着急,毕竟光线追踪这门课程本身的体系就比较庞大,下文也会给出答案并做进一步深入。