GAMES101知识梳理：光线追踪

本节说明

本节的第一部分阴影贴图严格来说是属于光栅化的内容，阴影、全局光照等内容是光栅化较难做的，GAMES101课程也是由阴影贴图讲起引出后续的光线追踪内容。

Shadow Mapping-阴影贴图

阴影贴图是一种图像空间的做法，在生成阴影的过程中不需要知道场景中的几何信息，其阴影有走样现象。

基本思想：对于不在阴影中的点，相机可以看到该点，光源也可以看到该点。

处理步骤

将光源视为相机看向场景，也就是说从光源位置做一次光栅化，但是不进行着色，仅记录看到的点的深度；
从真正的相机位置看向场景。对于看到的点，将其投影回光源位置的相机，计算其深度。如果该点计算出的深度与上一步中记录的深度一致，那么说明该点是同时被相机与光源看到的，不在阴影中，反之则是光源看不到该点，在阴影中。

存在问题

深度是浮点数，因此深度相等的判断较困难；
从光源看向场景的光栅化记录Shadow Map也具有分辨率，在其分辨率相比渲染场景的分辨率较小的时候，阴影会出现走样问题；
经典的阴影贴图只能处理点光源，且其生成的阴影是硬阴影(Hard Shadow)。

硬阴影与软阴影
产生软阴影一定是考虑光源有一定体积或是有多个光源的情况
本影Umbra：完全看不到光源
伴影Penumbra：部分看到光源
HardAndSoftShadows

Ray Tracing-光线追踪

光栅化方法不适合计算软阴影、Glossy Reflection以及非直接光照，因此需要光线追踪。
光栅化是比较快但是比较近似的方法，而光线追踪是一种准确但是较慢的方法，因此光栅化适合实时渲染，而光线追踪常用于离线渲染。

光线追踪方法对光线的前提假设：

光沿直线传播(尽管这并不对)；
光线间不会发生碰撞(尽管这也不对)；
光路可逆，即光线从光源发射出经过反射折射等进入人的眼睛，那么从人眼出发一定可以逆向找回光源位置。

Ray Casting-光线投射

从眼睛位置发射经过屏幕像素的Eye Rays，这些射线会碰撞到场景中最近的点(这一步也同时解决了深度测试问题)，对于每个点，发出从该点到光源位置的射线(Shadow Ray)，若无遮挡说明该点可被光源照亮，结合眼睛位置和光源位置以及该点法向量，就可以进行着色。
RayCasting

Recursive(Whitted-Style) Ray Tracing

在光线投射的基础上，Whitted-Style Ray Tracing进一步假设光线在与场景物体的碰撞点都发生镜面反射，由此不断碰撞、反射，然后对每一个碰撞点进行着色计算，最终对着色结果进行加权就是该光线发射对应的像素的值。
实际中，会考虑反射折射的占比、能量消耗、最大反射次数等，防止能量一直累加。
RecursiveRayTracing

Ray-Object Intersection-光线与物体的交点

光线由其起点与方向定义

\bm{r}(t)=\bm{o}+t\bm{d}, 0 \leq t \leq \infty

光线与球面求交

IntersectionSphere
球面定义

\bm{p}:(\bm{p}-\bm{c})^2-R^2=0

球面上的交点一定在光线 $\bm{r}(t)$ 上，即

(\bm{o}-t\bm{d}-\bm{c})^2-R^2=0

解得

\begin{aligned} a&=\bm{d}\cdot\bm{d} \\ b&=2(\bm{o}-\bm{c})\cdot\bm{d} \\ c&=(\bm{o}-\bm{c})\cdot(\bm{o}-\bm{c})-R^2 \\ t&=(-b\pm \sqrt{b^2-4ac})/2a \end{aligned}

此处 $t$ 必须是正实数根。

光线与隐式表面求交

隐式表面定义

\bm{p}:f(\bm{p})=0

隐式表面上的交点一定在光线上，因此

\bm{p}:f(\bm{o}+t\bm{d})=0

解出正实数根 $t$ 。

光线与三角形网格求交

IntersectionTri
光线与三角形的交点数量只可能为0个或1个，不考虑平行的情况。

求交过程分两步：

