闫令琪：Games101 现代计算机图形学-光线追踪(三)：渲染方程和路径追踪path ray tracing & 作业Assignment07解析

路径追踪技术（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'，分母除了和球的立体角相关，还和$r^2$相关，而$r^2$可以表示为面积。从这两个位置上的能量可以发现，这个球上任意一点的能量，和立体角、面积两个变量相关。
那么可以很自然地推出一条光线的数学表示：

这个分母上的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这篇文章里提供了两种解释，为什么定义成这两个量的比值，而不是$\frac{dL}{L}$等之值，从数学上讲，当出射方向立体角趋于0时，这个比值的极限也是0，所以没有意义，从实验的角度，入射方向的radiance却非常难测，所以采用这个比值。
公式中为什么是$d{L}_r({\omega }_r)$，因为定义里就是指定方向的入射光和出射光，所以这里的$\omega$是被微分的量。

1.4 渲染方程

渲染方程要解决的问题就是，怎么样描述看到的物体表面一个点上的光亮。这个点进入观察点的那个物理量，啰嗦一点，其实就是那个点漫反射到观察点的光线，而光线就是radiance。
由上面的BRDF公式可以很容易写出：
$$d{L}_r({\omega }_r)=f_r{L}_i({\omega }_i)cos{\theta }_{i}d{\omega }_i$$
写成积分形式就是该光线完整的表达：
$$L_r({\omega }_r)={\int }_{}^{}{f}_r{L}_i({\omega }_i)cos{\theta }_{i}d{\omega }_i$$
接下来用闫老师的课上PPT表示：

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

2 渲染方程的求解

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

对于这个漫反射来说，很自然会有疑问：绿色墙上的颜色会漫反射到立方块上，立方块上的也会漫反射到墙上，这就有点鸡生蛋、蛋生鸡的循环，这个过程到底是怎么样？我们求解渲染方程的话，终止光线弹射的条件是什么？
假如我们把时间放慢到可以观察光线传播的程度，以物体表面一个点Point点为例，观察：
①可以看到会有一根光线打到这个点，这个时候物体有一次亮度；
②所有物体表面都有一根光线被打量后，光线开始第一次弹射，Point点会被周围弹射的光线打中，叠加刚才的直接光照，会有第二次光亮；
③开始第二次弹射……
④第三次弹射等等……
闫老师的课上给了一个例子：

然而，到目前为止，跟渲染方程还没什么关系？？是的，哈哈哈……
闫老师课上通过一系列神操作，通过渲染方程推导出下面的数学表达式，即某点的光照可以分为自己发的光$E$，直接光照$KE$，二次弹射$K^{2}E$等。上面的思想实验，其实是可以通过数学推导出来，虽然我真的没听懂……[捂脸]

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

2.2 蒙特卡洛积分

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

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

$$F={\int }_a^{b}f(x)dx$$

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

$$F^N=\frac{1}{N}\sum _{i=1}^N\frac{f(X_i)}{pdf(X_i)}$$

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

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

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

pdf和pmf

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

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

$$E[X]=\sum _{x_i}pmf(x_i)x_i$$

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

$$E[X]={\int }_{-\infty }^{\infty }pdf(x)xdx$$

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

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

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

pdf
借用http://www.scratchapixel.com/的一个很好的例子来说明：

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

$$F^N=\frac{1}{N}\sum _{i=1}^N\frac{f(X_i)}{pdf(X_i)}$$

在此图中，做了四次随机采样，得到了四个随机样本$x_i$：$x_1、x_2、x_3、x_4$，并且进而得到了这四个样本的$f(x_i)$值：$f(x_1)、f(x_2)、f(x_3)、f(x_4)$。（原文没有提及如何得到$f(x_i)$。函数f是奇形怪状的，不太可能有表达式存在，难道是用尺子量的？暂且忽略这个事吧。）

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

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

$$Area=\frac{1}{4}(f(x_1)(b-a)+f(x_2)(b-a)+f(x_3)(b-a)+f(x_4)(b-a))$$

$$=\frac{1}{4}(b-a)(f(x_1)+f(x_2)+f(x_3)+f(x_4))$$

$$=\frac{1}{4}(b-a)\sum _{i=1}^{4}f(x_i)$$

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

$$F^N=\frac{1}{N}\sum _{i=1}^N\frac{f(X_i)}{pdf(X_i)}$$

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

$$Area=\frac{1}{4}\sum _{i=1}^4\frac{f(x_i)}{\frac{1}{b-a}}$$

于是可以知道$pdf(x_i)$等于：

