Cornell Box 里的两个现象

Cornell Box（康奈尔盒子）是测试全局光照的标准场景。画面里有两个直接光照产生不了的现象：

• 全局光照：光线在墙壁间反复弹射，角落和间接被照亮的区域也有柔和的光，不是死黑。
• 颜色溢出（Color Bleeding）：红墙旁的白色球体泛着红色——红墙反射的红光打到了白球上。这是间接光照携带颜色的直接证据。

这两个现象是判断一个渲染器是否真正实现了全局光照的标志。

参考

• Cornell box（维基百科）
• Color bleeding（计算机图形学）

路径追踪 vs 直接光照

区别

直接光照：光线打到表面后，连到光源算一次着色就结束。只算光源直接照到的部分。

路径追踪：光线打到表面后随机选一个方向继续弹射，弹到的下一个表面可能再弹，直到碰到光源或被概率终止。整条路径收集的光照贡献叠加起来，就是全局光照（Global Illumination, GI）。

一张表

维度	直接光照	路径追踪
光线弹射	1 次（到光源）	多次
间接光照	无	有
颜色溢出	无	有
计算量	小	大（需多次采样）
结果	确定性	随机（采样少有噪点）

容易看错的点

• 路径追踪 ≠ 经典光线追踪的简单延伸。Whitted 光线追踪只在镜面/折射方向追踪，漫反射处就停了；路径追踪在每个表面随机采样所有方向，这是它能算全局光照的关键。
• “随机”是必须的，不是偷懒。要积的方向是连续半球（B04 的积分维度太高），确定性遍历不可能，只能靠蒙特卡洛采样估计（B05）。
• 一条路径不够。单个像素只发一条路径方差极大（全噪点，B06），每个像素要发很多条路径求平均。

参考

• Path tracing（维基百科）
• Global illumination（维基百科）

渲染方程（The Rendering Equation）

方程

$L_o(x, \omega_o) = L_e(x, \omega_o) + \int_{\Omega} f_r(x, \omega_i, \omega_o)\, L_i(x, \omega_i)\, \cos\theta_i\, d\omega_i$

一个表面点 $x$ 沿方向 $\omega_o$ 发出的光 $L_o$，等于它自己发的光 $L_e$，加上来自整个半球 $\Omega$ 的入射光经表面反射后的贡献总和。这是几乎所有渲染算法都在近似的”大一统”公式，由 Kajiya 于 1986 年提出。

各项含义

项	含义	smallpt 里
$L_o(x, \omega_o)$	$x$ 沿 $\omega_o$ 的出射辐射度	radiance() 返回值
$L_e(x, \omega_o)$	自发光（只有光源非零）	球的 e 字段
$f_r(x,\omega_i,\omega_o)$	BRDF，入射 $\omega_i$ 散射到 $\omega_o$ 的比例	由 refl 材质类型决定
$L_i(x, \omega_i)$	从 $\omega_i$ 方向来的入射光	别的点的 $L_o$（递归）
$\cos\theta_i$	入射角余弦（Lambert 定律）	漫反射的余弦权重
$\int_\Omega d\omega_i$	对整个上半球积分	蒙特卡洛采样近似

为什么解不动

两个难点：

1. 高维积分：半球方向连续且无穷多，无法穷举。
2. 递归结构：被积项里的 $L_i$ 本质上是另一个点 $x'$ 的 $L_o$，方程套方程到无穷。

解析求解不可能，所以用蒙特卡洛方法随机采样来近似（见 B05）。

容易看错的点

• $L_e$ 只有光源非零。普通表面 $L_e = 0$，其 $L_o$ 完全来自反射积分。smallpt 里 obj.e 只有光源球不为 0。
• $\cos\theta_i$ 是 Lambert 余弦定律：同样的光斜着照，单位面积分到的能量按余弦减少。不是任意加权。
• BRDF 取决于材质：漫反射、镜面、玻璃的 $f_r$ 完全不同——这是 radiance 三个分支处理的（下集）。

参考

• Rendering equation（维基百科）
• Kajiya 1986 原论文 The Rendering Equation

蒙特卡洛积分（Monte Carlo Integration）

思想

要算积分 $\int f(x)\,dx$，不必遍历定义域每个点：随机采样 $N$ 个点，取函数值平均来估计积分。均匀采样时：

$\int_a^b f(x)\,dx \approx \frac{b-a}{N} \sum_{i=1}^{N} f(x_i)$

采样越多，估计越接近真实值。

为什么对路径追踪管用

渲染方程（B04）是对整个半球方向的积分。方向连续、无穷多，蒙特卡洛说：随机弹射几个方向，把那些方向的光照贡献平均，就是对积分的无偏估计。一次随机弹射 = 一次蒙特卡洛采样。

