8.2 辐射度方法
站内 /graphics-real-time-gi/ 那篇做的是实时 GI：用探针与光照缓存把间接光近似出来，关心的是每帧能不能塞进预算。本节换一条路——不近似，也不采样：把腔体切成有限块，让每一块的辐出度成为未知量，于是渲染方程变成一个线性方程组，一次解到底。上一节已经看到，按弹射次数展开的级数在 Q16.16 里被定点精度冻死在 89.7% 上，而「解一个线性系统」不存在这个问题。
从积分方程到线性方程组
把场景表面切成 N 块（patch），假设每块上的辐出度 B 与反照率 ρ 都是常数。第 i 块收到的能量来自所有能看见它的块：
B_i = E_i + ρ_i · Σ_j F_ij · B_j
E_i 是自发光，F_ij 是形状因子（form factor）——从第 i 块出发、最终落在第 j 块上的能量比例。写成矩阵就是 B = E + ρ⊙(F·B)，或者 (I − ρF)B = E。第 i 行所有形状因子之和必须等于 1：第 i 块朝外辐射的能量总得落到某个地方。
这一节的成败全压在 F 上。F 一旦不准，后面迭代再稳也只是在稳稳地解一个错方程。
最小复现：正方形腔的四元系统
先用一个能逐项核对的最小实例把方法立起来：正方形腔，每面墙一块 patch，四块两两可见。正方形腔的形状因子有闭式解——相邻 (2−√2)/2 = 0.292893，正对 √2−1 = 0.414214，每行和恰好 1。
#!/usr/bin/env python3
# file: lantern/sim/cg_trace_08_2_four_patch.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_trace_08_2_four_patch.py
"""8.2 最小复现：正方形腔、每面墙一块 patch 的四元辐射度系统。
四块 patch 两两可见，形状因子可以手工验证：正方形腔里 F(相邻) = (2-√2)/2 = 0.292893、
F(正对) = √2-1 = 0.414214，行和恰好 1。整套方法（形状因子矩阵 -> 迭代求解 -> 收敛
判据）在这个 4 元系统上全部看得见，而它小到可以逐项核对。
三个对照：Hottel 精确形状因子 vs 主脚本用的质心近似 cosθi·cosθj·Aj/(πr²)；行和是否
为 1，以及不做行归一化时解系统性偏低多少；Jacobi 与 Gauss-Seidel 各要多少轮。
全程 Q16.16 整数，除法一律 cglib.fdiv（向零截断），不用裸 `/` 或 `//`。
"""
import os
import sys
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cglib # noqa: E402
ONE = cglib.ONE
PI = 205887 # π 的 Q16.16
RHO = cglib.fdiv(cglib.f_from_int(217), cglib.f_from_int(255))
TOL = 2 # 收敛判据：相邻两轮最大变化 <= 2 ulp
F32 = cglib.f_from_int
HALF = cglib.fdiv(ONE, F32(2))
# 正方形腔四条边，逆时针：下 -> 右 -> 上 -> 左
SEGS = [((0, 0), (ONE, 0)), ((ONE, 0), (ONE, ONE)), ((ONE, ONE), (0, ONE)), ((0, ONE), (0, 0))]
CENT = [(HALF, 0, 0), (ONE, HALF, 0), (HALF, ONE, 0), (0, HALF, 0)]
NORM = [(0, ONE, 0), (-ONE, 0, 0), (0, -ONE, 0), (ONE, 0, 0)]
def dist(p, q):
return cglib.fsqrt(cglib.fmul(q[0] - p[0], q[0] - p[0])
+ cglib.fmul(q[1] - p[1], q[1] - p[1]))
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, F32(2))
return v if v > 0 else 0
def exact_matrix():
return [[0 if i == j else hottel(SEGS[i], SEGS[j]) for j in range(4)] for i in range(4)]
def centroid_matrix():
"""主脚本的核：cosθi·cosθj·Aj/(π r²)，r 取两质心距离。"""
out = [[0] * 4 for _ in range(4)]
for i in range(4):
for j in range(4):
if i == j:
continue
d = (CENT[j][0] - CENT[i][0], CENT[j][1] - CENT[i][1], 0)
r = cglib.v_len(d)
ci = cglib.fdiv(cglib.v_dot(NORM[i], d), r)
cj = cglib.fdiv(cglib.v_dot(NORM[j], (-d[0], -d[1], 0)), r)
if ci > 0 and cj > 0:
out[i][j] = cglib.fdiv(cglib.fmul(cglib.fmul(ci, cj), ONE),
cglib.fmul(PI, cglib.fmul(r, r)))
return out
def iterate(F, E, in_place):
"""Gauss-Seidel（in_place=True）或 Jacobi（False），返回 (解, 轮数)。"""
B = [0] * 4
for used in range(1, 401):
old = list(B)
delta = 0
for i in range(4):
src = B if in_place else old
nb = E[i] + cglib.fmul(RHO, sum(cglib.fmul(F[i][j], src[j]) for j in range(4)))
delta = max(delta, abs(nb - old[i]))
B[i] = nb
if delta <= TOL:
return B, used
return B, 400
def row_normalize(F):
return [[cglib.fdiv(F[i][j], sum(F[i])) for j in range(4)] for i in range(4)]
def dump(tag, M):
for i in range(4):
print(f" {tag} " + " ".join(f"{M[i][j] / ONE:.6f}" for j in range(4))
+ f" 行和 {sum(M[i]) / ONE:.6f}")
def main():
Fx, Fc = exact_matrix(), centroid_matrix()
Fn = row_normalize(Fc)
E = [ONE, 0, 0, 0]
print("== 8.2 四元辐射度系统（正方形腔，每面一块 patch）==")
dump("精确 ", Fx)
dump("质心 ", Fc)
sq2 = cglib.fsqrt(F32(2))
adj = cglib.fdiv(F32(2) - sq2, F32(2))
opp = cglib.fdiv(cglib.fmul(F32(2), sq2) - F32(2), F32(2))
print(f" 解析核对：(2-√2)/2 = {adj / ONE:.6f}（0.292893）；√2-1 = {opp / ONE:.6f}"
f"（0.414214）；与矩阵的差 {Fx[0][1] - adj} / {Fx[0][2] - opp} ulp")
print(f" 互易性 max|A_i F_ij - A_j F_ji|：精确 "
f"{max(abs(Fx[i][j] - Fx[j][i]) for i in range(4) for j in range(4))}，"
f"质心 {max(abs(Fc[i][j] - Fc[j][i]) for i in range(4) for j in range(4))}")
print()
print("同一个源项（patch 0 发射 1.0，反照率 0.8510）在三种矩阵下的解：")
Be, ie = iterate(Fx, E, True)
Bj, ij = iterate(Fx, E, False)
Bn, _ = iterate(Fn, E, True)
Bc, ic = iterate(Fc, E, True)
for tag, B in (("精确 ", Be), ("质心+行归一 ", Bn), ("质心不归一 ", Bc)):
print(f" {tag}: " + " ".join(f"{B[i] / ONE:.6f}" for i in range(4)))
print(f" 不归一 / 归一：逐项 {min(1000.0 * Bc[i] / Bn[i] for i in range(4)):.2f}"
f" ~ {max(1000.0 * Bc[i] / Bn[i] for i in range(4)):.2f}‰（系统性偏低）")
print(f" 归一 / 精确：逐项 {min(1000.0 * Bn[i] / Be[i] for i in range(4)):.2f}"
f" ~ {max(1000.0 * Bn[i] / Be[i] for i in range(4)):.2f}‰（把行和补回来就接近了）")
tot = sum(Be)
ex = cglib.fdiv(ONE, ONE - RHO)
print(f" 能量守恒：Σ B_i = {tot / ONE:.6f}，解析 1/(1-ρ) = {ex / ONE:.6f}，"
f"相对差 {1000.0 * (tot - ex) / ex:.4f}‰")
print()
print(f"收敛（容差 {TOL} ulp）：精确矩阵 Jacobi {ij} 轮 / Gauss-Seidel {ie} 轮"
f"（比值 {ij / ie:.2f}）；两解逐项差 "
f"{max(abs(Be[i] - Bj[i]) for i in range(4))} ulp")
print(f" 质心不归一 Jacobi {ic} 轮 / 质心归一 Jacobi {iterate(Fn, E, False)[1]} 轮："
f"行和缺口把谱半径压小，迭代反而更快——但它收敛到的是错解")
if __name__ == "__main__":
main()
跑一遍：
== 8.2 四元辐射度系统（正方形腔，每面一块 patch）==
精确 0.000000 0.292892 0.414200 0.292892 行和 0.999985
精确 0.292892 0.000000 0.292892 0.414200 行和 0.999985
精确 0.414200 0.292892 0.000000 0.292892 行和 0.999985
精确 0.292892 0.414200 0.292892 0.000000 行和 0.999985
质心 0.000000 0.318329 0.318298 0.318329 行和 0.954956
质心 0.318329 0.000000 0.318329 0.318298 行和 0.954956
质心 0.318298 0.318329 0.000000 0.318329 行和 0.954956
质心 0.318329 0.318298 0.318329 0.000000 行和 0.954956
解析核对：(2-√2)/2 = 0.292892（0.292893）；√2-1 = 0.414200（0.414214）；与矩阵的差 0 / 0 ulp
互易性 max|A_i F_ij - A_j F_ji|：精确 0，质心 0
同一个源项（patch 0 发射 1.0，反照率 0.8510）在三种矩阵下的解：
精确 : 2.264954 1.458969 1.525574 1.458969
质心+行归一 : 2.261322 1.482300 1.482300 1.482315
质心不归一 : 1.924240 1.137390 1.137375 1.137390
不归一 / 归一：逐项 767.30 ~ 850.94‰（系统性偏低）
归一 / 精确：逐项 971.63 ~ 1016.00‰（把行和补回来就接近了）
能量守恒：Σ B_i = 6.708466，解析 1/(1-ρ) = 6.709930，相对差 -0.2183‰
收敛（容差 2 ulp）：精确矩阵 Jacobi 57 轮 / Gauss-Seidel 32 轮（比值 1.78）；两解逐项差 9 ulp
质心不归一 Jacobi 24 轮 / 质心归一 Jacobi 58 轮：行和缺口把谱半径压小，迭代反而更快——但它收敛到的是错解
三个可直接用的判据：
- 互易性
A_i F_ij = A_j F_ji在两种核下都精确为 0——这是构造性的，不能用来挑核。 - 行和才是分水岭。精确核每行 0.999985（只差 1 ulp）；质心核每行 0.954956，缺口 4.5%。形状因子的「总量守恒」在质心近似下被打破了。
- 行和缺口把谱半径压小，迭代反而更快（24 轮 vs 58 轮），但收敛到的是错解——比精确核低 15%~23%。收敛快不是好消息。
真的场景：384 块、六面墙
回到三色渗色盒子。辐射度版本把六面墙按 patchRes 细分，天花板中央一块自发光，用质心核 F_ij ≈ cosθi·cosθj·Aj/(πr²) 建矩阵：
n patch 数 行和最大缺口(ulp) 行和平均缺口(ulp) 互易性 Jacobi 轮 GS 轮 GS 真残差 raw GS 轮
2 24 17880 17880 0 23 13 1 25
3 54 16255 11842 0 24 14 1 20
4 96 15868 8840 0 24 14 1 18
6 216 15631 5813 0 25 14 2 17
8 384 15488 4253 0 25 14 1 17
最大缺口几乎不随细分变化（17880 → 15488，约 27.3% → 23.6%），平均缺口却从 17880 掉到 4253。原因很直白：角上的小块几乎收不到正对方向的能量，fdiv 的向零截断把它们的形状因子吃成 0；块越细，这种小块越多，但它们在行里占的比例越小。缺口是分布不均的，均值会骗人。
收敛本身很干脆：容差 2 ulp 下 Jacobi 2325 轮、Gauss-Seidel 1314 轮，Gauss-Seidel 的首末残差从 262144 掉到 2。GS 比 Jacobi 少一半轮数，因为它是就地更新、每轮用上刚算出来的新值。
行归一化：补回来，但只是补回来
把每行除以行和（normalizeRows）之后：
n patch 数 归一化后能量平衡% 不归一能量平衡%
2 24 -18.3249 +34.0499
3 54 -18.2113 +11.1056
4 96 -18.1870 +2.2818
6 216 -18.2848 -5.5097
8 384 -18.4500 -9.2242
不归一的那一列在 n=2 时多出 34.0% 的能量（远处小块的形状因子被截断吃成 0，行和偏大），细到 n=8 又变成少 9.2%——误差换了号。归一化之后恒定在 −18.2% 附近，符号与量级都稳定了。
那个 −18.2% 不是归一化的错：它是辐射度模型自身的固有偏差——能量被强制留在腔体里，而路径追踪允许一部分能量从相机那一面（front 面）离开。两种口径本来就不该相等。真正要看的是它稳不稳。
质心核 vs 蒙特卡洛核
质心核是一阶近似。同一套 patch 用蒙特卡洛积分（每对 2048 个样本，rowSumMcMin/Max = 2048/2048）算形状因子，逐对比：
patch 数 96
样本/对 2048
逐对绝对差 均值 216.9245（Q16，≈0.00331）
逐对绝对差 最大 9693
非零对 7680
相对差 > 1% 的对 5640（均值相对差 18.0709%）
二项噪声本身 132.2775
二项分布给出的噪声底是 132.28；质心核与 MC 的平均绝对差是 216.92，只比噪声底大 64%。在 96 块这个分辨率下，质心核的误差和蒙特卡洛自己的噪声是一个量级——想再准就得加样本，而不是换公式。
细分能收敛吗
把一块 96×96 的板从 2×2 一路切到 128×128，用一行固定的探针形状因子（解析值 1.0）看行和：
细分 2 4 8 16 32 64 128
Q16 37820 63954 81024 68274 63192 54382 36863
浮点 0.5773 0.9765 1.2391 1.0537 1.0120 1.0028 1.0007
浮点那一行从 0.5773 单调爬到 1.0007，是教科书里的收敛曲线。Q16 那一行却在 8×8 处冲到 1.239（比真值大 24%），然后一路掉回 0.5625。原因是分子分母都在变小：块越小，cosθi·cosθj·Aj 这个乘积越快掉到 1 ulp 以下，fdiv 把它截成 0，行和跟着塌。细分不是免费的：Q16.16 能表达的形状因子有一个分辨率下限，切过那个下限，精度不升反降。
P 层：54 块的解向量逐位复现
8.2 的 P 层用 BigInt 重算 patch → 形状因子 → Gauss-Seidel 全链，n=3（54 块）的解向量 sha256 必须与 Python 侧相同。
// snippet: core/test/cg_trace_08_2_radiosity.test.ts（节选）
// verify: node --test 'core/test/cg_trace_08_2_radiosity.test.ts'
test("【核心】辐射度解向量的 sha256 与 Python 侧相同", () => {
const { values, iters } = solve();
assert.equal(values.length, report.parity.patches * 3);
assert.equal(iters, report.parity.gsIters);
const s = values.map((v) => v.toString()).join(",");
assert.equal(createHash("sha256").update(s).digest("hex"), report.parity.sha256);
});
test("解向量逐项相同（不只是哈希）", () => {
const { values } = solve();
assert.deepEqual(values.map(Number), report.parity.values);
});
test("互易性 A_i F_ij = A_j F_ji 在质心近似下精确成立", () => {
下一节换一种完全不同的做法：不再解方程，而是把能量当成一颗颗可以搬运的粒子，扔进场景里弹，最后在着色点附近数一数——光子映射。
阅读导航：上一节：8.1 间接光照 · 下一节：8.3 光子映射 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。