光线追踪Ray Tracing

Ray Tracing

之前的Blinn-Phone模型是用来计算直接光源对物体的作用，从而可以使我们看到物体所反射的光，但是，在某种情况下，不仅仅只有直接光源，可能物体所反射的光会会在一个空间中进行多次反射或者折射才可以进入到我们的眼睛当中，而blinn-phong模型就不再适合了，所以就有了更加全局光照的概念

以上模型，比如软阴影和间接光源，都是通过折射或者反射来计算的，运用全局光照模型才能取得更好的效果

光线追踪的原理

该图就很好的解释了该原理，虽然这个在物理上是不对的，因为我们可以看见光线是因为光线在物体上进行反射的光进入了我们的眼睛，我们眼睛是不会发光的，但是，我们知道，光路是可逆的，也就是说，光线路径可逆，从A发出的到B的光线，一定也可以从B发出到A（中途可发生反射和折射），这样可以很方便让我们计算我们是否可以看到某个物体：

也就是只需要从我们眼睛发出光，看这束光经过反射和折射之后是否能够到达光源，如果能，那么我们就可以认为这个物体可以被我们看到，如果不能，那么我们就可以认为这个物体不可以被我们看到，这就是大致光线追踪的原理（仅仅是大概的意思）

First：Casting Ray 眼睛发光

从人眼或摄像机向近投影平面上的每一个像素点发射一条光线，判断与场景物体的交点

当然一条光线自然可能会与不止一个物体相交，但是考虑遮挡关系，只去找最近的交点。接着连接该交点和光源，只要判断这条连线之间是否有物体存在就可以知道该交点是否在阴影之中，接下来就是计算颜色以及光照强度了，利用Blinn-Phong模型对这个点进行局部光照模型计算，得到该像素的颜色，那么遍历所有近投影平面上的像素就能得到一张完整的图像

但如果光线追踪仅仅是在第一步Ray Casting就停止的话，那么它的效果与局部光照模型是一样的，因此我们需要第二步，真正的考虑全局效果

Second：Recursive Ray Tracing 全局光照追踪 （Whitted-Style）

在第一步中所做的Ray Casting，该条光线第一个与圆球物体相交，假设该圆球是一个玻璃球，那么便会发生镜面反射

当然除了镜面反射之外，自然也存在折射，同时反射与折射出去的光线会可能与场景中的物体再次碰撞，发生第二次折射与反射，发生第三次第四次等等等等

从图中可以见到，不仅仅是与圆球相交的那一点可以贡献光到达眼睛，折射与反射之后再与物体相交的点也可以贡献光(光路可逆原理)。简而言之，除了直接从光源照射到圆球交点再沿着 eye rays(第一条发射的光线)到眼睛中，也可能存在这样一种情形，有光照射到其他物体，再沿着eye rays的反射或折射的光线方向传回人眼！

所以，每个交点的颜色都是来自这几种类型光的总和，即为：直接光源，反射间接光，折射间接光等

下一步将这些所有交点与光源连接，称这些线为shadow rays（用来检测阴影），计算这些所有点的局部光照模型的结果，将其按照光线能量权重累加，最终得到近投影平面上该像素点的颜色！而这就是一个全局光照模型了，因为不仅仅考虑了直接光源的贡献，还考虑各种折射与反射光线的贡献。

由上述过程我们也可以得出，本身这个模型要计算各种反射和折射光线，而现实情况中，这些光线是会被无数次的反射和折射，那么我们不可能去模拟计算无数次的结果，因此有：

• 一定的递归终止条件，比如说允许的最大反射或折射次数为10。

• 光线在每次反射和折射之后都有能量损耗的，由系数决定，因此越往后的折射和反射光贡献的能量越小，最后，直接为0时也会退出递归

Ray-Surface Intersection 射线曲面相交

Ray Equation 表示一条光线

由原理我们知道，我们要从某一个像素点向外发射一条光线，那么，我们该如何用数学来表达该光线？---- Ray Equation 即为：将光线表示为条射线，那么，在给定射线的方向，那么我们就可以唯一的确定一条光线了，该条光线将会有两个属性，一个是出射点的坐标，一个是光线的方向向量
比如：

