2.3 容差与自交
站内 /graphics-raytracing-hybrid/
那篇是实时混合光追，阴影偏置靠手调常数。本卷做离线光追的完整链路，本节给的是求交的容差边界分析：EPSILON 不是一个可以试出来的魔数，它被两条可测量的边界夹住。
两条边界
命中点 p 重发阴影光线时，起点不能正好是 p——定点求交的命中点带着几百分之一世界单位的误差，从 p 出发的阴影光线有概率重新打回自己所在的三角形，把每个面都判成背光（acne）。所以起点要沿法线偏移一个 EPSILON。这个量被两侧夹住：
- 下界（必须大于）：命中点的数值不确定度。小于它，自交清不掉。
- 上界（必须小于）：场景里最小的可见几何间隙。大于它，偏移会把起点推过另一个表面，撕出可见缝隙。
完整脚本
#!/usr/bin/env python3
# file: lantern/sim/cg_trace_02_3_epsilon.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_trace_02_3_epsilon.py
"""2.3 自交与容差：EPSILON 的两条边界，以及为什么沿法线而不是沿光线偏移。
1. 几何间隙：球顶点到地形网格的最小距离（实测 174 Q16.16 = 0.002655 世界单位，
出现在中心沉入地形的球 1 的接触环上）。EPSILON 必须**小于**它，否则沿法线的
偏移会把阴影光线起点推过另一个表面，在球-地形接触处撕出可见缝隙。
2. 自交：命中点重发阴影光线时打回自己所在的三角形。EPSILON 必须**大于**命中点的
数值不确定度，否则每个面都判成背光（acne）。地形自交在 eps >= 16 清零。
3. 沿法线 vs 沿光线：掠射命中处 |n·d| 很小，沿光线偏移的法向分量趋零，清不掉表面。
产物：out/cg_trace_02_3_epsilon.json
"""
import json
import os
import sys
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cglib # noqa: E402
import cgtrace # noqa: E402
ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
OUT = os.path.join(ROOT, "out")
ONE, FMUL, FDIV = cglib.ONE, cglib.fmul, cglib.fdiv
EPS_SWEEP = [0, 1, 8, 16, 32, 64, 256, 1024, 4096, 16384, 32768, 65536]
TERRAIN_TRIS = 512
def ptd2(p, a, b, c):
"""点到三角形的最小平方距离（Q16.16 进，Q32.32 平方出）。Ericson 的 Voronoi 区判定。
标度：`d` 把两个 Q16.16 乘成 Q32.32，故边参数必须显式左移 16 位成 Q16.16。
"""
s = lambda u, v: (u[0] - v[0], u[1] - v[1], u[2] - v[2])
d = lambda u, v: u[0] * v[0] + u[1] * v[1] + u[2] * v[2]
at = lambda o, e, t: tuple(o[i] + ((e[i] * t) >> 16) for i in range(3))
ab, ac, ap = s(b, a), s(c, a), s(p, a)
d1, d2 = d(ab, ap), d(ac, ap)
if d1 <= 0 and d2 <= 0:
return d(ap, ap)
bp = s(p, b)
d3, d4 = d(ab, bp), d(ac, bp)
if d3 >= 0 and d4 <= d3:
return d(bp, bp)
vc = d1 * d4 - d3 * d2
if vc <= 0 and d1 >= 0 and d3 <= 0:
q = at(a, ab, (d1 << 16) // (d1 - d3) if d1 != d3 else 0)
return d(s(p, q), s(p, q))
cp = s(p, c)
d5, d6 = d(ab, cp), d(ac, cp)
if d6 >= 0 and d5 <= d6:
return d(cp, cp)
vb = d5 * d2 - d1 * d6
if vb <= 0 and d2 >= 0 and d6 <= 0:
q = at(b, s(c, b), (d2 << 16) // (d2 - d6) if d2 != d6 else 0)
return d(s(p, q), s(p, q))
va = d3 * d6 - d5 * d4
if va <= 0 and d4 - d3 >= 0 and d5 - d6 >= 0:
q = at(b, s(c, b), ((d4 - d3) << 16) // ((d4 - d3) + (d5 - d6)))
return d(s(p, q), s(p, q))
den = va + vb + vc
if den == 0:
return d(ap, ap)
wq = (vc << 16) // den
q = at(a, ab, (vb << 16) // den)
q = (q[0] + ((ac[0] * wq) >> 16), q[1] + ((ac[1] * wq) >> 16), q[2] + ((ac[2] * wq) >> 16))
return d(s(p, q), s(p, q))
def contact_gap(spec, tris):
"""每个球的顶点到地形网格的最小距离（Q16.16），以及中心到地形的竖直间隙。"""
import math
per = []
for si, sp in enumerate(spec["spheres"]):
sv, _ = cglib.build_sphere(sp["centerQ16"], sp["radiusQ16"], sp["subdiv"])
md2 = min(ptd2(v, (t[0], t[1], t[2]), (t[3], t[4], t[5]), (t[6], t[7], t[8]))
for v in sv for t in tris[:TERRAIN_TRIS])
tr = spec["terrain"]
h = cglib.terrain_height(sp["centerQ16"][0], sp["centerQ16"][2],
tr["noiseSeed"], tr["heightScaleQ16"])
per.append({"sphere": si, "minDistQ16": math.isqrt(md2), "vertices": len(sv),
"verticalGapQ16": sp["centerQ16"][1] - sp["radiusQ16"] - h})
return per
def collect_primary(spec, tris, nodes, order):
"""跑一遍主光线：返回命中缓存 (p, n, ti, dirv) 与命中点 |n·d| 的分布。"""
cam = spec["camera"]
eye = tuple(cam["eyeQ16"])
view = cglib.mat_look_at(eye, tuple(cam["targetQ16"]), tuple(cam["upQ16"]))
right, up, back = ((view[0], view[1], view[2]), (view[4], view[5], view[6]),
(view[8], view[9], view[10]))
w, h = spec["width"], spec["height"]
th = cam["tanHalfFovYQ16"]
tx = FMUL(th, FDIV(cglib.f_from_int(w), cglib.f_from_int(h)))
tmax = cam["farQ16"]
buckets = {"<0.01": 0, "0.01-0.1": 0, "0.1-0.5": 0, "0.5-0.9": 0, ">=0.9": 0}
names, hits, mn = list(buckets), [], None
for py in range(h):
ay = FMUL(FDIV(cglib.f_from_int(2 * (h // 2 - py) - 1), cglib.f_from_int(h)), th)
for px in range(w):
ax = FMUL(FDIV(cglib.f_from_int(2 * (px - w // 2) + 1), cglib.f_from_int(w)), tx)
dirv = cglib.v_norm((
FMUL(ax, right[0]) + FMUL(ay, up[0]) - back[0],
FMUL(ax, right[1]) + FMUL(ay, up[1]) - back[1],
FMUL(ax, right[2]) + FMUL(ay, up[2]) - back[2]))
hit = cgtrace.bvh_intersect(nodes, order, tris, eye, dirv, tmax)
if hit is None:
continue
t, ti, u, v = hit
rec = tris[ti]
w0 = ONE - u - v
n = cglib.v_norm((FMUL(rec[9], w0) + FMUL(rec[12], u) + FMUL(rec[15], v),
FMUL(rec[10], w0) + FMUL(rec[13], u) + FMUL(rec[16], v),
FMUL(rec[11], w0) + FMUL(rec[14], u) + FMUL(rec[17], v)))
if cglib.v_dot(n, dirv) > 0:
n = (-n[0], -n[1], -n[2])
hits.append((cglib.v_add(eye, (FMUL(dirv[0], t), FMUL(dirv[1], t), FMUL(dirv[2], t))),
n, ti, dirv))
nd = abs(cglib.v_dot(n, dirv))
mn = nd if mn is None else min(mn, nd)
f = nd / ONE
buckets[names[0 if f < 0.01 else 1 if f < 0.1 else 2 if f < 0.5 else 3 if f < 0.9 else 4]] += 1
return hits, {"minAbsNd": mn, "buckets": buckets, "total": len(hits)}
def shade_scan(spec, tris, nodes, order, hits):
"""对每个 EPSILON 重发阴影光线，统计自交与「丢失的真实遮挡」。"""
dl, pl = spec["directionalLight"], spec["pointLight"]
ldir, ppos = tuple(dl["dirQ16"]), tuple(pl["posQ16"])
pr2, tmax = FMUL(pl["rangeQ16"], pl["rangeQ16"]), spec["camera"]["farQ16"]
rows, occ_by_eps, src = [], {}, {}
for eps in EPS_SWEEP:
st = ss = srays = shits = rid = 0
occ = {}
for p, n, ti, _d in hits:
so = cglib.v_add(p, (FMUL(n[0], eps), FMUL(n[1], eps), FMUL(n[2], eps)))
dirs = [ldir] if cglib.v_dot(n, ldir) > 0 else []
to_l = cglib.v_sub(ppos, p)
d2 = cglib.v_dot(to_l, to_l)
if 0 < d2 < pr2 and cglib.v_dot(n, lp := cglib.v_norm(to_l)) > 0:
dirs.append(lp)
for dv in dirs:
srays += 1
sh = cgtrace.bvh_intersect(nodes, order, tris, so, dv, tmax)
occ[rid], src[rid] = (-1 if sh is None else sh[1]), ti
rid += 1
if sh is None:
continue
shits += 1
if sh[1] == ti:
st += ti < TERRAIN_TRIS
ss += ti >= TERRAIN_TRIS
occ_by_eps[eps] = occ
rows.append({"eps": eps, "selfTerrain": st, "selfSphere": ss, "selfHits": st + ss,
"shadowRays": srays, "shadowHits": shits})
for r in rows:
occ = occ_by_eps[r["eps"]]
r["lostOccluders"] = sum(1 for rid, hh in occ_by_eps[64].items()
if hh != -1 and hh != src[rid]
and occ.get(rid, -1) in (-1, src[rid]))
return rows
def main():
spec = cglib.load_spec()
tris, _ = cgtrace.build_scene_tris(spec)
nodes, order = cgtrace.build_bvh(tris)
gaps = contact_gap(spec, tris)
gap = min(g["minDistQ16"] for g in gaps)
print("== 2.3 球-地形的最小几何间隙（顶点采样）==")
for g in gaps:
print(f"球 {g['sphere']}（{g['vertices']} 顶点）: 顶点-地形最小距离 {g['minDistQ16']} "
f"Q16.16（{g['minDistQ16'] / ONE:.6f}） 中心竖直间隙 {g['verticalGapQ16'] / ONE:+.4f}")
print(f"全局最小距离 : {gap} Q16.16（{gap / ONE:.6f}）；间隙 / EPSILON = {gap / 64:.3f}\n")
hits, nd = collect_primary(spec, tris, nodes, order)
rows = shade_scan(spec, tris, nodes, order, hits)
print(f"主光线命中像素 : {len(hits)}")
print("== 2.3 EPSILON 扫描：自交（地形/球分列）与丢失的真实遮挡 ==")
print(f"{'EPSILON':>8} {'自交-地形':>9} {'自交-球':>8} {'自交合计':>8} "
f"{'阴影光线':>8} {'遮挡':>6} {'丢遮挡':>7}")
for r in rows:
print(f"{r['eps']:>8} {r['selfTerrain']:>9} {r['selfSphere']:>8} {r['selfHits']:>8} "
f"{r['shadowRays']:>8} {r['shadowHits']:>6} {r['lostOccluders']:>7}")
print("\n== 2.3 命中点 |n·d| 分布（沿光线偏移为何无效）==")
print(f"最小 |n·d| : {nd['minAbsNd']} Q16.16（{nd['minAbsNd'] / ONE:.6f}）")
for k, v in nd["buckets"].items():
print(f" |n·d| ∈ {k:<10} : {v}")
os.makedirs(OUT, exist_ok=True)
payload = {"contactGapQ16": gap, "contactGapPerSphere": gaps, "epsilonDefault": 64,
"sweep": rows, "normalDot": nd}
with open(os.path.join(OUT, "cg_trace_02_3_epsilon.json"), "w", encoding="utf-8") as f:
json.dump(payload, f, ensure_ascii=False, indent=2)
f.write("\n")
print("落盘 out/cg_trace_02_3_epsilon.json")
return 0
if __name__ == "__main__":
sys.exit(main())
跑一次（本机约 24 秒）：
# snippet: 本机实跑（工作区 lantern/）
cd lantern
/opt/homebrew/bin/python3 sim/cg_trace_02_3_epsilon.py
== 2.3 球-地形的最小几何间隙（顶点采样）==
球 0（66 顶点）: 顶点-地形最小距离 1612 Q16.16（0.024597） 中心竖直间隙 +0.2585
球 1（66 顶点）: 顶点-地形最小距离 174 Q16.16（0.002655） 中心竖直间隙 -0.7993
球 2（66 顶点）: 顶点-地形最小距离 7527 Q16.16（0.114853） 中心竖直间隙 +0.3371
全局最小距离 : 174 Q16.16（0.002655）；间隙 / EPSILON = 2.719
主光线命中像素 : 43731
== 2.3 EPSILON 扫描：自交（地形/球分列）与丢失的真实遮挡 ==
EPSILON 自交-地形 自交-球 自交合计 阴影光线 遮挡 丢遮挡
0 1226 3617 4843 77211 5533 23
1 1173 3612 4785 77211 5475 22
8 12 3371 3383 77211 4092 0
16 0 3017 3017 77211 3711 0
32 0 2094 2094 77211 2767 0
64 0 717 717 77211 1357 0
256 0 53 53 77211 672 22
1024 0 8 8 77211 629 24
4096 0 3 3 77211 624 38
16384 0 0 0 77211 628 83
32768 0 0 0 77211 655 144
65536 0 0 0 77211 759 255
== 2.3 命中点 |n·d| 分布（沿光线偏移为何无效）==
最小 |n·d| : 2 Q16.16（0.000031）
|n·d| ∈ <0.01 : 50
|n·d| ∈ 0.01-0.1 : 727
|n·d| ∈ 0.1-0.5 : 23072
|n·d| ∈ 0.5-0.9 : 19385
|n·d| ∈ >=0.9 : 497
落盘 out/cg_trace_02_3_epsilon.json
上界：最小几何间隙 174 Q16.16
球是八面体细分网格，地形是 16×16 的高度场网格。把每个球的 66 个顶点与 512 个地形三角形两两求最小距离：球 0 是 1612 Q16.16，球 1 是 174，球 2 是 7527。全局最小值 174 Q16.16 = 0.002655 世界单位，出现在球 1——它的中心竖直间隙是 −0.7993，也就是球 1 的中心沉进了地形，174 是接触环上某个顶点到地形面的距离。
EPSILON = 64 比这个间隙小 2.719 倍。这条边界是真的紧：默认值只留了不到三倍余量。
下界：地形自交在 eps ≥ 16 清零
扫描 12 个 EPSILON，把自交按地形/球分列。地形自交从 1226（eps=0）掉到 12（eps=8）、0（eps=16）——干净地在 16 处清零，所以下界是 16 Q16.16，默认的 64 留了四倍余量。64 落在 [16, 174] 之间，两侧各有约 2.7~4 倍余量。
球自交是另一回事：3617（eps=0）只掉到 717（eps=64），要到 16384 才归零。原因不在容差——球面的法线是插值出来的（三个顶点法线按重心坐标加权再归一化），在面片内部与几何法线有夹角，沿插值法线偏移清不掉这个夹角带来的自交。这是法线插值的固有误差，不是 EPSILON 能修的，如实分列就是为了不把它混进「容差」的账里。
上界实证：丢遮挡在 eps ≥ 256 出现
自交是「该清的没清」，丢遮挡是「清过头了」：在 eps = 64 时被非本源三角形挡住的阴影光线，在别的 EPSILON 下漏了光。eps = 64 时是 0；256 时跳到 22，4096 时 38，65536 时 255。撕裂的起点 256 恰好越过了最小间隙 174——几何间隙这条边界被实测印证。（eps ≤ 1 时丢遮挡也是非零的，但那是相反的原因：eps 太小时阴影光线先打到自己，根本走不到真实遮挡物。）
为什么沿法线而不是沿光线
沿光线偏移 eps·d 的法向分量是 eps·|n·d|；掠射命中处 |n·d| 很小，这个分量趋零，清不掉表面。实测命中点的 |n·d| 分布：最小值只有 2 Q16.16（0.000031），|n·d| < 0.01 的有 50 条，0.01~0.1 有 727 条。沿光线偏移在这些掠射点上几乎无效；沿法线偏移的法向分量恒为 eps，与入射角无关。这就是 EPSILON 乘的是 n 而不是 d 的理由。
跨语言核对
P 层在 TS 侧用 BigInt 重建点-三角形距离与整条 EPSILON 扫描，逐行比对：
// snippet: 节选自 lantern/core/test/cg_trace_02_3_epsilon.test.ts
/** 原始内积：Q16.16 × Q16.16 -> Q32.32。不能用 vec.ts 的 vDot，后者走 fmul 会再 >>16。 */
const rdot = (u: Vec3, v: Vec3): bigint => u[0] * v[0] + u[1] * v[1] + u[2] * v[2];
// 起点沿**法线**偏移，而不是沿光线
const so = vAdd(p, [fmul(n[0], eps), fmul(n[1], eps), fmul(n[2], eps)]);
const sh = bvhIntersect(bvh, tris, so, dirv, tmax);
if (sh !== null && sh[1] === ti) { // 打回自己所在的三角形
if (ti < TERRAIN_TRIS) selfTerrain += 1;
else selfSphere += 1;
}
// 两条边界：地形自交在 eps>=16 清零（下界），丢失遮挡在 64 为 0、256 起出现（上界）
assert.equal(base.selfTerrain, 0, "默认 EPSILON 下地形自交应为 0");
assert.equal(rowOf(64).lostOccluders, 0, "默认 EPSILON 下不应丢失真实遮挡");
assert.equal(byEps.get(16)!.selfTerrain, 0, "eps=16 起地形自交应清零");
assert.ok(byEps.get(8)!.selfTerrain > 0, "eps=8 时地形自交应仍存在");
assert.ok(rowOf(256).lostOccluders > 0, "eps=256 起应出现丢失遮挡");
assert.ok(s.contactGapQ16 < 256, "最小间隙应小于撕裂起始的 eps");
两条测试都过：三个球的顶点-地形最小距离与全局最小值逐条相同，且断言「球 1 的中心应沉入地形」；12 行扫描的 selfTerrain / selfSphere / shadowRays / shadowHits / lostOccluders 逐行相同。扫描要重发 12 × 77211 条阴影光线，是本卷 P 层最重的一条测试（本机约 7 秒）。同一个 EPSILON、同一套 BVH 遍历，两侧各写一遍，计数仍逐行一致。
阅读导航：上一节：2.2 Möller–Trumbore 与数值稳定性 · 下一节：3.1 包围盒与 BVH 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。