8.3 光子映射
前两节都在解方程：8.1 把渲染方程展开成级数，8.2 把它离散成线性系统。实时 GI 那篇走的也是「解一个近似系统」的路子，只是把方程压到能实时求（见 实时全局光照 ）。光子映射换一种完全不同的表述——把能量当成粒子。从光源发射 N 个光子，每个光子带着一份功率在场景里随机弹射，落在哪个表面就在那里记一笔（位置 + 入射方向隐含的功率）。要着色时，在着色点附近找出最近的 k 个光子，用它们的总功率除以搜索球的体积，得到该点的辐照度。
代价从「解一个 N×N 的方程」变成「建一棵 kd-tree + 一次 k 近邻查询」，而且天然支持任意几何、任意材质，不需要形状因子，也不需要网格。代价是它有一个绕不开的旋钮：k。
光子怎么发、功率怎么带
自发光球光源的辐射功率是 P = A·π·L_e，其中球面积 A = 4πR²。本节场景里 R = 18350（Q16，约 0.28），L_e 逐通道 [262144, 258048, 253952]：
球面积 A(Q16) 64552（0.984985）
逐通道总功率 P(Q16) [811182, 798507, 785833]
一个关键决定：每个光子携带的是总功率 P，不是 P/N。Q16.16 里 P/N 在 N = 256000 时只剩 2 ulp，乘一次反照率（×0.667 再 >>16）就归零，整张光子图会静默变空。1/N 留到最后一步再除——这就是「一份功率走完全程，只在着色时归一化」。
发射方向用均匀面积 + 余弦加权半球采样：球面上的点按面积均匀取（z = 1 − 2u 均匀、方位角量化到 1024 档），出射方向按 cosθ/π 取。采样器是卷三底座的 SplitMix64（cgtrace.rng_next / rng_q16），确定性、可复现。
最小复现：一维上的 k 近邻密度估计
密度估计器的骨架是 λ̂ = k/(N·2r)——k 近邻半径 r 定带宽，估计值正比于 k/r²。三维里没法把「偏差」和「噪声」分开量，因为参考解自己就有噪声。降到一维、用一个已知解析密度采样就没这个问题了：取 p(x) = 2(1+x)/3，逆 CDF 是 x = −1 + sqrt(1+3u)，真值随取随算。
#!/usr/bin/env python3
# file: lantern/sim/cg_trace_08_3_knn_1d.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_trace_08_3_knn_1d.py
"""8.3 最小复现：一维 k 近邻密度估计——带宽、偏差，以及 Q16.16 的半径下限。
三维光子图的密度估计 L ≈ ρ·ΣΦ'/(N·π²r²) 的骨架就是「k 近邻半径 r 定带宽，估计值
∝ k/r²」。把它降到一维、用一个**已知解析密度**采样，就能把偏差与方差分开量——
三维里做不到，因为参考解本身有噪声。采样密度 p(x) = 2(1+x)/3，逆 CDF 是
x = -1 + sqrt(1+3u)，用 SplitMix64 采样；估计器 λ̂ = k/(N·2r)。
第二件事：k 小时 r 小到 Q16.16 表示不了——r ≈ 11 ulp 时 r² 只有 121 个 Q32.32 单位，
右移 16 位归零，估计器直接除零。对照「半径由 Q16.16 平方距离开方」与「完整 Q32.32」。
"""
import math
import os
import sys
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cglib # noqa: E402
import cgtrace # noqa: E402
ONE = cglib.ONE
F32 = cglib.f_from_int
N = 20000
MAPS = 24
K_LIST = (8, 32, 128, 512)
QUERIES = tuple(cglib.fdiv(F32(t), F32(10)) for t in (1, 3, 5, 7, 9))
SEED = 0x9E3779B97F4A7C15
def true_density(x):
"""Q16 的真概率密度 p(x) = 2(1+x)/3——估计器估的就是它，不是每单位长度的
光子数 N·p(x)。这一条错位会让整张表整体偏 -1/3，而且不报任何错。"""
return cglib.fdiv(cglib.fmul(F32(2), ONE + x), F32(3))
def sample(state):
"""逆 CDF 采样：u -> x = -1 + sqrt(1 + 3u)，全程 Q16。"""
state, u = cgtrace.rng_q16(state)
return state, cglib.fsqrt(ONE + cglib.fmul(F32(3), u)) - ONE
def maps(count):
"""MAPS 张独立光子图，每张 count 个 Q16 位置。"""
state, out = SEED, []
for _ in range(MAPS):
row = []
for _ in range(count):
state, x = sample(state)
row.append(x)
out.append(row)
return out
def ratio(dsq, k, q, count, wide):
"""第 k 近邻平方距离 -> λ̂/p(x)。`wide` 决定半径怎么开方。"""
r = math.isqrt(dsq) if wide else cglib.fsqrt(dsq >> 16)
if r <= 0:
return None
return cglib.fdiv(F32(k), cglib.fmul(F32(count), cglib.fmul(F32(2), r))) / true_density(q)
def stats(vals):
"""(均值偏差‰, 中位偏差‰, 相对标准差‰, 相对 RMSE‰)。均值与中位数分开报：带宽
小时比值分布严重右偏，均值被长尾拽上去，那部分不是偏差而是方差的形状。"""
mean = sum(vals) / len(vals)
var = sum((v - mean) ** 2 for v in vals) / len(vals)
sv = sorted(vals)
return (1000.0 * (mean - 1), 1000.0 * (sv[len(sv) // 2] - 1), 1000.0 * var ** 0.5,
1000.0 * (var + (mean - 1) ** 2) ** 0.5)
def collect(ms, k, count, wide):
out = []
for row in ms:
for q in QUERIES:
v = ratio(sorted((p - q) * (p - q) for p in row)[k - 1], k, q, count, wide)
if v is not None:
out.append(v)
return out
def main():
print("== 8.3 一维 k 近邻密度估计：带宽、偏差与 Q16.16 半径下限 ==")
print(f"采样密度 p(x) = 2(1+x)/3，逆 CDF x = -1 + sqrt(1+3u)，SplitMix64 采样；"
f"N = {N}，{MAPS} 张独立图，探测点 " + " ".join(f"{q / ONE:.1f}" for q in QUERIES)
+ f"；真值 p 从 {true_density(QUERIES[0]) / ONE:.4f} 到 "
f"{true_density(QUERIES[-1]) / ONE:.4f}")
ms = maps(N)
print()
print("半径由 Q16.16 平方距离开方（r² 先 >>16，与三维主脚本同口径）：")
print(f"{'k':>4} {'r²中位(ulp)':>12} {'r=0占比':>9} {'有效查询':>9} {'均值偏差‰':>11} {'标准差‰':>10}")
for k in K_LIST:
r2s, vals = [], collect(ms, k, N, False)
for row in ms:
for q in QUERIES:
r2s.append(sorted((p - q) * (p - q) for p in row)[k - 1] >> 16)
r2s.sort()
head = (f"{k:>4} {r2s[len(r2s) // 2]:>12d} "
f"{100.0 * sum(1 for v in r2s if v == 0) / len(r2s):>8.2f}% {len(vals):>9d} ")
print(head + (f"{stats(vals)[0]:>11.4f} {stats(vals)[2]:>10.4f}" if vals
else f"{'—— 全部除零':>11} {'——':>10}"))
print()
print("半径由完整 Q32.32 平方距离开方（不先 >>16）：")
print(f"{'k':>4} {'均值偏差‰':>11} {'中位偏差‰':>11} {'相对标准差‰':>12} {'相对RMSE‰':>11}")
for k in K_LIST:
b, m, s, e = stats(collect(ms, k, N, True))
print(f"{k:>4} {b:>11.4f} {m:>11.4f} {s:>12.4f} {e:>11.4f}")
print()
ms2 = maps(2000)
print(f"对照：光子数降到 1/10（N = 2000）。网格相对变粗、位置撞车变少，"
f"若 k 小时的偏差是网格造成的，它应当缩小。")
print(f"{'k':>4} {'N=20000 中位偏差‰':>18} {'N=2000 中位偏差‰':>17} {'半径(ulp, N=20000)':>18}")
for k in (8, 32, 128):
rr = sorted(math.isqrt(sorted((p - QUERIES[2]) * (p - QUERIES[2]) for p in row)[k - 1])
for row in ms)
print(f"{k:>4} {stats(collect(ms, k, N, True))[1]:>18.4f} "
f"{stats(collect(ms2, k, 2000, True))[1]:>17.4f} {rr[len(rr) // 2]:>18d}")
print()
print("标准差随 k 下降（每张图用到的光子更多），中位偏差随 k 下降（球越大，位置量化"
"那 1 ulp 的相对扰动越小）；而 k 小到半径掉进 1 ulp 时，估计器不是变差，"
"是直接除零失效。")
if __name__ == "__main__":
main()
跑一遍：
== 8.3 一维 k 近邻密度估计：带宽、偏差与 Q16.16 半径下限 ==
采样密度 p(x) = 2(1+x)/3，逆 CDF x = -1 + sqrt(1+3u)，SplitMix64 采样；N = 20000，24 张独立图，探测点 0.1 0.3 0.5 0.7 0.9；真值 p 从 0.7333 到 1.2667
半径由 Q16.16 平方距离开方（r² 先 >>16，与三维主脚本同口径）：
k r²中位(ulp) r=0占比 有效查询 均值偏差‰ 标准差‰
8 0 100.00% 0 —— 全部除零 ——
32 0 100.00% 0 —— 全部除零 ——
128 0 73.33% 32 68.7697 77.2709
512 10 0.00% 120 29.9253 51.1736
半径由完整 Q32.32 平方距离开方（不先 >>16）：
k 均值偏差‰ 中位偏差‰ 相对标准差‰ 相对RMSE‰
8 281.3178 80.2683 668.9306 725.6775
32 33.7245 9.5408 192.3769 195.3105
128 6.0047 3.4212 94.2595 94.4506
512 5.6176 0.2289 45.6892 46.0333
对照：光子数降到 1/10（N = 2000）。网格相对变粗、位置撞车变少，若 k 小时的偏差是网格造成的，它应当缩小。
k N=20000 中位偏差‰ N=2000 中位偏差‰ 半径(ulp, N=20000)
8 80.2683 35.8822 11
32 9.5408 -0.2168 51
128 3.4212 -6.6553 209
第一张表是 Q16.16 的硬墙：k = 8 时半径只有 11 ulp，r² 在 Q16.16 里直接归零，100% 的查询除零失效。不是「估计得差」，是估计器根本跑不动。第二张表换成完整 Q32.32 平方距离开方，同一批光子立刻能用。
第二张表里两个量朝相反方向走：标准差随 k 下降（669‰ → 192‰ → 94‰ → 46‰，接近 1/√k），中位偏差也随 k 下降（80.3‰ → 9.5‰ → 3.4‰ → 0.23‰）。第三张表解释了偏差的来源：把光子数降到 1/10（网格相对变粗、位置撞车变少），k=32 与 k=128 的中位偏差掉到 0 附近，k=8 也从 80.3‰ 降到 35.9‰。小带宽的偏差主要不是估计器的，是「把 20000 个光子塞进 65536 个整数格」这件事的。
均值偏差那一列（281‰ / 33.7‰ / 6.0‰ / 5.6‰）比中位偏差大得多，且随 k 变小反而变大。这不是矛盾：带宽小时比值分布严重右偏，均值被长尾拽上去——那部分不是偏差，是方差的形状。所以这张表里必须同时报两个数。
真的场景：光子图
回到三色渗色盒子，从自发光球发光子：
N 记录数 光线数 追踪 ms 建树 ms 最大深度 每查询访问节点数
4000 14412 15050 223 18 11 544.4
16000 57721 60247 869 77 13 470.9
64000 230518 240366 3533 352 15 430.0
256000 923260 962575 13738 1763 17 432.4
每条光子的记录数随 N 线性增长（15050/4000 = 3.76，962575/256000 = 3.76），最大深度随 N 缓慢增加——N 越大，随机游走越容易钻进深处，代价是超线性的（追踪时间 223 → 13738 ms，61.6 倍，而 N 只涨 64 倍）。
kd-tree 用中位数分割（每次切当前子区间最宽的轴，叶子上限 8 个光子），k 近邻查询走分支限界。「每查询访问节点数」那一列证明剪枝真的在起作用：N 从 4000 涨到 256000（64 倍），访问节点数只从 544.4 降到 432.4。没有剪枝的话这个数会跟着 N 线性涨。
记录条件写成 b >= 1：直接项被整个丢掉
第一版把记录条件写成了「第 1 次及以后的命中」，理由是「第 0 次是发射事件本身（落在光源上）」。这是错的。 第 0 轮的那条射线已经离开光源、打在场景的第一个表面上，它对应的正是路径追踪里「从光源直接飞到该点」的那部分能量——也就是整个直接项。跳过它，光子图里就只剩 2 阶以上的间接光了。
判据很干脆：光源悬在天花板下方 0.22 处，天花板正对光源那一块拿到的光子数占比应当约 28%。写成 b >= 1 时只剩 18%。这个判据被 P 层回归测试钉住（见下文），比「看起来暗了」可靠得多。
危险之处在于它不报任何错：光子图照建，kd-tree 照查，画面照出，只是整张图暗下去，而且暗的幅度还随场景变化——很难靠肉眼看出来。
精度下限：r² 掉到 3 ulp 之后
改对记录条件之后，出现了一个反直觉的现象：k 小的时候，光子越多反而越差。
N k=8 k=32 k=128 k=512 （相对参考解的均值偏差%）
4000 +17.45 +6.05 +6.02 +10.21
16000 +17.06 +3.90 +2.72 +4.82
64000 +23.84 +5.43 +2.11 +2.55
256000 +63.28 +10.54 +3.15 +1.65
k=512 那一列老实收敛（10.21 → 1.65），k=8 那一列却从 +17.45 一路涨到 +63.28。看 r² 的 ulp 表就明白了：
N k r² 中位(ulp) r² 最小(ulp) r² < 64 ulp 占比
4000 8 258 36 0.25%
16000 8 62 13 51.66%
64000 8 15 2 100.00%
256000 8 3 1 100.00%
256000 32 15 5 100.00%
256000 128 65 27 47.15%
光子越多，k 近邻半径越小，r² 在 Q16.16 里被压成个位数 ulp。而 1/r² 对 r² 是凸函数：把 r² 四舍五入到整数，E[1/r²] > 1/E[r²]，估计值系统性偏高。这不是随机噪声，是有方向的偏差，而且加光子只会让它更严重。
把半径改由完整 Q32.32 平方距离开方（gather_wide）之后：
N k Q16.16 半径% Q32.32 半径%
64000 8 +23.84 +13.52
64000 32 +5.43 +3.53
64000 128 +2.11 +1.64
256000 8 +63.28 +32.01
256000 32 +10.54 +2.92
256000 128 +3.15 +1.36
256000 512 +1.65 +1.19
k=8、N=256000 从 +63.28% 掉到 +32.01%；k=32 从 +10.54% 掉到 +2.92%。换精度比加光子便宜得多。
参考解自己的噪声底
剩下的偏差还有一个来源：参考解不是真值。同一套路径追踪，128 spp 与 512 spp 两次独立渲染逐像素比：
L1 相对差 14.1107%
逐点相对标准差 18.8985%
渲染时间 215244 ms（512 spp）
参考解自己的 L1 就有 14.11%，逐点相对标准差 18.90%。而光子映射在 N = 256000、k = 512 时的 L1 是 17.52%（MAE(Q8) 9.4413、PSNR 28.63 dB；最差的一组是 N=256000/k=8，L1 77.55%、PSNR 16.74 dB）。光子映射的 L1 已经贴着参考解的噪声底了——想再往下压，得先把参考解做得更干净，否则量到的全是参考解的噪声。
带宽的另一半：偏差与方差的实测分解
上面的表都是单张光子图。把同一组参数跑 8 张独立的光子图，在 64 个着色点上量偏差与方差（N = 16000，参考均值 14239.8906）：
k 偏差% 逐点相对标准差% MSE‰
8 +27.9710 49.5049 440.9910
32 +11.8847 22.1599 273.3738
128 +9.3037 10.9580 254.1554
512 +11.4452 5.3048 266.6122
标准差从 49.5% 单调掉到 5.3%（k 越大，每个估计用到的光子越多），偏差却在 9%~12% 之间几乎不动。MSE 的最小值出现在 k = 128（254.16‰）而不是两端——这就是带宽选择的全部内容：k 太小方差压不住，k 太大偏差压不住，最优在中间。
P 层：20000 个光子的光子图逐位复现
8.3 的 P 层把「发射 → 记录 → 建 kd-tree → k 近邻密度估计」整条链用 BigInt 重算：20000 个光子、4 次弹射、k = 32，8 个探测像素的 (r², ΣΦ, 估计值) 连同位置/功率向量的 sha256 必须与 Python 侧相同。
// snippet: core/test/cg_trace_08_3_photon.test.ts（节选）
// verify: node --test 'core/test/cg_trace_08_3_photon.test.ts'
test("【核心】光子图的 sha256 与 Python 侧相同", () => {
const flat = parityFlat();
assert.equal(flat.length, 8 * 8 + 72 + 72);
const s = flat.map((v) => v.toString()).join(",");
assert.equal(createHash("sha256").update(s).digest("hex"), report.parity.sha256);
});
这一条测试同时钉住了 kd-tree 的剪枝距离必须与 d² 同刻度：d² 走 fmul 已经 >>16，剪枝判据里若写成裸乘 d * d，两者差 2^16 倍，判据恒为「太远」，kd-tree 退化成一条根到叶的路径——结果全错，而且不报任何错。
下一章换一个维度：让光线在介质里走，每一步都可能散射或被吸收，渲染方程要重新写一遍。
阅读导航：上一节：8.2 辐射度方法 · 下一节：9.1 体积渲染方程 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。