其中，加上t是为了表示光线行进的距离方便后续的判断

Ray Intersection With Implicit Surface 射线与隐式曲面相交

$$\begin{gathered} Ray:\mathbf{r}(t)=\mathbf{o}+t\mathbf{d},0\le t<\infty \\ Sphere:\mathbf{p}:(\mathbf{p}-\mathbf{c}{)}^2-R^2=0 \end{gathered}$$

联立上述方程即可求得相交时间$t$，以及相交点的坐标

这里的$p$是球面的一点，$o$是射线的出射点，$d$是射线的方向向量，$R$是球面的半径。

由此，我们可以推广出更一般的式子，就是和任意曲线的交点和时间

Ray Intersection With Triangle Mesh 射线与三角网格相交

真正在图形学中大量运用的其实是显示曲面，更具体来说就是许许多多个三角形，因此如何判断一条光线与显示曲面的交点，其实也就是计算光线与三角形面的交点。

我们知道，如何去确定一个平面，也就是首先确定平面的法向量，然后再确定一个点的坐标，那么通过这个点的平面就是我们想要的平面,如上图所示,我们确定了平面的方程

然后,直接利用平面方程和射线方程求时间,然后在求交点然后判断交点是否在我们想要判断的三角形内即可,如下图所示:

进阶做法,直接将点的形式用重心坐标的形式表示，随后利用克莱姆法则求解线性方程组即可!

这个理解意思即可,运用的要记住公式,到时候直接看笔记就可

Accelerating Ray-Surface Intersection

Bounding Volumes 包围盒(AABB)

由上可知,我们要计算每个从一个像素点出射的光线是否和曲面有交点,但是,如果一个画面当中,有物体的部分仅仅只占画面的一小部分,那么,这样做显然会浪费很多时间,因此,我们可以将物体用一个包围体将其包裹起来,然后,只需要计算那部分会穿过包围体的光线即可

也就是说:有的光线显然不会与一个物体相交的时候，那么自然也没有必要去遍历该物体的所有三角形面，因此利用一个包围盒包住该物体，在与该物体的三角面计算求交之前先判断光线是否与包围盒相交，倘若连包围盒都与光线没有交点的话，那么显然不会与物体的三角面有交点。

也就有:AABB 轴对齐包围盒,由三对平面的交集构成,AABB的任意一对平面都与x-axis，y-axis或者z-axis垂直，所以称之为轴对齐包围盒

例如:

实现光线与包围盒求交

以2D为例,3D进行类推即可

如上图最左边所示，求出光线与x平面的交点，将先进入的交点(偏小的那个)记为 $t_{min}$, 后出去的交点(偏大的那个)记为 $t_{max}$，紧接着如中间图所示计算出光线与y平面的两个交点同样记为另外一组$t_{min}$ , $t_{max}$,当然计算的过程中要注意如果任意的 t < 0，那么这代表的是光线反向传播与对应平面的交点。求得时间之后,我们就可以进行判断是否光线与盒子有交点:我们可以知道,只有当光线进入了所有的平面才算是真正进入了盒子中,只要当光线离开了任意平面就算是真正离开了盒子,所以由下面的步骤:

求出进入盒子的最大时间$t_{enter}$,求出射出盒子的最小时间$t_{exit}$,然后可以判断,当射入时间小于射出时间时,则说明光线会在盒子里待一段时间,也就是会有交点出现

注意:经过求解的时间可能会出现负值,那么,又该如何判断呢,出现这个原因是因为我们计算的时候,将射线当成直线来计算,就很有可能出现负值,具体由以下几种情况:

1. if $t_{exit}$ < 0:
   • The box is “behind” the ray —no intersection!
2. if t_{exit} >= 0 and t_{enter} < 0:
   • The ray’s origin is inside the box — have intersection!