光线与三角形所在平面求交；
判断交点是否在三角形内。

第二步判断点是否在三角形内我们已经在前面的作业中做过，就是用三个叉乘是否同号来判断，不再赘述。
第一步的光线-平面求交，我们先定义平面点法式方程

\bm{p}:(\bm{p}-\bm{p^{\prime}})\cdot\bm{N}=0

同样，将光线方程代入

(\bm{p}-\bm{p^{\prime}})\cdot\bm{N}=(\bm{o}+t\bm{d}-\bm{p^{\prime}})\cdot\bm{N}=0

解得

t=((\bm{p^{\prime}}-\bm{o})\cdot\bm{N})/(\bm{d}\cdot\bm{N})

Möller Trumbore Algorithm

这是一种光线与三角形求交的快速算法，它能同时完成上述求交过程的两步。
三角形内的点同时也是光线上的点，并且该点可以用三角形重心坐标表示，因此可得下式

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

利用克拉默法则求解出 $t,b_1,b_2$ ,(其中 $b_1,b_2\ge0$ )

\begin{bmatrix} t \\ b_1 \\ b_2 \end{bmatrix} =1/(\vec{\bm{S}}_1\cdot\vec{\bm{E}}_1) \begin{bmatrix} \vec{\bm{S}}_2\cdot\vec{\bm{E}}_2 \\ \vec{\bm{S}}_1\cdot\vec{\bm{S}} \\ \vec{\bm{S}}_2\cdot\vec{\bm{D}} \end{bmatrix}

其中

\begin{aligned} &\vec{\bm{E}}_1=\vec{\bm{P}}_1-\vec{\bm{P}}_0 \\ &\vec{\bm{E}}_2=\vec{\bm{P}}_2-\vec{\bm{P}}_0 \\ &\vec{\bm{S}}=\vec{\bm{O}}-\vec{\bm{P}}_0 \\ &\vec{\bm{S}}_1=\vec{\bm{D}}\times\vec{\bm{E}}_2 \\ &\vec{\bm{S}}_2=\vec{\bm{S}}\times\vec{\bm{E}}_1 \\ \end{aligned}

Accelerating Ray-Surface Intersection-光线与表面求交加速

Axis-Aligned Bounding Box-轴对齐包围盒

所谓包围盒(或称包围体,Bounding Volume)就是用简单的几何形体将物体包围在其中，若光线与包围盒不相交，则更不可能与其中的物体相交。
BoundingVolume

轴对齐包围体(简称AABB)是其中最为常用的包围盒，它是一个长方体包围盒，在实际使用中将其看作三组无限大对面划分出的空间的交集。
AABB

之所以使用AABB，是因为光线与划分AABB的任意一个无限大平面求交的计算很方便
WhyAA
以垂直于x轴的对面为例，相交时间 $t$ 可表示为

t=(p_x^{\prime}-o_x)/d_x

光线与AABB求交

对于每一组对面，都可以求出光线与两个面相交的两个时间 $t_{min},t_{max}$ 。
IntersectionWithAABB
其求交的逻辑就是对上述三对时间间隔求交集，通俗理解就是光线要进入三组对面划定的空间才算进入了AABB包围盒，而光线只需离开任意一对对面划定的空间就算离开了AABB包围盒。
具体做法就是取三组中最大的 $t_{min}$ 和最小的 $t_{max}$ ，即

\begin{aligned} t_{enter}&=\max\{t_{min}\} \\ t_{exit}&=\min\{t_{max}\} \\ \end{aligned}

$t$ 可能为出现负值，对于负值需要分情况讨论：

$t_{exit}<0$ ，则包围盒在光线后方，无交点；
$t_{exit}\ge0,t_{enter}<0$ ，则光线起点在包围盒内，有一个交点。

总结来说，当且仅当 $t_{enter}<t_{exit}\&\&t_{exit}\ge0$ ，光线与包围盒有交点。

利用AABB加速光线追踪

Uniform Grids-均匀网格

