6.1 渲染方程
前五节一直在回答「一个方向上的光是哪里来的」：求交给交点，着色给颜色，递归把镜面反射串成一条链。这些答案都是局部的——它们假设别处的光已经知道了。本节把这个假设拆掉。实时混合光追那篇用屏幕空间的近似回避了这个问题（见 实时光线追踪的混合方案 ）；本卷要从渲染方程完整推到一条可跑的路径追踪，第一步就是把「一个点的出射光」写成「它自己发的光」加上「半球上所有方向射来的光」的积分。
出射辐亮度等于自发光加反射
记 $x$ 是表面点，$\omega$ 是出射方向（指向相机），$\omega’$ 是入射方向。渲染方程写成：
$$ L_o(x, \omega) = L_e(x, \omega) + \int_{\Omega} f_r(x, \omega, \omega’), L_i(x, \omega’), \cos\theta’, \mathrm{d}\omega’ $$
五项各有明确分工：
| 记号 | 含义 | 本卷里的取值 |
|---|---|---|
| $L_o(x,\omega)$ | 从 $x$ 沿 $\omega$ 出去的辐亮度 | 待求量，最终进像素 |
| $L_e(x,\omega)$ | $x$ 自己发出来的辐亮度 | 自发光球为 4.0，其余为 0 |
| $f_r(x,\omega,\omega’)$ | 双向反射分布函数（BRDF） | Lambertian：$\rho/\pi$，$\rho$ 是 albedo |
| $L_i(x,\omega’)$ | 从 $\omega’$ 方向射到 $x$ 的辐亮度 | 本身又是一个 $L_o$ |
| $\cos\theta'$ | 入射方向与法线的夹角余弦 | 投影因子，防止斜射的能量被高估 |
积分域 $\Omega$ 是 $x$ 处法线正侧的整个半球，立体角测度 $\mathrm{d}\omega’ = \sin\theta’,\mathrm{d}\theta’,\mathrm{d}\phi’$。半球而不是整个球面，因为背面的入射光到不了 $x$ 的正面。
这五个量的关键性质是：$L_i(x,\omega’)$ 就是另一点沿相反方向的 $L_o$。也就是说，方程左右两边出现的是同一个函数，只不过自变量换了个位置。
为什么它是第二类 Fredholm 积分方程
把 $L_i(x,\omega’)$ 写成 $L_o(x’, -\omega’)$（$x’$ 是从 $x$ 沿 $\omega’$ 出发撞到的下一个点），方程变成：
$$ L_o(x, \omega) = L_e(x, \omega) + \int_{\Omega} K(x, \omega; x’, \omega’), L_o(x’, -\omega’), \mathrm{d}\omega’ $$
这里 $K = f_r\cos\theta’$ 是积分核。未知函数 $L_o$ 同时出现在积分号内外——这正是第二类 Fredholm 积分方程的定义。第一类方程（未知量只出现在积分号内）通常是不适定的，第二类则只要核的范数小于 1 就唯一可解，而这里的范数就是表面的反照率 $\rho$，物理上恒有 $\rho < 1$（否则能量不守恒）。
两类方程的这个区别直接决定了算法形态：第一类要靠正则化去逼近，第二类可以直接迭代——迭代的每一步就是一次弹射。渲染器里那些 depth 参数、递归树、俄罗斯轮盘赌，全部是这一条性质的直接后果。
积分域的测度与那个余弦
积分写在半球上，测度是立体角：
$$ \int_{\Omega} \cdots \mathrm{d}\omega’ = \int_{0}^{2\pi}\int_{0}^{\pi/2} \cdots \sin\theta’,\mathrm{d}\theta’,\mathrm{d}\phi’ $$
$\sin\theta’$ 是球面坐标的雅可比因子，$\theta’$ 的上限取 $\pi/2$ 而不是 $\pi$——半球。这个上限不是偷懒：法线背侧的 $\omega’$ 在 $x$ 处 $\cos\theta’ < 0$，$L_i$ 也取不到值，两块都被截掉了。
$\cos\theta’$ 是投影立体角的因子。它只出现一次，但常常被误算两次——一次来自「立体角 → 投影立体角」的换算，一次来自 BRDF 里 Lambertian 的 $\rho/\pi$（$\pi$ 是归一化常数，不是余弦）。本卷的积分器统一按投影立体角写：
$$ \mathrm{d}\omega^\perp = \cos\theta’,\mathrm{d}\omega’ $$
写成投影立体角之后，Lambertian 的估计量有一个漂亮的化简：如果方向按 $p(\omega’) = \cos\theta’/\pi$ 采样，那么
$$ \frac{f_r \cos\theta’}{p} = \frac{(\rho/\pi)\cos\theta’}{\cos\theta’/\pi} = \rho $$
权重退化成 albedo 本身。这就是 6.2 里「一条路径的贡献 = 沿途 albedo 的连乘」的来历，也是下一节采样器的设计依据。
算子形式与收敛速度
把 $\mathcal{T}L(x,\omega) = \int_\Omega K(x,\omega;x’,\omega’)L(x’,-\omega’),\mathrm{d}\omega’$ 记成算子，方程是 $(I - \mathcal{T})L = L_e$，形式解 $L = (I-\mathcal{T})^{-1}L_e$。Neumann 级数就是这个逆算子的展开。
收敛性由 $\mathcal{T}$ 的谱半径决定。对 Lambertian 表面，$|\mathcal{T}| = \max\rho$：反照率越接近 1，级数收敛越慢，需要的弹射阶数越多。这解释了两件常被当成 bug 的事：白色房间里路径追踪特别慢（$\rho \approx 0.9$，第 10 阶还有 $0.9^{10} \approx 0.35$ 的贡献），以及为什么次表面散射材质（$\rho$ 极高）必须换用别的算法。
Neumann 级数：每次弹射就是一项
把方程写成算子形式 $L = L_e + \mathcal{T}L$，其中 $\mathcal{T}$ 是积分算子。形式解是几何级数（Neumann 级数）：
$$ L = \sum_{k=0}^{\infty} \mathcal{T}^k L_e = L_e + \mathcal{T}L_e + \mathcal{T}^2 L_e + \cdots $$
每一项都有直白的几何含义：$\mathcal{T}^k L_e$ 是「光从光源出发、恰好经过 $k$ 次弹射后到达相机」的那部分。$k=0$ 是光源直射相机（本卷场景里没有），$k=1$ 是直接光照，$k=2$ 是一次间接弹射，依此类推。级数收敛的速度由核的范数决定：$|\mathcal{T}| = \rho$，所以弹射阶数越高、贡献越小。
要看到这件事的可计算形态，需要一个能解析求解的构型。双板腔就是：两块无限大平行板，反照率 $a$、$b$，板间均匀发射 $E$。均匀腔里的几何项是常数，积分方程退化成代数方程 $L = aE + ab,L$，解是
$$ L = \frac{aE}{1 - ab} = aE\sum_{k=0}^{\infty}(ab)^k $$
第 $k$ 项 $T_k = aE(ab)^k$ 就是恰好 $k+1$ 次弹射的那部分能量，而截断到 $K$ 项的相对偏差恰为 $(ab)^K$。这条等式是本节的判据：它把「少算了多少能量」写成了一个可以直接算出来的数。
#!/usr/bin/env python3
# file: lantern/sim/cg_trace_06_1_neumann.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_trace_06_1_neumann.py
"""6.1 最小复现：双板腔的 Neumann 级数——渲染方程在可解析构型上的展开。
两块无限大平行板，反照率 a、b，板间均匀发射 E。几何项在均匀腔里是常数，
积分方程退化成代数方程 L = aE + ab·L，于是 L = aE/(1 − ab) = Σ_k aE(ab)^k。
第 k 项就是「恰好 k+1 次弹射」的那部分能量，截断到 K 项的相对误差恰为 (ab)^K。
全程 Q16.16 整数（`fmul` 右移 16 位、`fdiv` 向零截断），与 TypeScript 侧逐位一致。
"""
import os
import sys
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cglib # noqa: E402
ONE = cglib.ONE
C8 = lambda v: (v << 16) // 255 # noqa: E731 8 位通道 -> Q16.16
def neumann(a, b, e, kmax):
"""返回 `(ab, 真值, [(k, 第 k 项, 前 k+1 项之和, 相对误差)])`。"""
ab = cglib.fmul(a, b)
exact = cglib.fdiv(cglib.fmul(a, e), ONE - ab)
rows = []
acc = 0
term = cglib.fmul(a, e)
for k in range(kmax + 1):
acc += term
rows.append((k, term, acc, float(exact - acc) / float(exact)))
term = cglib.fmul(term, ab)
return ab, exact, rows
def main():
a = b = C8(170)
e = ONE
ab, exact, rows = neumann(a, b, e, 16)
abf = ab / ONE
print("== 6.1 渲染方程：Neumann 级数 ==")
print(f"a = b = 170/255 -> Q16 {a}；ab = {ab}（{abf:.6f}）；E = 1.0")
print(f"真值 L = aE/(1−ab) = {exact} / {ONE} = {exact / ONE:.6f}")
print(f"{'k':>3} {'第 k 项':>10} {'前 k+1 项和':>12} {'相对误差':>14} {'(ab)^(k+1)':>14}")
for k, term, acc, err in rows:
print(f"{k:>3} {term:>10} {acc:>12} {err:>14.10f} {abf ** (k + 1):>14.10f}")
print("截断到 K 项：")
print(f"{'K':>3} {'估计值':>12} {'相对偏差':>14} {'(ab)^K':>14}")
for k in (1, 2, 4, 8, 16):
acc = rows[k - 1][2]
print(f"{k:>3} {acc / ONE:>12.8f} {(exact - acc) / exact:>14.10f} {abf ** k:>14.10f}")
frozen = next(k for k, t, _, _ in rows if t == 0)
print(f"第 {frozen} 项起被定点截断为 0，部分和冻结在 {rows[frozen][2]}，"
f"偏差停在 {rows[frozen][3]:.10f} 不再下降。")
if __name__ == "__main__":
main()
实跑输出（$a = b = 170/255$，$E = 1.0$）：
== 6.1 渲染方程：Neumann 级数 ==
a = b = 170/255 -> Q16 43690；ab = 29126（0.444427）；E = 1.0
真值 L = aE/(1−ab) = 78639 / 65536 = 1.199936
k 第 k 项 前 k+1 项和 相对误差 (ab)^(k+1)
0 43690 43690 0.4444232505 0.4444274902
1 19417 63107 0.1975101413 0.1975157941
2 8629 71736 0.0877808721 0.0877814486
3 3834 75570 0.0390264373 0.0390124889
4 1703 77273 0.0173705159 0.0173382225
5 756 78029 0.0077569654 0.0077055827
6 335 78364 0.0034969926 0.0034245728
7 148 78512 0.0016149748 0.0015219743
8 65 78577 0.0007884129 0.0006764072
截断到 K 项：
K 估计值 相对偏差 (ab)^K
1 0.66665649 0.4444232505 0.4444274902
2 0.96293640 0.1975101413 0.1975157941
4 1.15310669 0.0390264373 0.0390124889
8 1.19799805 0.0016149748 0.0015219743
16 1.19970703 0.0001907451 0.0000023164
第 13 项起被定点截断为 0，部分和冻结在 78624，偏差停在 0.0001907451 不再下降。
三条读数值得逐条看：
- 前三项吃掉九成能量。 第 0 项 0.666656 加第 1 项 0.296280 已经到 0.962936，前 4 项到 1.153107，也就是真值的 96.1%。
- 偏差确实按 $(ab)^K$ 走，但只走到第 4 项。 $K=1$ 时实测 0.4444232505 对理论 0.4444274902，$K=2$ 时 0.1975101413 对 0.1975157941，$K=4$ 时 0.0390264373 对 0.0390124889——相对差都在 0.1% 以内。到 $K=8$ 就分开了：实测 0.0016149748，理论 0.0015219743，差了 6.1%。
- 分开的原因不是数学，是定点算术。
fmul每步右移 16 位、丢掉最低位，8 次递推累计的损失已经和 $T_8$ 本身（65 个单位）同量级。第 13 项起 $T_k$ 被截断成 0，部分和冻结在 78624，偏差停在 0.0001907451——这是一条地板，不是渐近线。真值 78639 与冻结值 78624 相差 15 个 Q16 单位，即 $1 - ab$ 的定点表示偏大带来的舍入。
这条地板是 Q16.16 的固有属性，与「弹射多少阶才够」无关。想要更深的级数就得换更高精度，而不是换算法。
灯箱：每一阶的实际占比
双板腔的几何项是常数，真实场景不是。把同一套分阶记账搬到 64×48 的灯箱上（左墙红、右墙绿、自发光球、漫反射球，见 6.2 的场景配置），统计每个像素的路径在第几次弹射命中自发光体，得到每阶的能量占比：
| 弹射阶数 $k$ | 占全部能量比例 | 弹射阶数 $k$ | 占全部能量比例 |
|---|---|---|---|
| 1 | 0.592314 | 5 | 0.027179 |
| 2 | 0.183228 | 6 | 0.017312 |
| 3 | 0.104208 | 7 | 0.009115 |
| 4 | 0.053430 | 8 | 0.005397 |
（32×24，16 spp，maxDepth=16，RR 关。）第 1 阶 0.592314、第 8 阶之后合计 0.00781676——间接光占了 40.8%，这正是 6.2 里 maxDepth=1 只有参考解 77.7% 亮度的原因。相邻两阶的比值 $k=2 \to k=3$ 是 0.568736，与双板腔的 $ab = 0.444427$ 同量级：盒子里有效反照率在 0.44～0.57 之间，取决于那一阶主要打在哪个面上。
P 层：级数的逐位复现
TypeScript 侧复用同一套 fmul / fdiv（core/src/lantern/fixed.ts），把双板腔的 Neumann 级数重算一遍。完整模块 lantern/core/src/lantern/cg_trace_06_1_rendering_equation.ts 提供 twoPlate(a, b, e, kmax)，返回真值与逐阶部分和；下面这份最小测试只挑本节正文引用过的数字：
// file: lantern/core/test/cg_trace_06_1_neumann.test.ts
// tier: P
// verify: node --test 'core/test/cg_trace_06_1_neumann.test.ts'
// 6.1 的最小 P 复现：Neumann 级数的四项硬断言——真值、逐阶项、截断误差、冻结点。
// 复用的是完整模块 cg_trace_06_1_rendering_equation.ts，本文件只挑本节正文引用的数字。
import assert from "node:assert/strict";
import { test } from "node:test";
import { ONE } from "../src/lantern/fixed.ts";
import { c8, truncationBias, twoPlate } from "../src/lantern/cg_trace_06_1_rendering_equation.ts";
const sol = twoPlate(c8(170), c8(170), ONE, 24);
test("真值 L = aE/(1−ab) = 78639/65536 = 1.199936", () => {
assert.equal(sol.ab, 29126n);
assert.equal(sol.exact, 78639n);
assert.ok(Math.abs(Number(sol.exact) / 65536 - 1.199936) < 1e-6);
// 定点真值比实数解略小：1 − ab 的定点表示偏大
assert.ok(Number(sol.exact) / 65536 < 1.199936);
});
test("第 0 项是直接光照 aE = 43690，其后每项乘 ab", () => {
assert.equal(sol.rows[0].term, 43690n);
assert.deepEqual(sol.rows.slice(0, 5).map((r) => Number(r.term)), [43690, 19417, 8629, 3834, 1703]);
assert.deepEqual(sol.rows.slice(0, 5).map((r) => Number(r.partial)), [43690, 63107, 71736, 75570, 77273]);
});
test("截断偏差前四项贴近 (ab)^K，第 8 项起被定点截断拖偏", () => {
const abf = Number(sol.ab) / 65536;
assert.ok(Math.abs(truncationBias(sol, 1).relBias - abf) < 1e-5);
assert.ok(Math.abs(truncationBias(sol, 4).relBias - abf ** 4) < 4e-4);
assert.equal(truncationBias(sol, 1).relBias.toFixed(10), "0.4444232505");
assert.equal(truncationBias(sol, 4).relBias.toFixed(10), "0.0390264373");
assert.equal(truncationBias(sol, 8).relBias.toFixed(10), "0.0016149748");
assert.equal(truncationBias(sol, 16).relBias.toFixed(10), "0.0001907451");
});
test("第 13 项起被截断为 0，部分和冻结在 78624", () => {
assert.equal(sol.rows[12].term, 2n);
assert.equal(sol.rows[13].term, 0n);
assert.equal(sol.rows[13].partial, 78624n);
assert.equal(sol.rows[24].partial, 78624n);
assert.equal(truncationBias(sol, 16).relBias, truncationBias(sol, 24).relBias);
});
$ node --test 'core/test/cg_trace_06_1_neumann.test.ts'
✔ 真值 L = aE/(1−ab) = 78639/65536 = 1.199936
✔ 第 0 项是直接光照 aE = 43690，其后每项乘 ab
✔ 截断偏差前四项贴近 (ab)^K，第 8 项起被定点截断拖偏
✔ 第 13 项起被截断为 0，部分和冻结在 78624
四条断言里最硬的一条是真值逐位相同：sol.exact === 78639n。这不是「接近」，是同一个整数。第二条把第 0～4 项的项值与部分和整段比对（[43690, 19417, 8629, 3834, 1703] / [43690, 63107, 71736, 75570, 77273]），任何一处递推次序写反都会立刻暴露。第三条用 toFixed(10) 把 $K = 1, 4, 8, 16$ 的偏差钉成字符串，避免浮点比较的模糊空间。
从方程到算法的三步
把 $L = \sum_k \mathcal{T}^k L_e$ 逐项落地，一项一项地估：
- 采样方向：按某个 $p(\omega’)$ 在半球上取一个方向。取 $p = \cos\theta’/\pi$ 时权重退化成 $\rho$，取均匀 $p = 1/(2\pi)$ 时要显式算 $2\cos\theta’$——两种写法都能用，方差差很多。
- 求交：沿方向找最近的正交点。这是第 1～3 节的全部内容，也是耗时的大头。
- 累加权重：命中自发光体就把 $\text{throughput} \times L_e$ 加进结果；命中普通表面就把 throughput 乘上 $\rho$ 并回到第 1 步。
三项都不需要知道「整个场景」——这正是第二类方程的迭代解法比第一类的正则化便宜的根源。第 1 步和第 3 步构成下一节的路径采样器，第 2 步是已经写完的求交底座。
本节的位置
渲染方程已经把「一个像素的颜色」写成了一个第二类积分方程，Neumann 级数也给出了「每次弹射一项」的展开——但级数是无穷的，计算机只能算有限项。下一节要回答的是：怎么把这个积分变成可估计的形式、方向从哪里采样、截断在哪儿截，以及截断之后那 40.8% 的间接光到底丢了多少。
阅读导航：上一节：5.3 分层与低差异序列 · 下一节：6.2 路径采样 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。