前言
本篇博客作为Games101关于Ray Tracing部分的复习总结。
蒙特卡洛积分
在此之前先复习一下概率密度函数(PDF, Probability Density Functions),可以理解为取某一随机变量x的概率是多少的函数 p ( x ) p(x) p(x)

PDF具有如下特点
- $ p(x) >= 0 \ and \int{p(x)}=1$
- E [ X ] = ∫ x p ( x ) d x E[X]=\int{xp(x)}dx E[X]=∫xp(x)dx
对于复杂函数 f ( x ) f(x) f(x)求定积分,用牛顿莱布尼茨公式十分困难,这时候需要用到另一种积分方法,也就是下面的蒙特卡洛积分(Monte Carlo Integration)
∫ f ( x ) d x = 1 N ∑ i = 1 N f ( X i ) p ( X i ) X i ∼ p ( x ) \int{f(x)}dx = \frac{1}{N}\sum_{i=1}^{N}{\frac{f(X_i)}{p(X_i)}}\ \ \ \ \ \ \ X_i \sim p(x) ∫f(x)dx=N1i=1∑Np(Xi)f(Xi) Xi∼p(x)
这里的做法可以理解为按 p ( x ) p(x) p(x)的概率取随机变量 X i X_i Xi,在图形学里也就是随机采样函数值再累加做平均。
一些需要注意的地方:
- 采样点N越多,方差越小
- 在x上采样,也是在x上求积分
路径追踪
路径追踪(Path Tracing)根据渲染方程实现全局光照,或者说基于物理的渲染。个人理解叫追踪的原因是因为光线本来是打到物体上再多次弹射到人眼所以人才能看见东西,追踪就是逆着这个路径,从眼睛“射出”光线,“所见即所得”。
Whitted-Style光线追踪的局限性
Whitted-style ray tracing是一种理想化模型,存在一些问题,光线只有反射或折射,在漫反射材质停止弹射,没有漫反射,因此Whitted-style ray tracing渲染出来的效果只有光源和直接光照,如下图左侧所示。我们理应看到的效果是下图右侧这种,就像是照片拍出来的一样,需要实现这样真实的效果就需要用到路径追踪。

渲染方程
根据辐射度量学,我们一通操作最终求得渲染方程,即在物理意义上得出的正确的能够描述我们现实生活所见的方程。
L o ( p , ω o ) = L e ( p , ω o ) + ∫ Ω L i ( p , ω i ) f r ( p , ω i , ω o ) , ( n ⋅ ω i ) d ω i L_o(p,\omega_o)=L_e(p,\omega_o)+\int_{\Omega}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o),(n\cdot\omega_i)d\omega_i Lo(p,ωo)=Le(p,ωo)+∫ΩLi(p,ωi)fr(p,ωi,ωo),(n⋅ωi)dωi
某着色点到相机的radiance由其本身发光产生的部分,和从所有立体角接受到的radiance并反射到相机这个一定角度的部分组成。
为了方便计算和统一,这里将着色点本身产生的 L e ( p , ω o ) L_e(p,\omega_o) Le(p,ωo)省略掉,我们更关注如下这个渲染方程(或者说反射方程)。
L o ( p , ω o ) = ∫ Ω L i ( p , ω i ) f r ( p , ω i , ω o ) , ( n ⋅ ω i ) d ω i L_o(p,\omega_o)=\int_{\Omega}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o),(n\cdot\omega_i)d\omega_i Lo(p,ωo)=∫ΩLi(p,ωi)fr(p,ωi,ωo),(n⋅ωi)dωi
根据蒙特卡洛积分可以得到
L o ( p , ω o ) ≈ 1 N ∑ i = 1 N L i ( p , ω i ) f r ( p , ω i , ω o ) ( n ⋅ ω i ) p ( ω i ) L_o(p,\omega_o)\approx\frac{1}{N}\sum_{i=1}^{N}\frac{L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n\cdot\omega_i)}{p(\omega_i)} Lo(p,ωo)≈N1i=1∑Np(ωi)Li(p,ωi)fr(p,ωi,ωo)(n⋅ωi)
这个方程可以理解为按 p ( ω i ) p(\omega_i) p(ωi)的概率随机选取N个方向 ω i \omega_i ωi采样,再累加求平均。
注意,这里的 ω o \omega_o ωo和 ω i \omega_i ωi都是由着色点向外的
直接光照
仅在直接光照下,考虑计算radiance的场景如下

