9.2 ray marching
站内 /graphics-volumetric-clouds-atmosphere/ 那篇讲的是体积云工程：为了帧率，它必须用尽量少的步数把云画得像；本节反过来问一个更基础的问题——步数到底该取多少，以及步长过大时画面会坏成什么样。9.1 已经把透射率与相函数钉在整数上，本节把那个沿路径的积分真正换成求和，然后把「离散化」这件事量成一张表。
把积分换成求和
体积渲染方程是连续积分，计算机只会求和。最直接的换法是中点法（midpoint rule）：把 [0, s] 等分成 N 段，每段长 Δs = s/N，用段中点的值代表整段：
∫_0^s f(t) dt ≈ Σ_{i=0}^{N-1} f(t_0 + i·Δs + Δs/2) · Δs
中点法的误差是 O(Δs²)——步长减半，误差降到四分之一。这是理论。实际呢？本节量三件事：中点法单次散射相对解析解收敛多快、透射率「逐段连乘」的偏差怎么漂、以及一阶近似什么时候塌成条带。
步长 vs 误差：二阶收敛撞上定点地板
路径取 8 个世界单位，介质 σ_t = 0.5，于是整条路径的光学深度 τ_span = 4.0。解析单次散射（9.1 的闭式）在 Q16.16 下是 4095。脚本从 N = 1 扫到 N = 512（每档翻倍）：
# snippet: sim/cg_trace_09_2_ray_marching.py（节选：中点法与两层介质）
# verify: /opt/homebrew/bin/python3 sim/cg_trace_09_2_ray_marching.py
def march_two_layer(n_steps, span, y_top, y_iface, sigma_hi, sigma_lo):
"""双层介质的透射率：从 `y_top` 往下走 `span`，`y_iface` 以下是另一层。
每步用**步中点**所在层的 `σ_t`（ray marching 的标准做法）。步长与界面不重合时，
跨界面那一步会用整段的单一 `σ_t`，于是误差随 `N` 震荡。
"""
ds = span // n_steps
if ds <= 0:
return 0, 0
acc = 0 # 累积光学深度 τ（Q16.16）
for i in range(n_steps):
y_mid = y_top - (i * ds + ds // 2) # 步中点高度（往下走，y 减小）
sigma = sigma_hi if y_mid >= y_iface else sigma_lo
acc += fmul(sigma, ds)
return med.fexp_neg(acc), acc
| N | Δs(Q16) | σ_t·Δs | 中点法 | 相对误差% | 比前档 |
|---|---|---|---|---|---|
| 1 | 524288 | 4.0000 | 2259 | −44.8352 | — |
| 2 | 262144 | 2.0000 | 3486 | −14.8718 | 3.01 |
| 4 | 131072 | 1.0000 | 3932 | −3.9805 | 3.74 |
| 8 | 65536 | 0.5000 | 4056 | −0.9524 | 4.18 |
| 16 | 32768 | 0.2500 | 4088 | −0.1709 | 5.57 |
| 32 | 16384 | 0.1250 | 4095 | 0.0000 | inf |
| 64 | 8192 | 0.0625 | 4096 | +0.0244 | — |
| 128 | 4096 | 0.0312 | 4094 | −0.0244 | −1.00 |
| 256 | 2048 | 0.0156 | 4089 | −0.1465 | 0.17 |
| 512 | 1024 | 0.0078 | 4082 | −0.3175 | 0.46 |
这张表里有两段完全不同的行为，值得分开读：
- N ≤ 16 是真正的二阶收敛。 误差从 −44.84% 降到 −0.17%，每翻一倍步数误差大约降到 1/4（比值列 3.01 → 3.74 → 4.18 → 5.57，在 4 附近浮动）。这一段是中点法本身的截断误差，符合
O(Δs²)的理论预期。 - N ≥ 32 之后误差不再单调下降，反而回升。
N = 32恰好命中参考值 4095（误差 0），再加密到 512，误差从 0 变成 −0.32%。这不是中点法变差了，而是定点累加的地板：9.1 已经量出exp(-x)本身的相对误差上限是 2262 ppm ≈ 0.226%，N = 512的 −0.3175% 已经落在这个量级附近。加密步长买不到比定点精度更准的结果。
结论很实用：本节这组系数下，N ≈ 32（σ_t·Δs = 0.125）就已经踩在定点地板上，再往上加步数纯属浪费。至于「画面看不看得出差别」，N = 8（误差 −0.95%）与 N = 32 之间的肉眼差异，得靠 9.3 的图像 MAE 去判，而不是靠感觉。
透射率链式累乘的偏差
实现里另一个常见写法是：不用解析的 exp(-σ_t·s)，而是把路径切成 N 段、每段各算 exp(-σ_t·Δs) 再连乘。数学上 exp(-σ_t·s) = ∏ exp(-σ_t·Δs) 恒等，但每段的定点取整会一路累积。参考值仍是解析的 1200：
| N | 链乘 T | 解析 T | 偏差% |
|---|---|---|---|
| 1 | 1200 | 1200 | 0.00000 |
| 2 | 1201 | 1200 | +0.08333 |
| 4 | 1201 | 1200 | +0.08333 |
| 8 | 1203 | 1200 | +0.25000 |
| 16 | 1203 | 1200 | +0.25000 |
| 32 | 1202 | 1200 | +0.16667 |
| 64 | 1194 | 1200 | −0.50000 |
| 128 | 1187 | 1200 | −1.08333 |
| 256 | 1172 | 1200 | −2.33333 |
| 512 | 1138 | 1200 | −5.16667 |
偏差不是单调的：N 小的时候偏高（0.08% → 0.25%），到 N = 64 之后转成偏低，并且越偏越多，N = 512 时已经 −5.17%。原因很直白——每段 exp(-σ_t·Δs) 的定点误差方向固定，段数一多，同号的取整误差就线性累积。这不是「步数越多越准」：它是一条 U 形曲线，最准的点在 N ≈ 4~32 之间。所以实现里宁可保留一个解析的 exp(-σ_t·s)，也不要为了「通用」把透射率拆成 N 段连乘。
一阶近似：条带出现的阈值
比连乘更省的写法是把 exp(-σ_t·Δs) 一阶展开成 1 - σ_t·Δs。它在一段内很准，但一旦 σ_t·Δs ≥ 1，这个因子就 ≤ 0——透射率塌成 0：
| N | σ_t·Δs | 线性 T | 坏步数 |
|---|---|---|---|
| 1 | 4.0000 | 0 | 1 |
| 2 | 2.0000 | 0 | 2 |
| 4 | 1.0000 | 0 | 4 |
| 8 | 0.5000 | 256 | 0 |
| 16 | 0.2500 | 656 | 0 |
| 32 | 0.1250 | 910 | 0 |
| 64 | 0.0625 | 1047 | 0 |
| 128 | 0.0312 | 1112 | 0 |
| 256 | 0.0156 | 1132 | 0 |
| 512 | 0.0078 | 1120 | 0 |
阈值就是 σ_t·Δs = 1：N ≤ 4 时每段的光学深度 ≥ 1，因子被钳到 0，整条链乘立刻归零。这就是条带（banding）的数学来源——不是「看起来有点假」，而是透射率在几何上直接断成 0。更隐蔽的是 N = 8~512 这一段：T 值本身在 256~1132 之间漂，且并不单调趋向解析值 1200（N = 256 的 1132 比 N = 64 的 1047 更接近，但 N = 512 的 1120 又退回去了）。一阶近似永远偏低，因为 1 - x ≤ e^{-x} 恒成立——这是它在任何步长下都消不掉的系统性偏差。
双层介质：误差会震荡
最后看一个把「步长与几何不重合」暴露得最清楚的例子：介质分两层，界面在 4 个单位处，上层 σ_hi = 0.5、下层 σ_lo = 1.0，光线往下走 6 个单位跨过界面。解析透射率是 exp(-(0.5·4 + 1.0·2)) = exp(-4) = 1200。
| N | Δs(Q16) | march T | 相对误差% |
|---|---|---|---|
| 1 | 393216 | 3263 | +171.92 |
| 2 | 196608 | 728 | −39.33 |
| 3 | 131072 | 1200 | 0.0000 |
| 4 | 98304 | 1541 | +28.42 |
| 5 | 78643 | 982 | −18.17 |
| 6 | 65536 | 1200 | 0.0000 |
| 7 | 56173 | 1385 | +15.42 |
| 8 | 49152 | 1060 | −11.67 |
| 9 | 43690 | 1201 | +0.0833 |
| 12 | 32768 | 1200 | 0.0000 |
这里误差不收敛，而是震荡：N 是 3 的倍数时界面恰好落在步边界上，误差为 0；不是 3 的倍数时，跨界面那一步用整段的单一 σ_t，误差符号随 N 在正负之间翻。N = 1 时误差高达 +172%（整条 6 单位路径被当成上层 0.5 算，τ 少了一半多）。只要介质有内界面，增大 N 并不能保证误差变小——必须让界面与步边界对齐，否则误差会一直在 ±10% 量级来回跳。这也是为什么工程实现里会把「介质边界」和「步长」绑在一起考虑，而不是无脑加密。
P 层把数钉死
P 层用 TypeScript 重写同一条链路，全程 BigInt，断言对象是上面三张表加两层介质的每一行：
// file: lantern/core/test/cg_trace_09_2_march.test.ts
// tier: P
// verify: node --test 'core/test/cg_trace_09_2_march.test.ts'
// 9.2 的 P 复现：BigInt 内核重跑步长扫描、透射率链式累乘、一阶近似条带、双层介质，
// 逐位对齐 Python 侧 out/cg_trace_09_2_march.json。
import assert from "node:assert/strict";
import { readFileSync } from "node:fs";
import { dirname, join } from "node:path";
import { test } from "node:test";
import { fileURLToPath } from "node:url";
import { ONE } from "../src/lantern/fixed.ts";
import {
INV_4PI_Q16, analyticSingleScatter, marchSingleScatter, marchTransmittance, marchTwoLayer,
} from "../src/lantern/cg_trace_09_medium.ts";
const root = join(dirname(fileURLToPath(import.meta.url)), "..", "..");
const rep = JSON.parse(
readFileSync(join(root, "out", "cg_trace_09_2_march.json"), "utf8"),
) as {
sigmaSQ16: number; sigmaTQ16: number; spanQ16: number; refSingleQ16: number;
refTransQ16: number;
convergence: { nSteps: number; marchQ16: number; refQ16: number; relErrPct: number; tEndQ16: number }[];
transmittanceChain: { nSteps: number; chainQ16: number; refQ16: number; badSteps: number }[];
linearBanding: { nSteps: number; sigmaTdsQ16: number; linearQ16: number; badSteps: number }[];
twoLayer: {
yTopQ16: number; yIfaceQ16: number; sigmaHiQ16: number; sigmaLoQ16: number; walkQ16: number;
refQ16: number;
rows: { nSteps: number; tauQ16: number; marchQ16: number; relErrPct: number }[];
};
};
const SIGMA_S = BigInt(rep.sigmaSQ16);
const SIGMA_T = BigInt(rep.sigmaTQ16);
const LIGHT = 65536n;
const SPAN = BigInt(rep.spanQ16);
test("【核心】中点法单次散射与解析解、与 Python 侧逐位相同", () => {
const ref = analyticSingleScatter(SIGMA_S, SIGMA_T, SPAN, LIGHT, INV_4PI_Q16);
assert.equal(ref, BigInt(rep.refSingleQ16));
for (const r of rep.convergence) {
const [val, tEnd] = marchSingleScatter(SIGMA_S, SIGMA_T, r.nSteps, SPAN, LIGHT, INV_4PI_Q16);
assert.equal(val, BigInt(r.marchQ16), `N=${r.nSteps}`);
assert.equal(tEnd, BigInt(r.tEndQ16), `N=${r.nSteps} tEnd`);
}
});
test("【核心】透射率链式累乘与一阶近似条带逐位相同", () => {
const refT = marchTransmittance(SIGMA_T, 1, SPAN, "exact")[0];
assert.equal(refT, BigInt(rep.refTransQ16));
for (const r of rep.transmittanceChain) {
const [t, bad] = marchTransmittance(SIGMA_T, r.nSteps, SPAN, "exact");
assert.equal(t, BigInt(r.chainQ16), `exact N=${r.nSteps}`);
assert.equal(bad, r.badSteps);
}
for (const r of rep.linearBanding) {
const [t, bad] = marchTransmittance(SIGMA_T, r.nSteps, SPAN, "linear");
assert.equal(t, BigInt(r.linearQ16), `linear N=${r.nSteps}`);
assert.equal(bad, r.badSteps);
}
});
test("条带阈值：σ_t·Δs ≥ 1 时一阶近似透射率塌成 0", () => {
for (const r of rep.linearBanding) {
if (BigInt(r.sigmaTdsQ16) >= ONE) {
assert.equal(r.linearQ16, 0, `N=${r.nSteps} 应塌成 0`);
assert.ok(r.badSteps > 0);
} else {
assert.ok(r.linearQ16 > 0, `N=${r.nSteps} 不应为 0`);
assert.equal(r.badSteps, 0);
}
}
});
test("中点法二阶收敛：N 翻倍误差约降到 1/4（N ≤ 16 段）", () => {
const conv = rep.convergence;
for (let i = 1; i < conv.length && conv[i].nSteps <= 16; i += 1) {
const ratio = Math.abs(conv[i - 1].relErrPct) / Math.abs(conv[i].relErrPct);
assert.ok(ratio > 2.5 && ratio < 6.0, `N=${conv[i].nSteps} ratio=${ratio}`);
}
// N = 32 时误差已落到 0（定点地板），继续加密反而回升
const at32 = conv.find((r) => r.nSteps === 32)!;
const at512 = conv.find((r) => r.nSteps === 512)!;
assert.equal(at32.relErrPct, 0);
assert.ok(Math.abs(at512.relErrPct) > 0);
});
test("【核心】双层介质（步长与界面不重合）逐位相同，且 N 为 3 的倍数时对齐", () => {
const tl = rep.twoLayer;
for (const r of tl.rows) {
const [t, tau] = marchTwoLayer(
r.nSteps, BigInt(tl.walkQ16), BigInt(tl.yTopQ16), BigInt(tl.yIfaceQ16),
BigInt(tl.sigmaHiQ16), BigInt(tl.sigmaLoQ16),
);
assert.equal(t, BigInt(r.marchQ16), `N=${r.nSteps}`);
assert.equal(tau, BigInt(r.tauQ16), `N=${r.nSteps} tau`);
}
// 界面在 4 单位、步长 6/N：N 为 3 的倍数时界面落在步边界上，误差为 0
for (const r of tl.rows) {
if (r.nSteps % 3 === 0 && r.nSteps <= 6) assert.equal(r.marchQ16, tl.refQ16, `N=${r.nSteps}`);
}
// 震荡：误差符号在相邻档之间翻转（N 从 1 到 6）
const head = tl.rows.slice(0, 6);
for (let i = 1; i < head.length; i += 1) {
if (head[i].relErrPct !== 0 && head[i - 1].relErrPct !== 0) {
assert.ok(Math.sign(head[i].relErrPct) !== Math.sign(head[i - 1].relErrPct));
}
}
});
跑一遍：
# snippet: 本机实跑（工作区 lantern/）
cd lantern
/opt/homebrew/bin/python3 sim/cg_trace_09_2_ray_marching.py
== 9.2 ray marching ==
介质 σ_t=0.5，路径长 8.0 单位，τ_span=4.0
解析单次散射参考 = 4095（Q16.16）
N Δs(Q16) σtΔs 中点法 相对误差% 比前档
1 524288 4.0000 2259 -44.8352
2 262144 2.0000 3486 -14.8718 3.01
4 131072 1.0000 3932 -3.9805 3.74
8 65536 0.5000 4056 -0.9524 4.18
16 32768 0.2500 4088 -0.1709 5.57
32 16384 0.1250 4095 0.0000 inf
64 8192 0.0625 4096 0.0244
128 4096 0.0312 4094 -0.0244 -1.00
256 2048 0.0156 4089 -0.1465 0.17
512 1024 0.0078 4082 -0.3175 0.46
透射率链式累乘偏差（精确指数逐段连乘 vs 解析）：
N 链乘 T 解析 T 偏差%
1 1200 1200 0.00000
2 1201 1200 0.08333
4 1201 1200 0.08333
8 1203 1200 0.25000
16 1203 1200 0.25000
32 1202 1200 0.16667
64 1194 1200 -0.50000
128 1187 1200 -1.08333
256 1172 1200 -2.33333
512 1138 1200 -5.16667
一阶近似条带（σ_t·Δs ≥ 1 时透射率塌成 0）：
N σtΔs 线性 T 坏步数
1 4.0000 0 1
2 2.0000 0 2
4 1.0000 0 4
8 0.5000 256 0
16 0.2500 656 0
32 0.1250 910 0
64 0.0625 1047 0
128 0.0312 1112 0
256 0.0156 1132 0
512 0.0078 1120 0
双层介质（界面在中间，步长不重合 -> 震荡）：
N Δs(Q16) march T 解析 T 相对误差%
1 393216 3263 1200 171.9167
2 196608 728 1200 -39.3333
3 131072 1200 1200 0.0000
4 98304 1541 1200 28.4167
5 78643 982 1200 -18.1667
6 65536 1200 1200 0.0000
7 56173 1385 1200 15.4167
8 49152 1060 1200 -11.6667
9 43690 1201 1200 0.0833
10 39321 1327 1200 10.5833
11 35746 1096 1200 -8.6667
12 32768 1200 1200 0.0000
13 30247 1296 1200 8.0000
14 28086 1118 1200 -6.8333
15 26214 1201 1200 0.0833
16 24576 1278 1200 6.5000
17 23130 1132 1200 -5.6667
18 21845 1201 1200 0.0833
19 20695 1266 1200 5.5000
20 19660 1142 1200 -4.8333
21 18724 1201 1200 0.0833
22 17873 1256 1200 4.6667
23 17096 1150 1200 -4.1667
24 16384 1200 1200 0.0000
小结
- 中点法是
O(Δs²)：N翻倍误差降到约 1/4，直到定点地板（≈2262 ppm）接手。 - 本节系数下
N ≈ 32（σ_t·Δs = 0.125）已到地板，再加密只会在 ±0.3% 之间来回。 - 透射率逐段连乘的偏差是 U 形：
N = 512时偏到 −5.17%，不如保留解析exp(-σ_t·s)。 - 一阶近似
1 - σ_t·Δs在σ_t·Δs ≥ 1时透射率塌成 0——这就是条带的阈值。 - 介质有内界面时误差震荡而非收敛，必须让界面与步边界对齐。
下一节把步进搬进真正带几何遮挡的场景：体积阴影、单次与二阶散射，以及各向异性相函数对整幅图的影响。
阅读导航：上一节：9.1 体积渲染方程 · 下一节：9.3 体积阴影与散射 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。