(中文详解篇)smallpt: 99行代码完成全局光照Path Tracing

0. 什么是SmallPT
- Features
1. 光线追踪需要了解知识
- 1.1 什么是全局光照？
- 1.2 渲染方程
2. SmallPT代码分析
- 2.1 代码块1
- 2.2 代码块2
- - 2.2.1 主函数第1步: 设置呈像平面
  - 2.2.2 主函数第2步: 设置相机位置
  - 2.2.3 主函数第3步: 创建图像
  - 2.2.4 如何计算radiance？
- 2.3 shading 着色
- - 2.3.1 diffuse reflection 漫反射
  - 2.3.2 Ideal Specular (Mirror) Reflection 高光反射
  - 2.3.3 玻璃材质(Dielectric)
4. 渲染效果
参考文献

0. 什么是SmallPT

smallpt是一个全局光照渲染器(global illumination renderer). 其核心代码是99行C++代码组成，渲染上述场景是使用unbiased Monte Carlo path tracing (无偏的蒙特卡洛path tracing算法)。

本次博客的目的是介绍smallPT的原理，以及代码对应的意义（主要参考奥克拉荷马州立大学David Cline教授的ppt），希望能够通过本文，读者能够理解代码层面的光线追踪算法：Path tracing，并对整个渲染流程有一个清晰的认识。

Features

基于unbiased Monte Carlo path tracing的全局光照（开源C++代码）
使用OpenMP进行多线程加速
使用漫反射光照产生软阴影(Soft shadows)
高光反射，漫反射以及玻璃的BRDFs（渲染流程中的材质属性）
基于重要性采样的超采样的antialiasing方法
光线和球面相交(Ray-sphere intersection)
余弦重要性采样(Cosine importance sampling of the hemisphere for diffuse reflection)
俄罗斯轮盘赌(用于光线停止)

1. 光线追踪需要了解知识

1.1 什么是全局光照？

全局光照就是相当于对场景进行一个虚拟的照相操作。

下图是特定光照(Ad-hoc Lighting)与全局光照(Global Illumination)渲染出来图像的区别，可以看出，全局光照的龙更加亮，而且前后的阴影也更加和谐。

那么，我们如何才能获得一张全局光照渲染的结果呢？
答案是”逆向工程“：即从渲染出的图像的每个像素位置向透镜发射光线，根据一系列的折射和反射，最终达到光源上的一点。这叫做光线追踪（最典型的代表就是Path Tracing）。

1.2 渲染方程

什么是渲染方程？渲染方程是长发飘飘的吉姆·卡吉雅于1986年发明的方程，从此以后，这个方程直接影响了影视，动画，游戏，视觉特效等方方面面的行业。

这里是只考虑点光源的渲染方程的例子(不考虑光的弹射, bounce), 那么渲染方程很简单:

其中，等式左边的L(x→wr)L(x \rightarrow w_r)L(x→wr)是从表面上一点x在方向wrw_rwr的radiance，这是我们要计算的内容。而等式右边分别是材质自发光项(ambient 或者Emission)和入射光与材质(BRDF)的入射面积。

但是对于有interreflection的情况, 求解渲染方程就遇到困难了…, 因为对各个面反射出去的radiance不太清楚:

咋整呢? 首先先简写这个渲染方程, 便于后面的骚操作(u是出射, v是入射, 红色表示unknown):

I(u) = e(u) + fI(v)K(u,v)dv

还是不够简化, 我们把u和v先省略, 再简化为红色框的内容: 一个简单的矩阵方程。可以看出, L是递归(recursive)的.

接着, 根据L去找一个解析解, 可以看出根据Taylor级数展开, 我们可以看出L的形式:
根据这种形式, 很容易代入物理意义: 即哪项代表的是自发光, 哪项代表的第几次bounce(光线弹射几次)的伪光源.
这样可以看出, 反射也是一个间接光照.

通过这个式子, 老闫介绍: 光栅化直接可做的只有E和KE.
即光栅化不能很容易的模拟间接光照的影响. 这很扎心…

所以对要求更高的效果, 还是使用Ray Tracing，所以引出了本文的全局光照~