在求交之前，对场景进行预处理以实现加速。先建立整个场景的包围盒，对包围盒划分均匀网格，存储物体在网格中的信息。
UniformGrids
求交时，先对光线与网格进行求交，若与光线相交的网格内有物体表面，则进一步对光线与网格中的物体表面求交。
确定光线下一个相交的网格的方法类似线的光栅化。以二维为例，如果光线向右上方传播，则下一个网格一定在当前网格的上方或右方。
Ray-SceneIntersection

关于网格分辨率，在三维情况下，根据经验有 $\#cells=C*objs,C\approx27$ ，此时求交的效率最高。

Spatial Partitions-空间划分

对于物体分布比较不均匀的场景，Uniform Grids方法的加速效果不佳，这种情况下更为常用的是空间划分方法。
SpatialPartition
其中较为常用的是KD-Tree。其划分过程是将包围盒沿着轴方向循环划分，直到子包围盒为空或子包围盒内物体数量小于设定的阈值。
KD-Tree数据结构特点:

中间节点存储：
- 沿哪一个轴划分，(x,y或z轴)
- 划分的位置
- 指向子节点的指针
- 中间节点不存储场景物体
叶子节点存储：
- 节点中物体的列表

光线追踪过程就是遍历KD-Tree的过程，对KD-Tree进行前序遍历，并对光线经过的叶子节点中的物体进行求交。
KD-Tree
KD-Tree可能存在一个物体存在于多个叶子节点中，而三角形与划分面的相交计算又较难，因此KD-Tree如今应用已经逐渐减少。

Object Partitions-物体划分

物体划分是光线追踪中更为广泛使用的方法，通过这种方法划分形成的加速结构称为Bounding Volume Hierarchy(BVH)。
BVH划分的是物体，而非包围盒。如下图，先将根节点包围盒内的三角形分成两部分，对两部分三角形求各自的包围盒作为两个子节点，对两个子节点进行同样的操作，循环划分直至子包围盒内三角形数量小于阈值。
这样划分可以保证一个物体只会出现在一个叶子节点中。
BVH

BVH的划分原则：

一次选择一个维度进行划分
Heuristic #1：总是选择当前节点最长的轴进行划分
Heuristic #2：选择中间的物体作为划分的界限(即找到当前划分维度的坐标的中位数，利用快速选择算法)

BVH数据结构特点：

中间节点存储：
- 包围盒
- 指向子节点的指针
叶子节点存储：
- 包围盒
- 节点中物体的列表
- 所有物体的信息都在叶子节点中

BVH遍历算法伪代码
BVHtraversal

Spatial Partitions vs Object Partitions(BVH)

Spatial Partition(e.g. KD-Tree)：

将空间划分成互不重叠的区域
一个物体可被划分到多个区域中
Object Partitions(e.g. BVH)：
将物体划分为两个子集
计算形成的子包围盒可能有空间重叠

Basic Radiometry-辐射度量学

辐射度量学做的事就是描述光照，它准确地描述了光照的空间属性。
辐射度量学中定义的光照物理量：Radiant Flux(辐射通量), Intensity(辐射强度), Irradiance(辐照度，辐射通量密度), Radiance(辐射率)。

辐射度量学物理量

Radiant Energy and Flux(Power)-辐射能量和通量(功率)

Radiant Energy是电磁辐射的能量

Q[J=Joule(焦耳)]

Radiant Flux(Power)是单位时间发射、反射、传递、接收的能量

\begin{aligned} &\Phi \equiv \mathrm{d}Q/\mathrm{d}t \\ &[W=Watt][lm=lumen(流明)] \end{aligned}

Radiant Intensity-辐射强度

Radiant Intensity是单位立体角上的辐射功率

\begin{aligned} &I(\omega)\equiv \mathrm{d}\Phi/\mathrm{d}\omega \\ &[W/sr][lm/sr=cd=candela] \end{aligned}

其中 $sr$ 为立体角单位。
均匀辐射能量的点光源可表示如下

\begin{aligned} \Phi&=\int_{S^2}I\mathrm{d}\omega=4\pi I \\ I&=\Omega/4\pi \end{aligned}

RadiantIntensity

Solid Angle-立体角

二维中，角度等于圆的弧长除以半径，推至三维，立体角等于对应球面上的弧面面积除以半径的平方，即

\Omega=A/r^2

SolidAngle1