3. In summary, ray and AABB intersect iff(当且仅当) $t_{enter}<t_{exit}\&\&t_{exit}>=0$

虽然,经过AABB的方法可以使得效率有了不少的提升,那么考虑以下几种情况:

1. 整个场景只有一个极其复杂的单一人物模型，那么只对这一个物体做包围盒的话，相当于对效率没有任何提升
2. 整个场景充斥着大量的细小模型，如草，花之类的，每个模型可能只有很少的面，如果此时对每个物体求包围盒，得到的包围盒数量会相当之多，对于光线追踪效率来说效率提升有限

基于以上两点考虑，AABB并不应只局限于以物体模型为单位，可以更加精细的考虑到以三角面为单位。因此如何更好的划分场景形成不同的AABB，使得划分之后的AABB能够更好的加速光线追踪，这就是接下来要考虑的问题关键！

Uniform Spatial Partitions (Grids) 均匀空间划分

首先,就是最简单的均匀空间划分,具体步骤为:

1. 在要求的场景中,找到一个包围盒
   

2. 在这个包围盒中,按照一定的规则,将场景划分成一定的网格
   

3. 在每个重叠的包围盒上存储物体模型的信息
   

   关于此处对于重叠的理解,应该是以物体为单位建立包围盒,然后再将此包围盒和划分的包围盒相交,从而得出有物体的包围盒的位置,并记录下来

4. 根据光线的方向与判断出所有相交的方格，倘若方格中存储有物体，再进一步与方格中的物体模型或是三角形面求交。
   

此种方法的一个缺点考虑以下情况:

包围盒数量太少,会使得准确率下降,使得白白浪费掉判断的机会,如果太多,计算一条光线穿过的包围盒的数量也会变多,从而也是浪费时间,因此划分方格数的多少也是性能的关键，方格太少，没有加速效果，方格太多，判断与方格的求交可能会拖累效率）

KD-Tree空间划分

一些常用的空间划分方法:

• 第一种Oct-Tree，也就是八叉树:

  每次将空间分为8个相等的部分，再递归的对子空间进行划分，因为图中是2维例子，所以只划分了4部分。当划分的子空间足够小或是空间中三角形面的数量很少的时候会停止划分。这种方法的显著缺点是，随着维度的上升划分的空间数量会呈指数级增长。

• 第二种KD-Tree

  其 每次将空间划分为两部分，且划分依次沿着x-axis，y-axis，z-axis ，即如图中所示，第一次横着将2维空间分为上下，第二次再竖着将上下两个子空间分别划分为左右部分，依次递归划分，终止条件与八叉树类似

  

  上图中的动态过程就很明显的展示了动态划分的过程,每次都要沿着一个维度的坐标轴去进行划分,直到划分的空间太小或者是子空间内只有少量的三角形面的时候停止划分。

  

  注意！因图中做了简化，最大包围盒的左半边并没继续进行划分(实际上应该要划分的)，所以左半部分对应的1号空间是叶子节点，如果光线与之相交，进一步判断与存储与叶子节点的物体信息求交。左半边判断完之后，接着判断右半部分,同样如果对于有半部分存在相交情况，则对于右半部分的所有子空间，递归的执行这个步骤即可

  优点:
  利用KD-Tree的结构来构建AABB的好处是倘若光线与哪一部分空间不相交，那么则可以省略该部分空间所有子空间的判断过程，在原始的纯粹的AABB之上更进一步提升了加速效率。
  缺点：
  缺点是判断包围盒与三角面的是否相交较难，因此划分的过程不是那么想象的简单，其次同一个三角面可能被不同的包围盒同时占有，这两个不同包围盒内的叶节点会同时存储这一个三角形面

• 第三种BSP-Tree

  其与KD-Tree类似，唯一不同的是划分不再沿着固定一轴，可以任意方向划分，缺点自然是划分的空间没有规则性，求交困难。

Bounding Volume Hierarchy

BVH与前几种方法最显著的区别就是，不再以空间作为划分依据，而是从对象的角度考虑，即三角形面

