700字范文 > 闫令琪：Games101 现代计算机图形学-光线追踪(三)：渲染方程和路径追踪path ray tr

闫令琪：Games101 现代计算机图形学-光线追踪(三)：渲染方程和路径追踪path ray tr

时间：2022-07-09 16:04:15

文章目录

0 whitted光线追踪的局限1 辐射度量学1.1 光线的表示 Radiance1.2 物体表面上一个点的亮度 Irradiance1.3 BRDF(Bidirectional ReflectanceDistribution Function)和材质1.4 渲染方程2 渲染方程的求解2.1 渲染方程视角下的光线传播2.2 蒙特卡洛积分2.3 渲染方程的蒙特卡洛形式2.4 递归求解的构建和终止条件2.4.1 采样指数发散的问题2.4.2 采样终止条件：俄罗斯轮盘赌2.5 求解效率的提高3 渲染方程的代码实现3.1 path tracing伪代码3.2 作业代码4 path ray tracing的优化方向

路径追踪技术（Path tracing，PT）已经是当下工业中离线渲染使用的主流技术，不管是商业渲染器如皮克斯的RenderMan，Solid Angle的Arnold等，还是迪士尼的in-house渲染器Hyperion以及Weta Digital的Manuka都是基于路径追踪技术。

但是搜遍全网也没有一篇能够通俗易懂的讲清、讲全一个path tracing，本文希望达到这个目标！！！

0 whitted光线追踪的局限

闫老师课上举了下图这两个例子说明whitted ray tracing的局限性。其实根本原因在于whitted光线追踪是非物理的，首先对于材质的处理，没有办法模拟真实材质那种漫反射，折射和反射都是理想情况，漫反射直接渲染颜色，这样就无法表现下图右边那种情况；其次一个点的颜色其实是受直接光照和间接光照影响的，但是whitted光线追踪是无法反映间接光照的影响，路径追踪是求解渲染方程，在渲染方程中，某个点的颜色其实该点周围立体角上所有光线对该点的积分。

1 辐射度量学

1.1 光线的表示 Radiance

这个概念是渲染方程和路径追踪中最关键的概念。因为渲染方程要用数学的方式去描述光的传播，那很自然的一个问题就是：

光线用数学怎么描述？

上图是闫老师课件中给的光线的描述radiance。

光的能量功率Φ非常好理解，那么怎么由Φ定义出radiance呢?

光的能量功率Φ可以认为由光源发出的这一圈一圈的能量，那么一条光线便是这圈上的一个点，对于半径为1的球，上面一个点能量表示为上图中的E，对于外圈的球上一个点的能量，表示为E'，分母除了和球的立体角相关，还和r2r^2r2相关，而r2r^2r2可以表示为面积。从这两个位置上的能量可以发现，这个球上任意一点的能量，和立体角、面积两个变量相关。

那么可以很自然地推出一条光线的数学表示：

这个分母上的cosθ怎么理解？

1.2 物体表面上一个点的亮度 Irradiance

搞定了光线的表示，还有一个定义需要解决，物体表面一个点是需要接受四面八方的漫反射光线，这个点的最终亮度怎么表示？

一图胜千言：

对于上图中的那个点的亮度，显然只需要把立体角内所有的光线进行积分即可。

从上面光线的表示，可以很容易写出微分形式：

则某个点的亮度，写成积分即可：

这个变量，叫irradiance：是表示一个点接收的周围光线的亮度。

1.3 BRDF(Bidirectional ReflectanceDistribution Function)和材质

BRDF的中译名叫双向反射分布函数。

这里是查阅资料时发现的比较全面的综述和细致的讲解：

《RTR3》提炼总结】(六) 第七章 · 高级着色：BRDF及相关技术。

基于物理着色：BRDF

我们看到一个表面，实际上是周围环境的光照射到表面上，然后表面将一部分光反射到我们眼睛里。双向反射分布函数BRDF（Bidirectional Reflectance Distribution Function）就是描述表面入射光和反射光关系的。入射光就是上面讲的Irradiance，出射光就是radiance。

BRDF定义为：
对于一个方向的入射光，表面会将光反射到表面上半球的各个方向，不同方向反射的比例是不同的，我们用BRDF来表示指定方向的反射光和入射光的比例关系，数学公式为

基于物理着色：BRDF这篇文章里提供了两种解释，为什么定义成这两个量的比值，而不是dLL\frac{dL}{ L}LdL等之值，从数学上讲，当出射方向立体角趋于0时，这个比值的极限也是0，所以没有意义，从实验的角度，入射方向的radiance却非常难测，所以采用这个比值。