收敛速度：1/√N

估计的误差随采样数下降的速度是：

$\text{误差} \propto \frac{1}{\sqrt{N}}$

这是 B06 噪点的根源，也是最反直觉的地方：要让噪点减半，采样量要翻 4 倍。

目标	需要的采样量
噪点减半	$N \times 4$
噪点降到 1/10	$N \times 100$

1000 spp 看着不错，要做到 2 倍干净要 4000 spp。

无偏性

蒙特卡洛估计是无偏的：当 $N \to \infty$ 时收敛到真实积分值。路径追踪理论上没有系统误差，只要采得够多结果就是对的。smallpt 是无偏路径追踪器。

容易看错的点

• 采样少 ≠ 错，只是方差大。4 spp 全噪点不是 bug，是高方差；多采样就会消失。
• 收敛慢不是算法差，是数学本质。$1/\sqrt N$ 是蒙特卡洛的固有性质。要更快得用重要性采样。

参考

• Monte Carlo integration（维基百科）
• Monte Carlo method（维基百科）

蒙特卡洛收敛：噪点从哪来

每个像素独立地随机发射光线路径。采样少时，不同像素的估计值波动大，表现为噪点。噪点不是 bug，是尚未收敛的方差。 采样越多，对渲染方程积分的估计越准，噪点越少。

4 / 100 / 1000 采样对比

采样数 (spp)	状态
4	几乎全噪点，每像素只发 4 条路径
100	可辨认，仍有明显噪点
1000	较干净，接近最终效果

为什么不能”多采几次就立刻干净”

收敛速度是 $O(1/\sqrt{N})$（见 B05）：噪点减半要 ×4 采样。从 1000 spp 想做到 2 倍干净，需要 4000 spp；再干净一倍要 16000 spp。这是路径追踪慢的根本原因。

参考

• Monte Carlo integration — 收敛性

99 行代码的结构地图

行号	模块	本集 block
1–13	Vec 向量数学	B08
14	Ray 射线定义	B09
15–31	Sphere + 求交	B10/B11/B12
32–40	场景：9 个球	B14/B15
41–43	工具函数	B16
44–73	radiance 核心（下集）	—
74–99	main：相机+采样（下集）	—

本集讲基础设施（数据结构 + 求交 + 场景），下集讲灵魂（radiance 材质物理 + main 相机采样）。

参考

• smallpt 源码

Vec 结构体：三合一

Vec 只有 double x, y, z，但同时表示三种东西：

1. 空间位置（点坐标）
2. 方向（射线方向、法线）
3. 颜色（RGB 三通道，正好三个数）

这种复用让代码极度精简。代价是可读性（看到 Vec 要猜它代表什么）。

两种乘法要分清

Vec operator*(double b) const        { return Vec(x*b, y*b, z*b); }       // 标量乘
Vec mult(const Vec &b) const         { return Vec(x*b.x, y*b.y, z*b.z); } // 分量乘

运算	用途	场景
operator*(double)	标量乘（向量缩放）	射线 $o + t \cdot d$
mult(Vec)	分量相乘（Hadamard 积）	颜色混合：光打到红墙，反射 = 入射光 × 红色

mult 专门给颜色用：RGB 三通道各自相乘，模拟”表面吸收/反射各通道多少光”。

容易看错的点

• % 被重载为叉积，不是取模。这在正经项目里是大忌（违反直觉），99 行为了省掉 cross() 函数名这么干。读到 % 要反应过来是叉积。
• norm() 原地修改并返回引用：把自身除以长度，调用后原向量被改，别指望保留原值。

参考

• 向量点积 / 叉积（维基百科）
• Hadamard 积（分量乘）

Ray 结构体

struct Ray {
  Vec o, d;   // 起点(origin)、方向(direction)
};

参数方程

射线上任意一点用参数 $t$ 表示：

$P(t) = o + t \cdot d \quad (t \ge 0)$

• $t = 0$：起点 $o$
• $t$ 越大：离起点越远
• $t < 0$：射线反方向，不用

求交的本质

找最小的正 $t$ 使 $P(t)$ 落在物体表面。对球体代入球面方程解出 $t$（B11）；对场景遍历所有物体取最近 $t$（B16）。

容易看错的点

• 方向 d 必须归一化。否则求交公式 $a = d \cdot d = 1$ 的简化不成立（B11）。smallpt 在 main 里保证 r.d.norm()。
• $t$ 是距离当且仅当 d 归一化。若 d 非单位，$t$ 只是参数不是真实距离，影响 eps 判断和深度比较。

参考

• 射线（计算机图形学）

Sphere 结构体

enum Refl_t { DIFF, SPEC, REFR };  // 漫反射 / 镜面 / 折射

