8.1 间接光照
站内 /graphics-real-time-gi/ 那篇讲的是实时全局光照的工程做法：屏幕空间反射、光照探针、辐照度体积，目标是在一帧十几毫秒里把「看起来对」的间接光凑出来。本节反过来，不加任何近似，把渲染方程的积分项老老实实展开成按弹射次数排列的级数，然后问一个很朴素的问题：砍到第 K 项，画面差多少。
把渲染方程按弹射次数展开
渲染方程 L_o = L_e + ∫ f L_i cosθ dω 里，L_e 是自发光，积分项是反射。把积分项里的 L_i 用同一个方程再展开一次，就得到一条无穷级数：
- 第 0 项：光源自己发出的能量，也就是直接光；
- 第 1 项：被光源直接照亮的表面反射出来的能量；
- 第 k 项：被 k−1 次弹射照亮的表面反射出来的能量。
有限弹射（finite-bounce）渲染就是把这条级数砍到第 K 项。截断误差完全由场景决定——反照率越高、腔体越封闭，衰减越慢。能回答「K 取多少够」的只有实测，而且得先有一个不截断的参考解可比。
最小复现：封闭腔里的空间剖面
参考解本身要先可信，所以先在一个能手算的几何上把方法立起来：二维矩形腔（z 方向无限），四壁各切 8 段，段与段之间的形状因子用 Hottel 交叉线法精确求出——凸腔里端点连线不被遮挡，公式没有近似。天花板正中偏右一段自发光，其余反照率 217/255。
一个必须钉住的细节：Hottel 只在所有线段沿腔体同向绕行时给出正号。段的方向反了，交叉线与非交叉线互换，形状因子整体变负，被 v > 0 的钳位吞掉之后表现成「这一段收不到光」，而且不报任何错。第一版就是天花板与左墙的方向反了，地板上有一半的段完全黑掉——比错更危险的是它安静地错。
#!/usr/bin/env python3
# file: lantern/sim/cg_trace_08_1_ambient_profile.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_trace_08_1_ambient_profile.py
"""8.1 最小复现：封闭矩形腔里「面积光源」与「环境常数」两种源项的空间剖面。
二维矩形腔（z 方向无限），四壁各切 8 段，段间形状因子用 Hottel 交叉线法精确求出
（凸腔里端点连线不被遮挡，公式没有近似），于是整腔的辐射度方程可以精确解出、不需要
任何采样。天花板正中偏右一段自发光，其余反照率 217/255。
这个可解析的模型上量两件事：**弹射次数截断**留下的残余，以及**把间接项换成一个与
位置无关的常数**之后空间剖面被抹平了多少。全程 Q16.16 整数。
一个必须钉住的细节：Hottel 的交叉线法只在**所有线段沿腔体同向绕行**时才给出正号。
段的方向反了，交叉线与非交叉线互换，形状因子整体变负——被 `v > 0` 的钳位吞掉之后
表现成「这一段收不到光」，而且不报任何错。
"""
import os
import sys
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cglib # noqa: E402
ONE = cglib.ONE
W = H = ONE # 腔是 1 x 1
M = 8 # 每面墙切 8 段
SEGLEN = cglib.fdiv(ONE, cglib.f_from_int(M)) # 每段长 1/8
RHO = cglib.fdiv(cglib.f_from_int(217), cglib.f_from_int(255)) # 217/255 ≈ 0.8510
F32 = cglib.f_from_int
def segs():
"""四壁线段，**逆时针**绕行：地板左->右、右墙下->上、天花板右->左、左墙上->下。"""
out = []
for t in range(M): # 地板
out.append(((cglib.fmul(W, cglib.fmul(SEGLEN, F32(t))), 0),
(cglib.fmul(W, cglib.fmul(SEGLEN, F32(t + 1))), 0)))
for t in range(M): # 右墙
out.append(((W, cglib.fmul(H, cglib.fmul(SEGLEN, F32(t)))),
(W, cglib.fmul(H, cglib.fmul(SEGLEN, F32(t + 1))))))
for t in range(M): # 天花板，从右往左
out.append(((cglib.fmul(W, ONE - cglib.fmul(SEGLEN, F32(t))), H),
(cglib.fmul(W, ONE - cglib.fmul(SEGLEN, F32(t + 1))), H)))
for t in range(M): # 左墙，从上往下
out.append(((0, cglib.fmul(H, ONE - cglib.fmul(SEGLEN, F32(t)))),
(0, cglib.fmul(H, ONE - cglib.fmul(SEGLEN, F32(t + 1))))))
return out, 2 * M + (M // 2 - 1) # 天花板第 3 段，x ∈ [0.5, 0.625]
def dist(p, q):
dx = q[0] - p[0]
dy = q[1] - p[1]
return cglib.fsqrt(cglib.fmul(dx, dx) + cglib.fmul(dy, dy))
def hottel(a, b):
"""同向绕行下的 Hottel：F(a->b) = (同行线长和 - 交叉线长和) / (2 * L_a)。"""
along = dist(a[0], b[0]) + dist(a[1], b[1])
cross = dist(a[0], b[1]) + dist(a[1], b[0])
v = cglib.fdiv(along - cross, 2 * SEGLEN)
return v if v > 0 else 0
def matrix(ss):
n = len(ss)
return [[0 if i == j else hottel(ss[i], ss[j]) for j in range(n)] for i in range(n)]
def incident(F, B, i):
"""第 i 段收到的入射辐照度 Σ_j F_ij B_j（Q16）。"""
return sum(cglib.fmul(F[i][j], B[j]) for j in range(len(F)))
def solve(F, E, iters):
n = len(F)
B = [0] * n
for _ in range(iters):
delta = 0
for i in range(n):
nb = E[i] + cglib.fmul(RHO, incident(F, B, i))
delta = max(delta, abs(nb - B[i]))
B[i] = nb
return B, delta
def neumann(F, E, K):
"""前 K+1 次弹射的部分和：T_0 = E，T_k = rho * (F T_{k-1})。"""
n = len(F)
T = list(E)
acc = list(T)
for _ in range(K):
T = [cglib.fmul(RHO, sum(cglib.fmul(F[i][j], T[j]) for j in range(n))) for i in range(n)]
acc = [acc[i] + T[i] for i in range(n)]
return acc
def main():
ss, em = segs()
n = len(ss)
F = matrix(ss)
E = [0] * n
E[em] = ONE
B, delta = solve(F, E, 400)
floors = list(range(M))
rowsum = sum(F[em])
print("== 8.1 封闭腔：面积光源 vs 环境常数 ==")
print(f"段数 {n}，每段长 {SEGLEN}（Q16，即 1/8 = {SEGLEN / ONE:.4f}），"
f"自发光段 #{em}，反照率 {RHO} = {RHO / ONE:.4f}")
print(f"Gauss-Seidel 400 步后最大残差 {delta} ulp；互易性抽查 "
f"A_i F_ij - A_j F_ji = {SEGLEN * (F[0][M] - F[M][0])}；"
f"自发光段行和 {rowsum} = {rowsum / ONE:.4f}")
print()
print("弹射截断：前 K+1 次弹射的部分和对全解的相对偏差（‰，取地板 8 段最大值）")
print(f"{'K':>3} {'最大偏差‰':>12} {'第 K 项总量':>14}")
full = max(B[i] for i in floors)
for K in (0, 1, 2, 3, 4, 8, 12, 16, 20, 24, 32, 64):
part = neumann(F, E, K)
worst = max(abs(B[i] - part[i]) for i in floors)
tk = part if K == 0 else [
part[i] - neumann(F, E, K - 1)[i] for i in range(n)
]
print(f"{K:>3} {1000.0 * worst / full:>12.4f} {sum(tk) / ONE:>14.6f}")
fz = next(K for K in range(1, 200)
if sum(neumann(F, E, K)) == sum(neumann(F, E, K - 1)) != 0)
fzr = 1000.0 * max(abs(B[i] - neumann(F, E, fz)[i]) for i in floors) / full
print(f" 级数在第 {fz} 项起被定点截断冻死（T_k 恒为 0），残余停在 {fzr:.4f}‰——"
"这不是物理弹射次数的极限，是 Q16.16 的极限。")
direct = neumann(F, E, 1)
share = [1.0 - direct[i] / B[i] for i in floors]
print(f" 地板 8 段里「2 次以上弹射」的占比："
f"{min(share):.4f} ~ {max(share):.4f}（单次弹射直接照到的只有 "
f"{1 - max(share):.4f} ~ {1 - min(share):.4f}）")
print()
lo = min(B[i] for i in floors)
hi = max(B[i] for i in floors)
print("地板剖面（8 段，x 从小到大；除以最亮段）")
print(" 面积光源：", " ".join(f"{1000.0 * B[i] / hi:7.2f}" for i in floors))
print(f" 最亮/最暗 = {hi / lo:.4f}")
mean_inc = cglib.fdiv(sum(incident(F, B, i) for i in floors), F32(len(floors)))
print(f" 环境常数按地板均值拟合：A = {mean_inc} = {mean_inc / ONE:.6f}")
amb = [E[i] + cglib.fmul(RHO, mean_inc) for i in range(n)]
print(" 环境常数：", " ".join(f"{1000.0 * amb[i] / hi:7.2f}" for i in floors))
errs = [abs(amb[i] - B[i]) / B[i] for i in range(n) if B[i] > 0]
print(f" 逐段相对偏差：均值 {100.0 * sum(errs) / len(errs):.4f}% "
f"最大 {100.0 * max(errs):.4f}%")
dl = sum(amb[i] for i in floors) - sum(B[i] for i in floors)
print(f" 地板均值的偏差 {100.0 * dl / sum(B[i] for i in floors):.4f}%"
f"（均值对了，剖面没了：最亮/最暗从 {hi / lo:.4f} 变成 1.0000）")
fv = cglib.fdiv(cglib.fsqrt(cglib.fmul(W, W) + cglib.fmul(H, H)) - H, W)
print(f"解析核验：正对等宽平行段 F = (sqrt(a²+h²)-h)/a，a = h = 1 时 "
f"{fv} = {fv / ONE:.6f}（应为 0.414214）")
if __name__ == "__main__":
main()
跑一遍：
== 8.1 封闭腔：面积光源 vs 环境常数 ==
段数 32，每段长 8192（Q16，即 1/8 = 0.1250），自发光段 #19，反照率 55769 = 0.8510
Gauss-Seidel 400 步后最大残差 0 ulp；互易性抽查 A_i F_ij - A_j F_ji = 0；自发光段行和 65536 = 1.0000
弹射截断：前 K+1 次弹射的部分和对全解的相对偏差（‰，取地板 8 段最大值）
K 最大偏差‰ 第 K 项总量
0 1000.0000 1.000000
1 733.9852 0.850769
2 662.2267 0.719971
4 473.2830 0.511688
8 258.7003 0.251877
16 115.5682 0.043457
20 103.5571 0.005981
24 103.4031 0.000000
级数在第 22 项起被定点截断冻死（T_k 恒为 0），残余停在 103.4031‰——这不是物理弹射次数的极限，是 Q16.16 的极限。
地板 8 段里「2 次以上弹射」的占比：0.7328 ~ 0.7919（单次弹射直接照到的只有 0.2081 ~ 0.2672）
两件事。第一，地板 73% 以上的亮度来自 2 次以上的弹射，只算直接光会把地板压暗四分之三。第二，级数在这个腔里到第 22 项就被 Q16.16 冻死了：T_k 恒为 0，部分和停在 Gauss-Seidel 全解的 89.7% 上。这 103.4‰ 的缺口不是物理截断，是定点精度的下限；要真正看收敛得解线性系统，那正是 8.2 的主题。
真的场景：64×48 上的弹射扫描
config/cg_trace_08_scene.json 描述一个 64×48 的封闭盒子：自发光球悬在天花板下，一颗漫反射球落在地上，左右两面墙分别偏红、偏绿——一个标准的三色渗色（color bleeding）盒子。参考解用 256 spp、最多 32 次弹射、俄罗斯轮盘赌在第 2 次弹射后开启，共 2 178 100 条光线，26 434 毫秒（Apple M1 Pro、单线程 Python）。
扫描组固定 64 spp，只改 maxBounces，与参考解逐像素比：
K 光线数 光线/路径 均值辐射度(Q16) 均值偏差% L1 相对误差% MAE(Q8)
0 196608 1.0000 19656.0000 -41.4226 41.4226 54.0117
1 378240 1.9238 26023.8280 -22.4457 23.2960 30.3895
2 465258 2.3664 29398.1175 -12.3899 16.7578 21.8192
3 506607 2.5767 31229.0822 -6.9334 14.7672 19.2130
4 526723 2.6791 32031.8449 -4.5411 14.0816 18.3154
8 543789 2.7659 33233.1051 -0.9612 13.4933 17.5248
16 544662 2.7703 33392.0155 -0.4876 13.5791 17.6234
参考 2178100 — 33555.6251 — — —
「光线/路径」这一列解释了代价：从 1.0000 涨到 2.7703 就基本封顶——K 从 4 加到 16，一条路径平均只多打 0.0912 条光线，也就是只有 3.3% 的光线发生在第 5 次及以后的弹射上。物理上这很合理：绝大多数路径早早撞回光源或飞出盒子被终止了。
均值偏差随 K 单调收窄（−41.4% → −0.49%），但 L1 相对误差在 K=8 之后不再下降（13.4933 → 13.5791，甚至略微回升）。原因不在弹射次数：扫描组只有 64 spp，L1 的噪声底就在 13% 上下，K 再大也压不下去。这就是为什么后面每一节都要先量参考解自己的噪声——不然会把采样噪声读成截断误差。
P 层：16×12 的逐位复现
// snippet: core/test/cg_trace_08_1_indirect.test.ts（节选）
// verify: node --test 'core/test/cg_trace_08_1_indirect.test.ts'
test("【核心】16x12 对照图的像素缓冲 sha256 与 Python 侧相同", () => {
const buf = totalPng(par.width, par.height, par.maxBounces, par.spp);
assert.equal(createHash("sha256").update(buf).digest("hex"), par.sha256);
assert.deepEqual(Buffer.from(buf), readFileSync(join(root, "out", "cg_trace_08_1_parity.rgb")));
});
test("maxBounces = 0 只剩直接项：间接累加必须恒为 0", () => {
const { indirect } = renderSplit(cfg, planes, spheres, 8, 6, 0, 2, 0n, par.rrStart);
assert.ok(indirect.every((row) => row.every((v) => v === 0n)));
});
test("分解是精确的：direct + indirect 恰好等于总量", () => {
assert.equal(split.directMeanQ16 + split.indirectMeanQ16, split.totalMeanQ16);
assert.equal(split.decomposeMaxDiffQ16, 0);
assert.equal(split.indirectPermille.toFixed(4), "413.2661");
});
这个测试文件里还有一条负例值得单独说：定点下 |x| < 256 时 fmul(x, x) 恒为 0，所以圆盘判据里的 r2 == 0 不是罕见事件——20 万对采样里出现 5 次。给判据加一个下界 0 < r2 会在这 5 处多消耗一对随机数，此后同一像素的整条路径全部错位。在确定性 RNG 上，一个「看起来更严谨」的边界检查就是一颗定时炸弹。
下一节换一个思路：既然截断到第 K 项要等那么久，不如把整个腔体离散成有限块，一次把线性系统解到底。
阅读导航：上一节：7.3 多重重要性采样 · 下一节：8.2 辐射度方法 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。