本文最后更新于 2024年9月17日 上午
难！好难！超级难！
之前我们了解到，光栅化方法有一些缺点，如无法很好地支持软阴影，对于间接光的处理也过于简单。此外，对于Glossy类材质，光栅化也很难实现。
光线追踪（Ray Tracing）是一种更加精准的渲染方式，但它非常缓慢，（基本）只能在离线渲染中使用。
对于光线，我们认为它有三个基本性质：
光线始终以直线传播。
光线与光线交叉时不会发生碰撞。
光线始终从光源出发，最终到达眼睛。
由于光路具有可逆性，所以在实际处理光线追踪时，会从相机处发射光线，将弹射到光源的光路进行回溯，得到最终结果。
光线投射（Ray Casting）的核心思想是：
对于每一个像素，生成一条穿过该像素，由相机发射的光线。
当光线和场景物体相交时，将交点记录下来，并作一条交点-光源的连线。
判断交点-光源连线中间是否存在阻挡物。若不存在，则形成了有效光路；否则交点处于阴影内。
如图，相机发出的光线被称为视线（Eye Ray），交点-光源连线被称为阴影射线（Shadow Ray）。判定可见后，我们便有了ViewDir、LightDir与Normal，便可以使用光照模型进行着色了。
对于Whitted风格光线追踪，光线可以反射/折射不止一次：
对于所有光线（无论是否折射/反射光线）与场景物体的交点，都会进行阴影判断。若判断通过，则进行着色，并加到起初发射光线的像素的颜色上。
每次反射/折射都会有一定的能量损失。我们认为Eye Ray为主要光线（Primary Ray）、其他光线为次要光线（Secondary Rays）。
光线是射线，由起点和方向向量所定义。
光线公式：r ( t ) = o ⃗ + t d ⃗ , 0 ≤ t ＜ ∞ r(t) = \vec{o} + t\vec{d}, 0≤t＜\infty r ( t ) = o + t d , 0 ≤ t ＜ ∞
其中，t t t o o o d d d
对于光线-球体交点p p p r ( t ) = o ⃗ + t d ⃗ r(t) = \vec{o}+t\vec{d} r ( t ) = o + t d ( p ⃗ − c ⃗ ) 2 − R 2 = 0 (\vec{p}-\vec{c})^2-R^2=0 ( p − c ) 2 − R 2 = 0
其中，r(t)和p都指的是点坐标，因此联立可得：( o ⃗ + t d ⃗ − c ) 2 − R 2 = 0 (\vec{o}+t\vec{d}-c)^2-R^2=0 ( o + t d − c ) 2 − R 2 = 0
可解得t = − b ± b 2 − 4 a c 2 a t = \frac{-b\pm\sqrt{b^2-4ac}}{2a} t = 2 a − b ± b 2 − 4 a c a = d ⃗ ⋅ d ⃗ a = \vec{d} · \vec{d} a = d ⋅ d b = 2 ( o ⃗ − c ⃗ ) ⋅ d ⃗ b=2(\vec{o}-\vec{c})·\vec{d} b = 2 ( o − c ) ⋅ d c = ( o ⃗ − c ⃗ ) ⋅ ( o ⃗ − c ⃗ ) − R 2 c=(\vec{o}-\vec{c})·(\vec{o}-\vec{c})-R^2 c = ( o − c ) ⋅ ( o − c ) − R 2
b 2 > = 4 a c b^2>=4ac b 2 > = 4 a c
对于一般的隐式表面，有f ( p ) = 0 f(p)=0 f ( p ) = 0
类似于球体相交，联立得f ( o ⃗ + t d ⃗ ) = 0 f(\vec{o}+t\vec{d})=0 f ( o + t d ) = 0
解方程应当得到正实数解t。
我们知道，一个三角形必定可以确定一个平面。我们首先求出光线和平面的交点，然后判断交点是否在三角形内部即可。
平面可以被一个法线N和一个不在法线方向上的点p’所定义：
对于点P，只要满足( p ⃗ − p ’ ⃗ ) ⋅ N ⃗ = 0 (\vec{p}-\vec{p’})·\vec{N}=0 ( p − p ’ ) ⋅ N = 0
联立光线方程可得：t = p ’ ⃗ − o ⃗ ⋅ N ⃗ d ⃗ ⋅ N ⃗ t=\frac{\vec{p’}-\vec{o}·\vec{N}}{\vec{d}·\vec{N}} t = d ⋅ N p ’ − o ⋅ N
需要检查t是否大于0。
Mollder Trumbore算法用于快速判断交点是否在三角形内。
对于顶点P 0 ⃗ 、 P 1 ⃗ 、 P 2 ⃗ \vec{P_0}、\vec{P_1}、\vec{P_2} P 0 、 P 1 、 P 2
因此，我们可以联立方程：O ⃗ + t D ⃗ = ( 1 − b 1 − b 2 ) P 0 ⃗ + b 1 P 1 ⃗ + b 2 P 2 ⃗ \vec{O}+t\vec{D} = (1-b_1-b_2)\vec{P_0}+b_1\vec{P_1}+b_2\vec{P_2} O + t D = ( 1 − b 1 − b 2 ) P 0 + b 1 P 1 + b 2 P 2
我们可以解得：
[ t b 1 b 2 ] \begin{bmatrix}t\\b_1\\b_2\end{bmatrix} ⎣ ⎢ ⎡ t b 1 b 2 ⎦ ⎥ ⎤ 1 S 1 ⃗ ⋅ E 1 ⃗ [ S 2 ⃗ ⋅ E 2 ⃗ S 1 ⃗ ⋅ S ⃗ S 2 ⃗ ⋅ D ⃗ ] \frac{1}{\vec{S_1}·\vec{E_1}}\begin{bmatrix}\vec{S_2}·\vec{E_2}\\\vec{S_1}·\vec{S}\\\vec{S_2}·\vec{D}\end{bmatrix} S 1 ⋅ E 1 1 ⎣ ⎢ ⎡ S 2 ⋅ E 2 S 1 ⋅ S S 2 ⋅ D ⎦ ⎥ ⎤
其中，E 1 ⃗ = P 1 ⃗ − P 0 ⃗ \vec{E_1}=\vec{P_1}-\vec{P_0} E 1 = P 1 − P 0 E 2 ⃗ = P 2 ⃗ − P 0 ⃗ \vec{E_2}=\vec{P_2}-\vec{P_0} E 2 = P 2 − P 0 S ⃗ = O ⃗ − P 0 ⃗ \vec{S} = \vec{O}-\vec{P_0} S = O − P 0 S 1 ⃗ = D ⃗ × E 2 \vec{S_1}=\vec{D}\times{E_2} S 1 = D × E 2 S 2 ⃗ = S ⃗ × E 1 ⃗ \vec{S_2}=\vec{S}\times\vec{E_1} S 2 = S × E 1
只要解出的t、b1、b2、1-b1b2均大于0，则交点位于三角形内部。
如果不优化，对于每个光线，我们都需要遍历所有的三角形，判断是否与光线相交，性能开销非常大。
我们可以使用包围容器（Bounding Volumn）的方式加速这一计算过程。
包围容器将模型包围在一个简单的几何体中（如长方体、球体等）。其关键思想为：若一个光线不与包围盒相交，则必然不与其中的几何体相交。通过判断光线是否与包围盒相交，可以过滤掉大量的“无用”光线。
对于三维的情况，长方体形包围盒最常用。
长方体是三个不同的“对面”（Slab）的交集。下图展示了这一定义的含义：
如图，将长方体的每个面都视为一个无限扩展的平面，一个“对面”就代表着两个平行的面。
一般，我们使用轴对齐包围盒（AABB），这个包围盒的每条边都必定与世界空间的x、y、z轴平行。
如图，对于二维的情况，我们分别求出光线与x、y轴上对面的交点对应的时间t。随后，将得到的两个线段求交集，便可得到最终光线进入和出去盒子的时间差。若光线未与包围盒相交，则交集为空，有t e n t e r > = t e x i t t_{enter}>=t_{exit} t e n t e r > = t e x i t
简而言之，核心为：对于三个对面，第一次与光线相交最晚的时间t1才是光线真正进入包围盒（倘如相交）的时间；第二次与光线相交最早的时间t2是光线真正出去包围盒的时间。
t e n t e r = m a x ( t m i n ) t_{enter}=max({t_{min})} t e n t e r = m a x ( t m i n ) t e x i t = m i n ( t m a x ) t_{exit}=min(t_{max}) t e x i t = m i n ( t m a x )
同时，当t e x i t < 0 t_{exit}<0 t e x i t < 0
当t e x i t > = 0 t_{exit}>=0 t e x i t > = 0 t e n t e r < 0 t_{enter}<0 t e n t e r < 0
总之，当且仅当t e n t e r < t e x i t t_{enter}<t_{exit} t e n t e r < t e x i t t e x i t > = 0 t_{exit}>=0 t e x i t > = 0
为什么使用AABB？
以x轴为例，要想求得t，只需要用p’（定义平面的点）的x轴坐标减去光线原点o的x坐标，除以方向向量d的x分量，即可得到t。
完成包围盒的构建后，可以对其进行预处理（Preprocess）。
首先，将包围盒划分为若干网格，将其中与物体相交的网格打上标记。
随后，当光线射入时，我们首先专注于打标记的网格。对于这部分网格，光线有可能与其中的物体相交，其他网格则完全没可能。
同时，当光线步入一个网格时，我们也不需要对周围的全部四个网格做判断。我们只需要根据光线的位置，判断其行进位置上的两个网格即可。如：光线向右上方行进，我们判断当前网格上方和右方的网格是否与光线相交即可。
使用该技术需要仔细调整格子数量。过于稀疏会导致加速效率低，过于密集会导致性能消耗高。
在三维空间中，划分的格子数应当为：27*物体数量。
当物体数量多且分布均匀时，该技术效果较好，如草地、灌木等。而对于人造建筑场景则效果欠佳。
实际的三维场景中，一般会存在某些地方物体密集、某些地方物体稀疏的情况。对于物体密集的地方，网格可以密集一些；稀疏的地方则反之。这类技术被称为空间划分（Spatial Partitioning）。
常见的空间划分技术有八叉树（Oct-Tree）、KD树和BSP树。
对于八叉树方法，它会不断地将大立方格切分为八个小立方格，直到切分后的小立方格中有一定比例不与物体相交为止。
对于KD树方法，它会交替地沿坐标轴划分。例如，在图中，是先做一次水平，对于划分后的小网格再做一次垂直划分，以此类推。
对于BSP树方法，它会按不同方向进行划分。但由于其划分方向不与坐标轴平行，所以对于它划分出的格子，计算量较大。
对于KD树，我们需要满足下列特性的数据结构：
内部节点需要存储切分的轴向
内部节点需要存储切分的位置
内部节点包含子节点的指针
内部节点不存储对象
子节点存储一系列对象。
逐节点遍历，首先判断光线是否与该节点代表的包围盒相交。若不相交，则到下一节点（忽视下级所有子节点，前往下一同级的节点）。若相交，且遍历到叶子结点时，求光线是否与其中的物体相交。若相交，且遍历到内部节点时，继续遍历子节点。
但KD树存在一个重大缺陷，就是将一个三角形纳入某个包围盒的判定范围内是非常困难的。我们很难知道一个包围盒把哪些三角形囊括进去了，因此现在进行空间划分预处理很少用KD树。同时，一个物体可以存在多个叶子节点，会导致重复计算。
空间划分以空间为抓手，而物体划分（Object Partition）则以物体为基准。
包围容器层级（Bounding Volume Hierarchy，BVH）是一种物体划分技术，该技术得到了非常广泛的应用，因为它解决了KD树的两个重大缺陷。
对于BVH，同样是构建二叉树，首先将场景内所有物体囊括入一个包围盒，然后将所有物体组织为两部分，分别构建包围盒：
以此类推，划分到一个节点里最多有N个（一般为5）三角形为止。需要注意的是，尽管图上显示的大蓝色三角形有一部分在黄色包围盒里，但实际存储时不会出现交叠的现象。
总结：BVH过程如下：
找到包围盒
递归地将一系列物体划分到两个子集中
重新计算子集的包围盒
在必要时停止。
对于BVH，划分方式十分重要。我们可以每次选择一个不同的维度划分来保证划分均匀；也可以选择当前最长的轴向（包围盒中最长的边）进行划分，用于处理长条状包围盒的情况；也可以取中间的物体的位置（快速划分算法，std::nth_element）进行划分。
BVH数据结构如下：
辐射度量学（Radiometry）在物理上准确定义了光照。它定义了光的空间性质，包括辐射通量（Radiant Flux）、辐射强度（Intensity）、辐照度（Irradiance）与辐亮度（Radiance）。
辐射能量（Radiant Energy）是电磁辐射的能量，单位为焦耳：Q [ J = J o u l e ] Q[J=Joule] Q [ J = J o u l e ]
辐射通量（Radiant Flux / Power）是物体在每单位时间内发射、反射、转移或接受的能量，单位为瓦特或流明：$\Phi\equiv\frac{dQ}{dt}[W = Watt][lm = lumen]^* $
辐射强度（Radiant Intensity）定义了一个光源发射出去的能量。
辐照度（Irradiance）定义了落在一个表面的光能量。
辐亮度（Radiance）定义了光线在传播过程中携带的能量。
辐射强度是一个点光源每单位立体角发射出的能量，单位为瓦特每立体角或流明每立体角（即坎德拉，cd）。
I ( ω ) ≡ d Φ d ω I(\omega)\equiv\frac{d\Phi}{d\omega} I ( ω ) ≡ d ω d Φ
[ W s r ] [ l m s r = c d = c a n d e l a ] [\frac{W}{sr}][\frac{lm}{sr}=cd=candela] [ s r W ] [ s r l m = c d = c a n d e l a ]
何为立体角（Solid Angle）？
对于二维平面，对角弧与半径的比率被称为角。角的默认单位为弧度（Radian）。
立体角是角度在三维空间中的延伸概念。对于一个球，对角面的面积A与半径平方的比率就是立体角。即：Ω = A r 2 \Omega = \frac{A}{r^2} Ω = r 2 A 4 π 4\pi 4 π
当我们已知球面坐标方向和球半径时，d A = r 2 s i n θ d θ d Φ dA = r^2sin\theta\ d\theta d\Phi d A = r 2 s i n θ d θ d Φ d ω = d A r 2 = s i n θ d θ d Φ d\omega = \frac{dA}{r^2} = sin\theta\ d\theta d\Phi d ω = r 2 d A = s i n θ d θ d Φ
其中，θ \theta θ ϕ \phi ϕ
辐射强度I与辐射通量Φ的关系如下：
Φ = ∫ S 2 I d ω = 4 π I \Phi = \int_{S^2}Id\omega=4\pi I Φ = ∫ S 2 I d ω = 4 π I
I = Φ 4 π I = \frac{\Phi}{4\pi} I = 4 π Φ
辐照度（Irradiance）是一个表面点上每单位区域的功率，即：单位面积上的辐射通量。
辐照度以字母E定义，公式：E ( x ) ≡ d Φ ( x ) d A E(x)\equiv \frac{d\Phi(x)}{dA} E ( x ) ≡ d A d Φ ( x )
需要注意，这里的单位面积所在的面必须与光线入射方向垂直。若不垂直则需要对面做投影操作（Lambert Cosine Law）。
辐照度随距离以r^2衰减。
辐亮度（Radiance）是单位立体角、单位面积上由表面发射、反射、传播或接收的功率（即单位立体角单位面积上的辐射通量，即单位立体角上的辐照度或单位面积上的辐射强度），它定义了光在传播过程中携带的能量。
辐亮度以L表示，公式：L ( p , ω ) ≡ d 2 Φ ( p , ω ) d ω d A c o s θ L(p,\omega)\equiv \frac{d^2\Phi(p,\omega)}{d\omega dA\ cos\theta} L ( p , ω ) ≡ d ω d A c o s θ d 2 Φ ( p , ω )
简单理解辐亮度：辐照度指的是单位面积上的辐射通量。把单位面积想象成一个小圆片，这个小圆片会向其上方各个方向发射光线。如果我们专注于其中一个方向，计算其辐射通量，那么最终得到的就是辐亮度。
辐亮度可以分为入射辐亮度（Incident Radiance），公式：L ( p , ω ) = d E ( p ) d ω c o s θ L(p,\omega) = \frac{dE(p)}{d\omega\ cos\theta} L ( p , ω ) = d ω c o s θ d E ( p ) L ( p , ω ) = d I ( p , ω ) d A c o s θ L(p,\omega) = \frac{dI(p,\omega)}{dA\ cos\theta} L ( p , ω ) = d A c o s θ d I ( p , ω )
辐照度和辐亮度的关系如下：
辐照度Irradiance是被单位面积dA接收到的总辐射通量；辐亮度Radiance是从单位立体角（可以理解为特定方向）被单位面积dA接收到的辐射通量。
下面的公式描述了二者的关系：
d E ( p , ω ) = L i ( p , ω ) c o s θ d ω dE(p,\omega) = L_i(p,\omega)cos\theta\ d\omega d E ( p , ω ) = L i ( p , ω ) c o s θ d ω
E ( p ) = ∫ H 2 L i ( p , ω ) c o s θ d ω E(p)=\int_{H^2}L_i(p,\omega)cos\theta\ d\omega E ( p ) = ∫ H 2 L i ( p , ω ) c o s θ d ω
其中，c o s θ cos\theta c o s θ
双向反射分布函数（Bidirectional Reflectance Distribution Function，BRDF）描述了反射光线的角度于能量大小的分布关系。
我们可以将反射的本质理解为：一个小表面吸收了来自立体角方向ω i \omega_i ω i d E dE d E
但实际上，表面吸收了特定方向的能量后，散发出的能量是朝向四面八方的。为了了解特定方向散发出的能量的信息，我们引入BRDF函数。它描述了如何把表面接收到的能量反射到另一个方向。
公式如下：
f r ( ω i → ω r ) = d L r ( ω r ) d E i ( ω i ) = d L r ( ω r ) L i ( ω i ) c o s θ i d ω i f_r(\omega_i \rightarrow \omega_r) = \frac{dL_r(\omega_r)}{dE_i(\omega_i)} = \frac{dL_r(\omega_r)}{L_i(\omega_i)cos\theta_id\omega_i} f r ( ω i → ω r ) = d E i ( ω i ) d L r ( ω r ) = L i ( ω i ) c o s θ i d ω i d L r ( ω r )
单位为1/sr，即球面度的倒数。
我们可以将该公式直观地解释为：
微表面从入射方向吸收能量，并向出射方向反射能量。反射出来的能量相当于（出射方向的辐亮度的微分）除以（入射方向的辐亮度与入射方向立体角，以及入射方向与上向量夹角的乘积）
它描述了给定出射角占入射角总能量的比例。
对于给定的反射点p p p r r r H 2 H^2 H 2 L r ( p , ω r ) = ∫ H 2 f r ( p , ω i → ω r ) L i ( p , ω i ) c o s θ i d ω i L_r(p,\omega_r) = \int_{H^2}f_r(p,\omega_i\rightarrow\omega_r)L_i(p,\omega_i)cos\theta_i\ d\omega_i L r ( p , ω r ) = ∫ H 2 f r ( p , ω i → ω r ) L i ( p , ω i ) c o s θ i d ω i
其中，L i ( p , ω i ) c o s θ i d ω i L_i(p,\omega_i)cos\theta_id\omega_i L i ( p , ω i ) c o s θ i d ω i
一次BRDF仅考虑了直接由光源发射的入射光。但实际上，还会有反射和折射光打到微表面，即：任意点的出射辐亮度也可能称为其他点的入射辐亮度。因此，实际上的渲染方程（Rendering Equation）是一个递归的过程。
当我们考虑物体自发光的情况时，我们只需要将反射辐亮度项加上一个自发光辐亮度项，即：
L o ( p , ω o ) = L e ( p , ω o ) + ∫ Ω + L i ( p , ω i ) f r ( p , ω i , ω o ) ( n ⋅ ω i ) d ω i L_o(p,\omega_o)=L_e(p,\omega_o)+\int_{\Omega+}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n·\omega_i)d\omega_i L o ( p , ω o ) = L e ( p , ω o ) + ∫ Ω + L i ( p , ω i ) f r ( p , ω i , ω o ) ( n ⋅ ω i ) d ω i
该方程便被称为渲染方程。其中，Ω + \Omega+ Ω + H 2 H^2 H 2 n ⋅ ω i n·\omega_i n ⋅ ω i c o s θ i cos\theta_i c o s θ i
需要注意，渲染方程中所有的方向都默认是由“球心”指向“球外”的。
我们考虑一个点光源的情况。
对于点光源，只会有一条光线直接辐射到微表面。因此，渲染方程退化为：
L r ( x , ω r ) = L e ( x , ω r ) + L i ( x , ω i ) f ( x , ω i , ω r ) ( ω i , n ) L_r(x,\omega_r) = L_e(x,\omega_r)+L_i(x,\omega_i)f(x,\omega_i,\omega_r)(\omega_i,n) L r ( x , ω r ) = L e ( x , ω r ) + L i ( x , ω i ) f ( x , ω i , ω r ) ( ω i , n )
可翻译为：反射光 = 自发光 + 来自光源的入射光 * BRDF * 入射角余弦值
考虑到多个点光源的情况，只需要把反射辐亮度项叠加即可：
L r ( x , ω r ) = L e ( x , ω r ) + ∑ L i ( x , ω i ) f ( x , ω i , ω r ) ( ω i , n ) L_r(x,\omega_r) = L_e(x,\omega_r)+\sum L_i(x,\omega_i)f(x,\omega_i,\omega_r)(\omega_i,n) L r ( x , ω r ) = L e ( x , ω r ) + ∑ L i ( x , ω i ) f ( x , ω i , ω r ) ( ω i , n )
考虑面光源的情况，面光源可以看作是无数微小点光源的集合。只需要把这些点光源对微表面辐亮度的贡献积分，即可得到面光源对微表面辐亮度的贡献：
L r ( x , ω r ) = L e ( x , ω r ) + ∫ Ω L i ( x , ω i ) f ( x , ω i , ω r ) c o s θ i d ω i L_r(x,\omega_r) = L_e(x,\omega_r)+\int_\Omega L_i(x,\omega_i)f(x,\omega_i,\omega_r)cos\theta_id\omega_i L r ( x , ω r ) = L e ( x , ω r ) + ∫ Ω L i ( x , ω i ) f ( x , ω i , ω r ) c o s θ i d ω i
当我们考虑其他物体反射到微表面的光线时，我们可以把其他物体的反射面当成面光源。
L r ( x , ω r ) = L e ( x , ω r ) + ∫ Ω L r ( x ’ , − ω i ) f ( x , ω i , ω r ) c o s θ i d ω i L_r(x,\omega_r) = L_e(x,\omega_r)+\int_\Omega L_r(x’,-\omega_i)f(x,\omega_i,\omega_r)cos\theta_id\omega_i L r ( x , ω r ) = L e ( x , ω r ) + ∫ Ω L r ( x ’ , − ω i ) f ( x , ω i , ω r ) c o s θ i d ω i
可翻译为：反射光 = 自发光 + 被其他物体反射并入射到微表面的光 * BRDF * 入射角余弦值
需要注意，该方程没有考虑光源的直接光线。其中，x ’ x’ x ’
该形态的渲染方程可被简写为：I ( u ) = e ( u ) + ∫ I ( v ) K ( u , v ) d v I(u) = e(u) +\int I(v)K(u,v)dv I ( u ) = e ( u ) + ∫ I ( v ) K ( u , v ) d v L = E + K L L=E+KL L = E + K L
简化方程的目的是解出L的值。简化后我们得到L = ( I − K ) − 1 E L=(I-K)^{-1}E L = ( I − K ) − 1 E
对于( I − K ) − 1 (I-K)^{-1} ( I − K ) − 1
L = ( I − K ) − 1 E = ( I + K + K 2 + K 3 + . . . ) E = E + K E + K 2 E + K 3 E + . . . L = (I-K)^{-1}E = (I+K+K^2+K^3+...)E = E+KE+K^2E+K^3E+... L = ( I − K ) − 1 E = ( I + K + K 2 + K 3 + . . . ) E = E + K E + K 2 E + K 3 E + . . .
由此，我们便可以认为：反射光 = 自发光 + 直接光照 + 一次反射光照 + 两次反射光照 + …，以此类推。
对于连续的随机变量，概率密度函数（Probability Distribution Function, PDF）描述了概率在随机变量变化时的变化。连续随机变量的期望为：E [ X ] = ∫ x p ( x ) d x E[X] = \int xp(x)dx E [ X ] = ∫ x p ( x ) d x p ( x ) p(x) p ( x ) p ( x ) > = 0 a n d ∫ p ( x ) d x = 1 p(x)>=0\ and \int p(x)dx=1 p ( x ) > = 0 a n d ∫ p ( x ) d x = 1
假设随机变量X X X X ∼ p ( x ) X \sim p(x) X ∼ p ( x ) Y Y Y X X X Y = f ( X ) Y = f(X) Y = f ( X ) E [ Y ] = E [ f ( x ) ] = ∫ f ( x ) p ( x ) d x E[Y] = E[f(x)] = \int f(x)p(x)dx E [ Y ] = E [ f ( x ) ] = ∫ f ( x ) p ( x ) d x
当我们想计算一个复杂函数的定积分 时，可以使用蒙特卡洛积分法（Monte Carlo Integration）。
该方法的核心思想是：通过在积分域内随机采样若干函数值并将他们平均，来估算函数积分。
给定定积分∫ a b f ( x ) d x \int_a^bf(x)dx ∫ a b f ( x ) d x X i ∼ p ( x ) X_i\sim p(x) X i ∼ p ( x ) F N = ∫ f ( x ) d x = 1 N ∑ i = 1 N f ( X i ) p ( X i ) F_N=\int f(x)dx=\frac{1}{N}\sum\limits_{i=1}\limits^N\frac{f(X_i)}{p(X_i)} F N = ∫ f ( x ) d x = N 1 i = 1 ∑ N p ( X i ) f ( X i )
其中，f ( X ) p ( X ) \frac{f(X)}{p(X)} p ( X ) f ( X )
以均匀随机变量X i ∼ p ( x ) = C X_i\sim p(x) = C X i ∼ p ( x ) = C ∫ a b p ( x ) d x = 1 ⇒ ∫ a b C d x = 1 ⇒ C = 1 b − a \int_a^bp(x)dx = 1\Rightarrow\int^b_aCdx=1\Rightarrow C = \frac{1}{b-a} ∫ a b p ( x ) d x = 1 ⇒ ∫ a b C d x = 1 ⇒ C = b − a 1
则有蒙特卡洛近似：F N = b − a N ∑ i = 1 N f ( X i ) F_N = \frac{b-a}{N}\sum\limits_{i=1}\limits^{N}f(X_i) F N = N b − a i = 1 ∑ N f ( X i )
对于Whitted风格的光线追踪，其遵循两个原则：1. 在光滑物体上总是执行反射或折射 2. 在漫反射物体上停止弹射。但实际上的光线行为完全不是这么草率的。路径追踪（Path Tracing）是更实际的、更符合物理的光线追踪方法。
能几乎完全反射光线的材质被称为镜面材质（Specular Material），而略带磨砂质感的金属材质被称为光泽材质（Glossy Material），这类材质无法被Whitted风格光线追踪处理。
渗色（Color Bleeding）指无法接收到直接光照的表面由于接收到其他表面反射的简介光照而呈现出其他表面颜色的情况。如下图，左边黑暗的表面再右边会呈现出不属于自身表面颜色的颜色。
我们重新观察渲染方程：
L o ( p , ω o ) = L e ( p , ω o ) + ∫ Ω + L i ( p , ω i ) f r ( p , ω i , ω o ) ( n ⋅ ω i ) d ω i L_o(p,\omega_o) = L_e(p,\omega_o)+\int_{\Omega +}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n·\omega_i)d\omega_i L o ( p , ω o ) = L e ( p , ω o ) + ∫ Ω + L i ( p , ω i ) f r ( p , ω i , ω o ) ( n ⋅ ω i ) d ω i
解渲染方程有两个关键：1. 在半球上求积分 2. 递归求解
假设我们想在具备一个面光源、一个盒子和一个地板的简单场景中，对单个无自发光点进行直接光照渲染：
L o ( p , ω o ) = ∫ Ω + L i ( p , ω i ) f r ( p , ω i , ω o ) ( n ⋅ ω i ) d ω i L_o(p,\omega_o) = \int_{\Omega +}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n·\omega_i)d\omega_i L o ( p , ω o ) = ∫ Ω + L i ( p , ω i ) f r ( p , ω i , ω o ) ( n ⋅ ω i ) d ω i
可以看出，我们最终要求解的是一个积分域为半球的积分，可以使用蒙特卡罗方法：∫ b a f ( x ) d x = 1 N ∑ i = 1 N f ( X i ) p ( X i ) \int_b^a f(x)dx=\frac{1}{N}\sum\limits_{i=1}\limits^N\frac{f(X_i)}{p(X_i)} ∫ b a f ( x ) d x = N 1 i = 1 ∑ N p ( X i ) f ( X i )
在这里，f ( x ) f(x) f ( x ) L i ( p , ω i ) f r ( p , ω i , ω o ) ( n ⋅ ω i ) L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n·\omega_i) L i ( p , ω i ) f r ( p , ω i , ω o ) ( n ⋅ ω i ) p ( ω i ) p(\omega_i) p ( ω i ) p ( ω i ) = 1 / 2 π p(\omega_i) = 1/2\pi p ( ω i ) = 1 / 2 π
如此，蒙特卡罗方法中的f ( x ) 、 p ( x ) f(x)、p(x) f ( x ) 、 p ( x )
L o ( p , ω o ) ≈ 1 N ∑ i = 1 N L i ( p , ω i ) f r ( p , ω i , ω o ) ( n ⋅ ω i ) d ω i p ( ω i ) L_o(p,\omega_o) \approx \frac{1}{N}\sum\limits_{i=1}\limits^N\frac{L_i(p,\omega_i)f_r(p,\omega{i},\omega{o})(n·\omega_i)d\omega_i}{p(\omega_i)} L o ( p , ω o ) ≈ N 1 i = 1 ∑ N p ( ω i ) L i ( p , ω i ) f r ( p , ω i , ω o ) ( n ⋅ ω i ) d ω i
每次均匀地在半球上取一个入射方向ω i \omega_i ω i
伪代码如下：
1 2 3 4 5 6 7 8 shade(p,wo)
我们无需将其他表面反射的辐亮度与光源发射的辐亮度区别对待。我们可以把反射面当作一个光源。
如图，要想计算Q点到P点发射的辐亮度，我们只需要把该过程视作：相机处于P位置，计算Q点的直接光照向QP方向发射的辐亮度。如此我们便完成了一次递归过程。
写作伪代码：
1 2 3 4 5 6 7 8 9 10 shade(p,wo)
然而，目前的算法存在指数爆炸 的问题。假设每一根光线会引出N个出射光线，那么M次弹射后就会共存M^N个光线。要解决这一问题，N，即出射光线的数量只能为1。伪代码如下:
1 2 3 4 5 6 7 shade(p , wo)direction wi~pdf(w)a ray r(p , wi)Li * f_r * cosine / pdf(wi)object at q q , -wi))* f_r * cosine / pdf(wi)
该伪代码表示的算法便是路径追踪 ，其显著特点为一条入射光对应一条出射光，该出射光的方向是随机采样的。
由于唯一的出射光是随机采样的，因此单词采样必定会产生极大的噪声。为此，我们可以从一个像素发射多条光线，对最终结果进行平均。伪代码如下：
1 2 3 4 5 6 7 8 ray_generation(camPos,, pixel)within the pixel0.0 each sample in the pixela ray r(camPos, cam_to_sample)the scene at p1 / N * shade(p, sample_to_cam)
之前shade函数可能会发生无限递归的情况。为此，我们需要设置终止条件。
我们知道，在自然界中光线是无限次弹射的，每次弹射都会损失掉一定的能量。但我们无法在计算机中模拟无限次弹射，只能在弹射到一定次数时将其终止。而这种操作就会导致能量不守恒。
为了解决这一问题，我们引入俄罗斯轮盘赌技术（Russian Roulette，RR）：
RR方法用于决定何时停止弹射。我们手动定义一个概率P (0<P<1)，以该概率发射光线（即：P概率发射光线，1-P概率不发射光线），并将着色结果除以概率P。如此这般，最终结果的期望就是Lo：
E = P ∗ ( L o / P ) + ( 1 − P ) ∗ 0 = L o E = P * (Lo/P) + (1 - P) * 0 = Lo E = P ∗ ( L o / P ) + ( 1 − P ) ∗ 0 = L o
伪代码：
1 2 3 4 5 6 7 8 9 10 shade(p , wo)a probability P_RRa uniform dist. in [0,1] 0.0 ;direction wi~pdf(w)a ray r(p , wi)object at q q , -wi)*fr*cosine/pdf(wi)/P_RR
对于均匀随机分布，当面光源的面积非常小时，打到光源的概率就会变得非常小，会有很多光线被“浪费掉”。为此，我们可以优化采样方式。
我们可以直接从光源采样，这样所有光线都不会被浪费了。我们假定面光源的面积为A，则在光源表面采样的PDF等于1/A。由于原始渲染方程是对dω积分而非dA，所以我们可以将dA转化为dω：
d ω = d A c o s θ ’ ∣ ∣ x ’ − x ∣ ∣ 2 d\omega = \frac{dAcos\theta ’}{||x’-x||^2} d ω = ∣ ∣ x ’ − x ∣ ∣ 2 d A c o s θ ’
该操作实际上是将面光源的表面投影到半球上。
随后重写渲染方程：
L o ( p , ω o ) = ∫ Ω + L i ( p , ω i ) f r ( p , ω i , ω o ) ( n ⋅ ω i ) d ω i = ∫ A L i ( x , ω i ) f r ( x , ω i , ω o ) c o s θ c o s θ ′ ∣ ∣ x ′ − x ∣ ∣ 2 d A L_o(p,\omega_o) = \int_{\Omega +}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n·\omega_i)d\omega_i = \int_AL_i(x,\omega_i)f_r(x,\omega_i,\omega_o)\frac{cos\theta cos\theta'}{||x'-x||^2}dA L o ( p , ω o ) = ∫ Ω + L i ( p , ω i ) f r ( p , ω i , ω o ) ( n ⋅ ω i ) d ω i = ∫ A L i ( x , ω i ) f r ( x , ω i , ω o ) ∣ ∣ x ′ − x ∣ ∣ 2 c o s θ c o s θ ′ d A
对于该方程，f(x)是积分符号和dA之间的式子，p(x)为1/A，二者已知，可以使用蒙特卡罗方法处理。
我们可以将着色点受到的辐亮度分为两部分：
来自光源（直接光照，无需RR）
其他反射光（简介光照，使用RR）
与此同时，我们需要判定直接光照与着色点之间是否有物体遮挡。取连线并判断是否与其他物体相交即可。
最终伪代码：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 shade(p, wo)the light at x' (pdf_light = 1 / A)from p to x'the ray is not blocked in the middle 2 / pdf_light0.0 with probability P_RRthe hemisphere toward wi (pdf_hemi = 1 / 2 pi )at q
