路径追踪（Path Tracing）与渲染方程（Render Equation）

简介

利用路径追踪我们可以实现比whitted-style ray tracing更好的全局光照（GI）效果。它的理论基础是渲染方程，最开始由吉姆卡吉亚（Jim Kajiya）提出。渲染方程是在物理的基础上定义的，利用它我们可以实现基于物理的渲染(PBR），效果会比以往的blinn-phong等模型看上去更加真实。

可以说在图形学中，要实现真实的渲染效果就是求解这个方程。那么渲染方程到底是个啥呢？公式如下：

接下来我们一点点的进行剖析，首先先从辐射度量学开始入手，即理解渲染方程中的

这一部分，也是渲染方程的物理基础。

辐射度量学（Radiometry）

在图形学中，不管是blinn-phong着色模型还是whitted-style的光线追踪，对光线的定义都进行了一定的简化，因此得到的效果并不是最真实的。因此若我们想要模拟出物体受光照后的真实效果，那么我们必然要通过物理量来精确的描述光，其中包括光的能量守恒，传播方法，与物体表面的作用等等。而辐射度量学就可以为我们提供精准的光的一系列物理量。

辐射度量学是一种用来度量电磁场辐射（包括可见光）的手段，其中有很多种辐射度量（radiometric quantities）可以用来测量表面或者方向上的光。接下来会介绍到其中的辐射能（radiant energy），辐射通量（radiant flux），辐射强度（Radiant Intensity），辐照度（Irradiance）以及辐射率（Radiance）。

注：讨论的依旧是几何光学，不考虑光的波动性，光线的碰撞等，上述与光有关的属性都是在空间上定义的。

辐射能（Radiant energy）

辐射能是指电磁波中电场能量和磁场能量的总和，这里我们理解为光线辐射出来的总能量，例如太阳就以辐射形式不断的向周围空间释放能量，这些能量能够被物体表面所吸收。

辐射能我们常用符号Q来表示，而它的单位即是能量的单位：焦耳（Joule），单位符号为J。

我们知道光源会以球的形式不断的往外辐射能量，那么一定时间内，光源所辐射出来的总能量（辐射能量）就是一个球内所有的能量。如下图：

辐射能量（Radient Flux/power）

由于我们的光总是无时无刻的产生/发射辐射能(Q)，这些辐射能会在空气中传播，被物体吸收和反射，因此我们要研究它时汪汪会选择单位时间内的辐射能。**而辐射通量值的就是单位时间内的发射、反射、传播、或吸收的辐射能。**这就好比速度指的是单位时间内运动的距离。那么我们只需要取得一个极短时间（dt）内接收到的辐射能（dQ），即可算出辐射通量的值。

辐射通量我们常用符号
来表示，根据定义我们就可得到如下公式：

前面我们说过辐射能的单位是焦耳(J)，而单位时间指的就是秒（S）。也就是说上面式子我们得到的单位正是J/s，这也就是功率的单位瓦特（瓦特的定义是1焦耳/s）。因此辐射通量的单位即是功率的单位：瓦特（Watt），单位符号为W。

在一个极短的时间内，光源在某个时刻所辐射的能量自然会集中在一个球的表面上（如下图），辐射通量指的就是在这个极短时间内球面上所有能量。

辐射强度（Radiant Intensity）

为了方便理解，我们先来看一个2D的光源，如下图：

在2D的情况下，光源会一直以一个圆的方式向外辐射能量，也就是辐射能（Q）。那么在某一个极短的时间内，这些能量都会分布在一个圆上面（如上图中的圆）。此时该圆上面的所有能量就代表着这个极端事件内光源传播的能量，也就是辐射通量。
那么此时，如果我们要求弧度为θ对应的圆弧（图中红色圆弧）上的能量是多少，应该怎么求？
首先我们要理解弧度（radian）的概念。弧度是国际单位制（官方）中对角的度量单位，它的定义为弧长除以半径，也就是说1弧度对应的弧长为r，那么整个圆的弧度即为2π。

既然整个圆的能量为辐射通量，整个圆的弧度为2π，如果我们光源均匀的向外辐射能量，那么只需要把它们相除，即，就能得到单位弧度所对应的辐射通量，而它就是2D下的辐射强度。利用它就可以求出任意弧度θ在极短时间内对应的能量，也就是任意弧度θ的辐射通量。

搞清楚2D后，我们来扩展到3D的情况，前面在介绍辐射通量的时候说过，在某一个极短的时间内，光源辐射的能量都会分布在一个球面上，例如下图：

如上图，对于球面上的任意一点P，我们都可以使用θ和φ两个弧度把它表示出来。此时若我们将θ的φ分别增加dθ和dφ，在球面上我们就会得到一块由四个点围成的面积，如下图：

其中P，P1，P2, P3四个点和光源（球心）会形成一个类似于四棱锥的角，这个角的度量单位我们称之为立体角（solid angle），常用w来表示，单位为sr（steradian，球面度）。其值为该角对应的面积（P,P1,P2,P3四个点形成的面积，标记为A）除以半径的平方，那么即可得到公式：

在2D的时候我们说辐射强度为单位弧度所对应的辐射通量，那么对应到3D下，辐射强度即为单位立体角所对应的辐射通量，用I来表示，单位为W/sr。我们可以用一个极小的立体角（微分立体角，dw）上的辐射通量算出辐射强度：

接下来我们看看立体角对应的面积怎么计算，因为我们可以考虑的是极小的一个立体角dw对应的一块极小的面积dA，那么就可以把该它看作一个长方形，面积自然就是PP1的长度乘以PP3的长度。

其中PP3的长度非常好求，它在半径为r的圆上，对应的弧度为dθ因此弧长（长度）即为dθ*r。而PP1所在的圆的半径并不是r，而是根据θ的变换而变换，即图中的PO1的长度。根据三角形的正弦定理，我们可以算出PO1的长度为 rsinθ，因此PP1的弧长为dφrsinθ，所以整个立体角对应的面积为：

从式子中我们可以发现dA除了和dθ与dφ有关外，还和sinθ有关。也就是说当dθ与dφ相同时，θ的值越小，对应的面积就越小，立体角也越小。
将dA的值带入带立体角的公式中，可得微分立体角的值:

由于微分立体角足够的小，因此我们也可以把它当做是原点在球心的一个三维空间中的方向，那么辐射强度我们也可以当做是光源在某个方向上的辐射通量。
我们设整个球面的能量（辐射通量）为Φ，球面上所有的微分立体角对应的面积积分起来自然就是整个球的面积4πr²：

注：其中θ的取值范围是0到π，φ的取值范围是0到2π。
因此整个球对应的立体角大小即为：4π。如果我们的光源均匀的往每个方向辐射能量，那么该光源的辐射强度即：

我们还可以发现辐射强度与球的半径无关，也就是说不管我们光线传播多久，辐射强度都不会衰减。

辐照度

如果此时我们使用一个面积极小的平面（dA）放在光源附近，那么该平面自然会接受到来自光源的能量，如下示意图：

此时在极短时间内整个面积接收到的能量即为dΦ(图中红色弧线)，那么单位面积接受到的辐射通量即为辐照度，常用E来表示，单位为：W/m²：

此外辐照度的计算还遵循Lambert的余弦定率，即物体表面吸收的辐射通量与光线入射角的余弦成比例。例如下图，若我们把该平面移到光源正下方，可以发现此时的dΦ会明显大于之前的，也就是说辐照度更大。

我们换个角度来看，如下图：

两束相同的光分别打在了两个不同的平面上，此时它们的辐射通量Φ是相同的，但是它们的面积并不相同，其中A=A’cosθ，即A’=A/cosθ，因此他们对应的辐照度为：

因此可以得出结论：
假设入射光线与平面的法线夹角为θ，那么：

但是通常情况下，辐照度的公式里不写cosθ，我们默认指的就是投影后的光线能量。
如果我们的点光源均匀的往各个方向辐射能量，那么球面上一点的辐照度即为：。可以看出辐照度和半径的平方成反比关系，也就是说点光源随着光线的传播，辐照度成半径平方衰减。例如上图中，我们把平面放的越来越远，其能接受到的辐射通量就越来越小，因为面积dA不变，那么辐照能也会越来越小。
因此我们会把辐照度用来定义光的强度，光的强度越大即辐照度越大，那么物体表面（dA相同）在一定时间内（dt相同）接受到的辐射能也越大，也就会越亮。对于平行光（例如阳光）而言，我们认为他不会随着距离衰减（并不需要设置光源的坐标），因此光的强度不会发生变化，但是对于点光源（例如电灯）以及聚光灯（例如手电筒）而言。光照的强度会随着物体表面到光源的距离的增大而衰减。因此通常情况下，针对点光源和聚光灯的光强度我们会设置r=1时的辐照度，即，**此时的辐照度正好与辐射强度相同。**然后根据表面距离光源的实际距离，按照半径平方对光的强度进行衰减。

上面的例子中我们只有一个光源，且不考虑间接光的情况下，那么这个极小的面积dA只会接收到这个光源的直接光，那么辐照度就是这个光源的在该面积上的辐射通量除以面积，但是如果考虑到间接光，那么dA就还可能接收到别的物体反射过来的间接光的辐射通量，辐照度的值也会随之增加。如果场景中有多个光源，那么同样的就会接收到更多的光的辐射通量，辐照度的值同样增大。

也就是说辐照度值得是一个微笑表面可以接受到的所有的光的辐射通量，这个范围我们认为是一个半球（从下面往上打的光是无效的），这些光可以是别的物体反射出来的光，也可能是其他光源的光。示意图如下 ：

辐射率（Radiance）

辐射率可以帮助我们计算微笑的平面上接受到的来自某一个方向的辐射通量，如下图：

前面理解辐射强度的时候，我们知道可以用微分立体角来代表一个方向，那么某个方向上的辐射通量我们就可以认为是某个微分立体角上的辐射通量。辐射率的值得就是单位面积下接受到的来自单位立体角的辐射通量，通常用L来表示，其单位为：W/sr▪m²。
我们假设一个微笑的平面dA接收到来自某个微分立体角dw的辐射通量为dΦ，那么辐射率即为：

注1：式子中辐射通量中的2并不是平方的意思，而是二阶导的意思，因为它和立体角与面积分别求了一次导数。
注2：公式中我们为什么要除以一个cosθ呢？因为式子中的dΦ我们认为是投影过后垂直与平面的值，也就是说这个值被乘以过cosθ，参考辐照度的介绍。因此当我们要使其变为投影前的值时（原本方向），就需要除以cosθ。

当我们的立体角足够小，那么他就可以代表一个三维方向，即式子中的w，而当我们的面积足够的小，那么它就可以额代表一个点（坐标），即式子中的p，那么我们就可以用辐射率来代表单数光线的从某个方向射向某个位置的能量了。

辐照度与辐射率的关系

我们假设一开始只有一个光源的直接光打到了p点上，如下图：

我们设该光线的辐射率为L(p，w)，那么w（单位向量）方向上对P点的贡献的辐照度dE(p,w）等于多少呢？

首先我们需要做一个投影，即乘以cosθ，其中cosθ可以写作w▪n。然后由于辐射率值得是单位立体角对应的值，而我们现在要求这个微分立体角对应的值，因此还需要再乘以一个微分立体角dw。最终得到：

由于我们此时只有一根光线，所以P点的辐照度就等于该方向提供的辐照度，即E（p） = dE(p,w).
接着如果我们在增加几个光源，或者说加上间接光（间接光本身可以当做从物体表面上发射出来的光，至于间接光的辐射率怎么算，后面再介绍），如下图：

那么此时P点的辐照度应该等于所有光线提供的
总和，即

那么最终P点的辐照度就应该是半球内所有光线提供的辐照度的一个积分，示意图如下：

公式即为（Ω代表整个半球域）：

等号右边不就是我们渲染方程里其中的一部分内容嘛，我们再来回顾一下之前的渲染方程：

可以发现加号右边的式子和我们的辐照度的计算公式非常的相似了，唯独多了一项，那么它是什么呢？这里要引入BRDF的概念了。

双向反射分布函数（BRDF,Bidirectional Reflectance Distribution Function）

我们知道之所以能够看见物体的颜色那是因为光线被反射到了人眼当中，通过前面的辐射度量学我们已经可以定义出一个点所接受到的所有来自各种光线的辐射通量，也就是辐照度。那么如果我们能够计算出有多少辐射通量被反射到人眼/camera中，不就可以计算出该点的颜色了嘛，而BRDF正是一个可以帮我们做这个计算的函数。
我们先从一个简单的只被一束光照射的例子来理解，示意图如下：

我们先来定义一下什么是反射camera的光，通过前面我们知道，P点其实就是面积为dA的一个微小平面，而camera我们也认为是一个点，因此反射到camera的光其实就是往camera方向（w0,微分立体角）的辐射通量。这不就是辐射率的定义嘛，因此我们可以将反射到camera的光定义为L（p，w0），我们要求的也就是它的值。
我们设入射光为L(p，wi)那么P点的辐照度为，那么我们只需要知道有多少比例的辐射度量会被反射到w0方向上，就可以求得L（p，wi）的值了。

如果是镜面反射，那么只有在w0▪n = wi▪n（与法线夹角相等）时，入射光会被反射到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函数就是一个双向反射分布函数（无限简化版），即告诉我们某个点上有多少被反射的辐射通量会被分布到w0方向上，同时也可以通过反射的辐射率逆推出入射的辐射率，因此叫做双向反射分布函数。

此外我们可以发现如果我们简单的修改下Function里的计算，就可以模拟出不同的材质，这里要注意能量守恒，即反射光的总量不能大于入射光的总量。

例如当wo▪n和wi▪n比较接近时才有大于0的返回值，就可以模拟出glossy的材质，如下图：

如果不管wo▪n和wi▪n的值为多少都有相同的返回值，那么就是漫反射的材质，如下图：

也就是说BRDF决定了物体的材质。
BRDF我们常用来表示，同时我要对P点接收到的所有光都要计算往camera方向的分布值，然后把它们积分起来，公式如下：

这个方程被称之为反射方程。
此时他只比渲染方程少了项，这一项的意义是：如果P点自己会发光（Emission），那么自然还要加上自发光往camera方向的辐射通量。那么反射方程在加沙昂自发光就是kajiya提出的渲染方程。

想要实现好的效果，BRDF远远没有前面所说的那么简单，如今已有很多的BRDF模型，最常用的是Cook-Torrance BRDF模型，具体介绍可以参考learnopengl里的文章。
此外P点接受到的Li（p,wi）也可能是来自其他点反射出来的间接光，即：

因此我们还应该去计算P’点往P点的反射光，当然了P’里也可能有P’'反射过来的光，因此整个计算过程是一个递归的过程。
这样我们就可以通过渲染方程来计算出camera通过某一个像素看到的极小面积点反射到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.10.2.

我们把所有的概率积分起来，自然就是PDF下面的面积，其值为1：
总结：对于任何一个连续性随机变量，我们假设其PDF为f(x)，那么它满足如下两个性质：


然后我们来看蒙特卡洛积分是啥，对于一些简单的函数，例如f(x)=x²，我们可以利用它的不定积分来求x在[0,1]之间的定积分，也就是函数下的面积，如下图：

但是对于一些复杂的函数（如下图），我们很难写出它的不定积分，那么这么求它的定积分（蓝色区域面积）呢？蒙特卡洛可以帮我们来求。

我们先在[a，b]间随机采样一个值x1，可求得f（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点往wi方向打一根光线，若没有打到光源上，那么就代表没有光照，。若打到光源上，那么说明受到了面光源对应位置的直接光照，的值也就是面光源上被打到的点往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’上的渲染方程来计算，即上面的。那么整个过程就是一个递归的过程，之前的代码我们可以再进行优化，实现全局的光照：

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根光线，即有200根光线，而且还会递归下去变得更多，这个运算量明显扛不住。由于是指数增长，即是每次只打2根光线，次数多了之后也都可能扛不住，那么难道每次只随机打1根光线么，别说，还真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方法存在的第二个问题，即我们的递归函数可能造成无限递归的问题。最简单的方法即是加一个次数的限制，但是这样弹射到限定次数后的能量被浪费掉了，不符合能量守恒。并且弹射的次数少的话，就可能造成场景偏暗，达不到预期的效果，因为真实的环境下光线就是可以无限次弹射的。那么应该怎么处理呢？针对这个问题，人们引入了俄罗斯轮盘赌（Russiann Roulette）的方法。

在很多警匪片里我们经常看到这样的画面，在左轮手枪里放入若干数量的子弹，然后对自己开枪来赌自己是否会被打死，如下图：

这个玩法就是俄罗斯轮盘赌。
左轮手枪一共可以放6颗子弹，如果放了一颗子弹，那么对自己开枪一次，生存概率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)相加，即：

很神奇吧？也就是说这样处理，最终的期望依旧是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 sequency）.
这样毅力啊就得到了一个正确的path tracing的算法了。

从光源采样

上面的算法已经正确了，但是它的效率并不高。因为shade方法中只有最终打到光源的路径/光线才能计算出对应的Lo，当我们光源越小，那么路径打到光源的概率就越低，也就是说需要往每个像素发射更多的路径才能得到比较好的效果，否则得到的效果噪声会非常的大，例如下图：

这也说明了路径追踪很难实现渲染点光源的场景，通常情况下会用一个较小的面光源来替代点光源。
那么之前的算法会造成很多路径的浪费，如下图，如果我们只考虑直接光照，那么只有紫色的路径是有效的。

如图，我们的camera会随机往像素内的某个位置发射光线，会打在场景中的某个区域上，假设交点为P。那么如果我们能够在P点的随机方向采样时限定在如图的虚线范围内的话，那就可以保证P点射出的光线能够打到光源上了。那么怎么可以实现这样的限制呢？很简单，我们可以每次在面光源上采样一个点（Q点，面积dA的区域），然后不再是随机采样一个方向射出光线，而是直接将光线射向Q点，即
这样就可以保证该光线肯定能够打到光源上，从而在不浪费任何光线的情况下计算出该方向上直接光照的结果。示意图如下：

上图为光源上采样一个区域，往该区域发射光线。
也就是说我们要把反射方程里的随机反向dwi用光源上的dA来表示，在学习立体角的时候，说立体角w的值就是球面上对应的一个面积除以半径平方。半径很好求，即为|Q-P|，那么dA怎么对应到以P为圆心，半径为|Q-P|的球面上呢？因为球面上的面积其法线方向肯定指向P，但是我们光源不一定，因此我们只需要做一个余弦变化即可。设dA的法线为n’，那么夹角的余弦值即为n’*wi，因此可以得到：

那么我们的渲染方程就可以写成如下形式：

这样我们就又可以用蒙特卡洛方法来解这个积分，可以用一个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;
}

渲染方程的一些理解

不考虑自发光的情况下，渲染方程的定义为：

假设只有一个光源，且不考虑间接光的情况下，公式为:

为反射光的辐射率(单位：W/sr•m²)，是我们要求的结果。
为入射光对表面产生的辐照度（单位为：W/m²），我们用辐照度来表示光照强度，并将其颜色化。也就是说：

则为BRDF，Microfacet Cook - Torrance BRDF的公式为：

其中F代表菲涅尔项和ks含义相同（可忽略ks），在不考虑漫反射的情况下，镜面反射光的计算公式即为:

D为法线分布函数，G为几何函数，此时公式中的所有值都将是已知的。

法线分布函数与几何函数的相关知识，以及分母的的由来推导，可参考文章，
BRDF中的法线分布函数。

，