4.2 镜面反射与折射
站内 /graphics-water-ocean-rendering/ 那篇做的是水面渲染工程：解析波形叠加、菲涅尔近似、屏幕空间反射，目标是把一片海画得像海，关心帧时间与观感。本节换一个尺度，把「水面为什么在某些角度变成镜子、在某些角度又透明见底」拆到方向公式一级，并在 Q16.16 定点下逐位算清楚。上一节的可见性是离散的（挡 / 不挡）；这一节的三个量——反射方向、折射方向、菲涅尔权重——都是连续的，误差会累积，所以每一个都要给出可断言的整数判据。
镜面反射：r = d - 2(d·n)n
镜子把入射方向 d 关于法线 n 做一次镜像，公式只有一行：r = d - 2(d·n)n。d 指向表面（与 n 反向），n 是单位法线。因为 d·n < 0，减掉两倍的法向分量正好把光线「翻」到另一侧。在定点里，2(d·n) 是整数加倍，没有精度损失；两个减法各只有一次 fmul。
判据是单位性：d 与 n 都是单位向量时，r 也必须是单位向量。S 层用 SplitMix64 造4096 组随机单位向量，两两配对算 r，量 | |r|² - 1 |：
- 最大偏差 55（Q16.16 的一个 LSB 是 1/65536）
- 平均偏差 4.6565
偏差不为零，全部来自 fmul 的向下取整——d[i] - fmul(k, n[i]) 里那个乘法被截断了一次。4096 组里最大 55 说明误差不会失控，但它确实存在，而且逐位可复现：把同样的随机序列喂给 TypeScript 的 BigInt 版本，得到的是同一个 32 位打包摘要。
折射与 Snell 定律
折射是光线穿过界面时的方向偏折，由 Snell 定律给出：sinθi / sinθt = eta = n1 / n2。写成向量形式（eta = n1/n2，cosθi = -d·n）就是 k = 1 - eta²(1 - cos²θi) 与 t = eta·d + (eta·cosθi - sqrt(k))·n。
k 是判别式：k < 0 时无实数解，光线根本折射不出去——这就是全内反射（TIR）。所以 refract() 的返回类型是「方向或 None」，调用方必须处理 None，否则就是漏了一条分支。
定点实现里有两个必须走对的函数：
fsqrt(k)求sqrt(k)，实现是isqrt(k << 16)——把k左移 16 位再开方，等价于对 Q16.16 的量开方。eta = fdiv(n1·ONE, n2·ONE)：fdiv向零截断，而 Python 的//向下取整。凡是与 TypeScript 对齐的整数除法都必须走fdiv，否则负数会分叉。
同样的 4096 组随机方向，这次先把 n 翻到与 d 反向（保证正面入射），再算 t：
- 成功折射 1035 条，全内反射 3061 条（eta = 1.5）
- 成功那部分的单位性偏差：最大 46、平均 5.4444
全内反射占了 75%——因为随机方向各向均匀，绝大多数入射角都超过了玻璃的临界角。这恰好说明 k < 0 这条分支不是罕见情况，而是主路径。
全内反射与临界角
令 k = 0 解出临界角：cosθc = sqrt(eta² - 1) / eta。超过这个角度，光线全部被反射回介质内部。S 层用两种互相独立的方法算它：一是公式 fdiv(fsqrt(eta² - ONE), eta)；二是整数扫描，从 cosθ = 0 往上逐格增加，找 k 首次变成非负的那一格。
| 材料 | eta (Q16.16) | cosθc 公式 | cosθc 扫描 | 差 | θc |
|---|---|---|---|---|---|
| 水 | 87381 | 43347 | 43348 | 1 | 48.5915° |
| 玻璃 | 98304 | 48847 | 48848 | 1 | 41.8112° |
| 金刚石 | 158597 | 59678 | 59680 | 2 | 24.4097° |
两种方法在 1–2 个 LSB 内一致（扫描落在整数格点上，天然带一个格点的量化误差），θc 与教科书值吻合：水约 48.6°、玻璃约 41.8°、金刚石约 24.4°。金刚石差 2，是因为它的 eta 更大，fsqrt 的截断误差在 eta² 上被放大了一档。
顺带一个工程结论：折射率越大，临界角越小——金刚石只要偏离法线 24.4° 以上就全反射，所以它「火彩」强。这不是光学奇观，是 cosθc = sqrt(eta²-1)/eta 的单调性。
菲涅尔项：那个「不逐位相同」的结果
真实界面既反射又折射，能量怎么分由菲涅尔项决定。本节用最简的 Schlick 近似 (1 - V·H)^5，其中五次幂有两种写法：
- 逐次相乘：
x*x -> *x -> *x -> *x，4 次fmul。 - 平方链：
x*x -> (x²)² -> *x，3 次fmul。
教科书上它们是同一个多项式，实数下当然相等。在 Q16.16 下不相等：4096 个采样里有 429 个差 1 个 LSB（最大差就是 1）。
归因可以坐实：fmul 的交换律零反例（4096 组三元组里 0 组 fmul(a,b) != fmul(b,a)），而结合律有 979 组反例。根因是 fmul(a, b) = (a*b) >> 16 的向下取整：同一个乘积在不同分组点被截断，少截断一次就落回另一侧。这不是 bug，是定点运算的固有性质——但它提醒我们：在定点里，「代数上等价」不等于「逐位相同」。
本节因此把两件事分开声明：
- 平方链与逐次相乘不是逐位相同（实测 429/4096）。
- 平方链本身在 Python 与 TypeScript 之间是逐位相同的（同一个 32 位打包摘要）。
前一句是数学事实，后一句才是跨语言判据。
完整脚本
#!/usr/bin/env python3
# file: lantern/sim/cg_trace_04_2_reflection_refraction.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_trace_04_2_reflection_refraction.py
"""4.2《镜面反射与折射》：反射向量、Snell 定律、全内反射与菲涅尔项。
三个方向公式在 Q16.16 下逐位算清楚，并给出可断言的量：
1. **镜面反射** `r = d - 2(d·n)n`：对 4096 组随机单位向量验证 `|r|² = 1`。
2. **折射**（Snell）：`eta = n1/n2`，`k = 1 - eta²(1-cos²θi)`；`k < 0` 即全内反射。
用 `fsqrt` 求 `sqrt(k)`、用向零截断的 `fdiv` 求 `eta`，统计全内反射条数。
3. **临界角** `cosθc = sqrt(eta²-1)/eta`：公式与整数扫描（从 cosθ=0 往上找 `k` 首次
变号）互证，给出水 / 玻璃 / 金刚石三档。
4. **菲涅尔** `(1-V·H)^5`：平方链（3 次乘法）与逐次相乘（4 次乘法）各算一遍。实测
**不逐位相同**——4096 个样本里 429 个差 1 个 LSB。归因用两个反例计数器坐实：`fmul`
的**交换律**零反例（整数乘法可交换、位移相同），**结合律** 4096 组里有 979 组反例
（`(a*b)>>16` 在不同分组点向下取整，少截断一次就落回另一侧）。
5. **一次镜面弹射**：反射率 0.5，主光线命中后加一条反射光线，与基线对照。
产出：
out/cg_trace_04_2_optics.json 计数、摘要与探针（供正文引用与 P 层对拍）
out/cg_trace_04_2_snell.csv 三档折射率的临界角表
out/cg_trace_04_2_mirror.png 一次弹射的镜面成片
"""
import json
import math
import os
import sys
import numpy as np
HERE = os.path.dirname(os.path.abspath(__file__))
ROOT = os.path.dirname(HERE)
sys.path.insert(0, HERE)
import cglib # noqa: E402
import cgtrace # noqa: E402
import cg_trace_04_2_mirror as mirror # noqa: E402
import cg_trace_04_optics as optics # noqa: E402
OUT = os.path.join(ROOT, "out")
ONE = cglib.ONE
fdiv, fsqrt, fmul = cglib.fdiv, cglib.fsqrt, cglib.fmul
v_dot = cglib.v_dot
SAMPLES = 4096
SEED = 20261015
def probe_snell():
"""三档折射率的临界角：公式 vs 整数扫描。"""
rows = []
for name, num, den in (("water", 4, 3), ("glass", 3, 2), ("diamond", 242, 100)):
eta = fdiv(num * ONE, den * ONE)
eta2 = fmul(eta, eta)
cosc_formula = fdiv(fsqrt(eta2 - ONE), eta)
cosc_sweep = ONE
for c in range(0, ONE + 1): # 从 cosθ=0 往上，找 k 首次 >= 0
if ONE - fmul(eta2, ONE - fmul(c, c)) >= 0:
cosc_sweep = c
break
rows.append({
"material": name, "eta_q16": eta,
"cos_theta_c_formula": cosc_formula, "cos_theta_c_sweep": cosc_sweep,
"agree_within_1lsb": abs(cosc_formula - cosc_sweep) <= 1,
"theta_c_deg": round(math.degrees(math.acos(cosc_formula / ONE)), 4),
})
return rows
def main():
os.makedirs(OUT, exist_ok=True)
# 1) 反射：单位性 + 逐位摘要
state, errs, vals, probe_r = SEED, [], [], []
for i in range(SAMPLES):
state, d = optics.rand_unit(state)
state, n = optics.rand_unit(state)
r = optics.reflect(d, n)
errs.append(optics.unit_err(r))
vals.extend(r)
if i < 4:
probe_r.append(list(d) + list(n) + list(r))
# 2) 折射：单位性 + 全内反射统计
eta_glass = fdiv(3 * ONE, 2 * ONE)
state, rerrs, rvals, probe_t = SEED, [], [], []
ok = tir = 0
for i in range(SAMPLES):
state, d = optics.rand_unit(state)
state, n = optics.rand_unit(state)
if v_dot(d, n) > 0: # 保证正面入射
n = (-n[0], -n[1], -n[2])
t = optics.refract(d, n, eta_glass)
if t is None:
tir += 1
rvals.extend([0, 0, 0])
else:
ok += 1
rerrs.append(optics.unit_err(t))
rvals.extend(t)
if i < 4:
probe_t.append(list(d) + list(n) + (list(t) if t else [None] * 3))
# 3) 临界角
snell = probe_snell()
# 4) 菲涅尔 (1-V·H)^5：平方链 vs 逐次相乘
state, diff, max_diff, pow_vals, probe_p = SEED, 0, 0, [], []
for i in range(SAMPLES):
state, x = cgtrace.rng_q16(state)
a, b = optics.pow5_sequential(x), optics.pow5_chain(x)
if abs(a - b):
diff += 1
max_diff = max(max_diff, abs(a - b))
pow_vals.append(b)
if i < 4:
probe_p.append([x, a, b])
# 4b) 归因：fmul 破坏的是结合律，不是交换律
state, comm_bad, assoc_bad = SEED, 0, 0
for _ in range(SAMPLES):
state, a = cgtrace.rng_q16(state)
state, b = cgtrace.rng_q16(state)
state, c = cgtrace.rng_q16(state)
comm_bad += fmul(a, b) != fmul(b, a)
assoc_bad += fmul(fmul(a, b), c) != fmul(a, fmul(b, c))
# 5) 一次镜面弹射
refl_q = ONE >> 1
spec, base, _, _, _ = mirror.render_mirror(0, 0)
_, mir, _, mir_s, mir_ms = mirror.render_mirror(refl_q, 1)
cglib.save_png(os.path.join(OUT, "cg_trace_04_2_mirror.png"), mir)
with open(os.path.join(OUT, "cg_trace_04_2_mirror.rgb"), "wb") as f:
f.write(np.ascontiguousarray(mir, dtype=np.uint8).tobytes())
cglib.write_csv(
os.path.join(OUT, "cg_trace_04_2_snell.csv"),
["material", "eta_q16", "cos_theta_c_formula", "cos_theta_c_sweep", "theta_c_deg"],
[[r["material"], r["eta_q16"], r["cos_theta_c_formula"], r["cos_theta_c_sweep"],
r["theta_c_deg"]] for r in snell])
report = {
"samples": SAMPLES,
"reflect_max_unit_err": max(errs),
"reflect_mean_unit_err": round(sum(errs) / len(errs), 4),
"reflect_digest": optics.pack32(vals),
"refract_eta_glass_q16": eta_glass,
"refract_ok": ok, "refract_tir": tir,
"refract_max_unit_err": max(rerrs),
"refract_mean_unit_err": round(sum(rerrs) / len(rerrs), 4),
"refract_digest": optics.pack32(rvals),
"snell": snell,
"pow5_chain_mults": 3, "pow5_sequential_mults": 4,
"pow5_diff_count": diff, "pow5_max_abs_diff": max_diff,
"pow5_chain_digest": optics.pack32(pow_vals),
"fmul_commutativity_diff": comm_bad,
"fmul_associativity_diff": assoc_bad,
"probes": {"reflect": probe_r, "refract": probe_t, "pow5": probe_p},
"reflectivity_q16": refl_q,
"baseline_sha256": cgtrace.digest(base),
"mirror_sha256": cgtrace.digest(mir),
"mae_baseline_vs_mirror": round(cgtrace.mae(base, mir), 4),
"psnr_baseline_vs_mirror": round(cgtrace.psnr(base, mir), 4),
"baseline_mean_channel": round(float(np.asarray(base, dtype=np.int64).mean()), 4),
"mirror_mean_channel": round(float(np.asarray(mir, dtype=np.int64).mean()), 4),
"mirror_reflect_rays": mir_s["reflect_rays"],
"mirror_shadow_rays": mir_s["dir_rays"] + mir_s["pt_rays"],
"ms_mirror": mir_ms,
}
with open(os.path.join(OUT, "cg_trace_04_2_optics.json"), "w", encoding="utf-8") as f:
json.dump(report, f, ensure_ascii=False, indent=2)
f.write("\n")
print("== 4.2 镜面反射与折射 ==")
print(f"反射单位性偏差 : max {max(errs)} / mean {report['reflect_mean_unit_err']}"
f"（{SAMPLES} 组随机单位向量）")
print(f"折射单位性偏差 : max {max(rerrs)} / mean {report['refract_mean_unit_err']}"
f"（成功 {ok} / 全内反射 {tir}）")
for r in snell:
print(f" {r['material']:8s} eta={r['eta_q16']} cosθc={r['cos_theta_c_formula']}"
f"（扫描 {r['cos_theta_c_sweep']}，差"
f" {abs(r['cos_theta_c_formula'] - r['cos_theta_c_sweep'])}）θc={r['theta_c_deg']}°")
print(f"菲涅尔 x^5 : 平方链 3 次乘法 / 逐次 4 次，差 {diff} / {SAMPLES}"
f"（最大 {max_diff} LSB）")
print(f"fmul 归因 : 交换律反例 {comm_bad} / 结合律反例 {assoc_bad}"
f"（{SAMPLES} 组三元组）")
print(f"镜面弹射 : 反射光线 {mir_s['reflect_rays']}，"
f"阴影光线 {report['mirror_shadow_rays']}")
print(f"基线 vs 镜面 MAE : {report['mae_baseline_vs_mirror']}"
f"（PSNR {report['psnr_baseline_vs_mirror']} dB）")
print(f"平均通道 : 基线 {report['baseline_mean_channel']}"
f" -> 镜面 {report['mirror_mean_channel']}")
return 0
if __name__ == "__main__":
sys.exit(main())
# snippet: 本机实跑（工作区 lantern/）
cd lantern
/opt/homebrew/bin/python3 sim/cg_trace_04_2_reflection_refraction.py
== 4.2 镜面反射与折射 ==
反射单位性偏差 : max 55 / mean 4.6565（4096 组随机单位向量）
折射单位性偏差 : max 46 / mean 5.4444（成功 1035 / 全内反射 3061）
water eta=87381 cosθc=43347（扫描 43348，差 1）θc=48.5915°
glass eta=98304 cosθc=48847（扫描 48848，差 1）θc=41.8112°
diamond eta=158597 cosθc=59678（扫描 59680，差 2）θc=24.4097°
菲涅尔 x^5 : 平方链 3 次乘法 / 逐次 4 次，差 429 / 4096（最大 1 LSB）
fmul 归因 : 交换律反例 0 / 结合律反例 979（4096 组三元组）
镜面弹射 : 反射光线 43731，阴影光线 82515
基线 vs 镜面 MAE : 20.8846（PSNR 21.7343 dB）
平均通道 : 基线 59.1939 -> 镜面 39.3528
一次镜面弹射的成片
公式正确之后，把它接回渲染链路：主光线命中后，沿反射方向再追一条光线，按反射率 0.5 与直接光照混合。
- 反射光线 43731——恰好等于命中像素数，因为每个命中点发且只发一条。
- 阴影光线 82515，比 4.1 的 77211 条多 5304 条——多出来的那部分来自反射光线自己的命中点：它们在
depth + 1处也各求了一次直接光照。 - 与基线的 MAE 20.8846、PSNR 21.7343 dB。
- 平均通道从 59.1939 降到 39.3528——反射把一半能量换成了「看到的那个东西」，而不是「本来的直接光照」。
逐字节判据：refl = 0 时得到的 sha256 与卷三基线完全相同（213b4130…），refl = 0.5 时是 202f9b18…。前者证明这条链路没有意外改动，后者证明弹射真的发生了。
P 层
# snippet: core/test/cg_trace_04_2_reflection_refraction.test.ts（节选：菲涅尔的诚实声明）
# verify: node --test 'core/test/cg_trace_04_2_reflection_refraction.test.ts'
assert.equal(packDigest(chainVals), golden.pow5_chain_digest);
assert.equal(diff, golden.pow5_diff_count);
assert.ok(diff > 0, "平方链与逐次相乘确实不是逐位相同");
最后一行断言是反向的：它要求两者不相等。如果哪天有人「优化」成让它们相等了（比如改成向上取整），这条断言会失败——因为那意味着定点语义变了。P 层在这里守的不是「正确」，而是「口径不变」。
小结
- 反射
r = d - 2(d·n)n：4096 组随机方向，单位性最大偏差 55、平均 4.6565。 - 折射
k = 1 - eta²(1 - cos²θi)，k < 0即全内反射；eta = 1.5 时 4096 组里 3061 组全反射。 - 临界角公式与整数扫描在 1–2 LSB 内互证：水 48.5915°、玻璃 41.8112°、金刚石 24.4097°。
- 菲涅尔
x^5的平方链与逐次相乘不逐位相同（429/4096 差 1 LSB）；归因是fmul的结合律失效（4096 组三元组 979 组反例），交换律则零反例。 - 一次镜面弹射：反射光线 43731、阴影光线 82515，与基线 MAE 20.8846。
一次弹射只用到了主光线的方向。要让反射光线自己再反射，就得把每一层的入射方向带下去——那是下一节的事。
阅读导航：上一节：4.1 硬阴影与阴影光线 · 下一节：4.3 递归的光线树 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。