公式中为什么是dLr(ωr)dL_r(ω_r)dLr(ωr)，因为定义里就是指定方向的入射光和出射光，所以这里的ω\omegaω是被微分的量。

1.4 渲染方程

渲染方程要解决的问题就是，怎么样描述看到的物体表面一个点上的光亮。这个点进入观察点的那个物理量，啰嗦一点，其实就是那个点漫反射到观察点的光线，而光线就是radiance。

由上面的BRDF公式可以很容易写出：

dLr(ωr)=frLi(ωi)cosθidωidL_r(ω_r) = f_r L_i(ω_i) cos\theta_id\omega_idLr(ωr)=frLi(ωi)cosθidωi

写成积分形式就是该光线完整的表达：

Lr(ωr)=∫frLi(ωi)cosθidωiL_r(ω_r) =\int^{}_{}f_r L_i(ω_i) cos\theta_i{d\omega_i}Lr(ωr)=∫frLi(ωi)cosθidωi

接下来用闫老师的课上PPT表示：

老师课上的公式中的符号说的并不清楚，网上发现了这张图，说的还挺不错的：

2 渲染方程的求解

2.1 渲染方程视角下的光线传播

对于这个漫反射来说，很自然会有疑问：绿色墙上的颜色会漫反射到立方块上，立方块上的也会漫反射到墙上，这就有点鸡生蛋、蛋生鸡的循环，这个过程到底是怎么样？我们求解渲染方程的话，终止光线弹射的条件是什么？

假如我们把时间放慢到可以观察光线传播的程度，以物体表面一个点Point点为例，观察：

①可以看到会有一根光线打到这个点，这个时候物体有一次亮度；

②所有物体表面都有一根光线被打量后，光线开始第一次弹射，Point点会被周围弹射的光线打中，叠加刚才的直接光照，会有第二次光亮；

③开始第二次弹射……

④第三次弹射等等……

闫老师的课上给了一个例子：

然而，到目前为止，跟渲染方程还没什么关系？？是的，哈哈哈……

闫老师课上通过一系列神操作，通过渲染方程推导出下面的数学表达式，即某点的光照可以分为自己发的光EEE，直接光照KEKEKE，二次弹射K2EK^2EK2E等。上面的思想实验，其实是可以通过数学推导出来，虽然我真的没听懂……[捂脸]

光栅化的局限就在于，它从原理上，是只能反映第一次光照，因此无法呈现全局光照

2.2 蒙特卡洛积分

渲染方程显然没有解析解，并且数值离散解也几乎不可能，所以这时候就得祭出蒙特卡洛大法。以下内容来源于：蒙特·卡罗(Monte Carlo)积分详解。

对于一个连续函数f，它的积分公式为：

F=∫abf(x)dxF = \int _{a}^{b}f(x)dx F=∫abf(x)dx

对应的，f的蒙特·卡罗积分公式如下：

FN=1N∑i=1Nf(Xi)pdf(Xi)F^{N} = \frac {1}{N}\sum _{i=1}^{N}\frac {f(X_{i})}{ pdf(X_{i}) } FN=N1i=1∑Npdf(Xi)f(Xi)

蒙特卡罗最关键的就是理解这条公式了。其他延伸探讨都可以暂时忽略。那么这条公式如何理解呢？首先第一点是，虽然这条公式没有积分符号∫\int∫，但是它认被称为积分，这是因为这公式的作用相当于在对f(x)做积分，只不过不那么“精确”，即蒙特·卡罗积分是对理想积分的近似。

那么这个近似是如何完成的？很简单，核心就是两个字：采样(Sampling)。对一个连续函数的采样方法是在该函数的定义域中随机挑N个值，并求出对应的N个f(Xi)f(X_{i})f(Xi)，就得到了样本集合。再对这些样本集合做一些换算，就可以得到一个近似的积分了。对于蒙特·卡罗积分，采样样本越多，就越逼近真实的积分结果，这是蒙特·卡罗积分的最核心特性。

继续观察上面的公式，里面还有一个极其重要的参数：pdf(probability distribution function，概率分布函数)。pdf还有个近亲pmf，下面小节详解pdf、pmf的由来。

pdf和pmf

pmf(probability mass function)，指的是离散的随机变量的概率分布函数pdf(probability distribution function)，指的是连续的随机变量的概率分布函数

离散的随机变量X的数学期望为：

E[X]=∑xipmf(xi)xiE[X] = \sum _{ x_{i} }pmf(x_{i})x_{i} E[X]=xi∑pmf(xi)xi

连续的随机变量X的数学期望为：

