9.3 体积阴影与散射
站内 /graphics-volumetric-clouds-atmosphere/
那篇为了帧率，把体积阴影做成一张低分辨率的三维贴图、再模糊；本节不追求实时，而是把「介质里的阴影到底是怎么算出来的」拆成可断言的数：一条阴影光线给出一整条透射率曲线，两次散射的贡献能量化成百分比，相函数的各向异性可以一次渲染重加权出任意 g 的整幅图。9.1 给了方程与误差地板，9.2 给了步长与误差的关系，本节把两者放进一个带几何遮挡体的场景里。
阴影光线在介质里不再回答「挡不挡」
第 4.1 节的阴影光线是二值的：挡 = 0，不挡 = 1。介质里没有这个二值答案——光从采样点射向光源，一路上都在被衰减。所以介质里的阴影光线返回的是一条透射率曲线 T_light ∈ (0, 1]，由沿光路的步进累积得到；只有当几何里真有个不透明体挡在中间时，它才截断为零。
场景沿用 9.1 的介质（σ_t = 0.5），水平板 y ∈ [0, 8]，方向光斜射，另在光路上放两个球当遮挡体。渲染 160×120、视线步 16、阴影步 8，共 13910 个像素的视线穿过介质。
# snippet: sim/cg_trace_09_3_volume_shadow_scattering.py（节选：阴影透射率与两次散射的累积）
# verify: /opt/homebrew/bin/python3 sim/cg_trace_09_3_volume_shadow_scattering.py
def light_transmittance(p, ldir, n_light, sigma_t, y_lo, y_hi, tmax):
"""从介质内一点 `p` 朝光源方向的透射率：沿阴影光线穿过介质累积 `exp(-σ_t Δs)`。
方向光在无穷远，故路径终点是介质板的出射面。返回 Q16.16 的 `T_light ∈ [0, ONE]`。
"""
span = med.ray_slab_span(p[1], ldir[1], y_lo, y_hi, tmax)
if span is None:
return ONE
length = span[1] - span[0]
if length <= 0:
return ONE
return med.march_transmittance(sigma_t, n_light, length, "exact")[0]
# …… 视线步进主循环（每个采样点）……
tl = light_transmittance(p, ldir, n_light, sigma_t, y_lo, y_hi, tmax)
if cgtrace.bvh_intersect(nodes, order, tris, p, ldir, tmax) is not None:
tl = 0 # 几何遮挡：透射率直接截断为零
# 单次散射：T_view · σ_s · T_light · L · Δs（5 个 Q16.16 因子 = Q80.64）
# 归一到 Q32.32 需 >> 48（80 − 32 = 48），最后统一 >> 16 得 Q16.16。
acc += (tv * sigma_s * tl * light * ds) >> 48
# 二阶散射：采一条均匀球面方向 ω'，估一阶散射光再散射回视线的量
li = (l2 * tq * sigma_s * ph1 * tlq * light) >> 80 # 6 因子 Q96.96 -> Q16.16
term += (tv * sigma_s * med.FOURPI_Q16 * ph2 * li * ds) >> 64 # -> Q32.32
跑一遍：
# snippet: 本机实跑（工作区 lantern/）
cd lantern
/opt/homebrew/bin/python3 sim/cg_trace_09_3_volume_shadow_scattering.py
== 9.3 体积阴影与散射 ==
介质 σ_t=0.5，分辨率 160x120，视线步 16 / 阴影步 8，命中介质 13910 像素
y 光程 解析 T 步2 步8 步32 被挡
1.0 12.12 152 152 152 150 False
2.0 10.39 363 363 364 361 False
3.0 8.66 862 862 865 863 False
4.0 6.93 2052 2052 2056 2054 False
6.0 3.46 11598 11598 11613 11598 True
7.0 1.73 27570 27587 27586 27581 True
单次散射均值 : 2372.3468（1354 ms） 阴影光线 222560，被挡 2369
单次+二阶均值 : 2802.6451（24952 ms） 二阶贡献 18.1381% 二阶样本 1780480
两图 MAE : 19.4097
各向异性（一次渲染重加权）：
g 均值辐射 前向 p̃ 后向 p̃ 前/后
0.000 2372 65536 65536 1.0
0.100 1915 89005 48746 1.8259
0.200 1523 122886 36410 3.3751
0.300 1190 173886 27145 6.4058
0.400 911 254879 20062 12.7046
0.500 677 393216 14563 27.001
0.600 484 655563 10240 64.0198
体积阴影：解析 vs 步进
在介质里沿方向光取六个高度，把「从该点射向光源」的透射率用解析式 exp(-σ_t·d) 与 2/8/32 步的步进各算一遍，并记录该点的光路是否被球挡住：
| y | 光程 | 解析 T | 步2 | 步8 | 步32 | 被挡 |
|---|---|---|---|---|---|---|
| 1.0 | 12.12 | 152 | 152 | 152 | 150 | False |
| 2.0 | 10.39 | 363 | 363 | 364 | 361 | False |
| 3.0 | 8.66 | 862 | 862 | 865 | 863 | False |
| 4.0 | 6.93 | 2052 | 2052 | 2056 | 2054 | False |
| 6.0 | 3.46 | 11598 | 11598 | 11613 | 11598 | True |
| 7.0 | 1.73 | 27570 | 27587 | 27586 | 27581 | True |
两点值得停一下。第一，步进与解析几乎逐位相同（差不超过 5），因为 σ_t·Δs 很小、9.1 的定点指数地板在这里远没触及。第二，y = 1.0（光程最长，12.12）解析 T 只有 152，而 y = 7.0 的 T 高达 27570——越靠近板的顶面，射向光源的路越短，透射率越高。被挡的两行（y = 6、y = 7）在几何上确实被球截住，但注意它们的 T_light 仍然被算成了非零的解析值——遮挡是「乘以 0」这一步施加的，不是透射率函数本身的责任。这个分工很重要：透射率只管介质，几何遮挡由 BVH 求交单独判定，两者在着色时相乘。
单次 vs 二阶散射
把上面的单次散射累积成整幅图，再额外算一次散射（在采样点沿均匀球面方向采一条再散射方向，估一阶散射光再散射回视线的量）。方向采样走 cgtrace.rng_next 的 SplitMix64，每像素 8 个样本（secondOrderSpp = 8），保证两侧逐位一致：
| 量 | 单次散射 | 单次 + 二阶 |
|---|---|---|
| 介质内平均辐射亮度(Q16) | 2372.3468 | 2802.6451 |
| 二阶贡献 | — | 18.1381% |
| 样本数 | 222560 视线样本 | 1780480 二阶样本 |
| 阴影光线 / 被挡 | 222560 / 2369 | — |
| 图 sha256 | 649c9823… | 647a8c82… |
| 两图 MAE | — | 19.4097 |
二阶散射贡献 18.14%——这个数不是随便来的，它要和 9.1 的泊松阶数模型对上：9.1 里 a = 0.8、τ = 2.0 的「截断误差」是 59.48%，那是「只算单次、漏掉 1..16 全部高阶」的总和；本节的 18.14% 是「单次 → 单次+二阶」再加一阶的增量。两者口径不同，但都指向同一结论：多次散射不可忽略。这里 τ_span = 4.0 比 9.1 的 2.0 更厚，单次占比自然更低，二阶的相对增量 18.14% 落在合理区间。
一张 160×120 的图，单次散射约 1354 ms，加上二阶采样涨到约 24952 ms——二阶散射的代价接近单次的 18 倍（每个视线样本要多采 8 条方向、每条再走一次阴影步进）。这是「质量换时间」最直白的一笔账：想拿到那 18% 的能量，就得付出近 20 倍的时间。
各向异性：一次渲染重加权
Henyey–Greenstein 相函数的 g 只通过 p(cosθ) 进入散射项，而 cosθ = d·l_dir 沿一条主光线是常数。所以一次「不含相函数」的渲染，就能重加权出任意 g 的整幅图——只要把每个像素的权重乘上 p(cosθ, g)：
| g | 均值辐射 | 前向 p̃ | 后向 p̃ | 前/后 |
|---|---|---|---|---|
| 0.0 | 2372 | 65536 | 65536 | 1.0 |
| 0.1 | 1915 | 89005 | 48746 | 1.83 |
| 0.2 | 1523 | 122886 | 36410 | 3.38 |
| 0.3 | 1190 | 173886 | 27145 | 6.41 |
| 0.4 | 911 | 254879 | 20062 | 12.70 |
| 0.5 | 677 | 393216 | 14563 | 27.00 |
| 0.6 | 484 | 655563 | 10240 | 64.02 |
注意均值辐射随 g 单调下降（2372 → 484）：这不是能量凭空消失，而是相函数把能量重新分配了。g 越大，越多的光被推到前向（cosθ 大）的方向，而本节这组相机-光源夹角下，视线恰好落在后向附近，于是整幅图变暗。前向/后向峰值比从 1.0 涨到 64.0，与 9.1 的相函数表逐位相同——两侧用的是同一个定点实现。这条「一次渲染重加权」的捷径成立的前提是 cosθ 沿光线恒定，这一点对方向光成立、对点光源不成立（点光源的 cosθ 逐采样点变化）。
P 层把数钉死
P 层用 TypeScript 重写整条链路，全程 BigInt，断言对象是阴影探针表、24×18 预览图的像素缓冲 sha256、二阶散射的整数和，以及各向异性表：
// file: lantern/core/test/cg_trace_09_3_shadow.test.ts
// tier: P
// verify: node --test 'core/test/cg_trace_09_3_shadow.test.ts'
// 9.3 的 P 复现：BigInt 内核重算体积阴影的透射率表、相函数各向异性表，并重画 24x18
// 预览图——像素缓冲 sha256 与二阶散射均值必须与 Python 侧逐位相同。
import assert from "node:assert/strict";
import { createHash } from "node:crypto";
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 { bvhIntersect } from "../src/lantern/cgtrace.ts";
import { hgPhaseTimes4pi, marchTransmittance, raySlabSpan, transmittance } from "../src/lantern/cg_trace_09_medium.ts";
import {
lightTransmittance, pngFrom, render, sceneBvh, type MediumCfg,
} from "../src/lantern/cg_trace_09_3_shadow.ts";
const root = join(dirname(fileURLToPath(import.meta.url)), "..", "..");
const cfg = JSON.parse(
readFileSync(join(root, "config", "cg_trace_09_medium.json"), "utf8"),
) as MediumCfg;
const rep = JSON.parse(
readFileSync(join(root, "out", "cg_trace_09_3_shadow.json"), "utf8"),
) as {
sigmaTQ16: number; viewSteps: number; lightSteps: number; insidePixels: number;
shadowProbe: {
yQ16: number; pathQ16: number; analyticQ16: number; occluded: boolean;
march2Q16: number; march8Q16: number; march32Q16: number;
}[];
singleOrder: { meanRadianceQ16: number; shadowRays: number; occluded: number; sha256: string };
secondOrder: { contributionPct: number };
anisotropy: { gQ16: number; meanRadianceQ16: number; forwardQ16: number; backwardQ16: number }[];
preview: { width: number; height: number; insidePixels: number; sha256: string; secondOrderSumQ16: number };
};
const SIGMA_T = BigInt(rep.sigmaTQ16);
const Y_LO = BigInt(cfg.slabYLoQ16);
const Y_HI = BigInt(cfg.slabYHiQ16);
const TMAX = BigInt(cfg.camera.farQ16);
test("【核心】体积阴影透射率表：解析 vs 步进，与 Python 侧逐位相同", () => {
const { tris, bvh } = sceneBvh(cfg);
for (const r of rep.shadowProbe) {
const p: [bigint, bigint, bigint] = [0n, BigInt(r.yQ16), 0n];
const span = raySlabSpan(p[1], BigInt(cfg.light.dirQ16[1]), Y_LO, Y_HI, TMAX)!;
const length = span[1] - span[0];
assert.equal(length, BigInt(r.pathQ16));
assert.equal(transmittance(SIGMA_T, length), BigInt(r.analyticQ16));
assert.equal(marchTransmittance(SIGMA_T, 2, length, "exact")[0], BigInt(r.march2Q16));
assert.equal(marchTransmittance(SIGMA_T, 8, length, "exact")[0], BigInt(r.march8Q16));
assert.equal(marchTransmittance(SIGMA_T, 32, length, "exact")[0], BigInt(r.march32Q16));
// lightTransmittance 走 8 步步进，故等于 march8Q16（不是解析值）
assert.equal(lightTransmittance(p, [BigInt(cfg.light.dirQ16[0]), BigInt(cfg.light.dirQ16[1]),
BigInt(cfg.light.dirQ16[2])], rep.lightSteps, SIGMA_T, Y_LO, Y_HI, TMAX), BigInt(r.march8Q16));
const hit = bvhIntersect(bvh, tris, p, [BigInt(cfg.light.dirQ16[0]),
BigInt(cfg.light.dirQ16[1]), BigInt(cfg.light.dirQ16[2])], TMAX) !== null;
assert.equal(hit, r.occluded, `y=${r.yQ16}`);
}
});
test("【核心】24x18 预览图的像素缓冲 sha256 与 Python 侧相同", () => {
const pcfg: MediumCfg = { ...cfg, width: rep.preview.width, height: rep.preview.height };
const r = render(pcfg, rep.viewSteps, rep.lightSteps, false, BigInt(cfg.hgGQ16));
const insideCount = r.inside.flat().filter(Boolean).length;
assert.equal(insideCount, rep.preview.insidePixels);
const buf = pngFrom(pcfg, r.base, r.cosmap, r.inside, BigInt(cfg.hgGQ16));
assert.equal(createHash("sha256").update(buf).digest("hex"), rep.preview.sha256);
});
test("【核心】24x18 二阶散射均值与 Python 侧逐位相同", () => {
const pcfg: MediumCfg = { ...cfg, width: rep.preview.width, height: rep.preview.height };
const r = render(pcfg, rep.viewSteps, rep.lightSteps, true, BigInt(cfg.hgGQ16));
let sum = 0n;
let n = 0;
for (let py = 0; py < pcfg.height; py += 1) {
for (let px = 0; px < pcfg.width; px += 1) {
if (r.inside[py][px]) { sum += r.s2[py][px]; n += 1; }
}
}
assert.equal(sum, BigInt(rep.preview.secondOrderSumQ16));
assert.equal(n, rep.preview.insidePixels);
assert.equal(r.stats.secondOrderSamples, rep.preview.insidePixels * rep.viewSteps * cfg.secondOrderSpp);
});
test("各向异性：p̃ 前向 / 后向峰值与 Python 侧逐位相同，前向峰值随 g 单调上升", () => {
for (const r of rep.anisotropy) {
const g = BigInt(r.gQ16);
assert.equal(hgPhaseTimes4pi(g, ONE), BigInt(r.forwardQ16));
assert.equal(hgPhaseTimes4pi(g, -ONE), BigInt(r.backwardQ16));
}
const fwd = rep.anisotropy.map((r) => r.forwardQ16);
for (let i = 1; i < fwd.length; i += 1) assert.ok(fwd[i] > fwd[i - 1]);
});
test("二阶散射贡献量级合理：介于 10% 与 30% 之间（与 9.1 阶数模型同量级）", () => {
assert.ok(rep.secondOrder.contributionPct > 10 && rep.secondOrder.contributionPct < 30);
});
P 层为逐位复现准备了一张 24×18 的预览图（cg_trace_09_3_preview.png）：介质内 312 个像素，单次散射图 sha256 = 31f1d477dd4d32682b7f9f8c04ad74b6a4ebd97da60f826d457d60c2239250e1，二阶散射的整数和 = 138083。小图的意义在于把「相机光线、遮挡 BVH、SplitMix64 流」三件事一起钉死——只要其中任何一处与 Python 侧分叉，这个 sha256 就对不上。
小结
- 介质里的阴影光线返回透射率曲线，几何遮挡是「乘以 0」单独施加的，两者在着色时相乘。
- 体积阴影的步进与解析几乎逐位相同（差 ≤ 5），因为
σ_t·Δs远未触及 9.1 的定点地板。 - 单次散射均值 2372.35，二阶散射贡献 18.14%，代价却是近 20 倍时间。
- 各向异性可用一次渲染重加权出任意
g；本节参数下g从 0 到 0.6，均值辐射从 2372 掉到 484。 - 24×18 预览图的 sha256 把相机、BVH、RNG 三件事一起钉死。
下一章离开介质，去看程序化世界的构造：噪声、分形地形与场景生成。
阅读导航：上一节：9.2 ray marching · 下一节：10.1 噪声的构造 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。