注意:包围盒会重叠，但一个三角形面只会被存储在唯一的包围盒内，而这也就解决了KD-Tree的缺点！

每次划分一般选择最长的那一轴划分，假设是x轴，那么划分点选择所有三角面的重心坐标在x坐标上的中位数进行划分，如此便能保证划分的三角形左右两边三角形数量尽可能差不多，当然也就使得树形结构建立的更加平衡，深度更小，平均搜索次数更少，提高了效率

中间节点不存储物体三角面信息，只在叶节点中存储，终止条件可设定为当前包围盒内三角形数量足够少

辐射度量学

辐射度量学其实是对光照的一套测量系统和单位，它能够准确的描述光线的物理性质。

有关概念解释

1. 辐射能量(Radiant energy)

   
   简单来讲就是辐射出来的电磁能量，单位为焦耳

2. 辐射通量 (Radiant flux)

   
   所谓辐射通量或者说辐射功率，其实就是在辐射能量的基础之上除以时间，也就是单位时间的能量。可以用物理当中的功率来进行类比。

3. 辐射强度(Radiant intensity)
   Radiant itensity其实就是指从一个光源出发某一方向上的亮度

   
   Radiant intensity一句话来说就是从光源发出的每单位立体角上的功率 有关立体角的概念如下：
   
   Radiant intensity的物理含义就为光源向某一方向所发射出的单位立体角的功率，简而言之就是光源在某个方向上的亮度如何

4. irradiance

   
   irradiance是指每单位照射面积所接收到的power

5. radiance

   
   radiance是指每单位立体角，每单位垂直面积的功率 直观来看的话，很像是Intensity和irradiance的结合。它同时指定了光的方向与照射到的表面所接受到的亮度。

   注：在irradiance中定义的每单位照射面积，而在radiance当中，为了更好的使其成为描述一条光线传播中的亮度，且在传播过程当中大小不随方向改变，所以在定义中关于接收面积的部分是每单位垂直面积

   

总结关系：
$$L(\mathrm{p},\omega )=\frac{\mathrm{d}E(\mathrm{p})}{\mathrm{d}\omega \cos \theta }$$

观察一下积分后的式子， $E(p)$ 就是点的irradiance，其物理含义是上文所提到过的点p上每单位照射面积的功率，而 $L_i(p,\omega )$ 指入射光每立体角，每垂直面积的功率，因此积分式子右边的 $\cos \theta$ 解释了面积上定义的差异，而对 $d\omega$ 积分，则 是相当于对所有不同角度的入射光线做一个求和，那么该积分式子的物理含义便是，一个点(微分面积元)所接收到的亮度 (irradiance)，由所有不同方向的入射光线亮度(radiance)共同贡献得到。

双向反射分布函数(BRDF)

一个点(微分面积元)在接受到一定方向上的亮度 $\left(dE\left({\omega }_i\right)\right)$ 之后，再向不同方向把能量辐射出去 $\left(d{L}_r\left({\omega }_r\right)\right)$

直观的理解，不同物体表面材质自然会把一定方向上的入射亮度 $\left(dE\left({\omega }_i\right)\right)$ 反射到不同的方向的光线上 $\left(d{L}_r\left({\omega }_r\right)\right)$ ，如理想 光滑表面会把入射光线完全反射到镜面反射方向，其它方向则完全没有。如理想粗粘表面会把入射光线均匀的反射到所有方向。 因此所谓BRDF就是描述这样一个从不同方向入射之后，反射光线分布情况的函数， 定义如下:

借助BRDF，可以定义出反射方程如下：

所以，摄像机所接受到的 ${\omega }_r$ 方向上的反射光，是由所有不同方向上入射光线的irradiance贡献得到的(图中式子的 $L_i\left(\mathrm{p},{\omega }_i\right)\cos {\theta }_i\mathrm{d}{\omega }_i)$ ，而不同方向入射光线的irradiance对反射方向 ${\omega }_r$ 的贡献程度则由物体表面材质决定，所以乘上了 一个BRDF函数。

渲染方程