E[X]=∫−∞∞pdf(x)xdxE[X] = \int ^{\infty }_{-\infty }pdf(x)xdxE[X]=∫−∞∞pdf(x)xdx

pdf和pmf名字接近，含义也是接近。pdf、pmf函数的参数都是样本值x，返回值是概率，即表示一个样本出现的概率，所有样本的出现概率之和(概率的积分)应等于1。要注意的是，pdf、pmf的存在说明有可能每个样本的出现概率都是各不相同的。

pmf

pmf的简单例子就是基于均匀分布的离散的随机变量X，此时pmf(Xi)pmf(X_{i})pmf(Xi)恒等于1N\frac{1}{N}N1，含义是每个随机样本的出现概率等于1样本总数\frac{1}{样本总数}样本总数1。

通过这个例子也印证了pmf的性质：pmf函数的所有结果值之和等于1。

pdf

借用/的一个很好的例子来说明：

这个例子中，目标问题是求出该函数[a,b]段曲线下方的面积(最后一幅图的黑色区域)，也就是要求该函数[a,b]段的积分。基于蒙特·卡罗积分的解法，就要用上面给出的公式：

FN=1N∑i=1Nf(Xi)pdf(Xi)F^{N} = \frac {1}{N}\sum _{i=1}^{N}\frac {f(X_{i})}{ pdf(X_{i}) } FN=N1i=1∑Npdf(Xi)f(Xi)

在此图中，做了四次随机采样，得到了四个随机样本xix_{i}xi：x1、x2、x3、x4x_{1}、x_{2}、x_{3}、x_{4}x1、x2、x3、x4，并且进而得到了这四个样本的f(xi)f(x_{i})f(xi)值：f(x1)、f(x2)、f(x3)、f(x4)f(x_{1})、f(x_{2})、f(x_{3})、f(x_{4})f(x1)、f(x2)、f(x3)、f(x4)。（原文没有提及如何得到f(xi)f(x_{i})f(xi)。函数f是奇形怪状的，不太可能有表达式存在，难道是用尺子量的？暂且忽略这个事吧。）

有了这4个样本后，可以针对每一个样本求一个近似面积值，这个面积值等于f(xi)(b−a)f(x_{i}) (b - a)f(xi)(b−a)。为什么可以这样做呢？是因为每一个单独的样本是对原函数f的近似，即在每个样本中，认为f(x)f(x)f(x)恒等于f(xi)f(x_{i})f(xi)，从而让原函数曲线简化成一个矩形区域，而矩形的面积显然就是长(b-a)乘以宽f(xi)f(x_{i})f(xi)。

得到4个近似面积值后，再求出它们的均值(数学期望)，就完成了蒙特·卡罗积分。把上述流程汇总得到：

Area=14(f(x1)(b−a)+f(x2)(b−a)+f(x3)(b−a)+f(x4)(b−a))Area = \frac {1}{4}(f(x_{1})(b - a) + f(x_{2})(b - a) + f(x_{3})(b - a) + f(x_{4})(b - a)) Area=41(f(x1)(b−a)+f(x2)(b−a)+f(x3)(b−a)+f(x4)(b−a))

=14(b−a)(f(x1)+f(x2)+f(x3)+f(x4))= \frac {1}{4}(b - a)( f(x_{1}) + f(x_{2}) + f(x_{3}) + f(x_{4}) ) =41(b−a)(f(x1)+f(x2)+f(x3)+f(x4))

=14(b−a)∑i=14f(xi)= \frac {1}{4}(b - a)\sum _{i=1}^{4}f(x_{i}) =41(b−a)i=1∑4f(xi)

此时，对比下蒙特·卡罗积分公式：

FN=1N∑i=1Nf(Xi)pdf(Xi)F^{N} = \frac {1}{N}\sum _{i=1}^{N}\frac {f(X_{i})}{ pdf(X_{i}) } FN=N1i=1∑Npdf(Xi)f(Xi)

发现两个式子非常相似，对式子做下转换得到：

Area=14∑i=14f(xi)1b−aArea = \frac {1}{4}\sum _{i=1}^{4}\frac {f(x_{i})}{\frac {1}{b - a} } Area=41i=1∑4b−a1f(xi)

于是可以知道pdf(xi)pdf(x_{i})pdf(xi)等于：

pdf(xi)=1b−apdf(x_{i}) = \frac {1}{b - a }pdf(xi)=b−a1

这意味着，对于连续函数f，f的每个可能取值x的出现概率等于x的取值范围[a,b]的倒数1b−a\frac{1}{b-a}b−a1。

在实际应用场合，随机变量X要写成F(X)，即可能需要对X做一个转换再使用。这时候要注意F(X)的pdf不等于X的pdf。