struct Sphere {
  double rad;       // 半径
  Vec p, e, c;      // 球心(p)、自发光(e)、表面颜色/反照率(c)
  Refl_t refl;      // 材质类型
};

字段	含义	备注
rad	半径	超大值(1e5)用来当墙(B15)
p	球心位置	—
e	自发光颜色	只有光源非零，对应渲染方程 $L_e$（B04）
c	表面颜色（反照率 albedo）	决定吸收/反射多少光，用 mult 混色（B08）
refl	材质类型	决定 BRDF 分支（下集核心）

容易看错的点

• e（自发光）与 c（反照率）别混淆。e 只有光源球非零，表示主动发光；c 是所有表面都有的反照率，表示反射光比例。渲染方程 $L_o = L_e + \int ...$ 里分别对应两项。
• 反照率取值 0–1：0 全吸收（黑），1 全反射。红墙 (0.75, 0.25, 0.25) 反射 75% 红光、25% 绿、25% 蓝。

参考

• Albedo（反照率，维基百科）
• BRDF（双向反射分布函数）

射线-球体求交

怎么求交

球面上任意点到球心 $C$ 的距离等于半径 $r$：

$\|P - C\| = r$

射线上的点用参数 $t$ 表示（$o$ 为起点，$d$ 为单位方向）：

$P(t) = o + t \cdot d \quad (t \ge 0)$

求交就是找最小的正 $t$ 使 $P(t)$ 落在球面上。把射线方程代入球面方程的平方形式 $\|P - C\|^2 = r^2$：

$\|o + t d - C\|^2 = r^2$

令 $oc = o - C$，展开点积 $(oc + td)\cdot(oc + td)$：

$(d \cdot d)\, t^2 + 2(oc \cdot d)\,t + (oc \cdot oc) - r^2 = 0$

这是关于 $t$ 的一元二次方程 $at^2 + bt + c = 0$。

判别式决定交点个数

判别式 $\Delta = b^2 - c$：

$\Delta$	几何意义	交点
$< 0$	射线偏离球	无
$= 0$	射线与球相切	1 个
$> 0$	射线穿过球	2 个（进/出）

两个解 $t_{1,2} = -b \mp \sqrt{\Delta}$，取较小的正根就是最近的交点。

smallpt 代码里 det 为什么是 Δ/4

这是读 smallpt 源码最容易卡住的地方。标准二次方程 $at^2 + Bt + C = 0$ 的判别式是 $B^2 - 4AC$。smallpt 做了两步简化：

1. $d$ 归一化使 $a = d \cdot d = 1$，消去 $a$
2. 令 $b = \tfrac{B}{2} = oc \cdot d$（取半系数）

于是代码里的 det = b*b - op%op + rad*rad 实际算的是 $\Delta/4$，sqrt(det) 是 $\sqrt{\Delta}/2$。求根公式相应变成 $t = b \mp \sqrt{det}$——省掉了因子 2 和 4。这是几何法射线-球求交的经典优化，数值与标准公式等价。

double intersect(const Ray &r) {
  Vec op = p - r.o;
  double b = op % r.d;                  // % 是点积; b = B/2
  double det = b*b - op%op + rad*rad;   // det = Δ/4
  if (det < 0) return 0;                // 判别式 < 0，未命中
  else det = sqrt(det);                 // 现在 det = √Δ/2
  double t = b - det;                   // 先试较近的解
  if (t > eps) return t;
  t = b + det;                          // 否则试较远的解
  if (t > eps) return t;
  return 0;
}

eps 防自相交

eps = 1e-4 丢弃过近的交点。从某个表面弹射出的新射线，因浮点误差，起点可能”略陷进”球内，导致算出极小的正 $t$，误判为”立刻又命中自己”。eps 把这些虚假交点排除掉。取值是权衡：太小挡不住自相交，太大漏掉真实近处交点，$10^{-4}$ 是渲染领域的常用经验值。

易错点

• 方向向量必须归一化。否则 $a = d \cdot d \ne 1$，上述 det = Δ/4 的简化不成立，要用完整求根公式。
• det 是 $\Delta/4$，不是判别式本身，也不是 $\sqrt\Delta$。代码先判 det < 0（等价判别式 < 0），再 sqrt(det) 转距离。读 smallpt 时这个命名最容易误解。

参考

• 射线-球求交（维基百科）
• smallpt 官方页面与源码

求交代码逐行解读

B11 推导了公式，这里是代码落地，重点在一个容易被忽略的小数 eps。

代码

double intersect(const Ray &r) {
  Vec op = p - r.o;
  double b = op % r.d;             // % 是点积；b = B/2
  double det = b*b - op%op + rad*rad;
  if (det < 0) return 0;           // 判别式 < 0 → 未命中
  else det = sqrt(det);
  double t = b - det;              // 先试近解
  if (t > eps) return t;
  t = b + det;                     // 再试远解
  if (t > eps) return t;
  return 0;
}

