2.1 球、平面与三角形
站内 /graphics-raytracing-hybrid/ 那篇讲的是实时混合光追：光追只承担阴影与一次反射，主通道仍是光栅化，几何基本只有三角形。本卷做离线光追的完整链路，本章给的是求交的容差边界分析——先把三种解析原语摆清楚，看清「无交」与「数值退化」是怎么被原语如实报出来的，后面几章才有共同的底座。
三条闭式解，一个返回约定
三种原语都是闭式解，代价与场景规模无关：
| 原语 | 方程 | 解 |
|---|---|---|
| 光线-球 | t² + 2bt + c = 0，b = d·(o−c)，c = (o−c)·(o−c) − r² | t = −b ∓ √(b² − c)，判别式 b² − c |
| 光线-平面 | n·(o + t·d) + w = 0 | t = −(n·o + w) / (n·d) |
| 光线-三角形 | Möller–Trumbore：e1 = v1 − v0，e2 = v2 − v0，p = d × e2 | det = e1·p，t = (e2·q)/det，重心 (u, v) 由 p、q 给出 |
三条式子都写着「除以某个量」。除数趋近零时会发生什么，就是本章的主线：球那里是判别式与 √，平面与三角形那里是 det。三条原语共用同一个返回约定——命中返回 t（三角形再带上重心 (u, v)），不命中返回 None。没有「近似命中」这个返回值，2.2 会讲为什么。
三个退化情形：原语一律返回 None
底座把三种「看起来像命中、其实是数值退化」的情形全部折叠进 None：
- 方向与平面平行：
n·d = 0，t的分母为零。几何上这条光线永远碰不到平面。 - 交点在光线背后：
t ≤ 0。光线的参数域是t ∈ (0, tmax]，负根落在反方向的延长线上。 - 方向落在三角形平面内：
det = 0。此时t与重心坐标都是0/0，解不存在。
原语不替调用方做决定：它只回答「这条光线按定义有没有交点」。容差、剔除、偏置全部留给调用方——第 4 章的阴影光线与 2.3 的 EPSILON 都在调用方一侧。
完整脚本
#!/usr/bin/env python3
# file: lantern/sim/cg_trace_02_1_primitives.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_trace_02_1_primitives.py
"""2.1 的三种求交原语：光线-球（二次方程）、光线-平面、光线-三角形（Möller–Trumbore）。
脚本分两半。前半把三种解析解摆在**手算可核对**的规范例子上跑一遍，包括三个退化
情形：方向与平面平行（`denom == 0`）、交点在光线背后（`t <= 0`）、方向与三角形
平面共面（`det == 0`）——原语对这三者一律返回 `None`，把「无交」和「数值退化」
留给调用方判断。后半拿真实的 896 个三角形，把**解析球**与**细分球网格**的命中
距离摆在一起：两者的差就是八面体细分带来的面化误差（sagitta），它是第 14 章
「两条路径对照」里几何这一项的来源。
全部是 Q16.16 定点整数，两侧（Python / TS）逐位一致，故正文数字可逐条复现。
"""
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 = cglib.ONE
# 场景三角形布局：地形 512 + 三个球各 128（八面体 subdiv=2 => 8 * 4^2）
TERRAIN_TRIS = 512
SPHERE_TRIS = 128
def q16_to_str(v):
"""Q16.16 -> 人类可读的定点串（整数部分 + 4 位小数）。"""
sign = "-" if v < 0 else ""
a = abs(v)
return f"{sign}{a >> 16}.{((a & 0xFFFF) * 10000) >> 16:04d}"
def canonical_cases():
"""三种原语各跑「命中 + 退化」两类规范例子，返回值供 JSON 落盘。"""
rows = []
# --- 光线-球：单位球在原点，光线从 -3 沿 +z，命中点 z = -1，t = 2.0
o = (0, 0, -3 * ONE)
d = (0, 0, ONE)
t = cgtrace.intersect_sphere(o, d, (0, 0, 0), ONE)
rows.append(("sphere/hit", t, "t = 2.0（球面 z = -1）"))
# 光线从侧面 2.0 处平行掠过 -> 无交
t = cgtrace.intersect_sphere((2 * ONE, 0, -3 * ONE), d, (0, 0, 0), ONE)
rows.append(("sphere/miss", t, "垂直距离 2.0 > 半径 1.0 -> None"))
# 光线背向球 -> 两个根都 <= 0 -> None
t = cgtrace.intersect_sphere(o, (0, 0, -ONE), (0, 0, 0), ONE)
rows.append(("sphere/behind", t, "两根皆 <= 0 -> None"))
# --- 光线-平面：y = 0 平面，法线朝上，光线自 (0,1,-1) 垂直下落，t = 1.0
n = (0, ONE, 0)
t = cgtrace.intersect_plane((0, ONE, -ONE), (0, -ONE, 0), n, 0)
rows.append(("plane/hit", t, "t = 1.0（落到 y = 0）"))
# 方向与平面平行 -> denom == 0 -> None
t = cgtrace.intersect_plane((0, ONE, -ONE), (ONE, 0, 0), n, 0)
rows.append(("plane/parallel", t, "d·n = 0 -> None"))
# 平面在光线背后 -> t <= 0 -> None
t = cgtrace.intersect_plane((0, ONE, -ONE), (0, ONE, 0), n, 0)
rows.append(("plane/behind", t, "t <= 0 -> None"))
# --- 光线-三角形：z = 0 平面上的直角三角形，命中 (0.25, 0.25, 0)
v0, v1, v2 = (0, 0, 0), (ONE, 0, 0), (0, ONE, 0)
r = cgtrace.intersect_triangle((ONE // 4, ONE // 4, -ONE), (0, 0, ONE), v0, v1, v2)
rows.append(("triangle/hit", r, "t = 1.0，重心 (u, v) = (0.25, 0.25)"))
# 方向落在三角形平面内 -> det == 0 -> None
r = cgtrace.intersect_triangle((ONE // 4, ONE // 4, -ONE), (ONE, 0, 0), v0, v1, v2)
rows.append(("triangle/coplanar", r, "det = 0 -> None"))
# 命中点落在三角形之外（u + v > 1）-> None
r = cgtrace.intersect_triangle((3 * ONE // 4, 3 * ONE // 4, -ONE), (0, 0, ONE), v0, v1, v2)
rows.append(("triangle/outside", r, "u + v > 1 -> None"))
return rows
def sphere_vs_mesh():
"""解析球 vs 细分球网格：逐光线比较命中距离 t。"""
spec = cglib.load_spec()
tris, _ = cgtrace.build_scene_tris(spec)
cam = spec["camera"]
eye = tuple(cam["eyeQ16"])
view = cglib.mat_look_at(eye, tuple(cam["targetQ16"]), tuple(cam["upQ16"]))
right = (view[0], view[1], view[2])
up = (view[4], view[5], view[6])
back = (view[8], view[9], view[10])
sp = spec["spheres"][0]
center = tuple(sp["centerQ16"])
radius = sp["radiusQ16"]
lo = TERRAIN_TRIS
hi = TERRAIN_TRIS + SPHERE_TRIS
width, height = spec["width"], spec["height"]
aspect = cglib.fdiv(cglib.f_from_int(width), cglib.f_from_int(height))
tan_half = cam["tanHalfFovYQ16"]
tan_x = cglib.fmul(tan_half, aspect)
half_w, half_h = width // 2, height // 2
n_rays = 0
n_hit = 0
sum_abs = 0
max_abs = 0
for py in range(height):
sy = cglib.fdiv(cglib.f_from_int(2 * (half_h - py) - 1), cglib.f_from_int(height))
ay = cglib.fmul(sy, tan_half)
for px in range(width):
sx = cglib.fdiv(cglib.f_from_int(2 * (px - half_w) + 1), cglib.f_from_int(width))
ax = cglib.fmul(sx, tan_x)
dirv = cglib.v_norm((
cglib.fmul(ax, right[0]) + cglib.fmul(ay, up[0]) - back[0],
cglib.fmul(ax, right[1]) + cglib.fmul(ay, up[1]) - back[1],
cglib.fmul(ax, right[2]) + cglib.fmul(ay, up[2]) - back[2],
))
ta = cgtrace.intersect_sphere(eye, dirv, center, radius)
if ta is None:
continue
n_rays += 1
# 只在球 1 的三角形里找最近的网格命中（地形与其它球不参与）
tm = None
for i in range(lo, hi):
t = tris[i]
r = cgtrace.intersect_triangle(
eye, dirv, (t[0], t[1], t[2]), (t[3], t[4], t[5]), (t[6], t[7], t[8])
)
if r is not None and (tm is None or r[0] < tm):
tm = r[0]
if tm is None:
continue
n_hit += 1
diff = abs(ta - tm)
sum_abs += diff
if diff > max_abs:
max_abs = diff
mean_abs = sum_abs // n_hit if n_hit else 0
return {
"sphereRays": n_rays,
"sphereHits": n_hit,
"meanAbsQ16": mean_abs,
"maxAbsQ16": max_abs,
"radiusQ16": radius,
"meanPermille": round(mean_abs * 1000.0 / radius, 4),
"maxPermille": round(max_abs * 1000.0 / radius, 4),
}
def main():
spec = cglib.load_spec()
cases = canonical_cases()
print("== 2.1 三种求交原语的规范例子 ==")
for name, val, note in cases:
if val is None:
shown = "None"
elif isinstance(val, tuple):
shown = "t=" + q16_to_str(val[0]) + f" u,v=({q16_to_str(val[1])},{q16_to_str(val[2])})"
else:
shown = q16_to_str(val)
print(f"{name:<20} {shown:<42} {note}")
print()
stats = sphere_vs_mesh()
print("== 解析球 vs 细分球网格（球 1，半径 2.0）==")
print(f"解析命中的光线 : {stats['sphereRays']}")
print(f"网格也命中的光线 : {stats['sphereHits']}")
print(f"命中距离平均绝对差 : {stats['meanAbsQ16']} Q16.16（{q16_to_str(stats['meanAbsQ16'])}，"
f"{stats['meanPermille']}‰ 半径）")
print(f"命中距离最大绝对差 : {stats['maxAbsQ16']} Q16.16（{q16_to_str(stats['maxAbsQ16'])}，"
f"{stats['maxPermille']}‰ 半径）")
print()
os.makedirs(OUT, exist_ok=True)
payload = {
"cases": [
{"name": n, "hit": (v is not None), "t": (v[0] if isinstance(v, tuple) else v) if v is not None else None,
"note": note}
for n, v, note in cases
],
"sphereVsMesh": stats,
}
with open(os.path.join(OUT, "cg_trace_02_1_primitives.json"), "w", encoding="utf-8") as f:
json.dump(payload, f, ensure_ascii=False, indent=2)
f.write("\n")
print("落盘 out/cg_trace_02_1_primitives.json")
return 0
if __name__ == "__main__":
sys.exit(main())
跑一次（本机约 1 秒）：
# snippet: 本机实跑（工作区 lantern/）
cd lantern
/opt/homebrew/bin/python3 sim/cg_trace_02_1_primitives.py
== 2.1 三种求交原语的规范例子 ==
sphere/hit 2.0000 t = 2.0（球面 z = -1）
sphere/miss None 垂直距离 2.0 > 半径 1.0 -> None
sphere/behind None 两根皆 <= 0 -> None
plane/hit 1.0000 t = 1.0（落到 y = 0）
plane/parallel None d·n = 0 -> None
plane/behind None t <= 0 -> None
triangle/hit t=1.0000 u,v=(0.2500,0.2500) t = 1.0，重心 (u, v) = (0.25, 0.25)
triangle/coplanar None det = 0 -> None
triangle/outside None u + v > 1 -> None
== 解析球 vs 细分球网格（球 1，半径 2.0）==
解析命中的光线 : 975
网格也命中的光线 : 920
命中距离平均绝对差 : 7611 Q16.16（0.1161，58.0673‰ 半径）
命中距离最大绝对差 : 121134 Q16.16（1.8483，924.1791‰ 半径）
落盘 out/cg_trace_02_1_primitives.json
九个例子里三个是退化情形，全部返回 None，且三个 None 的理由各不相同——这就是本节要钉住的区分力。如果原语把「平行」与「背后」都悄悄当成「命中在很远处」，第 4 章的阴影光线会立刻把整幅画面染黑。
- 光线-球：从
z = −3沿+z打原点单位球，球面在z = −1，t = 2.0；垂直距离2.0 > 半径 1.0的擦过、以及背向的两个根都≤ 0，都是None。 - 光线-平面：
y = 0平面，光线自(0, 1, −1)垂直下落，t = 1.0；d·n = 0与平面在光线背后各是None。 - 光线-三角形：
z = 0上的直角三角形，命中t = 1.0、重心(u, v) = (0.25, 0.25)；方向共面（det = 0）与u + v > 1落在三角形外各是None。
解析球 vs 细分球网格：面化误差
场景里的三个球在底座里是八面体细分网格（subdiv = 2，每球 128 个三角形），求交时也可以换成解析球。两者的差就是面化误差：网格是内接多面体，永远不比解析球更靠近相机，所以解析命中距离恒 ≥ 网格命中距离，差恒非负（实测 920 条全部满足）。
在 320×240 的主光线上量（球 1，半径 2.0 = 131072 Q16.16）：解析球命中 975 条光线，其中 920 条网格也命中；命中距离平均绝对差 7611 Q16.16（0.1161 世界单位，58.0673‰ 半径），最大绝对差 121134 Q16.16（1.8483 世界单位，924.1791‰ 半径）。
平均差约千分之五十八，最大差接近一整个半径——「用网格近似球」在离线光追里不是可以忽略的近似。最大值出现在轮廓附近的掠射光线上：那里球面几乎与光线平行，法向几十毫弧度的偏差会被放大成沿光线的巨大距离差（那条光线解析命中 t = 1563476、网格命中 t = 1684610，Δ = 121134）。命中集合也不同（975 vs 920）：轮廓附近有一批光线解析球打得到、网格打不到，差的是八面体在轮廓处的凹陷。
跨语言核对
P 层在 TS 侧用 BigInt 重建三条原语与整条网格比对循环，逐位复核：
// snippet: 节选自 lantern/core/test/cg_trace_02_1_primitives.test.ts
const got: Record<string, bigint | [bigint, bigint, bigint] | null> = {
"sphere/hit": intersectSphere([0n, 0n, -3n * ONE], [0n, 0n, ONE], Z, ONE),
"plane/hit": intersectPlane([0n, ONE, -ONE], [0n, -ONE, 0n], [0n, ONE, 0n], 0n),
"triangle/hit": intersectTriangle(
[ONE / 4n, ONE / 4n, -ONE], [0n, 0n, ONE], Z, [ONE, 0n, 0n], [0n, ONE, 0n],
),
"triangle/coplanar": intersectTriangle( // 方向共面 -> det == 0
[ONE / 4n, ONE / 4n, -ONE], [ONE, 0n, 0n], Z, [ONE, 0n, 0n], [0n, ONE, 0n],
),
};
assert.deepEqual(got["sphere/hit"], 2n * ONE); // t = 2.0
assert.equal(got["plane/hit"], ONE); // t = 1.0
assert.deepEqual(got["triangle/hit"], [ONE, ONE / 4n, ONE / 4n]); // (u, v) = (0.25, 0.25)
assert.equal(got["triangle/coplanar"], null); // det == 0 -> None
assert.equal(got["plane/parallel"], null); // d·n == 0 -> None
assert.equal(got["sphere/behind"], null); // 两根皆 <= 0 -> None
assert.ok(maxAbs < radius, "面化误差应小于半径");
两侧的 t、u、v 逐位相同，sphereRays / sphereHits / meanAbsQ16 / maxAbsQ16 四个数也逐条相同，并且断言「网格命中距离恒不小于解析命中距离」与「最大面化误差小于半径」。同一套定点语义在 Python 的 int 与 TypeScript 的 BigInt 里给出同一个答案——这是卷三后面所有章节的前提。
阅读导航：上一节：1.3 递归与深度限制 · 下一节：2.2 Möller–Trumbore 与数值稳定性 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。