蒙特·卡罗积分的数学期望等于理想积分？

对于下面的FFF和FNF^{N}FN：

F=∫abf(x)dxF = \int _{a}^{b}f(x)dx F=∫abf(x)dx

FN=1N∑i=1Nf(Xi)pdf(Xi)F^{N} = \frac {1}{N}\sum _{i=1}^{N}\frac {f(X_{i})}{ pdf(X_{i}) } FN=N1i=1∑Npdf(Xi)f(Xi)

是否随着N变大，FNF^{N}FN会逼近FFF？即FNF^{N}FN的数学期望是否等于FFF?Monte Carlo Methods in Practice文章中给出了推导过程：

第二行到第三行是最不好理解的。因为这里其实用到了新的知识点：Law of the unconscious statistician(简称：LOTUS)。LOTUS的应用情景是，已知随机变量X的概率分布，但不知道f(x)的分布，此时用LOTUS公式能计算出函数f(x)的数学期望。LOTUS的公式如下：

f(x)是离散函数时:

E[f(X)]=∑xif(xi)pmf(xi)E[f(X)] = \sum _{x_{i}}f(x_{i})pmf(x_{i}) E[f(X)]=xi∑f(xi)pmf(xi)

f(x)是连续函数时:

E[f(X)]=∫−∞∞f(x)pdf(x)dxE[f(X)] = \int _{-\infty }^{\infty}f(x)pdf(x)dx E[f(X)]=∫−∞∞f(x)pdf(x)dx

(建议对比第二小节开头的两条公式来理解)

有了LOTUS公式，再来看第二行到第三行的转换，就好理解了：

E[f(Xi)pdf(Xi)]=E[f(x)pdf(x)]=∫−∞∞f(x)pdf(x)pdf(x)dxE[ \frac {f(X_{i})}{pdf(X_{i})} ] = E[ \frac {f(x)}{pdf(x)} ] =\int _{-\infty }^{\infty}\frac {f(x)}{pdf(x)}pdf(x)dx E[pdf(Xi)f(Xi)]=E[pdf(x)f(x)]=∫−∞∞pdf(x)f(x)pdf(x)dx

=∫−∞∞f(x)dx=\int _{-\infty }^{\infty}f(x)dx =∫−∞∞f(x)dx

上面推导充分证明了，蒙特卡洛积分是无偏的。可以进行非均匀采样，结果依然是无偏

上面的内容已经足够继续进行下面的旅程，此处贴出两个知乎上的优质贴，后续继续嗑蒙特卡洛积分。

蒙特卡洛积分 - chopper的文章 - 知乎

Monte Carlo数学原理 - papalqi的文章 - 知乎

2.3 渲染方程的蒙特卡洛形式

渲染方程没法直接求解，求解的思路就是上面的蒙特卡洛积分。参考两个公式，可以很方便地写出来

不考虑自发光的情况的渲染方程形式：

这个方程理论上是可以求解的，pdf(ωi)pdf(\omega_i)pdf(ωi)可以假设就是半球面上均匀采样，frf_rfr可以初定一个固定的值，但是还有一个问题：

radiance该怎么表示？比如最初的灯光，它怎么表示Li(p,ωi)L_i(p,\omega_i)Li(p,ωi)？

2.4 递归求解的构建和终止条件

2.4.1 采样指数发散的问题

一图胜千言，先上图。

在P点的出射光线按照上面的方程形式，需要采样N条光线作为入射光线，这N条光线，每条假如与物体交点为Q，对每一个Q点又需要采样N条光线，这样光线弹射几次，N∗N∗N……N*N*N……N∗N∗N……这个计算量就爆炸了。

解决的思路也非常的简单粗暴：

N=1N=1N=1，每次只采样一条光线。

公式里的求和∑N\sum_{N}∑N里的NNN不是通过对像素光线相交点采样NNN，改为对每个像素点采样NNN条光线，每条光线每次碰撞只采样一条。

上图！

2.4.2 采样终止条件：俄罗斯轮盘赌

上面的求解还有一个关键问题，递归的终止条件是什么？如果没有递归的终止条件，这个采样会无穷进行下去。

在whitted ray tracing中，是通过直接规定弹射次数来作为终止条件，粗暴地讲，在path ray tracing中其实也可以这么做，但是这么做不是最优的，存在的问题有：①因为是随机采样，有的可能是不重要的路径，对于这些不重要的路径，其实并不需要那么多次采样，浪费计算资源；②有的采样可能是重要路径，其实需要多次采样，出现采样不足的情况。蒙特卡洛估计作为一种概率估计，很重要的就是采样，怎么样正确、高效地采样，是path ray tracing中很重要的一个方向。