渲染方程知识在反射方程的基础之上添加了一个自发光项(Emission term)

$$L_o\left(p,{\omega }_o\right)=L_e\left(p,{\omega }_o\right)+{\int }_{{\Omega }^{+}}{L}_i\left(p,{\omega }_i\right)f_r\left(p,{\omega }_i,{\omega }_o\right)\left(n\cdot {\omega }_i\right)\mathrm{d}{\omega }_i$$

其中 $L_e\left(p,{\omega }_o\right)$ 为自发光项，反射方程中的 $\cos \theta$ 用， $n\cdot {\omega }_i$ 代替。

解渲染方程就暂且不涉及了，只需要知道渲染方程的基础知识就可以了。

现在直接给出最后特别容易理解的解：

E为光源发出的光，KE则代表对光源反射一次的结果，即直接光照，那么前两项之和就是光栅化当中着色所考虑的结果
对于全局光照来说，还考虑了 $K^{2}E$，即一次弹射的间接照明， $K^{3}E$ 就是两次弹射的间接照明，依次类推。
所以物体的展示效果就是：光源发光加上直接光照与多次间接光照的结果!

结果：

蒙特卡洛路径追踪(Monte Carlo Path Tracing)

如何通过该渲染方程计算出屏幕坐标上每一个坐标的像素值？ -利用蒙特卡洛路径追踪

蒙特卡洛积分(Monte Carlo Integration)

蒙特卡洛积分的目的： 当一个积分很难通过解析的方式得到答案的时候可以通过蒙特卡洛的方式近似得到积分结果

显然对于这样一个函数，很难去用一个数学式子去表示，因此无法用一般解析的方法直接求得积分值，而这时候就可以采用蒙特卡洛的思想了。
蒙特卡洛积分的原理及做法:对函数值进行多次采样求均值作为积分值的近似

该做法十分容易理解，想象一下如果对上图这个函数值进行均匀采样的话，其实就相当于将整个积分面积切成了许许多多个长方形，然后将这些小长方形的面积全部加起来。没错，该做法其实就与黎曼积分的想法几乎一致。但蒙特卡洛积分更加的general，因为它可以指定一个分布来对被积分的值进行采样，定义如下：

如果都使用均匀采样的话，就有以下公式：

也就是说，为什么要引入这个呢，就是要因为我们之前得到的渲染方程很难直接通过正常的积分算出原函数来进行计算，所以才引入这个方法来实现近似计算

蒙特卡洛路径追踪(Monte Carlo Path Tracing)

我们的渲染方程如下：

$$L_o\left(p,{\omega }_o\right)=L_e\left(p,{\omega }_o\right)+{\int }_{{\Omega }^{+}}{L}_i\left(p,{\omega }_i\right)f_r\left(p,{\omega }_i,{\omega }_o\right)\left(n\cdot {\omega }_i\right)\mathrm{d}{\omega }_i$$

在进入具体计算之前，对渲染方程做出一点小修改，即舍弃一下自发光项(因为除了光源其他物体不会发光), 以方便进行计算推导

$$L_o\left(p,{\omega }_o\right)={\int }_{{\Omega }^{+}}{L}_i\left(p,{\omega }_i\right)f_r\left(p,{\omega }_i,{\omega }_o\right)\left(n\cdot {\omega }_i\right)\mathrm{d}{\omega }_i$$

回想第一章所提的，对于一个困难积分只要选定一个被积分变量的采样分布即可通过蒙特卡洛的方法得到积分结果的近似值，而此时的被积分值为 ${\omega }_i$ ，选定${\omega }_i\sim p\left({\omega }_i\right)$ ，不难得出积分近似结果如下:

$$L_o\left(p,{\omega }_o\right)\approx \frac{1}{N}\sum _{i=1}^N\frac{L_i\left(p,{\omega }_i\right)f_r\left(p,{\omega }_i,{\omega }_o\right)\left(n\cdot {\omega }_i\right)}{p\left({\omega }_i\right)}$$