由此可以得到着色方程的shade伪代码
shade(p, wo)
Randomly choose N directions wi~pdf
Lo = 0.0
For each wi
Trace a ray r(p, wi)
If ray r hit the light // 弹射的光线必须要能打到光源,这个过程中不能被其他物体阻挡
Lo += (1 / N) * L_i * f_r * cosine / pdf(wi)
Return Lo
其中p是着色点,wo是入射光线,wi是出射光线,r是由眼睛穿过需要渲染的像素到实际着色点的光线。
间接光照
如果路径追踪经过一次弹射后没有打到光源而是打到物体,就构成了间接光照。我们考虑如下这个场景

P接收到的radiance是由Q反射光源得到的,因此为了求点P反射的光,我们需要先求点Q反射出的光。显然,这里用到了递归的思路,将上述shade伪代码稍作修改即可完成
shade(p, wo)
Randomly choose N directions wi~pdf
Lo = 0.0
For each wi
Trace a ray r(p, wi)
If ray r hit the light
Lo += (1 / N) * L_i * f_r * cosine / pdf(wi)
Else If ray r hit an object at q
Lo += (1 / N) * shade(q, -wi) * f_r * cosine / pdf(wi)
Return Lo
可以看到只是增加了一个判断条件,判断从点p反射出的的光线wi是射到了光源还是其他物体,如果是光源直接加,如果是其他物体则递归先求L_i,也就是shade(q,-wi)。
注意,由于shade的参数要求光线方向向外,所以对于点Q,需要将wi取反。
两个问题
上述shade伪代码还存在两个问题
- 每次都会采样N个方向,最终导致指数爆炸
- 递归没有出口
对于第一个问题,可以用下图形象表示

在计算间接光照时,由于每次都要采样N个点,最终导致 N x N^x Nx的指数增长。但有一个数是可以避免这种情况,那就是N取1。
在采样点每次采样只取1个方向,我们的shade伪代码又有了如下改动
shade(p, wo)
Randomly choose ONE direction wi~pdf(w)
Trace a ray r(p, wi)
If ray r hit the light
Return L_i * f_r * cosine / pdf(wi)
Else If ray r hit an object at q
Return shade(q, -wi) * f_r * cosine / pdf(wi)
由前文的蒙特卡洛积分我们知道,N取越大越等接近原定积分的值,N取越小结果的噪声或者说偏差就会越大。我们这里取N为1,可想而知会有很大偏差。为了解决这个问题,可以采取对同一像素多次采样的方法,根据需要设置对每个像素采样的次数(spp, samples per pixel)。实现的效果如下图所示,可以看到间接光照每次采样的弹射结果虽然都不一样,但多次采样求平均后就可以有效减小噪声。