微分立体角
SolidAngle2

\begin{aligned} \mathrm{d}A&=(r\mathrm{d}\theta)(r\sin\theta \mathrm{d}\phi)=r^2\sin\theta \mathrm{d}\theta \mathrm{d}\phi \\ \mathrm{d}\omega&=\mathrm{d}A/r^2=\sin\theta \mathrm{d}\theta \mathrm{d}\phi \end{aligned}

对单位立体角 $\mathrm{d}\omega$ 在整个球面积分的结果为 $4\pi$ ，对应整个球面大小的立体角。

Irradiance-辐照度

Irradiance是单位面积上的辐射功率(需要投影至垂直方向，Blinn-Phong模型中的Lambert’s Cosine Law也是同理)

\begin{aligned} &E(x)=\mathrm{d}\Phi(x)/\mathrm{d}A \\ &[W/m^2][lm/m^2=lux] \end{aligned}

Irradiance

有了Irradiance的概念，我们可以更好地理解此前着色部分所述的光源到着色点的能量衰减。其衰减的实际上就是Irradiance，而Intensity实际上并无衰减,因为对于固定的立体角，随着传播距离增加其立体角大小不会变，但是立体角对于的表面积会增大。
IrradianceFalloff

Radiance-辐射率

Radiance是单位面积上，单位立体角发射、反射、传递或接受的辐射功率

\begin{aligned} &L(p,\omega)\equiv d^2\Phi(p,w)/(\mathrm{d}\omega \mathrm{d}A\cos\theta) \\ &[W/(sr\cdot m^2)][cd/m^2=lm/(sr\cdot m^2)=nit] \end{aligned}

这里的单位nit其实就是机圈常说的手机亮度多少尼特，手机亮度一般也用Radiance描述

Radiance的理解
Recall

Irradiance: Power per projected unit area.
Intensity: Power per solid angle.

Radiance: Irradiance per solid angle.
Radiance: Intensity per projected unit area.

前者可以理解成接收能量，从四面八方传递过来的能量在这个单位面积上的垂直投影后的辐射功率就是Irradiance，如果限定能量来源在单位立体角内，就变成了Radiance，可表示如下

L(p,\omega)=\mathrm{d}E(p)/(\mathrm{d}\omega\cos\theta)

IncidentIncidentRadiance

同理，后者可以理解成发射出能量，表面对特定立体角发射的垂直投影后的辐射功率可以描述为Intensity，如果限定发射能量的表面有单位面积，则变成了Radiance。

L(p,\omega)=\mathrm{d}I(\omega)/(\mathrm{d}A\cos\theta)

ExitingRadiance

Irradiance vs. Radiance

Irradiance：dA区域接收到的总辐射功率
Radiance：dA区域接收到的来自dω方向的辐射功率

则Irradiance就是将四面八方的Radiance积分起来

\begin{aligned} \mathrm{d}E(p,\omega)&=L_i(p,\omega)\cos\theta \mathrm{d}\omega \\ E(p)&=\int_{H^2}L_i(p,w)\cos\theta \mathrm{d}\omega \end{aligned}

IrradiancevsRadiance
Irradiance和Radiance的差别就在于方向性，Radiance具有方向性。

Bidirectional Reflectance Distribution Function(BRDF)-双向反射分布函数

PS：初看这部分可能很绕很晕，为了便于理解和防止我自己后续遗忘，我翻来覆去加了不少废话。总之，BRDF的内容不理解就翻来覆去地看，多参考不同大佬的解释说明，多看几遍多思考总会理解的，我第一次差不多看懂BRDF是三刷这部分内容的时候，别怕看不懂，再难理解也不要放弃思考。

BRDF描述的就是入射点收集到的能量如何被反射。换言之，BRDF定义了某点接收到的Irradiance如何在反射中在各个出射方向进行分配。
BRDF

f_r(\omega_i \rarr \omega_r)=\mathrm{d}L_r(\omega_r)/\mathrm{d}E_i(\omega_i)=\mathrm{d}L_r(w_r)/(L_i(\omega_i)\cos\theta_i \mathrm{d}\omega_i) [1/sr]