先单独考虑直接光照，因此只有当采样的方向${\omega }_i$击中光源的时候，光源才会对该着色点有贡献，计算伪代如下：

显而易见的，单独仅仅考虑直接光照自然是不够的，还需要间接光照，即当采样的方向${\omega }_i$碰撞到了别的物体，如下图所示：

此时采样的光线碰撞到了另一个物体的Q点，那么该条路径对着色点P的贡献是多少呢？自然是在点Q的直接光照再乘上反射到该方向上的百分比了！

注：每次采样仅仅只能采样一个方向，防止计算递归而导致计算爆炸，如果选择多个方向，则会造成以下后果：

所以，伪代码为：

每次如果只采样一个方向那么所带来的问题也是显而易见的，积分计算的结果会非常的noisy，虽然蒙特卡洛积分是无偏估计，但样本越少显然偏差越大。

但该问题很好解决，如果每次只去寻找一条路径结果不好，那么重复多次寻找到多条路径，让这多条路径每次选择一个采样方向，那不就是近似的有多条采样方向了嘛，所以，最后将多条路径的结果求得平均即可！如下图所示：

伪代码如下：

通过对经过像素的光线重复采样，每次在反射的时候只按分布随机选取一个方向，解决了只对经过像素的光线采样一次，而对反射光线按分布采样多次所导致的光线爆炸问题。

但是，目前来说，对于一个递归函数，但没有终止条件是十分不妥的，因此：这里采用了俄罗斯轮盘赌的方法来实现递归的结束同时还可以使得结果无偏

首先设定一个概率 $P$, 有 $\mathrm{P}$ 的概率光线会继续递归并设置返回值为 $L_o/P$ ，有 $1-P$ 的概率光线 停止递归，并返回0。这样巧妙的设定之下光线一定会在某次反射之后停止递归，并且计算的结果依然是无偏的，因为 Radiance的期望不变，证明如下:
$$E=P^{\star }(Lo/P)+(1-P{)}^{\star }0=L_0$$

伪代码如下：

最后一个问题：在做路径追踪时，我们直接是对所有的区域进行的采样，而有的区域根本就没有光源，显然这一部分的时间资源消耗是白白浪费掉了

因此在计算直接光照的时候改进为直接对光源所在区域进行积分!这样所有采样的光线都一定会击中光源(如果中间没有别的物体)，没有光线再会被浪费了。

由于原始的方程是对w进行积分的，所以现在要将其变换到光源的范围A上去进行积分：也就是要找到两者之间的关系，然后再做个变换即可：

关系式中的 $\cos {\theta }^{\prime }$ 是为了计算出光源上微分面积元正对半球的面积，之后再按照立体角的定义 $\mathrm{d}\omega =\frac{\mathrm{d}A}{r^2}$ ，除以着色点 $\mathrm{x}$ 与 光源采样点x’距离的平方即可。于是根据图中二者的关系可将洹染方程改写如下:
$$\begin{array}{rl}{L}_o\left(x,{\omega }_o\right)& ={\int }_{{\Omega }^{+}}{L}_i\left(x,{\omega }_i\right)f_r\left(x,{\omega }_i,{\omega }_o\right)\cos \theta \mathrm{d}{\omega }_i\\ & ={\int }_A{L}_i\left(x,{\omega }_i\right)f_r\left(x,{\omega }_i,{\omega }_o\right)\frac{\cos \theta \cos {\theta }^{\prime }}{{\Vert x^{\prime }-x\Vert }^2}\mathrm{d}A\end{array}$$
这样便成功从 ${\omega }_i$ 积分转到了对光源面积A的积分，就可以利用蒙特卡洛的方法对光源进行采样从而计算直接光照的积分值 了，对于间接光照，依然采用先前的方法进行光线方向的均匀采样。最终伪代码如下，分直接光照和间接光照两部分计算:

计算直接光照的时候还需要判断光源与着色点之间是否有物体遮挡，该做法也很简单，只需从着色点x向光源采样点x’发出一条检测光线判断是否与光源之外的物体相交即可，如图所示: