6 Ray tracing¶

文本统计：约 4648 个字 • 98 行代码

6.0 Introduction¶

Rasterization couldn’t handle global effects well, Such as soft shadows, glossy reflection, indirect illumination. 光栅化很难处理一些全局的效果

与光栅化相比：光栅化快但是不准确，使用在实时的场景中；光线追踪则更准确但是比较慢，适用于一些离线的场景中。

6.1 Whitted-Style Ray Tracing¶

Light Rays:

Light travels in straight lines 光沿着直线传播
Light rays do not “collide” with each other if they cross 光之间不会相互影响的
Light rays travel from the light sources to the eye, 并且光路具有可逆性

Ray Casting¶

Generate an image by casting one ray per pixel 通过每个像素投射一条光线来生成图像
Check for shadows by sending a ray to the light 通过向光源发送光线来检查阴影

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

当然一条光线自然可能会与不止一个物体相交，但是考虑遮挡关系，只去找最近的交点。接着连接该交点和光源，只要判断这条连线之间是否有物体存在就可以知道该交点是否在阴影之中（怎么样，是不是比shadow mapping那一套简单了许多）：

紧接着，自然可以利用Blinn-Phong模型对这个点进行局部光照模型计算，得到该像素的颜色，那么遍历所有近投影平面上的像素就能得到一张完整的图像。但如果光线追踪仅仅是在第一步Ray Casting就停止的话，那么它的效果与局部光照模型是一样的，因此我们需要第二步，真正的考虑全局效果

Recursive (Whitted-Style) Ray Tracing¶

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

因此每一个交点的颜色贡献来自这样种几类型直接光照，反射方向间接光，折射方向间接光（如果有折射的话）

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

总的来说，实现 Whitted-Style Ray Tracing 需注意

（1）整体过程是一个递归的过程，因此需要一定的递归终止条件，比如说允许的最大反射或折射次数为10。

（2）光线在每次反射和折射之后都有能量损耗的，由系数决定，因此越往后的折射和反射光贡献的能量越小，这也是为什么在上文中提到根据光线能量权重求和。e.g. 反射系数为0.7，那么第一次反射折损 \(30\%\)，第二次反射折损 \(1-（70\%\times70\%）\)，依次类推。

（3）如果反射或折射光线没有碰撞到物体，一般直接返回一个背景色。

（4）有一些关于光线表示，及如何求交点的实现细节在下一节里讨论。

RayTracing(original_point, ray_direction, objects, depth)
{
    if (depth > maxDepth)
        return color(0,0,0);

    if (IsHitObject(original_point, ray_direction, objects))
    {
        hitPoint = GetHitPoint();
        normal = GetNormal();

        reflectionDirection = reflect(ray_direction, normal);
        refractionDirection = refract(ray_direction, normal); // 如果没有折射，舍弃折射项

        local_color = BlinnPhongShader(original_point, normal, light_position)

        return local_color
            + k1 * RayTracing(hitPoint, reflectionDirection, objects, depth + 1)
            + k2 * RayTracing(hitPoint, refractionDirection, objects, depth + 1);
    }
    else
        return background_color;
}

6.2 Ray-Surface Intersection¶

Ray Equation : Ray is defined by its origin and a direction vector

for example : \(\pmb r (t)=\pmb{o}+ t \pmb d \quad \ 0\le t < \infty\)

6.2.1 Ray Insertion With Sphere¶

\[ \pmb p:(\pmb p-\pmb c)^2-R^2=0 \]

只需要解一个关于t的二次方程：

\[ (\pmb o+t\pmb d-\pmb c)^2-R^2=0 \]

6.2.2 Ray Intersection With Implicit Surface¶

General implicit surface

\[ \pmb p:f(\pmb p)=0 \]

ray equation

\[ f(\pmb o+t\pmb d)=0 \]

solve for real and positive roots

6.2.3 Ray Intersection With Triangle Mesh¶

我们先求光线与平面的交点，然后判断交点是否在三角形内

Plane equation :

Plane is defined by normal vector and a point on plane

