路径追踪（Path Tracing）与渲染方程(Render Equation) agile Posted on Oct 2 2021 优秀博文 > 本文由 [简悦 SimpRead](http://ksria.com/simpread/) 转码， 原文地址 [zhuanlan.zhihu.com](https://zhuanlan.zhihu.com/p/370162390) 简介 -- 利用路径追踪我们可以实现比 Whitted-style ray tracing 更好的全局光照（Global Illumination）效果。它的理论基础则是渲染方程，最开始由吉姆 · 卡吉雅（Jim Kajiya）提出。渲染方程是在物理的基础上定义的，利用它我们可以实现基于物理的渲染（PBR，Physically Based Rendering），效果会比以往的 Blinn-Phong 等模型看上去更加真实。 可以说在图形学中，要实现真实的渲染效果就是求解这个方程。那么渲染方程到底是个啥呢？公式如下： >  what？懵逼三连。  接下来让我们来一点点的剖析它，首先我们从辐射度量学开始入手，即理解渲染方程中的  这一部分，它也是渲染方程的物理基础。 辐射度量学（Radiometry） ----------------- 在图形学中，不管是 Blinn-Phong 着色模型还是 Whitted-Style 的光线追踪，对光线的定义都进行了一定的简化，因此得到的效果并不是最真实的。因此若我们想要模拟出物体受光照后的真实效果，那么我们必然要通过**物理量**来精确的描述光，其中包括光的能量守恒，传播方法，与物体表面的作用等等。而辐射度量学就可以为我们提供精准的光的一系列物理量。 辐射度量学是一种用来度量电磁场辐射（包括可见光）的手段，其中有很多种辐射度量 (radiometric quantities) 可以用来测量表面或者方向上的光。接下来会介绍到其中的 辐射能（radiant energy），辐射通量（radiant flux），辐射强度（Radiant Intensity），辐照度（Irradiance）以及 辐射率（Radiance）。 注：讨论的依旧是几何光学，不考虑光的波动性，光线的碰撞等，上述与光有关的属性都是在空间上定义的。 ### 辐射能（radiant energy） 辐射能是指电磁波中电场能量和磁场能量的总和，这里我们理解为**光线辐射出来的总能量**，例如太阳就以辐射形式不断向周围空间释放能量，这些能量能够被物体表面所吸收。 辐射能我们常用符号 Q 来表示，而它的单位即是能量的单位：焦耳（Joule），单位符号为 J 。 我们知道光源会以球的形式不断的往外辐射能量，那么一定时间内，光源所辐射出来的总能量（辐射能）就是一个球内所有的能量，如下图：  ### 辐射通量（radiant flux/power） 由于我们的光总是无时无刻的产生 / 发射辐射能（Q），这些辐射能会在空气中传播，被物体吸收和反射，因此我们要研究它时往往会选择单位时间内的辐射能。而**辐射通量指的就是单位时间内的发射、反射、传播或吸收的辐射能**。这就好比速度指的是单位时间内运动的距离。那么我们只需要取得一个极短时间（dt）内接收到的辐射能（dQ），即可算出辐射通量的值。 辐射通量我们常用符号  来表示，根据定义我们就可得到如下公式： >  前面我们说过辐射能的单位是焦耳（J） ，而单位时间指的就是秒（s），也就是说上面式子我们得到的单位正是 J/s，这也就是功率的单位瓦特（瓦特的定义是 1 焦耳 / 秒）。因此辐射通量的单位即是功率的单位：瓦特（Watt），单位符号为 W 。 在一个极短的时间内，光源在某个时刻所辐射的能量自然会集中在一个球的表面上（如下图），辐射通量指的就是在这个极短时间内球面上的所有能量。  ### 辐射强度（Radiant Intensity） 为了方便理解，我们先来看一个 2D 的光源，如下图：  在 2D 的情况下，光源会一直以一个圆的方式向外辐射能量，也就是辐射能（Q）。那么在某一个极短的时间内，这些能量都会分布在一个圆上面（如上图中的圆）。此时该圆上面的所有能量就代表着这个极短时间内光源传播的能量，也就是辐射通量（  ）。 那么此时，如果我们要求弧度为θ对应的圆弧（图中红色圆弧）上的能量是多少，应该怎么求？ 首先我们要理解**弧度（radian）**的概念。弧度是国际单位制（官方）中对角的度量单位，它的定义为**弧长除以半径**，也就是说 1 弧度 对应的弧长为 r ，那么整个圆的弧度即为 2π 。 既然整个圆的能量为辐射通量（  ），整个圆的弧度为 2π，如果我们光源均匀的向外辐射能量，那么只需要把它们相除，即  ，就能得到**单位弧度所对应的辐射通量，而它就是 2D 下的辐射强度。**利用它就可以求出任意弧度 θ 在极短时间内对应的能量，也就是任意弧度 θ 的辐射通量。 搞清楚 2D 后，我们来扩展到 3D 的情况，前面在介绍辐射通量的时候说过，在某一个极短的时间内，光源辐射的能量都会分布在一个球面上，例如下图：  如上图，对于球面上的任意一点 P，我们都可以使用 θ 和 φ 两个弧度把它描述出来。此时若我们将 θ 和 φ 分别增加 dθ 和 dφ，在球面上我们就会得到一块由四个点围成的面积，如下图：  其中 P，P1，P2，P3 四个点和光源（球心）会形成一个类似于四棱锥的角，这个角的度量单位我们就称之为**立体角（solid angle），常用 ω 来表示，单位为 sr（steradian，球面度）。**其值为该角对应的面积（P，P1，P2，P3 四个点形成的面积，标记为 A）除以半径的平方，那么即可得到公式： >  在 2D 的时候我们说辐射强度为单位弧度所对应的辐射通量，那么对应到 3D 下，**辐射强度即为单位立体角所对应的辐射通量，用 I 来表示，单位为 W/sr。**我们可以用一个极小的立体角（**微分立体角，**dω）上的辐射通量（  ）算出辐射强度： >  接下里我们来看看立体角对应的面积怎么计算，因为我们可以考虑的是极小的一个立体角 dω 对应的一块极小的面积 dA，那么就可以把该它**看作是一个长方形**，面积自然就是 PP1 的长度乘以 PP3 的长度。 其中 PP3 的长度非常的好求，它在半径为 r 的圆上，对应的弧度为 dθ 因此弧长（长度）即为 **dθ*r** 。而 PP1 所在的圆的半径并不是 r，而是根据 θ 的变换而变换，即图中的 PO1 的长度。根据三角形的正弦定理，我们可以算出 PO1 的长度为 r*cosθ ，因此 PP1 的弧长为 **dφ*r*cosθ**。所以整个立体角对应的面积即为： > dA = (dθ*r)*(dφ*r*cosθ) 从式子中我们可以发现 dA 除了和 dθ 与 dφ 有关外，还和 cosθ 有关。也就是说当 dθ 与 dφ 相同时，θ 的值越小，对应面积就越小，立体角也越小。 将 dA 的值带入到立体角的公式中，可得微分立体角的值： >  由于微分立体角足够的小，因此我们也可以把它当做是原点在球心的一个三维空间中的方向，那么**辐射强度我们也可以当做是光源在某个方向上的辐射通量**。 我们设整个球面的能量（辐射通量）为  ，球面上所有的微分立体角对应的面积积分起来自然就是整个球的面积  ： >  注：其中θ的取值范围是 0 到π，φ的取值范围是 0 到 2π。 因此整个球对应的立体角大小即为：  ，如果我们的光源**均匀的往每个方向**辐射能量，那么该光源的辐射强度即为： >  我们还可以发现辐射强度与球的半径无关，也就是说**不管我们光线传播多久，辐射强度都不会衰减。** ### 辐照度（Irradiance） 如果此时我们使用一个面积极小的平面（dA）放在光源附近，那么该平面自然会接收到来自光源的能量，如下示意图：  此时在极短时间内整个面积接收到的能量即为  （图中红色弧线），那么**单位面积接收到的辐射通量即为辐照度，常用 E 表示，单位为**  ： >  此外辐照度的计算还**遵循 Lambert 的余弦定率，即物体表面吸收的辐射通量与光线入射角的余弦成比例**。例如下图，若我们把该平面移到光源的正下方，可以发现此时的  会明显大于之前的，也就是说辐照度更大。  我们换个角度来看，如下图：  两束相同的光分别打在了两个不同的平面上，此时它们的辐射通量（  ）是相同的，但是它们的面积并不相同，其中 A=A'cosθ，即 A'=A/cosθ ，因此它们对应的辐照度为： > 左边：  > 右边：  因此可以得出结论：假设入射光线与平面的法线夹角为θ，那么： >  但是通常情况下，辐照度的公式里不写 cosθ，我们默认指的就是投影后的光线能量。 如果我们的点光源均匀的往各个方向辐射能量，那么球面上一点的辐照能即为：  。可以看出辐照能和半径的平方成反比关系，也就是说**点光源随着光线的传播，辐照能成半径平方衰减。**例如上图中，我们把平面放得越来越远，其能接收到的辐射通量就越来越小，因为面积 dA 不变，那么辐照能也越来越小。 上面的例子中我们只有一个光源，不考虑间接光的情况下，那么这个极小的面积 dA 只会接受到这个光源的直接光，那么辐照度就是这个光源的在该面积上的辐射通量除以面积。但是如果考虑间接光，那么 dA 就还可能接收到别的物体反射过来的间接光的辐射通量，辐照度的值也会随之增加。如果场景中有多个光源，那么同样的就会接收到更多光的辐射通量，辐照度的值同样增大。 也就是说**辐照度指的是一个微小的表面可以接收到的所有的光的辐射通量，这个范围我们认为是一个半球（认为从下面往上打的光是无效的）**，这些光可能是别的物体反射出来的光，也可能是其他光源的光。示意图如下：  ### 辐射率（Radiance） 辐射率可以帮助我们计算微小的平面上接收到来自某一个方向的辐射通量，如下图：  前面理解辐射强度的时候，我们知道可以用微分立体角来代表一个方向，那么某个方向上的辐射通量我们就可以认为是某个微分立体角上的辐射通量。**辐射率的指的就是单位面积下接收到的来自单位立体角的辐射通量，**常用 L 来表示，其单位为 。 我们假设一个微小的平面 dA 接收到来自某个微分立体角 dω 的辐射通量为  ，那么辐射率即为： >  注 1：式子中**辐射通量中的 2 并不是平方的意思，而是二阶导的意思**，因为它和立体角与面积分别求了一次导数。 注 2：公式中我们为什么要**除以一个 cosθ**呢？因为式子中的  我们认为是投影过后垂直于平面的值，也就是说这个值被乘以过 cosθ，参考辐照度的介绍。因此当我们要使其变为投影前的值时（原本方向）就需要除以 cosθ。 当我们的立体角足够的小，那么它就可以代表一个三维的方向，即式子中的 ω ，而当我们的面积足够的小，那么它就可以代表一个点（坐标），即式子中的 p ，那么我们就可以**用辐射率来代表单束光线的从某个方向射向某个位置的能量了。** ### 辐照度与辐射率的关系 我们假设一开始只有一个光源的直接光打到了 p 点上，如下图：  我们设该光线的辐射率为 L(p,ω)，那么 ω（单位向量）方向上对 P 点贡献的辐照度 dE(p,ω) 等于多少呢？ 首先我们需要做一个投影，即乘以 cosθ，其中 cosθ可以写作  。然后由于辐射率指的是单位立体角对应的值，而我们现在要求这个微分立体角对应的值，因此还需要再乘以一个微分立体角 dω。最终得到： >  由于我们此时只有一根光线，所以 P 点的辐照度就等于该方向提供的辐照度，即 E(p)=dE(p,ω) 。 接着如果我们在增加几个光源，或者说加上间接光（间接光本身可以当做从物体表面发射出来的光，至于间接光的辐射率怎么算，后面再介绍），如下图：  那么此时 P 点的辐照度应该等于所有光线提供的  总和，即 >  那么最终 P 点的辐照度就应该是半球内所有光线提供的辐照度的一个积分，示意图如下：  公式即为（ 代表整个半球域）： >  等号右边不就是我们渲染方程里其中的一部分内容嘛，我们再来回顾一下之前的渲染方程： >  可以发现加号右边的式子和我们的辐照度的计算公式非常的相似了，唯独多了一项  ，那么它是什么呢？这里就要引入 BRDF 的概念了。 双向反射分布函数（BRDF，Bidirectional Reflectance Distribution Function） -------------------------------------------------------------- 我们知道之所以能够看见物体的颜色那是因为光线被反射到了人眼当中，通过前面的辐射度量学我们已经可以定义出一个点所接收到的来自各种光线的辐射通量，也就是辐照度。那么如果我们能够**计算出有多少辐射通量被反射到人眼 / Camera 中**，不就可以计算出该点的颜色了嘛，而 BRDF 正是一个可以帮我们做这个计算的函数。 我们先从一个简单的只被一束光照射的例子来理解，示意图如下：  我们先来定义一下什么是反射到 Camera 的光。通过前面我们知道，P 点其实就是面积为 dA 的一个微小平面，而 Camera 我们也认为是一个点，因此反射到 Camera 的光其实就是往 Camera 方向（  ，微分立体角）的辐射通量。这不就是辐射率的定义嘛，因此我们可以将反射到 Camera 的光定义为  ，我们要求的也就是它的值。 我们设入射光为  那么 P 点的辐照度为  ，那么我们只需要知道有多少比例的辐射通量会被反射到  方向上，就可以求得  的值了。 如果是镜面反射，那么只有在  （与法线夹角相等）时，入射光会被反射到 Camera 上，那么我们就可以写出一个函数： ``` float Function(vec3 p, vec3 in, vec3 out) { vec3 n = ...;//通过p获取到法线 if(dot(in, n) == dot(out, n)) return 1.0; return 0.0; } ``` 那么就可以得到这样的一个公式： >  其中的 Function 函数就是一个双向反射分布函数（无限简化版），即告诉我们某个点上有多少被反射的辐射通量会被分布到  方向上，同时也可以通过反射的辐射率逆推出入射的辐射率，因此叫双向反射分布函数。 此外我们可以发现如果我们简单的修改下 Function 里的计算，就可以模拟出不同的材质，这里要**注意能量守恒，即反射光的总量不能大于入射光的总量**。 例如当  和  比较接近时才有大于 0 的返回值，就可以模拟出 glossy 的材质，如下图：  如果不管  和  的值为多少都有相同的返回值，那么就是漫反射的材质，如下图：  也就是说 **BRDF 决定了物体的材质**。 BRDF 我们常用  来表示，同时我要对 P 点接收到的所有光都要计算往 Camera 方向的分布值，然后把它们积分起来，公式如下： >  这个方程式被称为**反射方程**。 此时它只比渲染方程少了  项，这一项的意义是：如果 P 点自己会发光（Emission），那么自然还要加上自发光往 Camera 方向的辐射通量。那么反射方程再加上自发光就是 Kajiya 提出的渲染方程了。 >  想要实现好的效果，BRDF 远远没有前面所说的那么简单，如今已有很多的 BRDF 模型，最常用的是 **Cook-Torrance BRDF 模型**，具体介绍可参考下面文章： [理论 - LearnOpenGL CN](https://learnopengl-cn.github.io/07%20PBR/01%20Theory/) 此外 P 点接收到的  也可能是来自其他点反射出来的间接光，即： >  因此我们还应该去计算 P'点往 P 点的反射光，当然了 P'里也可能有 P''反射过来的光，因此整个计算过程是一个递归的过程。 这样我们就可以通过渲染方程来计算出 Camera 通过某一个像素看到的极小面积点 P 反射到 Camera 上的辐射通量了，也就是该像素最终应该显示的颜色。  ### 蒙特 · 卡洛（Monte Carlo）积分 知道渲染方程后，我们就要解出它的值，知道怎么求解后我们才能实现路径追踪的代码。而解渲染方程自然主要要求出里面积分的值，这里我们可以使用蒙特卡洛的方法来求积分。 我们先从概率论的角度来说起。掷骰子大家一定都玩过，从 1 到 6，每个数字出现的概率为 1/6，小学生都懂的事情。 接着我们引入**期望**（均值）的概念，期望指的就是每次可能的结果乘以它对应的概率，然后把它们相加，常用 E 表示。例如掷骰子的期望即为： >  我们骰子掷的次数越多，所有结果的算数平均值肯定会**收敛于期望值**。 掷骰子我们只能得到 1 到 6 之间的整数（即**离散**性随机变量），此时我们改进玩法，假设可以同等概率获得 1 到 6 之间的任何数，例如 1.5，3.333333 等等（即**连续**性随机变量）。那么此时的概率应该怎么算？这里就要引入**概率密度函数**（PDF，Probability Density Function）的概念了。 对于上面例子，它的概率密度函数为： f(x)=0.2 ，示意图如下：  注意：这里和以往的函数不一样，并不是指任意 x 的取值对应的 f(x)=0.2，这样我们随便取 6 个数，它们的总概率就超过 1 了。**概率密度函数 f(x) 指的是某个确定的取值点 x 附近的一小块区域 dx 对应的面积才是其概率**。例如 x=3，取它附近 dx = 0.1 的范围，对应概率为 f(3)=0.1*0.2 。  我们把所有的概率积分起来，自然就是 PDF 下面的面积，其值为 1： 。 总结：对于任何一个连续性随机变量，我们假设其 PDF 为 f(x)，那么它满足如下两个性质： > f(x)>0 （概率不可能为负的） >  （总概率为 1）  然后我们来看蒙特卡洛积分是啥，对于一些简单的函数，例如  ，我们可以利用它的**不定积分**来求 x 在 [0,1] 之间的**定积分**，也就是函数下的面积，如下图：  但是对于一些复杂的函数（如下图），我们很难写出它的不定积分，那么怎么求它的定积分（蓝色区域面积）呢？蒙特卡洛法可以帮我们来求。  我们先在 [a,b] 间随机采样一个值 x1，可求得 f(x1) 的值，然后我们就认为 f(x1)*(b-a) 的值就是函数 f(x) 在 [a,b] 区间的定积分。如下图，认为红色区域面积就是之前蓝色区域面积。  很明显误差很大，但是如果我们再随机采样 x2，x3，...，xn 的值算出它们对应的面积 f(xi)*(b-a) 并且求个平均呢？那么当我们采样的次数越多，得到的最终平均值一定越接近定积分的值，这就是蒙特卡洛法算定积分。根据描述我们可以得到下面一个式子： >  在蒙特卡洛法里我们提到了随机采样一个值，这不就是前面所说的连续性随机变量的问题了。即会有一个 PDF（这里写作 p(x) ）来影响我们的采样情况，例如 p(x)=1/(b-a) 那么就是均匀的采样。对于上面的式子，平均采样没有问题，但是若对应 p(x) 是一个使得采样到 x1 附近的值的概率为 0.8 的 PDF，那么用脚想也知道得到的结果肯定会更加的接近 f(x1)*(b-a)。为了避免这种情况的发生我们就需要根据采样的概率添加一个**权重**的计算，总权重为 b-a ，那么采样概率对应的权重即为 p(x)*(b-a)，因此最终蒙特卡洛法的公式为： >  路径追踪 ---- 了解了上面这些知识后，我们来正式的看下怎么解渲染方程。假设我们有如下一个场景，在场景中有一个面光源（点光源的集合），然后还有个小方块，此时 P 点能够接收到来自面光源的光以及小方块反射过来的光。  为了方便，我们先假设 P 点本身不会发光，这样的话渲染方程就等价于反射方程，即： >  接下来看看怎么用到蒙特卡洛积分法来求出  （后面简写为 Lo）。由于蒙特卡洛解的是  ，那么前面的渲染方程中间的一大串就是 f(x) ，即： >  因为 P 点能够接收到半球内来自四面八方的光，那么我们就可以**在半球上的各个方向进行随机的采样。**若我们做均匀的采样，由于半球对应的立体角是 2π，因此对应的 PDF 即为：p(x)=1/2π。但是一般我们会做**重要性采样**（importance sampling），这里就简单的假设其 pdf 为 p(x)。 那么根据蒙特卡洛的公式我们就可以得到： >  为了方便，我们先只考虑直接光的效果，我们可以从 p 点往  方向打一根光线，若没有打到光源上，那么就代表没有光照，  。若打到光源上，那么说明受到了面光源对应位置的直接光照，  的值也就是面光源上被打到的点往 p 方向的辐射率。 这样就可以写出路径追踪的伪代码了： ``` //输入p点和p到camera的方向 float shade(vec3 p, vec3 wo) { float lo = 0.0//初始化Lo for(int i=0; i<count; i++) { vec3 wi = random() by pdf;//根据pdf随机获取一个方向 ray r = ray(p, wi);//从p往wi方向发射射线 if(r hit light)//射线打到了光源上 { //打到了光源上，直接光 lo += (1.0 / count) * fr(p, wi, wo) * li * dot(n, wi) / pdf(wi); } } return lo; } ``` 当然了，我们做 Path tracing 可不能仅仅只考虑直接光，接下来我们间接光也想办法加入进去。  如图，P 点可能还会受到来自 P'的间接光照，那么我们可以认为 P 点往任意方向发射射线时，某一根打到了 P'上，而 P'到 P 的  不正是 P'点往 P 方向的  么，而这个值我们又可以用 P'上的渲染方程来计算，即上面的 shade(p', p-p') 。那么整个过程就是一个递归的过程，之前的代码我们可以再进行优化，实现全局的光照： ``` float shade(vec3 p, vec3 wo) { float lo = 0.0 for(int i=0; i<count; i++) { vec3 wi = random() by pdf; ray r = ray(p, wi); if(r hit light) lo += (1.0 / count) * fr(p, wi, wo) * li * dot(n, wi) / pdf(wi); else if(r hit object at o)//打到物体上的o点 { //递归计算o往p方向（-wi）的li lo += (1.0 / count) * fr(p, wi, wo) * shade(o, -wi) * dot(n, wi) / pdf(wi); } } return lo; } ``` 上面的代码看着好像很棒了，但实际上还存在两个问题，首先由于是递归的，那么当随机打出的光线数量很多时，递归几次后光线的数量就会爆炸。例如 P 点随机 100 根光线，打到了物体表面上 20 个点，这些点又要各自产生 100 根光线，即有 20000 根光线，而且还会递归下去变得更多，这个运算量明显扛不住。由于是指数增长，即是每次只打 2 根光线，次数多了之后也都可能扛不住，那么难道要**每次只随机打一根光线**么，别说，还 tm 真是。 代码可以简写成如下： ``` float shade(vec3 p, vec3 wo) { vec3 wi = random() by pdf; ray r = ray(p, wi); if(r hit light) return fr(p, wi, wo) * li * dot(n, wi) / pdf(wi); else if(r hit object at o) return fr(p, wi, wo) * shade(o, -wi) * dot(n, wi) / pdf(wi); } ``` 毋庸置疑，这样的话，虽然计算量不会爆炸，但是误差必然会大的离谱，怎么解决呢？之前是摄像机往某个像素打 1 根光线，打到的点上产生 n 根随机方向上的光线。现在被打到的点只能随机产生 1 根光线了，但是我们可以用 n 根光线打向像素啊。此外由于像素本身具有大小，我们还往像素中随机采样到的不同位置发射光线，最后把所有光线得到的结果求平均即可，这又是一个蒙特卡洛积分。 示意图如下：  注：由于每次打出的光线不会再变成多根，因此每根光线都会形成一个路径，如图中不同颜色的路径。 此时我们还需要一个方法来往像素不同方向发射射线，伪代码如下： ``` float generationRay(vec3 camera, vec3 pixel) { float pixel_radiance = 0.0;//初始化像素受到的辐射率 for(int i=0; i<count; i++) { vec3 pos = random() by pdf;//根据pdf随机像素内的一个位置 ray r = ray(camera, pos-camera);//从相机往pos发射射线 if(r hit object at p)//如果射线打到物体上，交点为p { //利用shade计算p点camera方向的radiance pixel_radiance += (1/count) * shade(p, camera-pos); } } return pixel_radiance; } ``` 这样我们第一个问题就解决，接下里来说说 shade 方法存在的第二个问题，即我们的递归函数可能造成**无限递归**的问题。最简单的方法即是加一个次数的限制，但是这样弹射到限定次数后的能量被浪费掉了，不符合能量守恒。并且弹射的次数少的话，就可能造成场景偏暗，达不到预期的效果，因为真实的环境下光线就是可以无限次弹射的。那么应该怎么处理呢？针对这个问题，人们引入了俄罗斯轮盘赌（Russian Roulette）的方法。 在很多警匪片里我们经常能看到这样的画面，在左轮手枪里放入若干数量的子弹，然后对自己开枪来赌自己是否会被打死，如下图：  这个玩法就是俄罗斯轮盘赌，**危险动作请勿模仿**。 左轮手枪一共可以放 6 颗子弹，如果放了 1 颗子弹，那么对自己开枪一次，生存概率 p 即为 5/6，死亡概率为 1-p = 1/6。 **我们用俄罗斯轮盘赌的概念来以一定概率来结束我们的递归，并且还能够保证我们得到的结果依旧是我们期望的结果。** 在 shade 方法里，我们要求的是 Lo，如果我们递归次数少必然会导致得到的值小于 Lo。解决方法是，我们每次以 P 的概率发射光线，P 的值在 (0,1) 之间，其值为多少可以我们自己设置，比如取随机值。那么我们就会**有 P 的概率发射光线，获得到的值我们把它除以 P，即 Lo/P，会有 1-P 的概率不发射光线，获得到的值，即为 0**。 回想一下掷骰子时提到的期望的概念，每次的取值（Lo/P 和 0）乘以它们对应的概率（P 和 1-P）相加，即： > E=Lo/P*P+0*(1-P)=Lo 很神奇吧，也就是说这样处理，最终的期望依旧是 Lo，完美的解决无限递归的问题。 这样我们就可以修改一下代码，如下： ``` float shade(vec3 p, vec3 wo) { float prob = 0.6;//随便指定一个概率值 float num = random(0,1);//随机获取一个0-1的数值 if(num > prob) { //1-prob 的概率不发射光线 return 0.0; } vec3 wi = random() by pdf; ray r = ray(p, wi); if(r hit light) return fr(p, wi, wo) * li * dot(n, wi) / pdf(wi) / p;//记得要除以p else if(r hit object at o) return fr(p, wi, wo) * shade(o, -wi) * dot(n, wi) / pdf(wi) / p; } ``` 注：更好的随机数生成方法，低差异化序列（low discrepancy sequences） 这样就得到了一个正确的 Path Tracing 的算法了。 ### 从光源采样 上面的算法已经正确了，但是它的**效率并不高**。因为 shade 方法中只有最终打到光源的路径 / 光线才能计算出对应的 Lo，当我们光源越小，那么路径打到光源的概率就越低，也就是说需要往每个像素发射更多的路径才能得到比较好的效果，否则得到效果噪声会非常的大，例如下图：  这也说明了路径追踪很难实现渲染点光源的场景，通常情况下会用一个较小的面光源来替代点光源。 那么之前的算法会造成很多路径的浪费，如下图，如果我们只考虑直接光照，那么只有紫色的路径是有效的。  如图，我们的 camera 会随机往像素内的某个位置发射光线，会打在场景中的某个区域上，假设交点为 P。那么如果我们能够在 P 点的随机方向采样时限定在如图的虚线范围内的话，那就可以保证 P 点射出的光线能够打到光源上了。那么怎么可以实现这样的限制呢？很简单，我们可以**每次在面光源上采样一个点**（Q 点，面积 dA 的区域），然后不再是随机采样一个方向射出光线，而是直接将光线射向 Q 点，即  这样就可以保证该光线肯定能够打到光源上，从而在不浪费任何光线的情况下计算出该方向上直接光照的结果。示意图如下：  也就是说我们要把反射方程  里的随机方向  用光源上的 dA 来表示。在学习立体角的时候，说立体角 ω 的值就是球面上对应的一个面积除以半径平方。半径很好求，即为 |Q-P|，那么 dA 怎么对应到以 P 为圆心半径为 | Q-P | 的球面上呢？因为球面上的面积其法线肯定指向 P，但是我们光源不一定，因此我们只需要做一个余弦变化即可。设 dA 的法线为 n' 那么夹角的余弦值即为  。因此可得到： >  那么我们的渲染方程就可以写成如下形式： >  这样我们就又可以用蒙特卡洛方法来解这个积分，可以用一个 PDF（pdf_light）在光源上随机采样 dA，假设光源总面积为 A，那么均匀采样的 pdf 即为 p(x)=1/A 。而 f(x) 同样是中间一大串，可得： >  这样对于直接光的部分我们就可以直接使用采样光源的方法来处理，这样不会导致路径的浪费。而间接光的部分依旧使用原来的逻辑，即从 P 点已某个概率随机往任意方向发射一条光线，若打到非光源的物体，则进入递归。如下图，P 随机打到 O 点，然后已 O 点进入递归，对于 O 的直接光部分依旧使用采样光源的方法。  当然采样光源的方法还会有个问题，若 Q 和 P 直接有障碍物（如下图），那么 P 点肯定就无法受到来着 Q 点的直接光照，因此我们还需要从 P 点再射一条光线来判断中间是否有障碍物。  这样我们就可以修改下之前的代码了，如下： ``` float shade(vec3 p, vec3 wo) { //直接光部分 vec3 q = random() by pdf_light;//光源上随机采样一个q点 vec3 wi = normalize (q-p); float l_dir = 0.0; ray r2light = ray(p, wi); if(r2light hit light at q)//没有障碍物 l_dir = fr(p, wi, wo) * li * dot(n, wi) * dot(n', wi) / pdf_light(q) / len(q-p)^2; //间接光部分 float l_indir = 0.0; float prob = 0.6; float num = random(0,1); if(num < prob) { vec3 wi = random() by pdf; ray r = ray(p, wi); if(r hit object at o) l_indir = fr(p, wi, wo) * shade(o, -wi) * dot(n, wi) / pdf(wi) / p; } return l_dir + l_indir; } ``` 本文很多个人理解，写的不好或者不对的地方望各位大佬们帮忙纠正~~ 【Unity笔记】ShaderLab与其底层原理浅谈 “《永劫无间》的动作与运动系统”笔记