由上面函数表达式的单位也不难看出，BRDF是利用该点的Irradiance计算反射中各方向的Radiance。

注：初见BRDF或许会很懵，为什么需要Irradiance进来而Radiance出去。其实可以这样理解，在实际渲染中，对于某一个着色点我们会考虑其光照/能量来自各个方向，但是我们仅考虑该着色点朝单一方向发射的光照/能量，这单一的特定方向就是我们相机所在的方向。这样理解就会发现BRDF是很有道理且很好用的。

下面我们只专注于一个出射方向(即相机方向)。
对于这个固定的出射方向，每个入射方向都会有对应的BRDF，即每个入射方向都可能对当前固定的出射方向有所贡献，因此如果将各入射方向贡献积分起来(做半球面积分)，就能得到出射方向的Radiance。
入射的Irradiance乘以BRDF就会变成来自当前单位立体角入射的在出射方向上的Radiance，再对其积分就变成了
ReflectionEquation
反射方程

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

但是，反射方程存在一个问题，其着色点的Irradiance除了来自光源，也可以来自于其他点反射的光，也就是说任何出射的Radiance都可以作为其他点入射Irradiance的一部分。那么这变成了一个递归的问题。

The 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)\mathrm{d}\omega_i

课程中渲染方程的理解部分推导过程比较长，具体见课程Lecture 15 ppt黑底部分。简而言之，就是将渲染方程理解为递归的过程，写成积分算子的形式最后可以化出如下表示

L=E+KE+K^2E+K^3E+...

即最终的Radiance被表示成直接来自光源的能量、光源的能量一次反射后的能量、光源的能量两次反射后的能量…
而其中的直接来自光源的能量和一次反射后的能量统称直接光照，这也是我们前述的光栅化完成的事情(光栅化也能计算间接光照，但是比较麻烦)；剩下各次反射后的能量就是间接光照。直接光照与间接光照加起来就是全局光照。

Path Tracing-路径追踪

Monte Carlo Integration-蒙特卡洛积分

蒙特卡洛积分可以解决不容易用解析法求解的定积分，是一种近似方法。具体做法是在积分域中进行随机采样，对于每次采样都把函数图像当成矩形，对各次矩形面积求平均就是函数积分的近似值。
MonteCarlo
假设定积分为

\int_a^bf(x)\mathrm{d}x

使用均匀分布的随机变量X进行采样

X_i \thicksim p(x)=1/(b-a)

则蒙特卡洛积分为

F_N=(1/N)\sum_{i=1}^N[f(X_i)/p(X_i)]=[(b-a)/N]\sum_{i=1}^Nf(X_i)

蒙特卡洛积分可以使用其他概率密度函数(Probability Density Function, PDF)进行采样。采样点越多与真实结果越接近。

Path Tracing的实现

Whitted-Style Ray Tracing存在一些问题：

无法实现Glossy材质
无法实现全局光照

我们可以认为Whitted-Style Ray Tracing的结果严格来说是错误的，而使用渲染方程可以得到正确结果。
但是渲染方程在实际实现时有两个主要问题：

需要求解半球积分
渲染方程是递归的

蒙特卡洛积分求解半球积分

解决半球积分的问题，需要用到前面所述的蒙特卡洛积分方法。
MonteCarloSolution
这里先只考虑直接光照，使用半球面上均匀分布的函数进行采样
直接光照Radiance

L_o(p,\omega_o)=\int_{\Omega+}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n \cdot \omega_i)\mathrm{d}\omega_i

采样函数PDF

p(\omega_i)=1/2\pi

于是直接光照Radiance的计算可写成蒙特卡洛积分的形式

\begin{aligned} L_o(p,\omega_o)&=\int_{\Omega+}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n \cdot \omega_i)\mathrm{d}\omega_i \\ &\approx(1/N)\sum_{i=1}^N[L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n \cdot \omega_i)/p(\omega_i)] \end{aligned}

上述求解过程的伪代码为：

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)

递归求解间接光照

对于多次反射的间接光照的求解，可以使用递归的方法。以下图为例，求解如图两次反射的间接光照相当于先假设相机在P点看向Q点求出Q点直接光照的情况，此时就求出了P点处来自Q点的Radiance。
IndirectIllumination
对上一部分的直接光照伪代码进行扩充得到支持全局光照的伪代码：

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