我们可以得到相应的平面方程：\(\pmb p:(\pmb p-\pmb p')\cdot \pmb N=0\)

解方程

\[ (\pmb o+t\pmb d-\pmb p')\cdot \pmb N=0 \]

得到

\[ t=\frac{(\pmb p'-\pmb o)\cdot \pmb N}{\pmb d\cdot \pmb N} \]

并验证t是不是正的。

有一个更快的算法：Moller Trumbore Algorithm

利用三角形的重心坐标，然后利用Crammer法则，解得三个变量。当三个量都是非负的时候，交点就在三角形内。

\[ \vec{O} + t\vec{D} = (1 - b_1 - b_2)\vec{P}_0 + b_1\vec{P}_1 + b_2\vec{P}_2 \]

\[ \begin{bmatrix} t \\ b_1 \\ b_2 \end{bmatrix} = \frac{1}{\vec{S}_1 \cdot \vec{E}_1} \begin{bmatrix} \vec{S}_2 \cdot \vec{E}_2 \\ \vec{S}_1 \cdot \vec{S} \\ \vec{S}_2 \cdot \vec{D} \end{bmatrix} \]

其中

\[ \left\{ \begin{aligned} &\vec{E}_1 = \vec{P}_1 - \vec{P}_0\\ &\vec{E}_2 = \vec{P}_2 - \vec{P}_0\\ &\vec{S} = \vec{O} - \vec{P}_0\\ &\vec{S}_1 = \vec{D} \times \vec{E}_2\\ &\vec{S}_2 = \vec{S} \times \vec{E}_1 \end{aligned} \right. \]

Note

具体推导过程:

\[ \vec{O} + t\vec{D} = (1 - b_1 - b_2)\vec{P}_0 + b_1\vec{P}_1 + b_2\vec{P}_2 \]

用矩阵来表达就是:

\[ \vec O - \vec P_2 = \begin{bmatrix} -\vec D & \vec P_1-\vec P_2 & \vec P_2 -\vec P_1 \end{bmatrix} \begin{bmatrix} t \\ b_1 \\ b_2 \end{bmatrix} \]

再进行换元:

\[ \vec O - \vec P_2 = \begin{bmatrix} -\vec D & \vec E_1 & \vec E_2 \end{bmatrix} \begin{bmatrix} t \\ b_1 \\ b_2 \end{bmatrix} \]

根据 Crammer 法则解得

\[ \left\{ \begin{aligned} t=\frac{|\begin{matrix} \vec S & \vec E_1 & \vec E_2 \end{matrix}|} {|\begin{matrix} -\vec D & \vec E_1 & \vec E_2 \end{matrix}|}=\frac{(\vec S_1\times\vec E_1)\cdot \vec E_2}{(\vec D\times \vec E_2)\cdot \vec E_1}\\ b_1=\frac{|\begin{matrix} -\vec D & \vec S & \vec E_2 \end{matrix}|} {|\begin{matrix} -\vec D & \vec E_1 & \vec E_2 \end{matrix}|}=\frac{(\vec D\times\vec E_2)\cdot \vec S}{(\vec D\times \vec E_2)\cdot \vec E_1}\\ b_2=\frac{|\begin{matrix} -\vec D & \vec E_1 & \vec S \end{matrix}|} {|\begin{matrix} -\vec D & \vec E_1 & \vec E_2 \end{matrix}|}=\frac{(\vec S\times\vec E_1)\cdot \vec D}{(\vec D\times \vec E_2)\cdot \vec E_1}\\ \end{aligned} \right. \]

再换元即得:

\[ \begin{bmatrix} t \\ b_1 \\ b_2 \end{bmatrix} = \frac{1}{\vec{S}_1 \cdot \vec{E}_1} \begin{bmatrix} \vec{S}_2 \cdot \vec{E}_2 \\ \vec{S}_1 \cdot \vec{S} \\ \vec{S}_2 \cdot \vec{D} \end{bmatrix} \]

6.3 Accelerating Ray-Surface Intersection¶

我们当然可以简单地就计算每个三角形与光线是否有交点，但这样计算量太大了，运行速度也比较慢。

于是我们引入一个方法：Bounding Volumes

6.3.1 Bounding Volumes¶

我们构建一个box，使其完全包裹住物体，于是如果光线与box并未相交，那么光线也不会与物体相交。

我们理解box的方法：把它看成三个对面的交。我们经常使用的box: Axis-Aligned Bounding Box (AABB)

那么如何判断光线和box是否相交？

我们由二维的情况引入：Compute intersections with slabs and take intersection of \(t_{min},t_{max}\) intervals

计算两个 \(x,y\) 的 \(t_{\min},t_{\max}\)，然后将两个区间求交集，即得光线在box中的部分（交集满足 \(x,y\) 均在长方体范围内）

在三维的情况就是光线在三个平面都在盒子中，于是只要分别计算\(t_{\min}\) 和 \(t_{\max}\)，然后 \(t_{enter}=\max\{t_{\min}\}\) , \(t_{exit}=\min\{t_{\max}\}\), 如果\(t_{enter}<t_{exit}\)，那么这个情况就是光线穿过盒子的情况。

然而光线并不是一条直线，而是一条射线，所以我们应该检查 \(t\) 是否小于0

\(t_{exit}<0\) \(\Rightarrow\) The box is “behind” the ray - no intersection
\(t_{exit}\geq 0\) and \(t_{enter}<0\ \Rightarrow\) The ray’s origin is inside the box - have intersection

In summary, ray and AABB intersect iff

\[ t_{enter}<t_{exit} \text{ and }t_{exit}\geq0 \]

6.3.2 Using AABBs to accelerate ray tracing¶

Pre-Processing，对画面进行预处理，让我们使用 bounding volumes 求交的时候会更快。

1 Uniform Spatial Partitions (Grids)¶

Preprocess

Find bounding boxCreate girdStore each object in overlapping cells

根据光线方向判断与那些盒子相交，如果这个盒子里面有物体，那么再与物体求交

Grid 的密度

如果只有一个 cell，那么就根本没有加速

如果格子太多了，那么光线传播过程冗余的太多了

于是我们可以找一个平衡的切分方式

Grid Resolution: #cells = C*#objs ( C 大约是 27 in 3D)

Grids work well on large collections of objects that are distributed evenly in size and space

但是在物体较为稀疏的情况下效果就不太好。

2 Spatial partition¶

Oct-Tree(八叉树)

在立体的情况下将一个大的box切成八块，根据某个标准切到某个大小为止，但是随着维度的增高，叉数越来越大，所以这个方法较少使用。

我们在这里着重介绍KD-Tree

Data structure for KD-trees

Internal nodes store

split axis: x-, y-, or z-axis
split position: coordinate of split plane along axis
children: pointers to child nodes
No objects are stored in internal nodes

Leaf nodes store

list of objects

我只需要在叶子节点中存储物体的信息，内部节点存储分割的轴以及分割的位置

Traversing a KD-tree

问题：较难判断这个box是否与三角形有交点；同一个物体会出现在多个叶子节点中

3 Objects Partition & Bounding Volume Hierarchy (BVH)¶

这种方法可以把保证一个物体仅可能出现在唯一一个叶子节点中，也可以很好地解决第一个问题

Build BVH

Step 1Step 2Step 3Step 4

按照物体进行划分，然后分别求包围盒

总结来说步骤如下

Find bounding box
Recursively split set of objects in two subsets
Recompute the bounding box of the subsets 重新计算子集的包围盒
Stop when necessary
Store objects in each leaf node

How to subdivide a node?

choose a dimension to split 选取一个维度去切分

Heuristic #1: Always choose the longest axis in node 选择最长的那一条轴
Heuristic #2 Split node at location of median object 我们找中间的那个物体作为划分（保证分裂的两边三角形数量差不多）

Termination criteria

Heuristic : Stop when node contains few elements

Data Structure for BVHs

Internal nodes store

Bounding box
Children: pointers to child nodes

Leaf nodes store

Bounding box
List of objects

Nodes represent subset of primitives in scene

All objects in subtree

BVH Traversal

Intersect(Ray ray, BVH node) {
    if (ray misses node.bbox) return;

    if (node is a leaf node)
       test intersection with all objs;
       return closest intersection;

    hit1 = Intersect(ray, node.child1);
    hit2 = Intersect(ray, node.child2);

    return the closer of hit1, hit2;
}

Spatial vs. Object Partition

6.4 Basic radiometry¶

——辐射度量学

Radiometry: Measurement system and units for illumination

Accurately measure the spatial properties of light

New terms

Radiant flux: 光通量，单位时间内通过的光总量

intensity: 光强，单位立体角光通量，单位：瓦特

irradiance: 辉度，单位面积光通量，单位：瓦特/平方米

radiance: 光亮度

Perform lighting calculations in a physically correct manner.

以更加准确（贴近物理）的方法定义光线，以获得更好的结果

6.4.1 Radiant Energy and Flux¶

Definition: Radiant energy is the energy of electromagnetic radiation. It is measured in units of joules, and denoted by the symbol:

\[ Q[\text{J=Joule}] \]

Definition: Radiant flux (power) is the energy emitted, reflected, transmitted or received, per unit time.

\[ \Phi \equiv \frac{\text{d}Q}{\text{d}t}[\text{W=Watt}][\text{lm=lumen}]^* \]

Flux - #photons flowing through through a sensor in unit time

6.4.2 Radiant Intensity (I)¶

Definition: The radiant (luminous) intensity is the power per unit solid angle (立体角) emitted by a point light source.

\[ I(\omega) \equiv\frac{\text{d} \Phi}{\text{d} \omega} \]

\[ [\frac{\text{W}}{\text{sr}}][\frac{\text{lm}}{\text{sr}}\text{=cd=candela}] \]

Angles and Solid Angles(角与立体角)

Angles : ratio of subtended arc length on circle to radius

\(\theta=\frac{l}{r}\)
Circle has \(2\pi\) radians

Solid Angles : ratio of subtended area on sphere to radius squared

\(\Omega=\frac{A}{r^2}\)
Sphere has \(4\pi\) steradians

\[ \text{d}A=(r\ \text{d}\theta)(r\sin\theta\ \text{d}\phi)=r^2\sin\theta \ \text{d}\theta\ \text{d}\phi \]

微分立体角

\[ \text{d}\omega=\frac{\text{d}A}{r^2}=\sin\theta\ \text{d}\theta\ \text{d} \phi \]

Isotropic Point Source

6.4.3 Irradiance (E)¶

Definition: The irradiance is the power per (perpendicular/projected) unit area incident on a surface point.

\[ E(\pmb x)\equiv \frac{\text{d}\Phi(\pmb x)}{\text{d}A} \]

\[ [\frac{\text{W}}{\text{m}^2}][\frac{\text{lm}}{\text{m}^2}=\text{lux}] \]

注意是投影方向上的光强，所以要乘上一个cos值（Lambert’s Cosine Law)，引入Irradiance这个概念后，我们可以更好地理解之前所说的光随着距离而衰弱的问题，Intensity并未减小，而是Irradiance减小了。

Review Lambert's Cosine Law

6.4.4 Radiance (L)¶

Radiance is the fundamental field quantity that describes the distribution of light in an environment （展示光在环境中辐射的特性）

Definition: The radiance (luminance) is the power emitted, reflected, transmitted or received by a surface, per unit solid angle, per projected unit area. (单位投影面积单位立体角辐射出去的能量)

根据之前Irradiance与Intensity的定义，我们也可以定义Radiance为单位立体角的Irradiance以及单位投影面积的Intensity

Incident Radiance （入射）¶

Incident radiance is the irradiance per unit solid angle arriving at the surface.

radiance相比于irradiance具有方向性，单位面积上某个角度光辐射的能量

Irradiance vs. Radiance

Exiting Radiance（出射）¶

Exiting surface radiance is the intensity per unit projected area leaving the surface.

6.4.5 Rendering equation¶

Differential irradiance incoming: \(\text{d}E(\omega_i)=L(\omega_i)cos\theta_i\ \text{d}\omega_i\)

Differential radiance exiting (due to \(\text{d}E(\omega_i)\)): \(\text{d}L_r(\omega_r)\)

Bidirectional Reflectance Distribution Function (BRDF)¶

The Bidirectional Reflectance Distribution Function (BRDF) represents how much light is reflected into each outgoing direction \(\omega_r\) from each incoming direction

\[ f_r(\omega_i \rightarrow \omega_r)=\frac{\text{d}L_r(\omega_r)}{\text{d}E_i(\omega_i)}=\frac{\text{d}L_r(\omega_r)}{L_i(\omega_i)cos\theta_i \ \text{d}\omega_i}[\frac{1}{sr}] \]

BRDF指的是从 \(\omega_i\) 入射的光的 irradiance 反射从 \(\omega_r\) 出去分配有有多少 radiance。那么相机在 \(\omega_r\) 处所观察到的 radiance 就是所有方向打到该面积反射后的反射光 radiance 的线性叠加，得到下面的式子。

\[ L_r(p, \omega_r) = \int_{H^2} f_r(p, \omega_i \rightarrow \omega_r) L_i(p, \omega_i) \cos \theta_i \, \mathrm{d}\omega_i \]

如果物体本身也发光，那么我们可以增加 \(L_e(p,\omega_o)\) 来代表物体发出的radiance。然后又因为cos为负数时光线会被遮挡，所以只考虑上半球面，接着用法向量与方向向量的乘积代表相应的余弦值，这样我们可以得到了Rendering equation。

\[ 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\cdot\omega_i) \ \text{d}\omega_i \]

同时我们也注意到 Reflected radiance depends on incoming radiance, But incoming radiance depends on reflected radiance (at another point in the scene)

\[ 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_i\text{ d}\omega_i \]

这样我们可以解一个积分方程（第二类 Fredholm 积分方程）

\[ I(u)=e(u)+\int I(v)K(u,v)dv \]

然后构建一个线性算子\(K\)意为反射一次的过程，上式就变为\(L=E+KL\)

\[ \begin{align*} L&=E+KL\\ IL-KL&=E\\ (I-K)L&=E\\ L&=(I-K)^{-1}E\\ L&=(I+K+K^2+K^3+\dots)E\\ L&=E+KE+K^2E+K^3E+\dots \end{align*} \]

对这个式子进行物理意义进行解释：\(E\) 为直接从光源发出的光，\(KE\) 就是光源直接照到物体表面进入眼睛的光，其余的都是间接光照，经过几次反射得到的结果。光栅化只处理了前两个部分，较难处理后面部分的内容。

6.5 Monte Carlo Path Tracing¶

6.5.1 Monte Carlo Integration¶

我们需要计算一个定积分：\(\int_a^bf(x)\text{d}x\)，但是我们无法求出该函数的原函数，我们便可以通过蒙特卡洛方法求出相应的数值。

怎么做呢？通过对函数值的随机抽样求平均值来估计函数的积分

\[ \begin{aligned} \text{Random variable}&: X_i\sim p(x)\\ \text{Monte Carlo estimator}&: F_N=\frac{1}{N}\sum_{i=1}^N\frac{f(X_i)}{p(X_i)} \end{aligned} \]

蒙特卡洛估计的证明

我们说明蒙特卡洛估计值的期望为所求定积分的值

\[ \begin{aligned} E(F_N)&=E(\frac{1}{N}\sum_{i=1}^N\frac{f(X_i)}{p(X_i)})\\ &=\frac{1}{N}E(\sum_{i=1}^N\frac{f(X_i)}{p(X_i)})\\ &=\frac{1}{N}\sum_{i=1}^N E(\frac{f(X_i)}{p(X_i)})\\ &=\frac{1}{N}\sum_{i=1}^N \int_a^b\frac{f(x)}{p(x)}\cdot p(x)\\ &=\frac{1}{N}\sum_{i=1}^N \int_a^b\frac{f(x)}{p(x)}\cdot p(x)\\ &=\int_a^b f(x) \end{aligned} \]

我们取得样本越多，得到的结果越准。我们要对x积分，那么我需要对x进行采样。

6.5.2 Path Tracing¶

我们之前讲了Whitted-Style Ray tracing，

Always perform specular reflections/refraction
Stop bouncing at diffuse surface 打到漫反射的物体直接停止

但是上面对光线传播的简化合理吗？显然不，所以我们要逐步改进，并引入路径追踪算法

Problem 1: 看下图中两个茶壶，左边的材质更像镜子，我们认为它是specular或者pure specular；右边的材质更像金属-接近磨砂的质感，我们认为它是glossy的材质。那如果应用Whitted-Style ray tracing,那么是不是只有左边的茶壶显示正确，右边的不对呢？右边的茶壶之所以看起来是雾面的，是因为光线打到茶壶表面，不仅仅只是沿着镜面反射方向反射，而是会反射到镜面反射方向周围的一个区域-它多少会有些粗糙程度，这是Whitted-Style 会出问题的地方

Problem 2: Whitted-Style Ray Tracing 的第二个问题是，光线打到漫反射物体上就停了，但这是不对的，光线在漫发射表面仍然会反射-向四面八方反射。在这种不考虑光线漫反射的情况下，漫反射物体和漫反射物体之前的反射的光线也就无法考虑进来，那么在表现上会有哪些问题呢？看下面的例子：

上面两幅图描述的场景中，只有一个位于天花板上的光源，光线直接往下打。对比左右两幅图可以看到，右边的才是我们想要的结果-没有被光源直接照射的地方不应该完全是黑的，这里还有一个有意思的现象叫 color bleeding-看右边这幅图，大盒子的左面被红色的漫发射墙面反射出的光照射，显示出红色，就好像颜色从墙面跑到了盒子的表面，这种现象被称为 Color bleeding，绿色墙面同理。Color bleeding 是全局光照带来的一种效果，是Whitted-Style ray-tracing 做不到的。

由于我们之前的渲染方程是基于物理原理推导出的。

\[ 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\cdot\omega_i) \ \text{d}\omega_i \]

我们现在需要解决上述方程的半球面积分以及递归的过程，但是我们只需要数值解，所以可以用到蒙特卡洛方法。

A Simple Monte Carlo Solution¶

我们考虑简单情况，即上图简单情形下的一个点的直接光照。这个场景中有一个较大的面光源，\(\omega_0\)为观测方向，\(\omega_i\)为各个不同的入射方向。

我们假设这个点本身不发光，那么它的直接光照结果来自四面八方的热射光，于是我们直接应用反射方程：

\[ L_o(p,\omega_o)=\int_{\Omega^+}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)\ (n\cdot\omega_i) \ \text{d}\omega_i \]

我们利用蒙特卡洛方法求解这个半球上的积分。

被积函数：\(L_i(p,\omega_i)\ f_r(p,\omega_i,\omega_o)(n\cdot\omega_i)\)

概率分布函数: \(p(\omega_i)=1/2\pi\)

于是积分式可以近似地写成：

\[ \begin{align*} L_o(p,\omega_r)&=\int_{\Omega^+}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)\ (n\cdot\omega_i) \ \text{d}\omega_i\\ &\approx\frac{1}{N}\sum_{k=1}^{N}\frac{L_i(p,\omega_k)\ f_r(p,\omega_k,\omega_o)(n\cdot\omega_k)}{p(\omega_k)} \end{align*} \]

这样我们就可以有一个算法，得到任意一个点出射的radiance是多少

shade(p, wo)
    Randomly choose N directions wi~pdf
    Lo = 0.0
    For each wi
        Trace a ray r(p, wi)
        If ray r hit the light
            Lo += (1 / N) * L_i * f_r * cosine / pdf(wi)
    Return Lo

Global Illumination (全局光照)¶

如果打到的不是光源（上述方程只考虑了打到光源的情况），打到了Q点，我们只需要利用Q在此方向上的radiance

shade(p, wo)
    Randomly choose N directions wi~pdf
    Lo = 0.0
    For each wi
        Trace a ray r(p, wi)
        If ray r hit the light
            Lo += (1 / N) * L_i * f_r * cosine / pdf(wi)
        Else If ray r hit an object at q
            Lo += (1 / N) * shade(q, -wi) * f_r * cosine / pdf(wi)
    Return Lo

这样我们就支持全局光照了，但是这个方法亦存在问题

Problem 1: Explosion of #rays as #bounces go up:

Key observation: #rays will not explode if N=1.

于是我们可以仅仅随机选取一个方向（N代表的是蒙特卡洛积分中的采样次数）

shade(p, wo)
    Randomly choose ONE direction wi~pdf(w)
    Trace a ray r(p, wi)
    If ray r hit the light
        Return L_i * f_r * cosine / pdf(wi)
    Else If ray r hit an object at q
        Return shade(q, -wi) * f_r * cosine / pdf(wi)

N=1 的时候就称为Path tracing（路径追踪）

但是N=1的时候，得到的结果会不精确，我们可以通过每个像素跟踪更多路径并平均它们的 radiance。

ray_generation(camPos, pixel)
    Uniformly choose N sample positions within the pixel
    pixel_radiance = 0.0
    For each sample in the pixel
        Shoot a ray r(camPos, cam_to_sample)
        If ray r hit the scene at p
            pixel_radiance += 1 / N * shade(p, sample_to_cam)
    Return pixel_radiance

Problem 2: The recursive algorithm will never stop!

由于在现实中光线确实会一直反射，如果减少弹射次数会减少到达像素的能量，而计算机并不能无限制地计算光线弹射，为了解决这一个问题，我们用俄罗斯轮盘赌的概率思想，以一定的概率来停止光线继续传播。

之前的做法：我们总是将光线打在着色点上，并计算出着色结果\(L_o\)，现在我们可以人为设定一个概率P（0<P<1）以一定的概率P，打出一根光线，并返回结果: \(\frac{L_o}{P}\), 以概率 \(1-P\)，不打出光线，那么此时便返回 0, 这样做，我们依然可以保证期望的着色结果是 \(L_o\)

\[ E=P∗(L_o/P)+(1−P)∗0=L_o \]

shade(p, wo)
    Manually specify a probability P_RR
    Randomly select ksi in a uniform dist. in [0, 1]
    If (ksi > P_RR) return 0.0;

    Randomly choose ONE direction wi~pdf(w)
    Trace a ray r(p, wi)
    If ray r hit the light
        Return L_i * f_r * cosine / pdf(wi) / P_RR
    Else If ray r hit an object at q
        Return shade(q, -wi) * f_r * cosine / pdf(wi) / P_RR

上述方法是一个正确的路径追踪算法，但是它并不高效。

SSP: samples per pixel (就是一个像素点打出多少光线)

Sampling the Light¶

我们之前的算法是均匀地打出光线的，是否碰到光源全凭运气。如果光源比较小的情况（比如第三幅图）我们每50000条光线才能打到光源，这样会造成很大的浪费。当然为避免这种浪费，我们可以改变蒙特卡洛方法中的pdf。

如果我们可以对光源进行均匀采样

\[ pdf=\frac{1}{A}(\text{because}\int pdf\ \text{d}A=1) \]

渲染方程的积分域是在立体角上，而蒙特卡洛积分要求采样是在积分域上进行，所以如果我们想要继续使用蒙特卡洛方法，就得重新写一下渲染方程，使其积分域是在光源上，这个也是可以办到的，只需要知道\(\text{d}A\)和 \(\text{d}w\) 的关系

\[ \text{d}\omega = \frac{\text{d}A cos\theta'}{||x'-x||^2} \]

这样我们可以重写渲染方程

\[ \begin{align*} L_o(x,\omega_o)&=\int_{\Omega^+}L_i(x,\omega_i)f_r(x,\omega_i,\omega_o)\cos\theta\ \text{d}\omega_i\\ &=\int_AL_i(x,\omega_i)f_r(x,\omega_i,\omega_o)\frac{\cos\theta \cos\theta'}{||x'-x||^2}\text{d}A \end{align*} \]

Previously, we assume the light is “accidentally” shot by uniform hemisphere sampling

Now we consider the radiance coming from two parts:

light source (direct, no need to have RR)
other reflectors (indirect, RR)

shade(p, wo)
    # Contribution from the light source.
    Uniformly sample the light at x’ (pdf_light = 1 / A)
    L_dir = L_i * f_r * cos θ * cos θ’ / |x’ - p|^2 / pdf_light 

    # Contribution from other reflectors.
    L_indir = 0.0
    Test Russian Roulette with probability P_RR
    Uniformly sample the hemisphere toward wi (pdf_hemi = 1 / 2pi)
    Trace a ray r(p, wi)
    If ray r hit a non-emitting object at q
        L_indir = shade(q, -wi) * f_r * cos θ / pdf_hemi / P_RR

    Return L_dir + L_indir

还有一个问题就是计算光源贡献的时候还需要计算光线是否被遮挡

# Contribution from the light source.
L_dir = 0.0
Uniformly sample the light at x’ (pdf_light = 1 / A)
Shoot a ray from p to x’
If the ray is not blocked in the middle
    L_dir = …

部分笔记摘自

计算机图形学十二：Whitted-Style光线追踪原理详解及实现细节_whitted style二级光线-CSDN博客

评论区

对你有帮助的话请给我个赞和 star =>

欢迎跟我探讨！！！