radiance() 是渲染方程的递归求解器

radiance 递归求解的是某条射线的出射辐射度（outgoing radiance）。它把光传输方程（Light Transport Equation, LTE）拆成五步代码：

Vec radiance(const Ray &r, int depth){
  double t;                               // 交点参数
  int id = 0;
  if (!intersect(r, t, id)) return Vec(); // 1. 未命中 → 环境黑
  const Object &obj = objects[id];
  Vec x = r.o + r.d*t;                    // 2. 交点
  Vec n = (x - obj.p).norm();             //    几何法线
  Vec nl = n.dot(r.d) < 0 ? n : n*-1;     //    朝向射线来源的定向法线
  Vec f = obj.c;                          //    反照率 albedo
  double p = f.x>f.y && f.x>f.z ? f.x : f.y>f.z ? f.y : f.z;
  if (++depth > 5) if (erand48(Xi) < p)  // 3. Russian Roulette
    f = f * (1/p); else return obj.e;
  if (obj.refl == DIFF) { /* 漫反射 */ }
  else if (obj.refl == SPEC) { /* 镜面 */ }
  else { /* 折射：Snell + Fresnel */ }
}

五步与光传输方程的对应

把场景表面的反射拆成”自发光 + 入射光经 BRDF 调制后的积分”，得到光传输方程（沿某方向 $\omega_o$ 的出射辐射度）：

$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$

代码步骤	方程中的角色
intersect 命中 → $x$	积分点 $x$
返回 obj.e	自发光项 $L_e$（光源才有）
if obj.refl == ...	选择 BRDF $f_r$
对半球采样方向 $\omega_i$	蒙特卡洛估计积分
递归 radiance(reflectionRay, depth)	$L_i = L_o(\text{邻居点})$

为什么不用迭代而用递归

渲染方程的积分里嵌着”别的点的出射辐射度” $L_i$，这是自引用结构。递归恰好是自引用的自然表达：每次 radiance 调用计算一个点的 $L_o$，其中需要邻居点的 $L_o$，于是再调一次 radiance。理论上无限递归，实践里靠 Russian Roulette 截断（见 B04）。

易错点

• return Vec() 不是”返回黑色光源”，而是”射向无穷远的天空（无环境光）“。smallpt 没有 ambient light，所有亮度只来自场景内的发射体（obj.e），这是物理正确的全局光照。
• depth 从 0 开始。相机射线打到地板是 depth 0；++depth>5 意味着第 6 次弹射才进 Russian Roulette，前 5 次强制继续。

参考

• Rendering equation (Kajiya 1986)
• Path tracing
• smallpt 源码 (Kevin Beason)

定向法线（oriented normal）：为什么要翻转

smallpt 命中球体后算两套法线：

Vec x = r.o + r.d*t;
Vec n  = (x - obj.p).norm();            // 几何法线：球心 → 交点
Vec nl = n.dot(r.d) < 0 ? n : n*-1;     // 定向法线：指向射线来源侧

几何法线 vs 定向法线

对球体，几何法线 $n = (x - C)/r$ 永远朝外。但射线可以从外面打进（空气→玻璃），也可以从内部射出（玻璃→空气）。这两种情形法线应该指向射线来源那一侧才方便后续 BRDF 计算，于是引入定向法线 $n_l$：

$n_l = \begin{cases} n, & n \cdot d < 0 \quad(\text{射线在外表面，迎面打来}) \\ -n, & n \cdot d > 0 \quad(\text{射线从内部射出}) \end{cases}$

判据是点积符号：$n \cdot d < 0$ 表示法线与射线方向相反（迎面），符合需要；否则翻转。

为什么这一步对折射致命

Snell 定律 $n_1 \sin\theta_1 = n_2 \sin\theta_2$ 里，$\theta_1$ 是入射方向与界面外法线的夹角。计算 $\cos\theta_1 = -d \cdot n_l$ 必须用朝向射线来源的法线，否则符号错乱：

• 用 $n$：射线从内部出射时 $d \cdot n > 0$，$\cos\theta_1 < 0$，没有物理意义。
• 用 $n_l$：恒有 $\cos\theta_1 = -d \cdot n_l \ge 0$，对应一个真实的入射角。

漫反射虽然对法线方向不敏感（余弦加权采样只需一个朝外的轴），但折射分支离开正确的 $n_l$ 就会得出荒谬结果（如算出 $\sin\theta_2 > 1$ 的假全反射）。

易错点

• 法线方向”朝外”是几何属性，“朝向射线来源”是渲染约定。很多引擎（如 Mitsuba、pbrt）显式区分 geometryNormal 和 shadingNormal，smallpt 用一个 nl 变量同时承担两个角色。
• 不要把 $n_l$ 误写为 n * sign(n·d)。当 $n \cdot d < 0$ 时应当保持 $n$（不翻转），是 n.dot(r.d) < 0 ? n : n*-1，符号容易写反。

参考

• Surface normal (维基百科)
• Orientation conventions in ray tracing

Russian Roulette：无偏的递归终止法

路径追踪的光线理论上无限弹射。最朴素的终止是限制最大深度 $D_{\max}$，但早于 $D_{\max}$ 终止的路径贡献被丢弃，期望值改变——这就是截断偏差（truncation bias），画面会偏暗、漏掉长路径的间接照明。

轮盘赌思想

每次弹射前掷”骰子”：以概率 $p$ 存活继续递归，以概率 $1-p$ 终止。存活时把贡献除以 $p$ 补偿。

为什么除以 p 保持无偏

设某条路径后续的真实贡献为 $L$。轮盘赌估计器 $\hat{L}$ 是随机变量：

$\hat{L} = \begin{cases} L / p, & \text{概率 } p \text{（存活）} \\ 0, & \text{概率 } 1-p \text{（终止）} \end{cases}$

取期望：

$E[\hat{L}] = p \cdot \frac{L}{p} + (1-p) \cdot 0 = L$

期望恰好等于真实贡献，故无偏。代价是方差增大：

$\mathrm{Var}(\hat{L}) = E[\hat{L}^2] - L^2 = \frac{L^2}{p} - L^2 = L^2 \cdot \frac{1-p}{p}$