好了，说的有点太多了，那么我还是回到本文关注的Path Tracing的逻辑流程本身，我们的叙述还是以这张图为例进行展开:

其流程如下:

① 对每个要计算颜色的像素点，设置samplesPerPixel个采样点(这个值越大，渲染的效果越好, 但是也会越慢)。
② 随机选择透镜上的一点PlensP_{lens}Plens和成像平面image plane上的一点PimageP_{image}Pimage，构造射线ray(从透镜上的随机点为原点出发，方向为normalize(PimageP_{image}Pimage - PlensP_{lens}Plens)。
③ 发射射线，如果射线直接接触到某个物体上（如上图的scatter point），那么就由这个scatter point作为新的射线的出发点，继续构造新的射线ray，直到其回到光源上的点。若没有接触到物体上，则使用背景颜色。
④ 最后，将C除以每个像素的采样点（求平均），即得到对应像素点渲染的结果。

2. SmallPT代码分析

原代码是Kevin Beason写的99行代码，但是有很多人抱怨说写的过于精简，难以看懂。因此，Oklahoma State University的David Cline教授将其扩展为218行，并将其主要部分列出，使用面向对象（OOP）的方式来完成SmallPT的易读版本~

2.1 代码块1

其作用是包含了常用的结构体，如Vec, Ray, Sphere (smallPT只包含球体渲染)，以及一些功能函数和球体的初始化。

以球体的初始化这步可以看出，其传入的5个参数按顺序分别是:

① 球体半径（radius): 有1e5，16.5和600.
② 球体的位置 (position): 渲染结果里面的2个球应该是对应半径为16.5的，位置分别为Vec(27, 16.5, 47)和Vec(73, 16.5, 78).
③ emission: 默认都是空，也就是环境光项都为0。
④ 颜色: 球体的颜色。
⑤ 材质: 一共有3种, DIFF, SPEC和REFR。

Vec的结构体的详细介绍，其目的是为POINTS， COLORS以及VECTORS提供基础的结构。

其中norm()的意义是为了求射线的方向，点积是为了求余弦角，叉积是为了求正弦角。

下面介绍射线的结构体，其形式如: P(t)=O+tDP(t) = O + tDP(t)=O+tD

再下面是球体结构, 通过球体的中心和半径，即可定义一个球，其方程和Vector形式都在下面给出。

判断球面和射线相交的解析解形式：

判断相交的代码如下，其中:

double b = op.dot(r.d) 就对应 b=2D⋅(O−C)b = 2D · (O - C)b=2D⋅(O−C)，DDD 就是代码里的r.d, OOO是射线的原点r.o，CCC的球面的中心p. 这里的b因此是0.5b。
det = (b2−4ac)4\frac{(b^2 - 4ac)}{4}4(b2−4ac), 因为a = (D⋅D)(D · D)(D⋅D), 其模为1，所以a=1. det = b2−cb^2 - cb2−c。
最后返回一个最小的正t，表示与射线相交的最近的平面上的点。

场景图如下所示，如果渲染正确，应该得到如下结果。

2.2 代码块2

其作用是计算radiance和主函数执行。
其中，主函数的作用包括：设置camera位置，对每个像素做超采样抗锯齿，并行指令等操作。

2.2.1 主函数第1步: 设置呈像平面

这里可以看出，是构造了一个高为384, 宽为512的image plane用于显示渲染出来的图像。显然，图像越小，渲染的计算代价也就越小，计算的也就越快。

2.2.2 主函数第2步: 设置相机位置

相机的位置和方向都是非常重要的，如果设置有误，很容易导致渲染出全黑或者部分可见的效果。
相机的设置包括cam(相机的位置和方向)，相机的水平方向和竖直方向。0.5135表示的是field of view的角度。

2.2.3 主函数第3步: 创建图像

这一步其实就是之前的Path Tracing的伪代码的具体实现，不过相较于伪代码版本，本版本加入了超采样antialiasing(抗锯齿)和tent filter(阴影贴图).

① 首先是OpenMP并行加速，这个的作用是可以并行执行for循环，因为每个像素的渲染是互不干扰的，所以可以使用OpenMP进行并行加速。