上述的采样方法用伪代码描述如下,与shade伪代码类似
ray_generation(camPos, pixel)
Uniformly choose N sample positions within the pixel
pixel_radiance = 0.0
For each sample in the pixel
Shoot a ray r(camPos, cam_to_sample)
If ray r hit the scene at p
pixel_radiance += 1 / N * shade(p, sample_to_cam)
Return pixel_radiance
现在我们来解决第二个问题,这个递归算法没有出口所以会一直持续下去,因此我们需要认为设置一个出口。
这里采用的方法是俄罗斯轮盘赌。设定一个值Russian Roulette (RR),光线每次弹射前需先与RR“赌博”,如果成功则继续弹射,否则停止。我们可以设在概率P_RR下,将射出光线返回的渲染方程结果Lo除以P_RR,那么将会有1-P_RR的概率停止弹射直接返回0。在此假设下,我们可以得到期望值E,是与我们期望得到的Lo相同的。
E = P ∗ ( L o / P ) + ( 1 − P ) ∗ 0 = L o E=P*(Lo/P)+(1-P)*0=Lo E=P∗(Lo/P)+(1−P)∗0=Lo
由此,我们的shade伪代码再次改动如下
shade(p, wo)
Manually specify a probability P_RR
Randomly select ksi in a uniform dist. in [0, 1]
If (ksi > P_RR) return 0.0;
Randomly choose ONE direction wi~pdf(w)
Trace a ray r(p, wi)
If ray r hit the light
Return L_i * f_r * cosine / pdf(wi) / P_RR
Else If ray r hit an object at q
Return shade(q, -wi) * f_r * cosine / pdf(wi) / P_RR
优化
到此,我们终于得到了一个正确的路径追踪算法,的的确确能正确渲染出一个像素,但是目前这种做法的效率其实并不高。在计算直接光照时,如果我们按照上面得伪代码计算直接光照部分,那么每次只会随机射出一根光线,实际上很难打到光源上,计算效率低下。
为了提高成功采样到直接光照的概率,我们可以对直接光照计算的方法进一步修改,在弹射点进行多次采样。与间接光照递归计算相比,因为直接光照只会弹射一次到光源,所以这样并不会造成指数爆炸。

如上图所示,随着着色点采样数的提高,成功采样的几率自然会提高。但是当光源很小时,如果我们继续采用平均采样,那么弹射的数跟光线中可能只有少数几根能打到光源上,其他采样光线可以说都是被浪费掉了,这同样是我们不希望看到的。
如果能对更加确定的东西采样,那么这个问题自然迎刃而解。考虑我们是想要找到弹射点能打到光源的有效光线,所以反过来考虑,我们对光源进行采样。并且实际上,蒙特卡洛积分也是允许我们自定义采样方法的。
L o ( p , ω o ) = ∫ Ω L i ( p , ω i ) f r ( p , ω i , ω o ) ( n ∗ ω i ) d ω i L_o(p,\omega_o)=\int_{\Omega}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n*\omega_i)d\omega_i Lo(p,ωo)=∫ΩLi(p,ωi)fr(p,ωi,ωo)(n∗ωi)dωi
上面是我们在前文我们得到的渲染方程,可以看到我们的积分目标是 d ω i d\omega_i dωi,可以解释为我们是从着色点出发考虑这个积分的。如果我们对光源进行采样,积分对象自然也需要变化,所以我们需要找到 d ω i d\omega_i dωi和 d A dA dA之前的关系。下图很好展现了 d ω i d\omega_i dωi和 d A dA dA之前的关系。