逐行对应 B11 的推导：op→$oc$、b→$B/2$、det→$\Delta/4$（sqrt(det) 后为 $\sqrt\Delta/2$）。算法本质：解二次方程，取较小的正根。

eps 与自相交

eps = 1e-4，解决一个麻烦：从某个表面弹射出的新射线，因浮点误差，起点可能”略陷进”球内，于是算出极小的正 $t$（如 $10^{-7}$），误判为”立刻又命中自己”。

eps 丢弃 $t < 10^{-4}$ 的解，只接受明显在前方的交点。这就是为什么先试 b - det（近解），不行才试 b + det（远解）。取值是权衡：太小挡不住自相交，太大漏掉真实近处交点，$10^{-4}$ 是渲染领域的常用经验值。

数值稳定性

smallpt 用 det = b*b - op%op + rad*rad 这种形式。当 rad 很大（如当墙的超大球，半径 1e5）时，b*b 和 op%op 都是巨数相减，会损失有效数字。教学/demo 级别够用；生产级渲染器会用更稳的求根公式。

参考

• 浮点精度（IEEE 754 灾难性抵消）

9 个球拼出 Cornell Box

角色	半径	说明
左/右/后/前/上/下墙	1e5	6 个超大球，左红右蓝余白色
镜面球	16.5	小球，SPEC 镜面反射
玻璃球	16.5	小球，REFR 折射
光源	—	自发光 e 非零

6 个超大球构成房间，2 个小球是主角，1 个是光源。

前墙为什么是黑色

前墙（朝相机那面）颜色设为黑色（不发光不反光），相当于这面墙不存在——房间”开口”朝向相机，相机才能看到里面。

为什么用 Cornell Box

1984 年 Cornell 大学提出的标准测试场景，专门检验全局光照算法。简单到能手工验证，又能暴露间接光照、颜色溢出等关键现象（B01）。几乎所有路径追踪教程都用它。

参考

• Cornell box（维基百科）

大球当墙：局部近似平面

技巧

球面是弯的，墙是平的，怎么用球当墙？半径趋近无穷大时，球面局部接近平面。

用左墙的数字算一遍

左墙球：球心 x = 100001，半径 = 100000。

球面上最靠近房间的点 = 球心朝房间方向退一个半径：

$x_{最近} = 100001 - 100000 = 1$

这个巨大的球从右侧”顶住”了 $x = 1$ 这个平面位置。

为什么看起来是平的

房间宽约 98，相对半径 $10^5$ 只占：

$\frac{98}{100000} \approx 0.1\%$

半径 10 万的球上取 0.1% 的小弧段，曲率几乎为 0，肉眼就是平面。类比：地球半径 6371 km，地表几公里范围感觉是平的——同理。

右墙对称：球心 x = -99901，半径 = 100000，$x_{最近} = -99901 + 100000 = 99$。两面墙形成宽 98 的房间。

工程意义

用大球当墙，避免了实现”平面求交”代码——所有几何统一用球体求交，99 行才装得下。

容易看错的点

• 球心位置不是墙的位置。左墙球心在 $x = 100001$（很远），墙实际在 $x = 1$（球面最近点）。
• 半径不是越大越好。1e5 够平且数值尚稳；用 1e10 会因巨数相减丢有效位（IEEE 754 灾难性抵消）。

参考

• Loss of significance（浮点灾难性抵消）

全局 intersect：找最近的交点

inline bool intersect(const Ray &r, double &t, int &id) {
  int n = sizeof(spheres)/sizeof(Sphere);
  double d, inf = t = 1e20;
  for (int i = n; i--;) {                 // 遍历所有球
    if ((d = spheres[i].intersect(r)) &&  // 该球有交点
        d < t) {                           // 且更近
      t = d; id = i;
    }
  }
  return t < inf;
}

算法：暴力遍历取最小

t 初始化为 $10^{20}$（当无穷），对每个球调它的 intersect（B12），命中且更近就更新 t 和 id。循环结束 t 是最近交点距离，id 是那个球的索引。

为什么暴力够用

只有 9 个球，遍历成本可忽略。工业级场景有上百万三角形，才需要加速结构（BVH / kd-tree）避免逐个测试。smallpt 场景小到暴力最优。

容易看错的点

• d && d<t 的短路：intersect 返回 0 表示未命中，d 为 0 即 false，&& 短路不比较距离。smallpt 把”命中判断”和”距离比较”压在一行。
• t 用引用传出：调用方拿到最近距离，用于算交点 x = r.o + r.d * t，id 用于查该球属性（颜色、材质）。

参考

• Bounding volume hierarchy（BVH 加速结构）