② 遍历image plane上的所有像素点

③ subpixel，用于抗锯齿，注意看c[i] = c[i] + Vec(…) * 0.25, 的0.25~

④ 像素索引 Pixel Index

计算像素点i的索引位置，注意为啥是i=(h−y−1)∗w+xi = (h-y-1) * w + xi=(h−y−1)∗w+x呢？这是因为OPENGL设备坐标系是左手坐标系，屏幕坐标系原点在左下角，向上向右增加。

⑤ Tent Filter

Tent Filter是啥？有老哥问了，很好~[3], 来自Nvidia的大佬Nathan Reed回答了这个问题:

sinc filter的波形图

”理论上最好的antialiasing方法是sinc filter[4]，因为sinc filter可以完美的去除所有高于Nyquist frequency的频率，并保留lower ones，所以，我们的目标是尽可能的接近sinc filter的波形，以期达到完美的antialiasing效果。“

至于为啥smallPT用了tent filter，Nathan老哥认为: tent filter代码少，几行就可以搞定~，而bicubic filter需要多行代码。

⑥ 使用cam.d, cx, cy来计算光线的方向，计算radiance，并进行subpixel估计的求和。

2.2.4 如何计算radiance？

① 进行相交判断，若不相交直接返回；相交的话取出最近的相交球体obj.

其中变量Xi是用来存储随机数的（由erand48()生成），is seeded using an arbitrary function of the row number to decorrelate (at least visually) the sequences from row-to-row.

② 获取表面的法线，颜色(BRDF modulator)等信息

注意看nl的方向问题，对于玻璃表面，ray tracer必须确定其是否是enter or exit glass，并以此来计算折射光(refraction ray)。我的理解是，需要根据射线的方向确定法线的方向，这对于某些材质（玻璃）的物体渲染非常重要。

③ 俄罗斯轮盘赌：用于终止recursive的radiance计算

2.3 shading 着色

2.3.1 diffuse reflection 漫反射

这里是随机构造角度r1和距离中心的距离r2，来计算标准正交的坐标(w, u, v).

接着是从计算表面的漫反射（光线的bounce，即实际世界中，物体并不一定只是被光源点亮，也会被其它物体发出的光影响），这里不解释具体的公式了，需要补的孩子们建议去看看冯乐乐小姐姐写的Unity的Shader那本书。

2.3.2 Ideal Specular (Mirror) Reflection 高光反射

这里用的是最简单的理想高光反射，即入射角==出射角。

2.3.3 玻璃材质(Dielectric)

这里David Cline教授分析了一下这个代码，首先，由于玻璃本身既可以反射(reflective)又可以折射(refractive)，所以这里计算是反射光线why?

并为glass设置nt为1.5. nnt 为1.5或者1/1.5.
（nt或者nnt是啥? 下面的内容解释了，1.5是玻璃的折射系数~）

当出现异常情况：比如light ray想要离开glass，但是其出射角度过小(shallow angle), 则将所有光反射。

使用菲涅尔项来计算计算折射光。

T是折射光的方向。

菲涅尔项是什么？

是基于入射角(θa\theta_aθa)的玻璃表面的光的反射/折射的百分比。

在法线入射方向上的反射比是Fo=(n−1)2(n+1)2F_o = \frac{(n-1)^2}{(n+1)^2}Fo=(n+1)2(n−1)2, 其它角度上的反射毕为Fr(θ)=Fo+(1−Fo)(1−cos⁡θ)5F_r(\theta) = F_o + (1 - F_o)(1 - \cos\theta)^5Fr(θ)=Fo+(1−Fo)(1−cosθ)5

4. 渲染效果

上面的每项的含义其实并没有完全解释清楚，有需要的同学还是得一行行debug来看[1, 2]. 我们可以看到，随着每个像素点采样数目的增加，path tracing的效果越好，显然，渲染效率也变得非常的底下了~

参考文献

David Cline教授介绍smallPT的slide
smallPT–kevin beason
Why use a tent filter in path tracing?—stackexchange
Sinc Filter