所以直接固定次数不是不行，但是显然不是很好的做法。

于是乎，俄罗斯轮盘赌算法应运而生。

其实原理很简单，就是设置一个概率值P∈[0,1]P\in[0,1]P∈[0,1]，每一次反射，生成一个相关概率，如果概率小于PPP就继续传播并且Lo=L/PL_o=L/PLo=L/P，否则就停止传播。这样既能实现对光线的终止，同时估计还是无偏的。

E期望=P∗(L/P)+(1−P)∗0=LE_{期望} = P*(L/P) + (1-P)*0 = LE期望=P∗(L/P)+(1−P)∗0=L

作业框架里，PPP就直接生成一个随机数，并没有和BRDF联系起来，是简化版本。

2.5 求解效率的提高

问题到此结束了吗？哈哈，当然没有！

上面的已经是一个正确的解决方案了，但是，存在一个很大的问题，就是效率太低了，举个栗子：

实际场景中，按照上述步骤，光线反射过程中能打中光源获得亮度的很少，这样绝大部分采样的光线返回都是0，这会造成大量的浪费，和噪点。

重新拿上面的传播图示意，对于P点，之前的思路是蓝色的采样线，如果他没有打中光源，这个采样就浪费了，但是实际上，对于P点，他肯定是有光照颜色的，因为有直接光照，因此挽救这个采样线的思路也很简单：加上直接光照。

所以每次采样的时候，返回的颜色，除了采样光线，再加上直接光照。其实跟2.1小节的分析也是一样的，第一轮光线传播是直接光照，所以光线漫反射过程中，至少第一轮的直接光照应该参与漫反射。

但是方程中是对立体角ω\omegaω进行积分，光源对应的是面积，需要把面积转换为立体角。其实非常简单，因为立体角的定义就是面积除以距离的平方。

因此方程变成下面的形式，对光源积分的话，pdf=1S面积pdf = \frac{1}{S_{面积}}pdf=S面积1

在实际写代码的时候，还要加上光源是否被遮挡的判断。

3 渲染方程的代码实现

3.1 path tracing伪代码

shade(p,wo){求光线相交的物体if(物体== 光源)：return 光源颜色//直接光对光源上采样一个点x',pdf_light=1/A；if(这个点跟P点不存在遮挡)L_dir = L_i * f_r * cos θ * cos θ’ / |x’ - p|^2 / pdf_light//间接光取一个随机概率P_RRif(P_RR满足俄罗斯轮盘赌){采样一个漫反射方向wi；确定漫反射相交点q；L_indir = shade(q, wi) * fr(wi, wo, N) * dotProduct(wi, N) / pdf(wi, wo, N)/ Scene::RussianRoulette);}Return L_dir + L_indir}

3.2 作业代码

// Implementation of Path TracingVector3f Scene::castRay(const Ray& ray, int depth) const{// TO DO Implement Path Tracing Algorithm here//求入射光线wo的交点Intersection intersection = intersect(ray);Vector3f hitcolor = Vector3f(0);//判断是不是直接采样到光源了if (intersection.emit.norm() > 0)hitcolor = Vector3f(1);else if (intersection.happened){//确定入射光线的方向、相交点坐标、交点法线Vector3f wo = normalize(-ray.direction);Vector3f p = intersection.coords;Vector3f N = normalize(intersection.normal);//对光源取一个采样点，确定采样点的方向、相交点坐标、交点法线float pdf_light = 0.0f;Intersection inter;sampleLight(inter, pdf_light);//光源内随机采样一个点Vector3f x = inter.coords;Vector3f ws = normalize(x - p);Vector3f NN = normalize(inter.normal);Vector3f L_dir = Vector3f(0.0f);//求直接光照if ((intersect(Ray(p, ws)).coords - x).norm() < 0.01) //判断灯光有没有遮挡{L_dir = inter.emit * intersection.m->eval(wo, ws, N) * dotProduct(ws, N) * dotProduct(-ws, NN)/ (((x - p).norm() * (x - p).norm()) * pdf_light);}//求间接光照Vector3f L_indir = Vector3f(0);float P_RR = get_random_float();if (P_RR < Scene::RussianRoulette){//采样一个漫反射方向Vector3f wi = intersection.m->sample(wo, N);L_indir = castRay(Ray(p, wi), depth) * intersection.m->eval(wi, wo, N) * dotProduct(wi, N) / (intersection.m->pdf(wi, wo, N) * Scene::RussianRoulette);}hitcolor = L_indir + L_dir;}return hitcolor;}