12.1 噪声与收敛的度量
站内 /graphics-gpu-profiling-debugging/ 那篇做的是 GPU 剖析与调试：抓一帧、看时间轴、定位瓶颈，量的是每帧多少毫秒。本章走一条完全不同的路——离线渲染不关心帧率，它关心的是再加一万个样本，画面还会不会变。而第一步得回答一个看起来很主观的问题：「这张图噪不噪」，能不能变成一个整数？
噪声是一个随机变量
蒙特卡洛渲染里，一个像素的颜色不是一个数，而是一个估计量：
I_N = (1/N) · Σ f(x_i) x_i 从光源面上独立均匀抽出
同一个像素跑两遍会得到两个不同的 I_N——这就是噪声。它没有确定取值，但有确定的方差。把样本记成 x_i，均值和样本方差是
I_N = Σx / N s² = (N·Σx² − (Σx)²) / (N(N−1))
第二式只用整数矩 Σx 与 Σx² 就能算出来。这是本节全部统计量的地基：整数矩可以跨语言逐位比对，浮点的方差不行。本章的三个脚本都只把这些矩落盘，方差、标准误、置信区间各自在需要的地方从矩还原。
于是「噪不噪」至少要分成三个层次来说：
| 层次 | 量 | 本节实测（N = 1） |
|---|---|---|
| 像素级典型值 | 相对偏差的中位数 p50 | 0.044403 |
| 像素级高位 | 相对偏差的 p90 | 0.134614 |
| 图像级 | 逐像素相对偏差的 RMS | 0.181985 |
| 分布尾部 | p99 / max | 0.260089 / 1.538384 |
注意 p90 已经是 RMS 的 0.74 倍，max 又是 RMS 的 8.45 倍。「噪声多大」这句话在没有指定统计量之前是不完整的——同一个像素集合，说 p50 还是说 max，差了 35 倍。这是本节要立的第一条纪律。
锚点：把噪声加在一条已验过的链路上
量噪声之前先得确认被加噪的那条链路是对的。做法是把面光源的两个半边向量 e1、e2 都取成 0——此时采样点恒等于光源中心，估计量退化成点光源，N = 1 的结果必须与卷三基线逐字节相同：
# snippet: sim/cg_trace_12_1_noise.py（节选：锚点这一段）
# verify: /opt/homebrew/bin/python3 sim/cg_trace_12_1_noise.py
flat = dict(C.make_light(spec, cfg), e1=(0, 0, 0), e2=(0, 0, 0))
w, h = spec["width"], spec["height"]
px0 = C.prepare_pixels(spec, tris, bvh, w, h, flat)
sums = []
for i, px in enumerate(px0):
if px is None:
sums.append((0, 0, 0))
continue
st = C.pixel_stream(0x5EED, ex["streamStride"], i)
st, lp = C.sample_light(st, flat["center"], flat["e1"], flat["e2"])
sums.append(C.light_term(px, lp, flat, tris, bvh))
img = C.image_from_pixels(px0, w, h, tuple(spec["backgroundQ8"]), sums, 1)
sha = cgtrace.digest(img)
assert sha == BASELINE_SHA, f"锚点失败：{sha}"
跑出来是 213b4130…，与卷三基线一致。这一步很关键：它把「噪声统计」钉在一条已经逐字节验证过的链路上，之后所有数字都不是「某个渲染器的输出」，而是同一条链路在随机化之后的行为。注意 e1 = e2 = 0 时采样偏移恒为 0，但代码照常消耗两次 rng_q16——流的位置必须与有半径时完全一致，否则锚点只是「碰巧一样」。
一条流跑到底，只记整数矩
扫掠脚本对每个命中像素只开一条流（起点 seed + index·stride，像素之间互不重叠），跑满 256 个样本，在 1, 2, 4, …, 256 共 9 个档位记下 Σx 与 Σx²。一条流而不是九条流，是为了让「N = 256 的估计」严格等于「N = 128 的估计再加上 128 个样本」——收敛曲线上的点必须是同一条路径上的点，否则量出来的斜率里混着不同路径之间的差异。
#!/usr/bin/env python3
# file: lantern/sim/cg_trace_12_1_noise.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_trace_12_1_noise.py
"""12.1 噪声与收敛的度量：把「噪不噪」变成一个可算的整数。
四步：锚点 -> 逐像素 256 样本扫掠 -> 8 批参考解 -> 误差表与 log-log 拟合。
锚点：面光源半径取 0、N = 1 时估计量退化成卷三基线（直接光照 + 硬阴影），320×240
像素缓冲的 sha256 必须等于 `213b4130…`。噪声是加在一条**已验过**的链路上的。
产物：out/cg_trace_12_1_noise.json、out/cg_trace_12_1_noise.csv、
out/cg_trace_12_1_{ref,n1,n256,noisemap}.png
"""
import json
import math
import os
import sys
import time
import numpy as np
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cglib # noqa: E402
import cgtrace # noqa: E402
import cg_trace_12_common as C # noqa: E402
ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
OUT = os.path.join(ROOT, "out")
BASELINE_SHA = "213b413069ca66c9f0cc091bce5328eac0f540d075afa3bab37788aea335ba71"
def anchor(spec, tris, bvh, cfg, ex):
"""半径 0 + N = 1：估计量必须逐字节退化成卷三基线。"""
flat = dict(C.make_light(spec, cfg), e1=(0, 0, 0), e2=(0, 0, 0))
w, h = spec["width"], spec["height"]
px0 = C.prepare_pixels(spec, tris, bvh, w, h, flat)
sums = []
for i, px in enumerate(px0):
if px is None:
sums.append((0, 0, 0))
continue
st = C.pixel_stream(0x5EED, ex["streamStride"], i)
st, lp = C.sample_light(st, flat["center"], flat["e1"], flat["e2"])
sums.append(C.light_term(px, lp, flat, tris, bvh))
img = C.image_from_pixels(px0, w, h, tuple(spec["backgroundQ8"]), sums, 1)
sha = cgtrace.digest(img)
assert sha == BASELINE_SHA, f"锚点失败：{sha}"
return sha
def main():
os.makedirs(OUT, exist_ok=True)
spec = cglib.load_spec()
cfg = C.load_area_cfg()
ex = cfg["experiment"]
levels, n_max = ex["levels"], max(ex["levels"])
w, h = ex["width"], ex["height"]
tris, _ = cgtrace.build_scene_tris(spec)
bvh = cgtrace.build_bvh(tris)
t0 = time.time()
anchor_sha = anchor(spec, tris, bvh, cfg, ex)
anchor_ms = int((time.time() - t0) * 1000)
light = C.make_light(spec, cfg)
pixels = C.prepare_pixels(spec, tris, bvh, w, h, light)
idx = [i for i, p in enumerate(pixels) if p is not None]
m = len(idx)
t0 = time.time()
stats = {"rays": 0, "occluded": 0}
S, Q = C.sweep(pixels, idx, light, tris, bvh, levels, n_max,
ex["seed"], ex["streamStride"], stats)
sweep_ms = int((time.time() - t0) * 1000)
t0 = time.time()
per, tot = C.batch_sums(pixels, idx, light, tris, bvh, ex["batches"],
ex["batchSamples"], ex["seed"], ex["batchStride"],
ex["streamStride"])
ref_ms = int((time.time() - t0) * 1000)
nb = ex["batches"]
ref = tot // (nb * ex["batchSamples"])
bmean = per / ex["batchSamples"]
ref_sem = bmean.std(axis=0, ddof=1) / math.sqrt(nb)
ref_lum = ref.mean()
rows, xs, ys = [], [], []
for n in levels:
mean = S[n] / n
# 样本方差 s² = (n·Σx² − (Σx)²) / (n(n−1))，全整数进、浮点出
var = ((Q[n] * n - S[n] * S[n]).astype(np.float64) / (n * (n - 1))
if n > 1 else np.zeros_like(mean, dtype=np.float64))
sem = np.sqrt(np.maximum(var, 0) / n) # 流内预测的「均值的标准差」
err = mean - ref
pred = float(np.sqrt((sem ** 2).mean()) / ref_lum)
act = float(np.sqrt((err ** 2).mean()) / ref_lum)
q = C.rel_stats(np.abs(err), ref)
pq = C.rel_stats(sem, ref)
rows.append({"n": n, "predictedRelRms": round(pred, 8),
"actualRelRms": round(act, 8), "actualP50": q["p50"],
"actualP90": q["p90"], "actualP99": q["p99"], "actualMax": q["max"],
"predictedP99": pq["p99"], "predictedPixels": pq["pixels"],
"mean8bitRmse": round(float(np.sqrt(
(((mean - ref) * 255.0 / 65536.0) ** 2).mean())), 6)})
if n >= 4:
xs.append(math.log10(n))
ys.append(math.log10(act))
slope, icept, resid = C.fit_slope(xs, ys)
ref_img = C.build_image(pixels, idx, w, h, tuple(spec["backgroundQ8"]), ref)
n1_img = C.build_image(pixels, idx, w, h, tuple(spec["backgroundQ8"]), S[1])
n256_img = C.build_image(pixels, idx, w, h, tuple(spec["backgroundQ8"]),
S[n_max] // n_max)
noise = np.zeros((h, w), dtype=np.int64)
for k, i in enumerate(idx):
y, x = divmod(i, w)
noise[y, x] = int(np.abs(S[1][k] - ref[k]).max())
cglib.save_png(os.path.join(OUT, "cg_trace_12_1_ref.png"), ref_img)
cglib.save_png(os.path.join(OUT, "cg_trace_12_1_n1.png"), n1_img)
cglib.save_png(os.path.join(OUT, "cg_trace_12_1_n256.png"), n256_img)
cglib.save_gray_png(os.path.join(OUT, "cg_trace_12_1_noisemap.png"), noise)
cglib.write_csv(os.path.join(OUT, "cg_trace_12_1_noise.csv"),
["n", "predictedRelRms", "actualRelRms", "actualP50", "actualP90",
"actualP99", "actualMax", "predictedP99", "mean8bitRmse"],
[[r["n"], r["predictedRelRms"], r["actualRelRms"], r["actualP50"],
r["actualP90"], r["actualP99"], r["actualMax"], r["predictedP99"],
r["mean8bitRmse"]] for r in rows])
report = {
"width": w, "height": h, "hitPixels": m, "levels": levels, "nMax": n_max,
"anchor": {"sha256": anchor_sha, "baseline": BASELINE_SHA,
"match": anchor_sha == BASELINE_SHA, "ms": anchor_ms},
"sweep": {"ms": sweep_ms, "samples": m * n_max, "shadowRays": stats["rays"],
"occluded": stats["occluded"],
"samplesPerSecond": round(m * n_max / (sweep_ms / 1000.0))},
"reference": {"batches": nb, "samplesPerBatch": ex["batchSamples"],
"totalSamples": nb * ex["batchSamples"], "ms": ref_ms,
"refLumaQ16": round(float(ref_lum), 4),
"relStandardError": round(float(ref_sem.mean() / ref_lum), 8),
"sha256": cgtrace.digest(ref_img)},
"rows": rows,
"fit": {"slope": round(slope, 6), "intercept": round(icept, 6),
"maxResidual": round(resid, 6), "over": "N >= 4"},
"noiseMap": {"min": int(noise.min()), "max": int(noise.max()),
"span": max(1, int(noise.max()) - int(noise.min()))},
}
with open(os.path.join(OUT, "cg_trace_12_1_noise.json"), "w", encoding="utf-8") as fh:
json.dump(report, fh, ensure_ascii=False, indent=2)
fh.write("\n")
print("== 12.1 噪声与收敛的度量 ==")
print(f"锚点（半径 0，N = 1，320x240）: {anchor_sha[:16]}… == 基线 "
f"{anchor_sha == BASELINE_SHA}（{anchor_ms} ms）")
print(f"命中像素 {m} / {w * h}；扫掠 {m * n_max} 样本，阴影光线 {stats['rays']}"
f"（遮挡 {stats['occluded']}），{sweep_ms} ms")
print(f"参考解：{nb} 批 × {ex['batchSamples']} 样本 = {nb * ex['batchSamples']} 样本/像素，"
f"{ref_ms} ms，平均亮度 {ref_lum:.2f} Q16.16，相对标准误 "
f"{ref_sem.mean() / ref_lum:.6f}")
print(f"{'N':>5} {'预测RMS':>12} {'实际RMS':>12} {'p99':>10} {'max':>10} {'8位RMSE':>10}")
for r in rows:
print(f"{r['n']:>5} {r['predictedRelRms']:>12.6f} {r['actualRelRms']:>12.6f} "
f"{r['actualP99']:>10.4f} {r['actualMax']:>10.4f} {r['mean8bitRmse']:>10.4f}")
print(f"log-log 斜率（N ≥ 4）= {slope:+.4f}（理论 -0.5），最大残差 {resid:.4f}")
print(f"噪声图跨度 {int(noise.min())} … {int(noise.max())}")
return 0
if __name__ == "__main__":
sys.exit(main())
跑一遍：
# snippet: 本机实跑（工作区 lantern/）
cd lantern
/opt/homebrew/bin/python3 sim/cg_trace_12_1_noise.py
== 12.1 噪声与收敛的度量 ==
锚点（半径 0，N = 1，320x240）: 213b413069ca66c9… == 基线 True（6595 ms）
命中像素 2734 / 4800；扫掠 699904 样本，阴影光线 540692（遮挡 7835），28621 ms
参考解：8 批 × 64 样本 = 512 样本/像素，57370 ms，平均亮度 4966.85 Q16.16，相对标准误 0.005313
N 预测RMS 实际RMS p99 max 8位RMSE
1 0.000000 0.181985 0.2601 1.5384 3.5170
2 0.129050 0.132530 0.1995 1.1853 2.5613
4 0.091551 0.093531 0.1332 1.1908 1.8076
8 0.063888 0.067887 0.1142 0.4855 1.3120
16 0.045857 0.047171 0.0838 0.3445 0.9116
32 0.032598 0.032953 0.0583 0.3489 0.6368
64 0.022906 0.021952 0.0417 0.3063 0.4243
128 0.016201 0.015823 0.0295 0.1626 0.3058
256 0.011483 0.012326 0.0246 0.1209 0.2382
log-log 斜率（N ≥ 4）= -0.5027（理论 -0.5），最大残差 0.0269
噪声图跨度 0 … 8088
实测收敛速率：-0.5027
把 N ≥ 4 的七个档位取 log-log 做最小二乘，斜率 -0.5027，理论值 -0.5，最大残差 0.0269（log10 单位）。蒙特卡洛的 O(N^(-1/2)) 在这里不是教科书上的一句话，而是用整数矩实测出来的一个数。
再看「预测 RMS」与「实际 RMS」两列——预测值由流内样本方差给出，实际值是与参考解比出来的：
| N | 预测 RMS | 实际 RMS | 实际/预测 |
|---|---|---|---|
| 2 | 0.129050 | 0.132530 | 1.027 |
| 16 | 0.045857 | 0.047171 | 1.029 |
| 64 | 0.022906 | 0.021952 | 0.958 |
| 256 | 0.011483 | 0.012326 | 1.073 |
注意输出里 N = 1 那一行的预测值是 0.000000——不是因为它准，而是因为一个样本估不出方差，s² 按定义为 0。这正是下一节 12.3 里「最小样本数」问题的根：方差估计本身是个随机变量，样本太少时它会骗人。
另一个值得注意的现象在最后一行：N = 256 时实际误差大于预测误差，比值 1.073。把参考解自身的相对标准误 0.005313 按平方和叠加进来，sqrt(0.011483² + 0.005313²) = 0.012653，与实际值 0.012326 同量级。也就是说参考解自己已经成了误差地板——再怎么加样本，这个比值也降不到 1 以下。这件事 12.2 会单独量。
图像级与像素级：为什么两个数差这么多
8 位通道上的 RMSE 从 N = 1 的 3.5170 降到 N = 256 的 0.2382，比值 14.76，与理想的 sqrt(256) = 16 同量级。这里不能要求它精确等于 16：逐像素的方差本身就不同，参考解在两端各贡献一份不确定度，8 位转换（乘 255 再降精度）又引入一点截断——所以「收敛速率」这个结论应该从 log-log 拟合去读（-0.5027），而不是从两个端点的比值去读。
而「噪声图跨度 0 … 8088」是另一回事：那是逐像素 |N=1 的估计 − 参考解| 的最大值，8088 Q16.16 ≈ 0.1234（线性光）。它和 8 位 RMSE 3.5170 描述的是同一批像素，但一个是最大值、一个是均方根，换算到同一量纲后差了近 9 倍。这两个数都不是错的，它们回答的是不同的问题：
- 「最坏的那个像素错多少」→ max；
- 「整体上平均错多少」→ RMS。
写质量报告时如果只报一个，就必须说清楚报的是哪个。本节把五个统计量（p50 / p90 / p99 / max / RMS）全部落盘，正是为了让「噪声多大」这句话每次都能落到具体的量上。
P 层把整数矩钉死
P 层用 TypeScript 重跑同一批 SplitMix64 子流，全程 BigInt，断言对象是整数矩本身，不是浮点的方差：
# snippet: core/test/cg_trace_12_1_noise.test.ts（节选：梯子与负例）
# verify: node --test 'core/test/cg_trace_12_1_noise.test.ts'
const got = ladder(px, light, tris, bvh, cfg.budget, cfg.checkEvery, seed, stride, i);
const want = cfg.ladder[String(i)];
for (let li = 0; li < want.length; li++) {
for (let k = 0; k < 3; k++) {
assert.equal(got.sum[li][k].toString(), String(want[li][k]));
assert.equal(got.sumSq[li][k].toString(), String(want[li][k + 3]));
}
}
// 负例：子流起点只差 1，但 SplitMix64 的雪崩性使两次结果毫无关系
assert.ok(a.sum[last][0] !== b.sum[last][0]);
8 个抽样像素 × 64 个检查点 × 3 通道 × 2 个矩 = 3072 个整数逐个相等，整个测试 1.6 秒跑完。最后那条负例是必要的：只测「同种子同结果」不够，还得测换种子必须换结果，否则一个返回常量的实现也能通过。
小结
- 噪声是随机变量的方差，不是观感形容词；本节把它拆成 p50 / p90 / p99 / max / RMS 五个可落盘的整数比。
- 锚点
213b4130…保证加噪的那条链路与卷三基线逐字节相同。 - 一条流跑满 256 样本，实测收敛斜率 -0.5027（理论 -0.5），最大残差 0.0269。
- N = 256 时实际误差是预测值的 1.073 倍——参考解自身的不确定度已经成了地板。
下一节把「参考解」这件事单独拿出来量：它自己有多少不确定度，凭什么当基准。
阅读导航：上一节：11.3 SDF 的法线与软阴影 · 下一节：12.2 参考解与误差 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。