$p$ 越小方差越大（少数存活路径被放大）。所以 $p$ 不能取太小，否则收敛慢、噪点多。

p 取颜色的最大分量

double p = f.x>f.y && f.x>f.z ? f.x : f.y>f.z ? f.y : f.z;
if (++depth > 5)
  if (erand48(Xi) < p) f = f * (1/p);   // 存活：补偿
  else return obj.e;                    // 终止：只返回自发光

$p = \max(f_r, f_g, f_b)$。颜色越暗（$p$ 小）越易终止——物理直觉是暗表面吸收更多光，反射出去的能量本来就少，提前终止的”相对误差”更可控。这也是一个路径重要性的近似：亮表面上的路径贡献大，应该多保留。

注意 f = f * (1/p) 补偿后，递归下一层用的 f 已经被放大了 $1/p$ 倍，整个路径的连乘积里就嵌入了这层补偿，不需要显式记账。

depth > 5 才启用

前 5 次弹射无条件执行，第 6 次起才进轮盘赌。原因：相机直接看到的前几跳贡献最大（直接照明 + 一次间接），如果一开始就轮盘赌会大量丢失这些”高权重”路径，导致整图偏暗。5 是 smallpt 的经验值；pbrt 等产品级渲染器通常配 3–10。

易错点

• erand48(Xi) < p 是”存活”，不是”终止”。erand48 返回 $[0,1)$ 均匀分布，落点小于 $p$ 的概率正是 $p$。读代码时这个比较方向容易看反。
• 补偿系数 $1/p$ 乘在 f 上，不是乘在最终结果上。因为整条路径贡献是各表面反照率的连乘积（见 B19），在某一层把 f 放大 $1/p$ 就等价于把该层补偿，递归自然传播。

参考

• Russian roulette (Monte Carlo) — 维基百科
• pbrt: Russian Roulette (Chapter 13)

Lambertian BRDF 与重要性采样抵消

漫反射（Lambertian）表面把入射光均匀散射到上半球，BRDF 是常数：

$f_r(x, \omega_i, \omega_o) = \frac{\rho}{\pi}$

其中 $\rho$ 是反照率（albedo，0–1 的 RGB）。代回光传输方程（仅反射项）：

$L_r(x, \omega_o) = \int_{\Omega^+} \frac{\rho}{\pi}\, L_i(x, \omega_i)\, \cos\theta_i\, d\omega_i$

把常数 $\rho/\pi$ 提出：

$L_r = \frac{\rho}{\pi} \int_{\Omega^+} L_i(\omega_i)\, \cos\theta_i\, d\omega_i$

蒙特卡洛估计与重要性采样

蒙特卡洛估计 $\int g(\omega)\,d\omega$ 需要一个采样概率密度 $p(\omega)$，估计量为 $g(\omega)/p(\omega)$ 的样本均值。关键选择：让 $p(\omega)$ 与被积函数形状一致，估计器就大幅简化。

这里被积函数核心是 $L_i \cdot \cos\theta$。选余弦加权分布：

$p(\omega_i) = \frac{\cos\theta_i}{\pi}$

（它确实是 $\Omega^+$ 上的合法 pdf：$\int_{\Omega^+}\frac{\cos\theta}{\pi}d\omega = \frac{1}{\pi}\int_0^{2\pi}\int_0^{\pi/2}\cos\theta\sin\theta\,d\theta\,d\phi = 1$。）

代入估计器：

$L_r \approx \frac{\rho}{\pi} \cdot \frac{L_i(\omega_i)\cos\theta_i}{\cos\theta_i / \pi} = \rho \cdot L_i(\omega_i)$

π 和 cosθ 全部抵消

分子里的 $\cos\theta_i$ 和分母里的 $\cos\theta_i$ 抵消，$\pi$ 和 $\pi$ 抵消，只剩 $\rho \cdot L_i$。这就是为什么漫反射分支代码里看不到除以 π：

return obj.e + f.mult(radiance(Ray(x, d), depth));
//            ^^^       ^^^^^^^^^^^^^^^^^^^^^^^^^^^
//          自发光     ρ × L_i(采样方向)

如果用均匀半球采样（$p = 1/2\pi$ 常数），就必须显式写 $\cos\theta/\pi$，且方差大得多——因为采样方向经常落到 $\cos\theta$ 小的”掠射”区域，贡献小却被同等权重采样，浪费样本。

易错点

• ”$\pi$ 消失” 不是因为公式写错，而是因为采样分布 $\cos\theta/\pi$ 与 BRDF 的 $1/\pi$ 配对抵消。换采样策略（均匀、BRDF 重要性等）公式形态会不同。
• BRDF 是 $\rho/\pi$ 不是 $\rho$。很多人误把反照率当 BRDF；少了 $1/\pi$ 会导致能量不守恒（射出多于射入）。$\pi$ 来自半球立体角积分 $\int \cos\theta\,d\omega = \pi$，它保证完美漫反射体对全方向均匀入射光反射 100%。

参考

• Lambertian reflectance (维基百科)
• Importance sampling (PBRT)

逆变换采样生成 cosθ 加权方向

目标是生成 $\omega \in \Omega^+$，pdf $p(\omega) = \cos\theta / \pi$。逆变换采样法（inverse transform sampling）的步骤：从极坐标 pdf 求边缘/条件 cdf，反演出 $\phi, \theta$ 作为均匀随机数的函数。

推导

球面微元 $d\omega = \sin\theta\,d\theta\,d\phi$，于是 pdf 用 $(\theta,\phi)$ 表示为：

$p(\theta,\phi) = \frac{\cos\theta}{\pi}\sin\theta$

它可分离：方位角 $\phi$ 在 $[0,2\pi)$ 均匀（$\cos\theta$ 不依赖 $\phi$），极角的边缘 pdf：

$p(\theta) = \int_0^{2\pi} \frac{\cos\theta\sin\theta}{\pi}\,d\phi = 2\cos\theta\sin\theta = \sin(2\theta), \quad \theta \in [0,\pi/2]$

求 cdf 并反演：

