9.1 体积渲染方程
站内 /graphics-volumetric-clouds-atmosphere/ 那篇做的是体积云工程：噪声形状、步进次数、光照近似，目标是在 GPU 上把一片云画得好看并压住帧时间。本节换一条路——把「光穿过一片参与介质时到底发生了什么」当成积分方程来推，再用定点整数把它算成一条条可断言的数。这一节不谈帧率，只谈方程长什么样，以及它落到整数上之后误差有多大。
光在介质里会被改道
表面渲染假设光只在几何表面与物质相互作用。参与介质（雾、烟、云、大气）推翻这个假设：光沿路径的每一个点都可能被吸收或散射。于是原本「只在表面取值」的渲染方程，变成沿路径的积分方程——这就是体积渲染方程（volume rendering equation），也叫辐射传输方程的一个稳态形式：
L(x, ω) = ∫_0^s T(t) · [ σ_a · L_e(t) + σ_s · ∫_{S²} p(ω, ω') L_i(t, ω') dω' ] dt
+ T(s) · L_surface
逐项读：
T(t) = exp(-σ_t · t)是透射率，从起点到t之间光没被挡住的概率；σ_t = σ_a + σ_s是消光系数。σ_a · L_e(t)是介质自身发射（火焰、发光雾）。σ_s · ∫ p(ω, ω') L_i(t, ω') dω'是散射进来的光：把来自所有方向ω'的入射辐射亮度，按相函数p加权重排到出射方向ω。T(s) · L_surface是穿过整片介质后没被碰过、直接来自背景表面的那部分。注意σ_a只出现在L_e一项里——吸收不产生光，它只通过T削弱别的光，这是「介质不是黑体」在方程里的落点。
三项里，先把透射率钉死
后面所有「离散化误差」都要跨过一个地板：透射率本身在定点下能算多准。exp(-x) 没有整数版本，我们把它拆成 2^(-k) · 2^(-f)，2^(-k) 用算术右移、2^(-f) 查一张由纯整数递推建出来的表。相函数用 Henyey–Greenstein 的闭式，同样全程整数。S 层脚本把这两件事各量一遍，再量单次与多次散射的取舍。
#!/usr/bin/env python3
# file: lantern/sim/cg_trace_09_1_volume_equation.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_trace_09_1_volume_equation.py
"""9.1 体积渲染方程：把「光子在介质里被吸收 / 被散射 / 自己发光」写成一项一项可断言的整数。
体积渲染方程把表面渲染方程里那条「只在表面取值」的积分，换成**沿路径连续取值**：
L(x, ω) = ∫_0^s T(t) · [ σ_a · L_e(t) + σ_s · ∫_{S²} p(ω, ω') L_i(t, ω') dω' ] dt
+ T(s) · L_surface
三项分别是**发射**、**吸收**（体现在 `T` 里）、**散射**。本脚本量三件事：
1. **透射率** `T = exp(-σ_t s)` 的定点实现——`fexp_neg`（整数 LUT）与浮点 `exp` 差多少，
这是后面一切「离散化误差」必须跨过的地板。
2. **相函数** `p(cosθ)` 的定点实现——Henyey–Greenstein，验证归一化 `∫p̃ dcosθ = 2`
与各向异性因子 `⟨cosθ⟩ = g`，并给出前向/后向峰值比。
3. **单次 vs 多次散射**——用能量守恒的泊松阶数模型，量「只算一次散射」相对
「所有阶求和」漏掉多少，随单次散射反照率 `a = σ_s/σ_t` 与光学深度 `τ` 变化。
产物：out/cg_trace_09_1_volume.json、out/cg_trace_09_1_volume.csv
"""
import json
import os
import sys
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cglib # noqa: E402
import cg_trace_09_medium as med # noqa: E402
ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
ONE = cglib.ONE
def poisson_orders(sigma_a, sigma_t, tau, light, phase, max_order):
"""能量守恒的泊松阶数模型：返回 `[L_0, L_1, ..., L_max]`（Q16.16）。
沿光学深度 `τ` 的路径，散射次数服从泊松分布 `P(k) = τ^k e^{-τ} / k!`；每次散射
保住能量比例 `a = σ_s/σ_t`。故第 `k` 阶贡献
L_k = L_in · p · a^k · τ^k · e^{-τ} / k!
所有阶求和 = `L_in · p · e^{-(1-a)τ}`——即「只有吸收才衰减」的极限。`L_1` 就是
**单次散射近似**，`Σ_k L_k` 是**多次散射的完整和**。注意这是一个**模型**（各向同性
再发射、只统计能量），不是体积渲染方程的完整解；它只用来量「阶数截断」这一项。
"""
a = cglib.fdiv(sigma_t - sigma_a, sigma_t)
exp_neg = med.fexp_neg(tau)
fact = 1
a_pow = ONE # a^k
tau_pow = ONE # τ^k
out = []
for k in range(max_order + 1):
if k > 0:
fact *= k
a_pow = cglib.fmul(a_pow, a)
tau_pow = cglib.fmul(tau_pow, tau)
if k == 0:
out.append(0) # 0 阶是未散射的直射光，不计入散射
continue
# L_k = light · phase · a^k · τ^k · exp_neg / k!
# 5 个 Q16.16 因子的原始乘积是 Q80.64，除 k! 后再 >> 64 得 Q16.16。
out.append((light * phase * a_pow * tau_pow * exp_neg) // (fact << 64))
return out
def poisson_total(sigma_a, sigma_t, tau, light, phase):
"""所有**散射**阶（k ≥ 1）之和 = `L_in · p · (e^{-(1-a)τ} − e^{-τ})`（Q16.16）。
减去 `e^{-τ}` 是为了扣掉 0 阶（未散射的直射光），使「阶数和」与 `poisson_orders`
的口径一致——两边都只统计**被散射过**的能量。
"""
a = cglib.fdiv(sigma_t - sigma_a, sigma_t)
total = med.fexp_neg(cglib.fmul(ONE - a, tau)) - med.fexp_neg(tau)
# 3 个 Q16.16 因子的乘积是 Q48.16，>> 32 得 Q16.16。
return (light * phase * total) >> 32
def main():
with open(os.path.join(ROOT, "config", "cg_trace_09_medium.json"), encoding="utf-8") as fh:
cfg = json.load(fh)
out_dir = os.path.join(ROOT, "out")
os.makedirs(out_dir, exist_ok=True)
sigma_a = cfg["sigmaAQ16"]
sigma_s = cfg["sigmaSQ16"]
sigma_t = cfg["sigmaTQ16"]
# ---- A. 定点指数 vs 浮点 exp：这是离散化误差的地板 ----
exp_rows = []
exp_max_rel = 0.0
for x in cfg["expProbeQ16"]:
q = med.fexp_neg(x)
ref = med.transmittance_float(ONE, x) # exp(-x)，x 已按 Q16.16 解释
refq = round(ref * ONE)
rel = abs(q - refq) / max(1, refq)
exp_max_rel = max(exp_max_rel, rel)
exp_rows.append((x, q, refq, q - refq, round(rel * 1e6, 4)))
# 把 x 当作光学深度 τ，再量一次「一条路径的透射率」的误差
tau_probe = cfg["tauSweepQ16"]
tau_rows = []
for tau in tau_probe:
q = med.fexp_neg(tau)
refq = round(med.transmittance_float(ONE, tau) * ONE)
tau_rows.append((tau, q, refq, q - refq))
# ---- B. HG 相函数的定点实现：归一化与各向异性 ----
g_rows = []
for g in cfg["gSweepQ16"]:
norm, asym = med.phase_quadrature(g, 4000)
fwd = med.hg_phase_times4pi(g, ONE)
bwd = med.hg_phase_times4pi(g, -ONE)
g_rows.append({
"gQ16": g, "g": round(g / ONE, 6),
"normQ16": norm, "normErrQ16": norm - 2 * ONE,
"asymQ16": asym, "asymErrQ16": asym - g,
"forwardQ16": fwd, "backwardQ16": bwd,
"fbRatio": round(fwd / bwd, 4) if bwd else None,
})
# ---- C. 单次 vs 多次散射：阶数截断的取舍 ----
light = cfg["light"]["intensityQ16"]
phase_iso = med.INV_4PI_Q16
order_rows = []
for a_q16 in cfg["albedoSweepQ16"]:
# 固定光学深度 τ = 2.0，扫反照率
tau = 2 * ONE
sa = cglib.fmul(ONE - a_q16, sigma_t) # σ_a = (1-a)·σ_t
orders = poisson_orders(sa, sigma_t, tau, light, phase_iso, cfg["maxOrder"])
total = poisson_total(sa, sigma_t, tau, light, phase_iso)
l1 = orders[1]
s_all = sum(orders)
order_rows.append({
"kind": "albedo", "aQ16": a_q16, "a": round(a_q16 / ONE, 4),
"tauQ16": tau, "tau": 2.0,
"L1Q16": l1, "sumQ16": s_all, "totalQ16": total,
"singleFracPct": round(100.0 * l1 / s_all, 4) if s_all else None,
"truncErrPct": round(100.0 * (s_all - l1) / s_all, 4) if s_all else None,
"ordersQ16": orders[1:],
})
for tau in cfg["tauSweepQ16"]:
a_q16 = cfg["albedoSweepQ16"][4] # 固定 a = 0.8，扫光学深度
sa = cglib.fmul(ONE - a_q16, sigma_t)
orders = poisson_orders(sa, sigma_t, tau, light, phase_iso, cfg["maxOrder"])
total = poisson_total(sa, sigma_t, tau, light, phase_iso)
l1 = orders[1]
s_all = sum(orders)
order_rows.append({
"kind": "tau", "aQ16": a_q16, "a": round(a_q16 / ONE, 4),
"tauQ16": tau, "tau": round(tau / ONE, 4),
"L1Q16": l1, "sumQ16": s_all, "totalQ16": total,
"singleFracPct": round(100.0 * l1 / s_all, 4) if s_all else None,
"truncErrPct": round(100.0 * (s_all - l1) / s_all, 4) if s_all else None,
"ordersQ16": orders[1:],
})
report = {
"sigmaAQ16": sigma_a, "sigmaSQ16": sigma_s, "sigmaTQ16": sigma_t,
"sigmaA": round(sigma_a / ONE, 4), "sigmaS": round(sigma_s / ONE, 4),
"sigmaT": round(sigma_t / ONE, 4),
"singleScatterAlbedo": round((sigma_t - sigma_a) / sigma_t, 4),
"expMaxRelErr": round(exp_max_rel, 8),
"expProbe": [{"xQ16": r[0], "fexpQ16": r[1], "refQ16": r[2], "diffQ16": r[3],
"relPpm": r[4]} for r in exp_rows],
"tauProbe": [{"tauQ16": r[0], "fexpQ16": r[1], "refQ16": r[2], "diffQ16": r[3]}
for r in tau_rows],
"hg": g_rows,
"orders": order_rows,
}
with open(os.path.join(out_dir, "cg_trace_09_1_volume.json"), "w", encoding="utf-8") as fh:
json.dump(report, fh, ensure_ascii=False, indent=2)
fh.write("\n")
cglib.write_csv(
os.path.join(out_dir, "cg_trace_09_1_volume.csv"),
["kind", "aQ16", "tauQ16", "L1Q16", "sumQ16", "totalQ16", "singleFracPct", "truncErrPct"],
[[r["kind"], r["aQ16"], r["tauQ16"], r["L1Q16"], r["sumQ16"], r["totalQ16"],
r["singleFracPct"], r["truncErrPct"]] for r in order_rows])
print("== 9.1 体积渲染方程 ==")
print(f"介质系数 : σ_a={report['sigmaA']} σ_s={report['sigmaS']} "
f"σ_t={report['sigmaT']}（单次散射反照率 {report['singleScatterAlbedo']}）")
print(f"定点 exp 最大相对误差 : {report['expMaxRelErr']:.3e}"
f"（{round(exp_max_rel * 1e6, 4)} ppm）")
print(f"{'x(Q16)':>9} {'fexp':>8} {'ref':>8} {'diff':>5} {'ppm':>8}")
for x, q, refq, d, ppm in exp_rows:
print(f"{x:>9} {q:>8} {refq:>8} {d:>5} {ppm:>8}")
print("HG 相函数（∫p̃ dcosθ 应 = 131072，⟨cosθ⟩ 应 = g）：")
print(f"{'g':>8} {'norm':>9} {'normErr':>8} {'asym':>9} {'asymErr':>8} "
f"{'前向 p̃':>10} {'后向 p̃':>10} {'前/后':>10}")
for r in g_rows:
print(f"{r['g']:>8.3f} {r['normQ16']:>9} {r['normErrQ16']:>8} {r['asymQ16']:>9} "
f"{r['asymErrQ16']:>8} {r['forwardQ16']:>10} {r['backwardQ16']:>10} "
f"{r['fbRatio']:>10}")
print("单次 vs 多次散射（泊松阶数模型，k ≤ 16）：")
print(f"{'kind':>7} {'a':>6} {'τ':>7} {'单次占比%':>11} {'截断误差%':>11} "
f"{'阶数和':>8} {'闭式和':>8}")
for r in order_rows:
print(f"{r['kind']:>7} {r['a']:>6} {r['tau']:>7} {r['singleFracPct']:>11} "
f"{r['truncErrPct']:>11} {r['sumQ16']:>8} {r['totalQ16']:>8}")
return 0
if __name__ == "__main__":
sys.exit(main())
跑一遍：
# snippet: 本机实跑（工作区 lantern/）
cd lantern
/opt/homebrew/bin/python3 sim/cg_trace_09_1_volume_equation.py
== 9.1 体积渲染方程 ==
介质系数 : σ_a=0.1 σ_s=0.4 σ_t=0.5（单次散射反照率 0.8）
定点 exp 最大相对误差 : 2.262e-03（2262.4434 ppm）
x(Q16) fexp ref diff ppm
32768 39768 39750 18 452.8302
65536 24117 24109 8 331.8263
98304 14624 14623 1 68.3854
131072 8874 8869 5 563.7614
196608 3263 3263 0 0.0
262144 1200 1200 0 0.0
327680 441 442 -1 2262.4434
393216 162 162 0 0.0
HG 相函数（∫p̃ dcosθ 应 = 131072，⟨cosθ⟩ 应 = g）：
g norm normErr asym asymErr 前向 p̃ 后向 p̃ 前/后
0.000 131072 0 -1 -1 65536 65536 1.0
0.100 131071 -1 6553 -1 89005 48746 1.8259
0.200 131071 -1 13106 -1 122886 36410 3.3751
0.300 131070 -2 19659 -2 173886 27145 6.4058
0.400 131071 -1 26213 -1 254879 20062 12.7046
0.500 131069 -3 32766 -2 393216 14563 27.001
0.600 131068 -4 39319 -3 655563 10240 64.0198
单次 vs 多次散射（泊松阶数模型，k ≤ 16）：
kind a τ 单次占比% 截断误差% 阶数和 闭式和
albedo 0.2 2.0 81.7391 18.2609 345 346
albedo 0.4 2.0 65.4292 34.5708 862 865
albedo 0.6 2.0 51.8043 48.1957 1635 1637
albedo 0.7 2.0 45.8682 54.1318 2154 2157
albedo 0.8 2.0 40.524 59.476 2786 2790
albedo 0.85 2.0 38.0349 61.9651 3155 3158
albedo 0.9 2.0 35.6922 64.3078 3561 3564
tau 0.8 0.5 81.4028 18.5972 1554 1554
tau 0.8 1.0 65.347 34.653 2349 2351
tau 0.8 2.0 40.524 59.476 2786 2790
tau 0.8 3.0 23.9892 76.0108 2597 2604
tau 0.8 4.0 13.6039 86.3961 2242 2248
2262 ppm 是这一章的地板
介质取 σ_a = 0.1、σ_s = 0.4、σ_t = 0.5（单次散射反照率 a = σ_s/σ_t = 0.8）。第一张表的八行是定点 exp(-x) 与浮点 exp(-x) 的对照：
最大相对误差 2262 ppm ≈ 0.226%，出现在 x = 5.0（exp(-5) ≈ 0.0067，定点值只剩 441，一个最低位就值 1/441 ≈ 2268 ppm）。这不是 LUT 建得差，而是定点数在深衰减处的分辨率极限：当结果本身只剩几百个最低位，任何实现都不可能更准。记住这条地板——9.2 里步长再密，误差也压不到 2262 ppm 以下。
相函数：HG 的不动点
Henyey–Greenstein 相函数把「散射有多偏前向」压成一个参数 g ∈ (-1, 1)：
p̃(cosθ) = (1 - g²) / (1 + g² - 2g·cosθ)^{3/2} （p̃ = 4π·p）
它有两个可以拿整数去核的性质：归一化 ∫_{-1}^{1} p̃ dcosθ = 2，以及各向异性 ⟨cosθ⟩ = g。脚本用 4000 段中点法对 cosθ 求积：
| g | ∫p̃ dcosθ | 误差 | ⟨cosθ⟩ | 误差 | 前向 p̃ | 后向 p̃ | 前/后 |
|---|---|---|---|---|---|---|---|
| 0.0 | 131072 | 0 | −1 | −1 | 65536 | 65536 | 1.0 |
| 0.1 | 131071 | −1 | 6553 | −1 | 89005 | 48746 | 1.83 |
| 0.2 | 131071 | −1 | 13106 | −1 | 122886 | 36410 | 3.38 |
| 0.3 | 131070 | −2 | 19659 | −2 | 173886 | 27145 | 6.41 |
| 0.4 | 131071 | −1 | 26213 | −1 | 254879 | 20062 | 12.70 |
| 0.5 | 131069 | −3 | 32766 | −2 | 393216 | 14563 | 27.00 |
| 0.6 | 131068 | −4 | 39319 | −3 | 655563 | 10240 | 64.02 |
目标值 2·ONE = 131072、g 本身的 Q16 值分别是 0/6554/13107/19661/26214/32768/39322。两列误差都不超过 4 个最低位——∫p̃ dcosθ 恒等于 131072±4，⟨cosθ⟩ 恒等于 g 自己±3。这说明相函数的定点实现是保归一化的：它不会凭空多出或少掉能量。前向/后向峰值比从 1.0 一路涨到 64.0，g 越大散射越偏前向——9.3 会把这个参数对整幅图的影响量出来。
单次 vs 多次散射
体积渲染方程里那个 ∫_{S²} 是递归的：散射进来的 L_i 自己又满足同一个方程。工程上最常见的近似是只算一次散射（single scattering），把 L_i 换成「直接来自光源的光」。这个近似漏掉多少？脚本用一个能量守恒的泊松阶数模型回答：沿光学深度 τ 的路径，散射 k 次的概率是 τ^k e^{-τ}/k!，每次保住能量比例 a，故第 k 阶贡献 L_k = L_in · p · a^k · τ^k · e^{-τ}/k!。L_1 就是单次散射，Σ_k L_k 是所有阶之和。
固定 τ = 2.0，反照率 a 从 0.2 涨到 0.9，单次占比从 81.74% 一路掉到 35.69%（截断误差从 18.26% 涨到 64.31%）；固定 a = 0.8，光学深度 τ 从 0.5 涨到 4.0，单次占比从 81.40% 掉到 13.60%（截断误差 18.60% → 86.40%）。两条曲线都说明同一件事：介质越「白」（a 越大）、越厚（τ 越大），单次散射越不够用。本节这组系数 a = 0.8、τ = 2.0 正好落在中间——单次散射只占 40.52%，漏掉 59.48% 的能量。所以 9.3 的散射图会把「单次」与「单次+二阶」并排放，让差额看得见。
另一件值得看的事是阶数和与闭式和贴得多紧：a = 0.8 那行是 2786（阶数和）对 2790（闭式和），差 4，0.14%；全部十二行都在 0.2% 内吻合。闭式和 L_in·p·(e^{-(1-a)τ} − e^{-τ}) 是无穷阶和的解析值，两边独立算出来却对得上，说明 k ≤ 16 的截断已经把能量收得差不多了，剩下的偏差来自「只数到 16 阶」和定点取整。
小结
- 体积渲染方程 = 沿路径的积分：
T(t)·[σ_a·L_e + σ_s·∫p·L_i]加上未碰撞的T(s)·L_surface。 - 定点
exp(-x)与浮点相差最大 2262 ppm（在深衰减处），这是本章一切误差的下界。 - HG 相函数在定点下保归一化：
∫p̃ dcosθ = 131072±4、⟨cosθ⟩ = g±3；g从 0 涨到 0.6，前向/后向峰值比从 1.0 涨到 64.0。本节系数下单次散射只占 40.5%，a或τ越大，多次散射越不可忽略。
下一节把这个积分真正离散化：沿光线取 N 个样点、把积分换成求和，然后量步长与误差的关系。
阅读导航：上一节：8.3 光子映射 · 下一节：9.2 ray marching 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。