$$pdf(x_i)=\frac{1}{b-a}$$

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

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

蒙特·卡罗积分的数学期望等于理想积分？
对于下面的$F$和$F^N$：

$$F={\int }_a^{b}f(x)dx$$
$$F^N=\frac{1}{N}\sum _{i=1}^N\frac{f(X_i)}{pdf(X_i)}$$

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

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

f(x)是离散函数时:

$$E[f(X)]=\sum _{x_i}f(x_i)pmf(x_i)$$

f(x)是连续函数时:

$$E[f(X)]={\int }_{-\infty }^{\infty }f(x)pdf(x)dx$$

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

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

$$E[\frac{f(X_i)}{pdf(X_i)}]=E[\frac{f(x)}{pdf(x)}]={\int }_{-\infty }^{\infty }\frac{f(x)}{pdf(x)}pdf(x)dx$$

$$={\int }_{-\infty }^{\infty }f(x)dx$$

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

上面的内容已经足够继续进行下面的旅程，此处贴出两个知乎上的优质贴，后续继续嗑蒙特卡洛积分。
蒙特卡洛积分 - chopper的文章 - 知乎
Monte Carlo数学原理 - papalqi的文章 - 知乎

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

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

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

这个方程理论上是可以求解的，$pdf({\omega }_i)$可以假设就是半球面上均匀采样，$f_r$可以初定一个固定的值，但是还有一个问题：

radiance该怎么表示？比如最初的灯光，它怎么表示$L_i(p,{\omega }_i)$？

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

2.4.1 采样指数发散的问题

一图胜千言，先上图。

在P点的出射光线按照上面的方程形式，需要采样N条光线作为入射光线，这N条光线，每条假如与物体交点为Q，对每一个Q点又需要采样N条光线，这样光线弹射几次，$N\ast N\ast N\dots \dots$这个计算量就爆炸了。
解决的思路也非常的简单粗暴：
$N=1$，每次只采样一条光线。
公式里的求和${\sum }_N$里的$N$不是通过对像素光线相交点采样$N$，改为对每个像素点采样$N$条光线，每条光线每次碰撞只采样一条。
上图！

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

上面的求解还有一个关键问题，递归的终止条件是什么？如果没有递归的终止条件，这个采样会无穷进行下去。
在whitted ray tracing中，是通过直接规定弹射次数来作为终止条件，粗暴地讲，在path ray tracing中其实也可以这么做，但是这么做不是最优的，存在的问题有：①因为是随机采样，有的可能是不重要的路径，对于这些不重要的路径，其实并不需要那么多次采样，浪费计算资源；②有的采样可能是重要路径，其实需要多次采样，出现采样不足的情况。蒙特卡洛估计作为一种概率估计，很重要的就是采样，怎么样正确、高效地采样，是path ray tracing中很重要的一个方向。
所以直接固定次数不是不行，但是显然不是很好的做法。
于是乎，俄罗斯轮盘赌算法应运而生。
其实原理很简单，就是设置一个概率值$P\in [0,1]$，每一次反射，生成一个相关概率，如果概率小于$P$就继续传播并且$L_o=L/P$，否则就停止传播。这样既能实现对光线的终止，同时估计还是无偏的。
$E_{期望}=P\ast (L/P)+(1-P)\ast 0=L$

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

2.5 求解效率的提高

问题到此结束了吗？哈哈， 当然没有！
上面的已经是一个正确的解决方案了，但是，存在一个很大的问题，就是效率太低了，举个栗子：

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

重新拿上面的传播图示意，对于P点，之前的思路是蓝色的采样线，如果他没有打中光源，这个采样就浪费了，但是实际上，对于P点，他肯定是有光照颜色的，因为有直接光照，因此挽救这个采样线的思路也很简单：加上直接光照。
所以每次采样的时候，返回的颜色，除了采样光线，再加上直接光照。其实跟2.1小节的分析也是一样的，第一轮光线传播是直接光照，所以光线漫反射过程中，至少第一轮的直接光照应该参与漫反射。

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

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

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

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_RR
	if(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 Tracing
Vector3f 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;
}

4 path ray tracing的优化方向

本文只是路径追踪的起点，其实读懂上面的过程就能发现，这里面可以优化的地方太多，比如
采样，闫老师课上包括作业框架只是介绍了平均采样，
还有蒙特卡洛积分数学角度更好的求解，
还有点光源怎么处理，
怎么提高速度，运算量这么大的算法怎么做到实时的？？

这篇文章，介绍了一些去除噪点的改进算法，非常棒
在路径追踪中分别采样直接与非直接光照