$P(\theta) = \int_0^\theta \sin(2\theta')\,d\theta' = \frac{1 - \cos(2\theta)}{2} = \sin^2\theta$

令 $P(\theta) = \xi_2$（均匀随机数），解 $\sin^2\theta = \xi_2$：

$\cos\theta = \sqrt{1 - \xi_2}$

方位角 $\phi = 2\pi\xi_1$。转笛卡尔坐标（局部系，z 朝法线）：

$\begin{aligned} x &= \cos\phi\,\sin\theta = \cos(2\pi\xi_1)\sqrt{1-\cos^2\theta} = \cos(2\pi\xi_1)\sqrt{\xi_2}\\ y &= \sin\phi\,\sin\theta = \sin(2\pi\xi_1)\sqrt{\xi_2}\\ z &= \cos\theta = \sqrt{1 - \xi_2} \end{aligned}$

对应到代码

double r1 = 2*M_PI*erand48(Xi);   // φ = 2πξ₁
double r2 = erand48(Xi);          // ξ₂
double r2s = sqrt(r2);            // sinθ = √ξ₂
double x = cos(r1)*r2s;           // 局部 x
double y = sin(r1)*r2s;           // 局部 y
double z = sqrt(1-r2);            // cosθ = √(1-ξ₂)

注意 z = sqrt(1-r2) 不是笔误：$\cos^2\theta + \sin^2\theta = 1$，$\sin^2\theta = \xi_2$ 故 $\cos\theta = \sqrt{1-\xi_2}$。r2s（即 $\sin\theta$）和 1-r2（即 $\cos^2\theta$）共用了同一个 $\xi_2$，所以两个分量由一个随机数生成——这是分离分布的必然结果。

为什么不直接均匀采样半球

均匀半球 pdf 是常数 $1/2\pi$，采到的方向各向同性。但渲染方程里被积函数带 $\cos\theta$，靠近法线（$\theta$ 小）的方向贡献大。均匀采样会把同样多的样本浪费在掠射角（$\theta\to\pi/2$，$\cos\theta\to 0$ 贡献几乎为零）方向。余弦加权采样让样本密度随贡献增长，方差显著降低——这就是 B05 里 $\cos\theta$ 能抵消的根源。

易错点

• $\xi_2$ 同时决定 $\sin\theta$ 和 $\cos\theta$，不能用两个独立随机数。pdf 是 $\sin(2\theta)$，分离后 $\theta$ 只剩一个自由度。
• z = sqrt(1-r2) 必须开方。常见错误是写 z = 1-r2，那对应 $\cos\theta = 1-\xi_2$（线性分布），整个 pdf 全错，重要性采样失效。

参考

• Inverse transform sampling (维基百科)
• Cosine-weighted hemisphere sampling (lighthouse3d)

ONB：以法线为 z 轴的正交基

B06 生成的方向 $(x,y,z)$ 是局部坐标系下的（z 朝法线）。要变成世界坐标，需要一组以法线为 z 轴的正交基 $\{u, v, w\}$（orthonormal basis, ONB），然后用线性组合：

$d_{\text{world}} = x\,u + y\,v + z\,w$

smallpt 的快速 ONB 构造（Frisvad 2012 风格的极简版）

Vec w = nl;                          // w = 法线
Vec u = ( (fabs(w.x) > .1 ? Vec(0,1,0) : Vec(1,0,0)) % w).norm();  // u = a × w
Vec v = w % u;                       // v = w × u

% 是 smallpt 重载的叉积。思路：

1. $w = n_l$（目标 z 轴）。
2. 取一个不与 $w$ 共线的辅助向量 $a$：若 $|w_x| > 0.1$ 取 $a=(0,1,0)$，否则取 $a=(1,0,0)$。
3. $u = \mathrm{normalize}(a \times w)$：得到一条与 $w$ 正交的切向。
4. $v = w \times u$：右手定则补出第三轴，自动正交且单位长（因为 $u, w$ 都已单位化且正交）。

为什么辅助向量要”挑剔”

叉积 $a \times w$ 当 $a \parallel w$ 时为 $\mathbf{0}$，归一化会得到 NaN。选 $a$ 时绕开 $w$ 的主分量：

• 若 $|w_x|$ 较大，说明 $w$ 接近 x 轴，取 $a=(0,1,0)$（沿 y 轴）保证不平行。
• 否则 $w$ 接近 yz 平面，取 $a=(1,0,0)$ 安全。

阈值 0.1 是经验值：足够小保证大多数情况选 $(1,0,0)$（数值稳定），又能在 $w$ 真的接近 x 轴时切换。正确性不依赖具体阈值，只依赖”$a$ 不与 $w$ 共线”——但若阈值过小，叉积结果长度趋近 0，归一化时浮点误差被放大。

为什么这一步等价于”旋转矩阵”

ONB $\{u, v, w\}$ 作为列组成的 $3\times3$ 矩阵 $R = [u\,|\,v\,|\,w]$ 是正交矩阵（$R^{-1} = R^\top$）。局部坐标 $(x,y,z)$ 到世界坐标的变换就是 $R \cdot (x,y,z)^\top = xu + yv + zw$。所以这套叉积构造等价于求一个把局部 z 轴对齐到法线 $n_l$ 的旋转矩阵，只是不开方、不调矩阵库，三步算术搞定。

易错点

• $v = w \times u$ 不是 $u \times w$。顺序决定右手系还是左手系；写反了 $v$ 会反向，采样方向镜像，漫反射”看起来对”但和真实几何不匹配（特定表面会出现偏色）。
• % 在 smallpt 是叉积不是取模。Vec 类把 % 重载为 cross，* 重载为点积——和数值含义完全相反，读源码极易混淆。

参考

• Orthonormal basis (维基百科)
• Building an orthonormal basis from a 3D vector (Frisvad 2012)

漫反射分支：6 行实现物理正确的间接光

// obj.refl == DIFF
double r1 = 2*M_PI*erand48(Xi), r2 = erand48(Xi), r2s = sqrt(r2);
Vec w = nl, u = ((fabs(w.x)>.1?Vec(0,1,0):Vec(1,0,0))%w).norm(), v = w%u;
Vec d = Vec(cos(r1)*r2s, sin(r1)*r2s, sqrt(1-r2));     // 局部余弦加权方向
d = (u*d.x + v*d.y + w*d.z).norm();                    // → 世界坐标
return obj.e + f.mult(radiance(Ray(x, d), depth));
//     ^^^^   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
//   自发光   ρ · L_i(ω_i)

每一步对应渲染方程的哪一项

回到 B05 推出的简化结果 $L_r \approx \rho \cdot L_i(\omega_i)$：

代码	对应
r1, r2 采样	生成 $\omega_i \sim p(\omega)=\cos\theta/\pi$
u, v, w ONB	局部→世界的正交变换 $R$
d = u*dx + v*dy + w*dz	$\omega_i$ 在世界系下的方向
f.mult(radiance(...))	$\rho \cdot L_i(\omega_i)$，f 即反照率
obj.e +	加上自发光 $L_e$（光源才有，否则为 0）

递归调用算的是 L_i

radiance(Ray(x, d), depth) 从交点 $x$ 沿新方向 $d$ 发出射线，命中下一个表面，返回那个点的 $L_o$——对当前点而言正是入射辐射度 $L_i(\omega_i)$。自引用的渲染方程被一行递归吃掉。

为什么这是无偏蒙特卡洛估计

由 B05：估计器 $\hat{L}_r = \frac{\rho}{\pi}\cdot\frac{L_i\cos\theta}{p(\omega)}$，代入 $p=\cos\theta/\pi$ 得 $\rho L_i$。这是单样本无偏估计：$E[\hat{L}_r] = L_r$。整张图每个像素发射 $N$ 条这样的路径取平均，由大数定律收敛到真实 $L_r$，没有系统偏差（除 Russian Roulette 引入的方差，但期望仍正确）。

易错点

• f.mult 是按分量相乘（RGB 各通道），不是点积。f 是反照率（颜色），radiance(...) 返回入射光颜色，两者逐通道相乘才是物理正确的”白光被红墙染红”。
• d 必须归一化。构造后 .norm() 不能省——递归下一层的射线求交假设方向是单位向量（见 smallpt1 B11 的 $a=d\cdot d=1$ 简化）。不归一化会让所有距离量级错乱。

参考

• smallpt 源码 (Kevin Beason)
• Path tracing & the rendering equation (scratchapixel)

镜面反射方向公式

理想的镜面（perfect specular）反射遵循反射定律：入射角 = 反射角，且入射光、反射光、法线共面。给入射方向 $d$（指向表面）和单位法线 $n$，反射方向：

$r = d - 2(d \cdot n)\, n$

推导：分解为法向 + 切向分量

把 $d$ 分解为沿法线的分量 $d_\perp = (d \cdot n)\,n$ 和切向分量 $d_\parallel = d - (d \cdot n)\,n$。反射时切向不变、法向翻转：

$r = d_\parallel - d_\perp = \bigl(d - (d\cdot n)n\bigr) - (d\cdot n)n = d - 2(d\cdot n)n$

直觉：法向分量”翻过去”，切向分量”保留”，相当于关于切平面做镜像。

smallpt 代码

// obj.refl == SPEC
return obj.e + f.mult(radiance(Ray(x, r.d - n*2*n.dot(r.d)), depth));
//                             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
//                             r = d - 2(d·n)n

这里用几何法线 n（不是 nl）也行，因为公式对 $n \to -n$ 对称：$(d\cdot(-n))(-n) = (d\cdot n)\,n$。但下一层的 nl 在命中新表面后仍要正确计算。

为什么镜面不需要随机采样

镜面 BRDF 是Dirac delta $f_r = \rho \cdot \delta(\omega_o - r(\omega_i))$——反射能量集中在唯一一个方向（反射方向 $r$），其余方向严格为零。蒙特卡洛采样这个 delta 函数最优策略就是直接采样反射方向（pdf 在该方向为 1，其余为 0），无需随机。

对比漫反射（B05）的常数 BRDF——能量散布整个半球，必须随机采样方向。镜面是”确定性”的极端：所有样本都走同一方向。

三种材质在数学上的统一

它们其实是同一族 BRDF 的不同特例：

材质	BRDF $f_r$	采样
漫反射	$\rho/\pi$（常数）	余弦加权随机
镜面	$\rho\,\delta(\omega_o - r)$	确定性反射方向
折射	透射 delta + Fresnel 反射	概率分支（B16）

易错点

• 公式里系数是 2，不是 1。常见笔误 $r = d - (d\cdot n)n$ 漏了 2，结果是”入射方向关于法线的投影”，不是反射。可用特例验证：$d$ 垂直入射（$d = -n$），应有 $r = -d = n$；代入 $-n - 2(-n\cdot n)n = -n + 2n = n$ ✓。
• 入射方向 $d$ 指向表面（从光源/相机出发），反射方向 $r$ 也指向”离开表面”的一侧。若读者习惯把入射方向定义成”从表面指向光源”，公式会差一个负号。

参考

• Reflection (mathematics) — 维基百科
• Specular reflection

Snell 定律与折射方向向量推导

Snell 定律给出折射角与入射角的关系：

$n_1 \sin\theta_1 = n_2 \sin\theta_2$

但渲染要的是完整的折射方向向量 $t$，不只是角度。下面从入射方向 $d$（单位向量，指向表面）和法线 $n$（指向入射侧介质）推出 $t$。

把 d 分解为切向 + 法向

设 $\eta = n_1 / n_2$。入射方向 $d = d_\parallel + d_\perp$：

• 切向分量 $d_\parallel = d + \cos\theta_1\, n$（注意 $d\cdot n = -\cos\theta_1$，因 $d$ 指向表面）
• 法向分量 $d_\perp = -\cos\theta_1\, n$

折射后切向方向不变但大小按 $\eta$ 缩放（因 $\sin\theta_2 = \eta \sin\theta_1$），法向重算为 $-\cos\theta_2\,n$：

$t = \eta\, d_\parallel - \cos\theta_2\, n$

用 cosθ₁ 表示 cosθ₂

由 $\sin^2\theta_2 = \eta^2 \sin^2\theta_1$ 和 $\sin^2 = 1 - \cos^2$：

$\cos^2\theta_2 = 1 - \eta^2(1 - \cos^2\theta_1)$

把 $\cos\theta_1 = -d\cdot n$ 代入，并展开 $d_\parallel = d + (d\cdot n)\,n \cdot(-1)\cdot(-1)$… 整理得经典向量折射公式：

$\boxed{\,t = \eta\,d + (\eta\cos\theta_1 - \cos\theta_2)\,n\,}$

smallpt 代码里的等价写法

Ray reflRay(x, r.d - n*2*n.dot(r.d));          // 反射方向（B10），Fresnel 时复用
bool into = n.dot(nl) > 0;                      // true = 从外向内
double nc=1, nt=1.5;                            // 空气 / 玻璃折射率
double nnt = into ? nc/nt : nt/nc;              // η = n₁/n₂
Vec tdir = (r.d*nnt - n*((into?1:-1)*(ddotn*nnt + sqrt(cos2t)))).norm();

代码把公式重排：r.d*nnt 对应 $\eta d$，括号里第二项是 $(\eta\cos\theta_1 - \cos\theta_2)\,n$，乘上 (into?1:-1) 是为了配合 smallpt 用几何法线 n（而非定向法线 $n_l$）计算——内外两种情形通过这一个符号分支统一。

cos2t 与全内反射

double ddotn = r.d.dot(nl);                     // = -cosθ₁（带符号）
double cos2t = 1 - nnt*nnt*(1 - ddotn*ddotn);   // cos²θ₂

验证：$\cos^2\theta_2 = 1 - \eta^2(1 - \cos^2\theta_1) = 1 - \eta^2(1 - \text{ddotn}^2)$ ✓。

若 cos2t < 0，即 $\sin^2\theta_2 > 1$，折射角无实数解——发生全内反射（见 B14）。

易错点

• $\eta$ 是 $n_1/n_2$（入射侧 / 折射侧），不是 $n_2/n_1$。光线从外进玻璃 $\eta = 1/1.5 < 1$，弯向法线（接近垂直）；从玻璃出空气 $\eta = 1.5 > 1$，弯离法线。搞反会让所有玻璃看起来”反向”。
• ddotn 带符号是 $-\cos\theta_1$，因为 $d$ 朝向表面而 $n_l$ 朝向入射侧。代码里 ddotn*ddotn 平方掉符号算 $\cos^2\theta_1$ 没问题，但重建 $t$ 时符号要小心（smallpt 用 (into?1:-1) 处理）。

参考

• Snell’s law (维基百科)
• Vector form of refraction (Bram de Greve)

全内反射：sin²θ₂ > 1 的几何含义

从光密介质（高折射率 $n_1$）射向光疏介质（低折射率 $n_2$）时，Snell 定律 $n_1\sin\theta_1 = n_2\sin\theta_2$ 给出：

$\sin\theta_2 = \frac{n_1}{n_2}\sin\theta_1$

由于 $n_1 > n_2$，比值 $n_1/n_2 > 1$。当入射角 $\theta_1$ 增大到使 $\sin\theta_2 = 1$，对应 $\theta_2 = 90°$（折射光沿界面方向）。再增大就要求 $\sin\theta_2 > 1$——数学上无实数解，物理上能量全部反射回光密介质。

临界角

令 $\sin\theta_2 = 1$（$\theta_2 = 90°$）求临界角 $\theta_c$：

$\sin\theta_c = \frac{n_2}{n_1}$

玻璃（$n=1.5$）→ 空气（$n=1$）：

$\theta_c = \arcsin\!\left(\frac{1}{1.5}\right) \approx 41.8°$

超过 $41.8°$ 入射的光全部反射——这就是潜望镜、双筒望远镜里全反射棱镜的原理，也是从水底抬头看、水面在某个角度以上变成镜子的原因（潜水员视野呈”Snell 窗”）。

在 smallpt 代码里

double cos2t = 1 - nnt*nnt*(1 - ddotn*ddotn);   // cos²θ₂
if (cos2t < 0) {                                 // sin²θ₂ = 1 - cos²θ₂ > 1
  // 全内反射：折射方向不存在
  return obj.e + f.mult(radiance(reflRay, depth));   // 走镜面反射分支
}

cos2t < 0 ⟺ $\cos^2\theta_2 < 0$ ⟺ $\sin^2\theta_2 > 1$。注意这只在 nnt > 1（从光密到光疏）时才可能发生：若 $\eta = n_1/n_2 < 1$，则 $\eta^2(1-\cos^2\theta_1) \le \eta^2 < 1$，cos2t 恒为正。所以全内反射是单向的——只有从玻璃里向外出射时才可能触发，从空气进玻璃永远不会。

物理后果：玻璃球的”内反射陷阱”

光线进入玻璃球后，在内壁反弹。如果忽略全内反射（很多简易渲染器图省事直接用 Schlick 近似），玻璃球内部会”漏光”，看起来通透得不真实。smallpt 在内表面也走完整的 Fresnel + 折射/反射分支，所以能渲染出真实玻璃那种内部反射带来的边缘变暗和焦点光斑。

易错点

• cos2t < 0 不是”折射失败”，而是”该界面情况下根本不存在透射光”。代码必须返回反射结果，而不是返回 0 或黑。漏掉这一步会让玻璃内部出现黑色空洞。
• 临界角对方向敏感。同一条射线在玻璃内多次弹跳，每次内壁命中都重新判 cos2t，可能在某次入射角超过临界角就不再透出——这是玻璃球的”陷光”效应，蒙特卡洛要靠 Russian Roulette 在深层弹射终止。

参考

• Total internal reflection (维基百科)
• Critical angle

Fresnel：界面反射率的物理根源

光打到不同折射率介质的界面时，能量按比例分为反射光和折射（透射）光。这个比例由Fresnel 方程精确描述，依赖入射角 $\theta_1$、折射角 $\theta_2$ 和两介质的折射率。

完整的 Fresnel 方程分 s 偏振和 p 偏振两套（电场垂直/平行于入射面）：

$r_s = \frac{n_1\cos\theta_1 - n_2\cos\theta_2}{n_1\cos\theta_1 + n_2\cos\theta_2}, \quad r_p = \frac{n_2\cos\theta_1 - n_1\cos\theta_2}{n_2\cos\theta_1 + n_1\cos\theta_2}$

非偏振光的反射率 $R = \tfrac{1}{2}(r_s^2 + r_p^2)$。两个方程、四个三角函数、还要算 $\theta_2$——渲染里每条射线都要算太贵。

Schlick 近似：一个 5 次多项式

Christophe Schlick (1994) 提出的近似，用入射角的余弦直接插值两个端点：

$R(\theta) \approx R_0 + (1 - R_0)(1 - \cos\theta)^5$

其中 $R_0$ 是正入射（$\theta=0$）时的反射率：

$R_0 = \left(\frac{n_1 - n_2}{n_1 + n_2}\right)^2$

对空气（$n_1=1$）→ 玻璃（$n_2=1.5$）：

$R_0 = \left(\frac{1 - 1.5}{1 + 1.5}\right)^2 = \left(\frac{-0.5}{2.5}\right)^2 = 0.04$

即正对着玻璃看只反射 4%，96% 透射。但 $\theta \to 90°$（掠射）时 $(1-\cos\theta)^5 \to 1$，$R \to 1$——这就是为什么斜着看玻璃窗或水面几乎全是反射。

smallpt 的实现

double a = nt - nc, b = nt + nc;
double R0 = a*a/(b*b);                          // 正入射反射率
double c = 1 - (into ? -ddotn : tdir.dot(n));   // 1 - cosθ₁（入射侧）
double Re = R0 + (1-R0)*c*c*c*c*c;              // Schlick: 5 次方
double Tr = 1 - Re;                             // 透射率

c*c*c*c*c 是 $(1-\cos\theta)^5$（编译器一般不会把 5 次幂优化成乘法链，手写更明确）。into 分支选 cosθ 的来源：从外进玻璃用入射角的余弦 -ddotn，从内出玻璃用折射方向的余弦 tdir.dot(n)——这是 Schlick 近似的入射侧约定。

为什么 5 次方

经验拟合：Fresnel 曲线在 $\theta \to 90°$ 附近迅速爬升，5 次多项式的”急转弯”形状与真实曲线吻合度最高（误差通常 < 1%）。3 次太平、7 次太陡。Schlick 论文用 5 次是基于大量材料实测数据的最小二乘拟合。

易错点

• Schlick 是近似，不是精确 Fresnel。对高折射率材料（如钻石 $n=2.4$）或带吸收的金属，误差会变大；产品级渲染器（Mitsuba、pbrt）对 conductor 用完整的 Fresnel-Conductor 方程。smallpt 只渲染电介质（玻璃）才用它。
• 全内反射时 Schlick 失效：$\theta_2$ 不存在，公式不能用折射角余弦。smallpt 在 B14 已经先判 cos2t < 0 直接走反射，不会走到 Schlick 这里。

参考

• Fresnel equations (维基百科)
• Schlick approximation (原论文, 1994)

分裂 vs 随机：路径数的指数爆炸与无偏终止

光打到玻璃界面，物理上同时有反射光（份额 $R_e$）和折射光（份额 $T_r = 1 - R_e$）。最忠实的做法是分裂成两条子射线分别递归：

$L = R_e \cdot L_{\text{refl}} + T_r \cdot L_{\text{refr}}$

问题：每次穿过界面路径数翻倍，10 个界面就是 $2^{10}=1024$ 条——指数爆炸，渲染不可行。

smallpt 的策略：分深度处理

double P=.25+.5*Re, RP=Re/P, TP=Tr/(1-P);   // RR 概率与补偿系数
return obj.e + f.mult(
  depth>2 ? (erand48(Xi)<P ?                              // 深层：RR 选一条
             radiance(reflRay,depth,Xi)*RP :
             radiance(Ray(x,tdir),depth,Xi)*TP)
          : radiance(reflRay,depth,Xi)*Re                  // 浅层：两条都算
            + radiance(Ray(x,tdir),depth,Xi)*Tr);

浅层（depth ≤ 2）：同时算反射和折射两条，按 $R_e, T_r$ 加权求和。这两条最初的弹射（相机直接看到的反射/折射）最影响画质，分裂保精度。

深层（depth > 2）：随机选一条，选中后乘以补偿系数 $RP$ 或 $TP$。注意 smallpt 选反射的概率不是 $R_e$，而是 $P = 0.25 + 0.5 R_e$——下面解释为什么。

为什么用 P = 0.25 + 0.5·Re 而不是 Re

直觉是”反射能量占 $R_e$，就按 $R_e$ 的概率选反射”（按贡献抽样）。但当 $R_e$ 很小（正入射玻璃只有 0.04）时，反射几乎永远不被选中，反射路径的样本极少，方差爆炸——尤其是玻璃球边缘那种强反射区域恰好对应少数样本。

smallpt 把选反射的概率下限抬高到 0.25：

$P = 0.25 + 0.5\,R_e$

• $R_e = 0$ 时仍以 $P=0.25$ 选反射（保证反射路径有 1/4 的样本）
• $R_e = 1$ 时 $P = 0.75$（全反射接近时几乎都选反射）

补偿系数 $RP = R_e / P$、$TP = T_r / (1-P)$ 保证期望不变：

$E[\hat{L}] = P \cdot \frac{R_e}{P} L_{\text{refl}} + (1-P)\cdot\frac{T_r}{1-P} L_{\text{refr}} = R_e L_{\text{refl}} + T_r L_{\text{refr}} = L \;\checkmark$

这种偏离贡献的抽样叫 splitting heuristic——牺牲一点理论最优性换取更均衡的样本分配，实测对玻璃场景降噪明显。

为什么乘 RP 而不是除以 P

很多人写成”选中除以概率”的形式 $L_i / P$。但 smallpt 这里返回的是带权重的贡献 $R_e \cdot L_{\text{refl}}$，所以补偿系数是”权重 ÷ 概率”：$RP = R_e / P$。两种写法等价：

• “返回 $L_{\text{refl}} / P$” ⟺ “返回 $R_e L_{\text{refl}} / R_e P$“——前者把 $R_e$ 视为路径的一部分，后者显式分离权重。
• smallpt 选后者，是因为浅层直接用 *Re，深层保持相同的”返回加权贡献”语义，便于读者对照。

易错点

• depth > 2 不是 >= 2。即 depth=0,1,2 都走分裂（前 3 次弹射），depth=3 起才进 RR。
• P 不是 $R_e$。读 RP = Re/P 时容易误以为 $P=R_e$ 导致 $RP=1$ 的错觉；实际上 $P$ 由启发式公式给出，$RP$ 通常介于 0.05–4 之间。

参考

• Russian roulette & splitting (PBRT 13.7)
• Sampling proportional to contribution (Veach thesis, Ch. 9)

折射分支完整代码逐行

把前面 B13/B14/B15/B16 的零件串起来，这是 smallpt 折射分支的全部：

Ray reflRay(x, r.d - n*2*n.dot(r.d));                  // 1. 反射方向（B10 公式）
bool into = n.dot(nl) > 0;                             // 2. 从外向内？
double nc=1, nt=1.5, nnt=into?nc/nt:nt/nc, ddn=r.d.dot(nl), cos2t;
if ((cos2t=1-nnt*nnt*(1-ddn*ddn))<0)                   // 3. cos²θ₂<0 → 全内反射
  return obj.e + f.mult(radiance(reflRay,depth,Xi));
Vec tdir = (r.d*nnt - n*((into?1:-1)*(ddn*nnt+sqrt(cos2t)))).norm();  // 4. 折射方向
double a=nt-nc, b=nt+nc, R0=a*a/(b*b), c=1-(into?-ddn:tdir.dot(n));
double Re=R0+(1-R0)*c*c*c*c*c, Tr=1-Re, P=.25+.5*Re,   // 5. Schlick + RR 概率
       RP=Re/P, TP=Tr/(1-P);
return obj.e + f.mult(                                  // 6. 概率分支返回
  depth>2 ? (erand48(Xi)<P ?
             radiance(reflRay,depth,Xi)*RP : radiance(Ray(x,tdir),depth,Xi)*TP)
          : radiance(reflRay,depth,Xi)*Re + radiance(Ray(x,tdir),depth,Xi)*Tr);

六步对应物理/数学

行	物理量	数学来源
1	反射方向 $r=d-2(d\cdot n)n$	B10 镜面公式
2	进入判断 $\eta=n_1/n_2$	B13 Snell 比值
3	$\cos^2\theta_2 = 1-\eta^2(1-\cos^2\theta_1)$	B13 推导，<0 即全反射
4	折射向量 $t = \eta d + (\eta\cos\theta_1 - \cos\theta_2)\,n$	B13 向量折射公式
5	$R_0=((n_1-n_2)/(n_1+n_2))^2$，$R_e=R_0+(1-R_0)(1-\cos\theta)^5$	B15 Schlick
6	$L = R_e L_{\text{refl}} + T_r L_{\text{refr}}$	B16 Fresnel 加权 + RR

关键技巧：复用 reflRay

第 1 行无条件先算 reflRay。哪怕最终走全反射（第 3 行）或选折射分支（第 6 行深层），反射方向至少在两个情形下要用：全反射时直接返回它；浅层分裂时也要乘 $R_e$。提前算避免重复。

(into?1:-1) 符号分支

第 4 行 n*((into?1:-1)*(ddn*nnt+sqrt(cos2t)))：into 决定法线方向系数。smallpt 用几何法线 n（不是定向法线 $n_l$），而 $t$ 必须指向折射侧介质。从外进玻璃，$t$ 应朝内（与 $n$ 反向）；从内出玻璃，$t$ 应朝外（与 $n$ 同向）。(into?1:-1) 正好补这个符号。

为什么这套代码”物理正确”

• 能量守恒：$R_e + T_r = 1$，反射和折射的能量加起来等于入射（无吸收的理想电介质）。
• 无偏估计：RR 补偿后期望等于真实加权贡献（B16）。
• 完整物理分支：全反射（B14）没有用近似糊弄，而是直接走反射；Fresnel 用 Schlick 但对电介质足够精确。

易错点

• ddn = r.d.dot(nl) 用的是定向法线 nl，不是 n。$ddn$ 在数值上是 $-\cos\theta_1$（$d$ 朝向表面、$n_l$ 朝向入射侧），后面 cos2t = 1 - nnt^2(1-ddn^2) 用 $ddn^2$ 故符号无关，但重建 $t$ 时用到的 ddn*nnt 必须是带符号的 $-\cos\theta_1$，所以 ddn 不能取绝对值。
• c = 1-(into?-ddn:tdir.dot(n)) 入射角余弦的来源：进入时用 $-ddn=\cos\theta_1$（入射侧），出去时用 tdir.dot(n)（折射侧角）。两边都对应”光离开介质的夹角”，是 Schlick 的入射侧约定。

参考

• smallpt 源码 (Kevin Beason)
• Reflections and Refractions in Ray Tracing (Bram de Greve, 2006)

三种材质在 BRDF 层面的统一比较

smallpt 的 DIFF / SPEC / REFR 三类材质，对应三种典型 BRDF，数学上可以放进同一框架：

材质	BRDF $f_r$	物理模型	采样	是否随机
DIFF 漫反射	$\rho/\pi$（常数）	Lambertian 散射	余弦加权 $\cos\theta/\pi$	是
SPEC 镜面	$\rho\,\delta(\omega_o - r)$	反射定律	反射方向 $r$（pdf=1）	否
REFR 折射	$\rho\bigl(R_e\delta_r + T_r\delta_t\bigr)$	Snell + Fresnel	概率分支（B16）	是

为什么 DIFF 需要 π、SPEC 和 REFR 不需要

DIFF 的 BRDF 是 $\rho/\pi$，这个 $\pi$ 是为了能量守恒：完美漫反射体要把入射能量均匀散到整个半球（立体角 $2\pi$，但带 $\cos\theta$ 权重积分是 $\pi$），所以 BRDF 必须有 $1/\pi$ 归一化。结合余弦加权采样，$\pi$ 抵消（见 B05）。

SPEC 和 REFR 的 BRDF 是 Dirac delta 函数——能量集中在单一方向（反射方向 $r$ 或折射方向 $t$）。delta 函数的”积分值”是该方向系数（$R_e$ 或 $T_r$），不需要 $\pi$ 归一化，蒙特卡洛直接采样该方向，pdf=1，估计器就是 $\rho \cdot L$。

为什么 DIFF 是”软”而 SPEC 是”硬”

• DIFF 把入射光散到所有方向，每个方向只分到很小一份——边缘平滑过渡、形成颜色溢出（color bleeding，红墙把白色地板染红）。
• SPEC 把所有能量集中在一个方向，完全不模糊——形成完美镜像。这就是为什么镜子里的像清晰、磨砂墙看不到像。
• REFR 介于两者：玻璃既透射又反射，且折射方向由 Snell 定律唯一决定（确定性），但哪条路径被采样是随机的（Fresnel 加权）。多次内部弹射后玻璃边缘出现复杂的反射纹路。

物理本质：散射 vs 镜面

这是材质建模的核心二分：

• 散射（scattering）：表面微观粗糙度远大于光波长，入射光被随机反弹到很多方向。Lambertian 是其最简模型。
• 镜面（specular）：表面光滑到原子尺度（如抛光金属、液体表面、玻璃），入射光按几何光学精确反弹。

真实材质（如塑料、皮肤）是两者的混合：一层镜面高光（clear coat）+ 底层漫反射。smallpt 为了教学清晰，把三种极端分开；产品级渲染器（如 Disney Principled BRDF）用一个参数 roughness 连续插值两者。

易错点

• REFR 的”反射”和 SPEC 的”反射”不一样。REFR 内部的反射分支（Fresnel 反射）受角度调制（掠射时 $R_e\to 1$），SPEC 是恒定全反射。前者是”玻璃表面偶尔反射一点”，后者是”镜子全反射”。
• 不要把 BRDF 当材质。BRDF 描述表面如何散射光，是 4D 函数 $f_r(x,\omega_i,\omega_o)$；材质（material）是 BRDF + 反照率 + 其他参数（折射率、粗糙度）的捆绑包。smallpt 的 obj.refl 枚举 + obj.c 颜色就是一个最小材质系统。

参考

• Bidirectional reflectance distribution function (维基百科)
• Disney Principled BRDF (Burley 2012)

一条路径 = 渲染方程积分的一个样本

每次 radiance 递归调用，物理上就是光弹射一次（一次 bounce）。从相机出发的一条路径：

相机 → 地板(depth 0) → 红墙(depth 1) → 天花板(depth 2) → 光源(depth 3)

最终返回的亮度 = 光源自发光 × 路径上每个表面反照率的连乘积：

$L = L_e \cdot f_{\text{地板}} \cdot f_{\text{红墙}} \cdot f_{\text{天花板}}$

为什么是连乘积

每次 radiance 返回 obj.e + f.mult(radiance(下一条射线, depth))。展开递归：

$L^{(0)} = f_0 \cdot L^{(1)} = f_0 \cdot f_1 \cdot L^{(2)} = \cdots = \prod_{k=0}^{n-1} f_k \cdot L_e^{(n)}$

(忽略第一次命中若是光源则直接 $L^{(0)} = L_e$，以及 RR 补偿因子 $1/p$)。连乘积正是几何衰减——每个表面都按反照率 $\rho$ 吸收一部分光，路径越长亮度越暗，这也是间接照明比直接照明弱的原因。

多条路径平均 = 蒙特卡洛积分

渲染方程的反射项是积分 $\int_\Omega f_r L_i \cos\theta\, d\omega$。一条随机路径对应一个采样方向 $\omega_i$，单条路径的估计 $\hat{L} = \prod f_k \cdot L_e$ 是这个积分的一个样本。每个像素发射 $N$ 条路径取平均：

$\bar{L} = \frac{1}{N}\sum_{i=1}^N \hat{L}_i \;\xrightarrow{N\to\infty}\; L$

由大数定律收敛到真实值，方差随 $1/N$ 下降（注意是 $\sqrt{N}$ 的标准差，所以 4 倍样本只降噪一半）。

路径追踪 vs 光线追踪的本质区别

• 光线追踪（Whitted-style）：从相机发射少量”重要”射线（镜面反射/折射方向），忽略散射表面（漫反射直接用常数着色）。快但物理不正确——间接照明缺失。
• 路径追踪（path tracing，smallpt）：从相机发射很多随机路径，覆盖所有散射方向，蒙特卡洛估计完整渲染方程。慢但物理正确。

smallpt 一个像素可能跑几百到几千条路径（如 1000 spp = samples per pixel），最终求平均消除噪声。这就是为什么小图渲染也要几十秒到几分钟。

间接照明的物理意义

红色墙把白光照地板，地板的红色”染色”就是一次间接照明（一次 bounce）。smallpt 自然产生这个效果，不需要任何特殊代码——因为递归路径会经过红墙（depth 1），把红墙的反照率 $f_{\text{红}}=(0.75, 0.25, 0.25)$ 乘进连乘积，地板接收到的光就被偏红了。

类似地，玻璃球下方出现的焦散（caustics）——光经球面聚焦形成的亮斑——也是路径追踪自动产生的（光线路径经过玻璃聚焦到桌面）。Whitted 光线追踪做不到，因为焦散需要追踪从光源出发的”折射方向”，不是镜面方向。

易错点

• depth 0 不是”光源”。depth 0 是相机射线第一次命中表面，可能是地板（吸收）、墙（反射）、或光源（直接返回自发光）。depth 与”是否光源”无关，只数弹射次数。
• 连乘积里的 f 已含 RR 补偿。B04 的 f = f * (1/p) 在某一层把 f 放大，所以连乘积展开后某层是 $f_k/p$，这是补偿该层提前终止的”幻影”路径——期望上等价于完整无限递归的结果。

参考

• Path tracing (维基百科)
• Monte Carlo integration
• Caustic (optics)