显然，至此问题还未解决，因为上述递归计算的过程没有跳出条件，因为一个点的接收很多个来自之前反射点的Radiance，这样一来Path Tracing过程的光线数量会爆炸。
RaysExplosion
如图所示，光线数量与每次反射新产生光线数量N是指数关系，但是如果每次反射只产生一条光线，就不会产生指数爆炸的问题。
再次修改伪代码：

shade(p, wo)
  Randomly choose one directions wi~pdf
  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也是可以得到接近真实的结果的。

首先，引入经过像素的光线生成的伪代码，这个与光线追踪的Ray Casting类似：

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

关于递归的终止，引入Russian Roulette(RR)方法，即每次反射都有一定的概率停止追踪该光路。
具体而言，我们设定以概率P继续追踪当前光路，然后将下一点计算的结果除以概率P；同时有概率(1-P)不再追踪当前光路，此时得到结果为0.
在这种方法下，我们计算得到的Lo的数学期望为

E=P*(Lo/P)+(1-P)*0=Lo

也就是说，只要经由像素发射的光路足够多，最终计算结果与真实情况就会足够接近。
这也是为什么说Path Tracing是一种无偏的计算全局光照的方法。
补：通过定义的概率P我们也能求出弹射次数的期望为P/(1-P)

由此我们得到最终的shade伪代码：

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 directions wi~pdf
  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

对光源采样

每个像素打出的光路数量(Samples Per Pixel, SSP)会影响结果，越多效果越好，但是使用过高的SPP效率较低。
SPP

Low-SPP情况效果不佳的主要原因是如果光源较小，我们可能需要很多个像素采样光路才能打到光源。由此我们也可以发现，实际上前面所述的Path Tracing方法浪费了许多光路，因为很多光路实际上是不会打到光源的。

如果我们对光源进行采样而不是半球面采样，则可以避免光路浪费
由于蒙特卡洛积分要求在积分域上采样，因此需要将对光源的采样表示为在半球面对于投影区域的采样，如下图
SamplingtheLight
重写渲染方程为

\begin{aligned} L_o(p,\omega_o)&=\int_{\Omega+}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)\cos\theta \mathrm{d}\omega_i \\ &=\int_AL_i(x,\omega_i)f_r(x,\omega_i,\omega_o)(\cos\theta\cos\theta^{\prime})/ \begin{Vmatrix} x^{\prime}-x \end{Vmatrix}^2 \mathrm{d}A \end{aligned}

于是我们再次改进算法，将radiance来源分成两部分：

直接来自光源(直接光照。在光源区域积分，对光源采样，无需使用RR)
来自其他物体反射(间接光照。在半球面积分，对半球面采样，使用RR)

修改shade伪代码：

shade(p, wo)
  # Contribution from the light source
  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 = 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

Path Tracing可以得到几乎完全正确的结果，照片级的真实感(Photo-Realistic)
PhotoRealistic

一些补充说明

早期Ray Tracing在概念上基本上就是说Whitted-Style Ray Tracing，而如今我们所说的Ray Tracing往往是指所有光线追踪方法的大集合，囊括了如Path Tracing, Photon Mapping, Metropolis light transport, VCM/UPBP等等各种光线追踪方法。
课程中未具体覆盖的内容：
1. Path Tracing中如何对半球均匀采样并未提及，这个做起来并不容易，作业代码框架中直接提供了所需函数;
2. 蒙特卡洛积分允许使用任意的PDF，选择最佳的PDF需要用到重要性采样理论(Importance Sampling);
3. 生成随机数的质量对效果有影响。低差异序列(Low Discrepancy Sequences)效果就比较好;
4. 将半球面采样和光源采样结合可以产生更佳的效果，称为Multiple Importance Sampling；
5. 每个像素打出多条光线，计算结果的平均就是像素的Radiance，这个问题由Pixel Reconstruction Filter相关理论说明；
6. Radiance与颜色不是一个东西，甚至两者不是线性关系，将计算出的Radiance转换成最终显示所需的颜色需要用到伽马矫正(Gamma Correction)。