回顾以前对立体角 ω \omega ω的定义
d ω = d A r 2 d\omega = \frac{dA}{r^2} dω=r2dA
单位立体角就是单位面积与距离平方的比值。因此在这里我们可以类比得出相似的关系,类似Bllin-Phong模型,由光源投射出的光线方向与半球面有一个夹角,所以这里的 A = A c o s θ ′ A=Acos\theta' A=Acosθ′,而r自然是两点之间的距离 r = x ′ − x r=x'-x r=x′−x。由此,我们得到 d ω i d\omega_i dωi和 d A dA dA之前的关系
d ω = d A c o s θ ′ ∣ ∣ x ′ − x ∣ ∣ 2 c o s θ ′ = n d\omega = \frac{dAcos\theta'}{||x'-x||^2} \ \ \ \ \ cos\theta'=n dω=∣∣x′−x∣∣2dAcosθ′ cosθ′=n
将上式代入渲染方程,我们就可以得到按光源采样的新蒙特卡洛积分形式
L o ( p , ω o ) = ∫ Ω L i ( p , ω i ) f r ( p , ω i , ω o ) c o s θ d ω i = ∫ Ω L i ( p , ω i ) f r ( p , ω i , ω o ) c o s θ c o s θ ′ ∣ ∣ x ′ − x ∣ ∣ 2 d A ≈ 1 N ∑ i = 1 N L i ( p , ω i ) f r ( p , ω i , ω o ) c o s θ c o s θ ′ ∣ ∣ x ′ − x ∣ ∣ 2 p ( A ) 其中 c o s θ = n ⋅ d ω c o s θ ′ = n ⋅ ( − d ω ) \begin{aligned} L_o(p,\omega_o) &=\int_{\Omega}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)cos\theta d\omega_i\\ &=\int_{\Omega}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)cos\theta\frac{cos\theta'}{||x'-x||^2}dA\\ &\approx\frac{1}{N}\sum_{i=1}^{N}\frac{L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)cos\theta\frac{cos\theta'}{||x'-x||^2}}{p(A)}\\ 其中&cos\theta=n\cdot d\omega \ \ \ \ \ cos\theta'=n\cdot(-d\omega)\\ \end{aligned} Lo(p,ωo)其中=∫ΩLi(p,ωi)fr(p,ωi,ωo)cosθdωi=∫ΩLi(p,ωi)fr(p,ωi,ωo)cosθ∣∣x′−x∣∣2cosθ′dA≈N1i=1∑Np(A)Li(p,ωi)fr(p,ωi,ωo)cosθ∣∣x′−x∣∣2cosθ′cosθ=n⋅dω cosθ′=n⋅(−dω)
对于上式,由pdf性质我们知道按 p ( A ) p(A) p(A)对面积积分的结果应该是1( ∫ p ( A ) d A = 1 \int p(A)dA=1 ∫p(A)dA=1),自然 p ( A ) = 1 A p(A)=\frac{1}{A} p(A)=A1。
针对直接光照部分,按面积对光源采样的shade伪代码如下所示
shade(p, wo)
# Contribution from the light source.
Uniformly sample the light at x’ (pdf_light = 1 / A)
L_dir = L_i * f_r * cos θ * cos θ’ / |x’ - p|^2 / pdf_light
# Contribution from other reflectors.
L_indir = 0.0
Test Russian Roulette with probability P_RR
Uniformly sample the hemisphere toward wi (pdf_hemi = 1 / 2pi)
Trace a ray r(p, wi)
If ray r hit a non-emitting object at q
L_indir = shade(q, -wi) * f_r * cos θ / pdf_hemi / P_RR
Return L_dir + L_indir
除此之外,还有一个遗漏了的问题,那就是在光源上采样的光线必须得能射到着色点上,也就是说着色点 p p p和光源 x ′ x' x′之间不能有物体遮挡。

对上述伪代码稍作修改,判断一下是否有物体遮挡
shade(p, wo)
# Contribution from the light source.
Uniformly sample the light at x’ (pdf_light = 1 / A)
Shoot a ray from p to x’
If the ray is not blocked in the middle // 判断是否有物体遮挡
L_dir = L_i * f_r * cos θ * cos θ’ / |x’ - p|^2 / pdf_light
# Contribution from other reflectors.
L_indir = 0.0
Test Russian Roulette with probability P_RR
Uniformly sample the hemisphere toward wi (pdf_hemi = 1 / 2pi)
Trace a ray r(p, wi)
If ray r hit a non-emitting object at q
L_indir = shade(q, -wi) * f_r * cos θ / pdf_hemi / P_RR
Return L_dir + L_indir
自此,我们终于得出了正确的,高效的渲染方式。
总结
为了尽可能展现我们看到的现实世界,实际上也就是渲染一张图片,再根本点就是渲染每一个像素。对每个像素进行多次采样,根据路径追踪得到的渲染方程计算每一次采样结果再累加平均,得到该像素理想的渲染结果。
对于渲染方程,是根据路径追踪得到的基于物理的渲染。与最初我们做的工作相比,之前是简单粗暴假设光线由半球均匀随机采样得到,存在递归指数爆炸,没有出口,效率低下等问题,而现在我们考虑radiance是有两部分组成,分别采用不同的方法计算。
- 直接光照,对光源按面积进行采样
- 间接光照,按半球随机采样一次,并且采用RR控制递归深度

参考
Games101:作业7(含提高部分)_Q_pril的博客-CSDN博客_games101 作业7