1.1 光线的表示与求交
站内 /graphics-raytracing-hybrid/ 那篇讲的是实时混合光追：光栅化主通道加上光追阴影与反射，跑在 GPU 上，关心帧时间与降噪。本卷换一个角度——做离线光追的完整链路，把每一步都落成可逐光线计数的整数。这条路线的第一个对象是光线：它怎么被表示，以及「命中一个物体」这件事在数学上到底是什么。
光线就是 o + t·d
一条光线写成 p(t) = o + t·d（t ≥ 0）。o 是起点，d 是方向，t 是沿光线的参数。这个式子里只有三个约定需要说清：
d必须归一化。只有当|d| = 1时，t才等于「从起点出发走了多远」；否则t只是一个无量纲的比例，不同像素的t之间没有可比性，后面所有关于「最近命中」和「距离衰减」的判据都会失去意义。- 主光线的
o是同一个点。卷三的 76800 条主光线全部从相机位置eye出发，逐像素不同的只有d。这跟阴影光线不同：阴影光线的起点是每个命中点，各不相同。 t的符号有意义。t < 0表示交点在相机背后，求交时必须丢掉。本卷的求交原语统一只接受t ∈ (ε, tmax]。
求交的一般形式：解一个隐式方程
「一条光线命中一个物体」在数学上只有一种写法：把 p(t) 代进物体的隐式方程，解出 t。三个原语只是三种不同的方程：
| 原语 | 隐式方程 F(p) = 0 | 代入 p = o + t·d 后 | 求解 |
|---|---|---|---|
球 |p-c|=r | |p-c|² - r² | |o + t·d - c|² - r² | 二次方程，取较小正根 |
平面 n·p + d = 0 | n·p + d | n·o + t(n·d) + d | 一次方程，t = -(n·o+d)/(n·d) |
| 三角形 | 所在平面的方程 + 重心内判 | 同上 | 先解平面 t，再判 u,v ≥ 0, u+v ≤ 1 |
三种原语解完 t 之后，都能做同一件事：把 t 代回 F(p)，看残差是不是 0。这个残差是本节的判据——它把「求交对不对」从一句断言变成了一个可以落盘的整数。球的 F 量纲是长度的平方（残差按平方放大），平面的 F 是线性量，三角形用未归一化的几何法线做 F，故残差的绝对值只能同类相比。
一条探针光线，三种原语
要在一帧里同时验证三种原语，得找一条同时命中球、地面和三角形的光线。本节用一份独立的探针配置：
{
"name": "lantern-ray-probe",
"samplePixels": [
[0, 0],
[159, 0],
[319, 0],
[0, 119],
[159, 119],
[319, 119],
[0, 239],
[159, 239],
[319, 239]
],
"groundPlane": {
"normalQ16": [0, 65536, 0],
"dQ16": 0
},
"probeSphere": {
"centerQ16": [-393216, 196608, 0],
"radiusQ16": 131072
},
"normResidualTolQ16": 64,
"pointResidualTolQ16": 8
}
samplePixels 是九个采样点：四角 + 四边中点 + 中心，覆盖整帧。probeSphere 与 groundPlane 是探针几何，不是场景几何——它们只用来验证求交公式，groundPlane 的 dQ16 = 0 表示平面过原点（y = 0 的地面）。
完整脚本
#!/usr/bin/env python3
# file: lantern/sim/cg_trace_01_1_ray.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_trace_01_1_ray.py
"""1.1 光线的表示与求交：把 `o + t·d` 与「求交 = 解一个隐式方程」摆成可跑的表。
三件事分开量：
1. 参数表示 `o + t·d`——主光线的 `o` 是同一个眼睛点，`d` 逐像素不同；
命中点由 `o + t·d` 复原，与三角形重心坐标复原的点逐分量比对。
2. 定点归一化——`v_norm` 之后 `|d|^2` 与 ONE 的残差有多大，全部 76800 条主光线都量。
3. 求交的一般形式——同一条光线对球（二次方程）、地面（一次方程）、
三角形（平面 + 重心内判）各求一次 `t`，再把 `t` 代回各自的隐式函数 `F(p)` 看残差。
产物：out/cg_trace_01_1_ray.csv、out/cg_trace_01_1_ray.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__)))
CONFIG = os.path.join(ROOT, "config", "cg_trace_01_1_ray.json")
def camera_basis(spec):
"""从 `mat_look_at` 反推相机基向量：右 / 上 / 后（世界空间，Q16.16）。
`mat_look_at` 返回的是**行主序**矩阵：第 0 行 right、第 1 行 up、第 2 行 back。
故取三**行**（`view[0:3] / view[4:6] / view[8:10]`），不是三列。
"""
cam = spec["camera"]
v = cglib.mat_look_at(tuple(cam["eyeQ16"]), tuple(cam["targetQ16"]), tuple(cam["upQ16"]))
return (v[0], v[1], v[2]), (v[4], v[5], v[6]), (v[8], v[9], v[10])
def raw_dir(px, py, width, height, right, up, back, tan_half, tan_x):
"""像素 -> **未归一化**方向 `ax·right + ay·up - back`（与底座逐字同式）。"""
sy = cglib.fdiv(cglib.f_from_int(2 * (height // 2 - py) - 1), cglib.f_from_int(height))
sx = cglib.fdiv(cglib.f_from_int(2 * (px - width // 2) + 1), cglib.f_from_int(width))
ax, ay = cglib.fmul(sx, tan_x), cglib.fmul(sy, tan_half)
return tuple(cglib.fmul(ax, right[k]) + cglib.fmul(ay, up[k]) - back[k] for k in range(3))
def point_on_ray(o, d, t):
"""`o + t·d`（Q16.16）。"""
return cglib.v_add(o, tuple(cglib.fmul(d[k], t) for k in range(3)))
def bary_point(v0, v1, v2, u, v):
"""重心坐标复原：`(1-u-v)·v0 + u·v1 + v·v2`。"""
w0 = cglib.ONE - u - v
return tuple(cglib.fmul(v0[k], w0) + cglib.fmul(v1[k], u) + cglib.fmul(v2[k], v)
for k in range(3))
def main():
spec = cglib.load_spec()
cfg = json.load(open(CONFIG, encoding="utf-8"))
w, h = spec["width"], spec["height"]
eye = tuple(spec["camera"]["eyeQ16"])
tan_half = spec["camera"]["tanHalfFovYQ16"]
tan_x = cglib.fmul(tan_half, cglib.fdiv(cglib.f_from_int(w), cglib.f_from_int(h)))
right, up, back = camera_basis(spec)
tris, _ = cgtrace.build_scene_tris(spec)
nodes, order = cgtrace.build_bvh(tris)
tmax = spec["camera"]["farQ16"]
def ray(px, py):
return cglib.v_norm(raw_dir(px, py, w, h, right, up, back, tan_half, tan_x))
# ---- 全图：归一化残差 + 命中统计 + 点复原残差 ----
norm_max = over = hits = 0
t_lo = t_hi = None
pres_max, worst = 0, None
for py in range(h):
for px in range(w):
d = ray(px, py)
res = abs(cglib.v_dot(d, d) - cglib.ONE)
norm_max = max(norm_max, res)
over += res > cfg["normResidualTolQ16"]
hit = cgtrace.bvh_intersect(nodes, order, tris, eye, d, tmax)
if hit is None:
continue
hits += 1
t, ti, u, v = hit
t_lo = t if t_lo is None else min(t_lo, t)
t_hi = t if t_hi is None else max(t_hi, t)
rec = tris[ti]
pr = point_on_ray(eye, d, t)
pb = bary_point(rec[0:3], rec[3:6], rec[6:9], u, v)
mx = max(abs(pr[k] - pb[k]) for k in range(3))
if mx > pres_max:
pres_max = mx
# 顺带记录最坏像素的 Möller–Trumbore 行列式 det：det 越小说明光线越掠射
e1 = cglib.v_sub(rec[3:6], rec[0:3])
e2 = cglib.v_sub(rec[6:9], rec[0:3])
worst = {"px": px, "py": py, "tri": ti, "tQ16": t, "uQ16": u, "vQ16": v,
"detQ16": cglib.v_dot(e1, cglib.v_cross(d, e2))}
# ---- 采样像素：逐条记录 ----
rows = []
for px, py in cfg["samplePixels"]:
u = raw_dir(px, py, w, h, right, up, back, tan_half, tan_x)
d = cglib.v_norm(u)
hit = cgtrace.bvh_intersect(nodes, order, tris, eye, d, tmax)
rows.append({"px": px, "py": py, "u": u, "lenU": cglib.v_len(u), "d": d,
"len2": cglib.v_dot(d, d), "t": hit[0] if hit else None})
# ---- 三种原语的「一般形式」：找一条同时命中球、地面、三角形的光线 ----
sp, gp = cfg["probeSphere"], cfg["groundPlane"]
sc, sr, pn = tuple(sp["centerQ16"]), sp["radiusQ16"], tuple(gp["normalQ16"])
probe = None
for py in range(h):
for px in range(w):
d = ray(px, py)
ts = cgtrace.intersect_sphere(eye, d, sc, sr)
hit = cgtrace.bvh_intersect(nodes, order, tris, eye, d, tmax)
if ts is None or hit is None:
continue
p = point_on_ray(eye, d, hit[0])
# 网格是球的内接多面体，最近命中点到球心的距离略小于半径
if abs(cglib.v_len(cglib.v_sub(p, sc)) - sr) > 3000:
continue
if cgtrace.intersect_plane(eye, d, pn, gp["dQ16"]) is None:
continue
probe = (px, py, d)
break
if probe:
break
probe_px, probe_py, pd = probe
t_sp = cgtrace.intersect_sphere(eye, pd, sc, sr)
oc = cglib.v_sub(point_on_ray(eye, pd, t_sp), sc)
f_sp = cglib.v_dot(oc, oc) - cglib.fmul(sr, sr) # F_sphere(p) = |p-c|^2 - r^2
t_pl = cgtrace.intersect_plane(eye, pd, pn, gp["dQ16"])
f_pl = cglib.v_dot(pn, point_on_ray(eye, pd, t_pl)) + gp["dQ16"] # F_plane(p) = n·p + d
th = cgtrace.bvh_intersect(nodes, order, tris, eye, pd, tmax)
rec = tris[th[1]]
geo_n = cglib.v_cross(cglib.v_sub(rec[3:6], rec[0:3]), cglib.v_sub(rec[6:9], rec[0:3]))
f_tr = cglib.v_dot(geo_n, cglib.v_sub(point_on_ray(eye, pd, th[0]), rec[0:3]))
out = os.path.join(ROOT, "out")
os.makedirs(out, exist_ok=True)
cglib.write_csv(
os.path.join(out, "cg_trace_01_1_ray.csv"),
["px", "py", "ux", "uy", "uz", "len_u_q16", "dx", "dy", "dz", "len2_q16",
"len2_residual", "t_q16"],
[[r["px"], r["py"], r["u"][0], r["u"][1], r["u"][2], r["lenU"],
r["d"][0], r["d"][1], r["d"][2], r["len2"], abs(r["len2"] - cglib.ONE),
r["t"] if r["t"] is not None else ""] for r in rows])
stats = {
"width": w, "height": h, "eyeQ16": list(eye),
"rightQ16": list(right), "upQ16": list(up), "backQ16": list(back),
"primaryRays": w * h, "normResidualTolQ16": cfg["normResidualTolQ16"],
"normResidualMaxQ16": norm_max, "normResidualOverTol": over,
"hitCount": hits, "hitPermille": round(hits * 1000.0 / (w * h), 4),
"tMinQ16": t_lo, "tMaxQ16": t_hi,
"pointResidualMaxQ16": pres_max, "pointResidualWorst": worst,
"probe": {"px": probe_px, "py": probe_py, "dirQ16": list(pd),
"sphere": {"tQ16": t_sp, "fQ16": f_sp},
"plane": {"tQ16": t_pl, "fQ16": f_pl},
"triangle": {"tQ16": th[0], "fQ16": f_tr, "tri": th[1],
"uQ16": th[2], "vQ16": th[3]}},
"sampleRows": [{"px": r["px"], "py": r["py"], "u": list(r["u"]), "lenU": r["lenU"],
"d": list(r["d"]), "len2": r["len2"], "t": r["t"]} for r in rows],
}
with open(os.path.join(out, "cg_trace_01_1_ray.json"), "w", encoding="utf-8") as fh:
json.dump(stats, fh, ensure_ascii=False, indent=2)
fh.write("\n")
print("== 1.1 光线的表示与求交 ==")
print(f"主光线总数 : {stats['primaryRays']}")
print(f"|d|^2 残差 max : {norm_max} Q16.16（超容差 {cfg['normResidualTolQ16']} 的条数 {over}）")
print(f"命中像素 : {hits}（{stats['hitPermille']}‰）")
print(f"t 范围 : [{t_lo}, {t_hi}] Q16.16")
print(f"点复原残差 max : {pres_max} Q16.16（o+t·d vs 重心复原）"
f" 最坏像素 {worst['px']},{worst['py']} det={worst['detQ16']}")
print(f"探针光线 {probe_px},{probe_py} : d = {pd}")
print(f" 球 t = {t_sp} F = {f_sp}")
print(f" 平面 t = {t_pl} F = {f_pl}")
print(f" 三角 t = {th[0]} F = {f_tr}（三角形 #{th[1]}）")
print("csv/json : out/cg_trace_01_1_ray.{csv,json}")
return 0
if __name__ == "__main__":
sys.exit(main())
跑一次：
# snippet: 本机实跑（工作区 lantern/）
cd lantern
/opt/homebrew/bin/python3 sim/cg_trace_01_1_ray.py
== 1.1 光线的表示与求交 ==
主光线总数 : 76800
|d|^2 残差 max : 5 Q16.16（超容差 64 的条数 0）
命中像素 : 43731（569.4141‰）
t 范围 : [451114, 3276183] Q16.16
点复原残差 max : 1437 Q16.16（o+t·d vs 重心复原） 最坏像素 211,112 det=1417
探针光线 106,94 : d = (-16223, -2735, -63438)
球 t = 1604219 F = 712
平面 t = 9422231 F = 0
三角 t = 1622517 F = -29（三角形 #528）
csv/json : out/cg_trace_01_1_ray.{csv,json}
五个数字怎么读
|d|² 残差最大 5 Q16.16，越界 0 条。 v_norm 之后每条主光线的方向长度平方理论上恰为 ONE = 65536，实测与 65536 的差最大是 5。这是定点开方与除法各截断一次留下的舍入量级（5/65536 ≈ 7.6e-5）。判据用 normResidualTolQ16 = 64 卡：比实测最坏值宽约 12 倍，又比「忘了归一化」的残差小得多——角点 (0,0) 的未归一化长度平方是 90796² ≈ 8.24e9，换算回 Q16.16 的残差是 60256（约 0.92 个 ONE），判据必须能拒绝它。这条实测数字顺带保证了后续一切「按 t 排序取最近命中」的判据成立。
命中 43731 个像素（569.4141‰）。 76800 条主光线里约 57% 打中了东西，其余穿出天空。这个数与卷三基线逐位相同，是相机与场景装配没有走样的第一道旁证。
t 落在 [451114, 3276183]。 命中点离相机最近约 6.9 个世界单位、最远约 50.0 个（除以 65536）。最近的那些是地形近端，最远的是球体背面——它们都远大于 EPSILON = 64（世界单位 1/1024），说明自交抑制偏移与场景尺度是匹配的。
点复原残差最大 1437，出现在像素 211,112。 这条量的是「o + t·d 复原的点」与「三角形重心坐标复原的点」差多少。绝大多数像素的残差是 0 或个位数，但最坏像素达到 1437（约 0.022 个世界单位）。原因是这个像素命中的三角形 #724 是一条掠射光线：它的 Möller–Trumbore 行列式 det = 1417，比典型值小两个数量级，定点除法在分母很小时把误差放大了。这不是 bug，是定点求交的固有边界，也是第 2 章《Möller–Trumbore 与数值稳定性》的主题——本节只负责把它记下来并留下坐标。
探针光线的三种 F(p) 残差：712、0、-29。 这是全节的核心。同一条光线 d = (-16223, -2735, -63438)：
- 球：
t = 1604219，F = |p-c|² - r² = 712。球的F量纲是长度平方，r² = 131072² ≈ 1.72e10，残差 712 相对它是 4.1e-8 量级——精度已经贴着 Q16.16 的上限。 - 平面：
t = 9422231，F = n·p + d = 0。平面的F是线性量，一次除法后残差精确为 0。这是三种原语里最干净的一条。 - 三角形：
t = 1622517，F = geo_n·(p - v0) = -29。用未归一化的几何法线，残差同样是舍入量级。
三条残差都落在 Q16.16 的舍入噪声里，说明「求交 = 解隐式方程」这个一般形式在三种原语上都成立，而且可以被一个整数判据逐条验证。
跨语言核对
P 层用 BigInt 重算同一份配置，与 S 侧 out/cg_trace_01_1_ray.json 的 golden 表逐项比对：
// snippet: 节选自 lantern/core/test/cg_trace_01_1_ray.test.ts
const d = vNorm(rawDir(row.px, row.py)); // 与 Python 的 cglib.v_norm 同序
const res = vDot(d, d) - ONE;
const absRes = res < 0n ? -res : res;
assert.ok(absRes <= 4n, `|d|^2 残差 ${absRes} > 4`); // 采样行的 |d|^2 残差必须在 ±4 内
assert.deepEqual(d.map(Number), row.d); // 方向逐分量与 S 侧一致
const ts = intersectSphere(eye, d, center, radius);
assert.equal(Number(ts), golden.probe.sphere.tQ16); // 球求交 t 逐位相同
两个负例把判据的区分力钉死：跳过 vNorm 后 |d|² 与 ONE 的残差会远大于 4（角点处 60256），判据必须拒绝；点复原残差最坏像素的 det 必须小于 4096（实测 1417），否则「掠射放大误差」这个结论就无从谈起。两者都被断言拒绝，所以「三种原语求交正确」这个结论才可信。
阅读导航：上一节：目录 · 下一节：1.2 相机与主光线 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。