6.3 能量守恒的验证
6.1 的 BRDF 和 6.2 的参数映射都建立在一个不成文的前提上：反射出去的能量不能超过入射的。BRDF 的公式本身给不出这个保证，只能实测。那篇讲实时 PBR 与阴影贴图的工程实现见 实时 PBR 与阴影贴图 ；本节把"能量不超 1"变成一个能跑、能断言的整数实验。
上界是什么
BRDF 的物理约束只有一条：对任意入射方向 $l$，反射到整个半球的总能量不超过入射能量，
$$ E(\mu_v)=\int_{\Omega} f_s(l,v),\cos\theta_l,d\omega_l\ \le\ 1,\qquad \mu_v=n\cdot v $$
这个量叫方向反射率。注意 $\le1$ 不是等式：真实材质会吸收，$E$ 通常小于 1，而差值就是被材质吃掉的部分。本节的全部工作就是把这个积分算出来，看它离 1 有多远。
为什么单次弹射必然丢能量
微表面 BRDF 的 $G$ 项只统计了"入射和出射各自不被挡住"的概率，它默认一次碰撞之后光就离开了。真实表面上光会在微面之间多次弹射——一次没反射出去的光可能撞第二次、第三次。所以单次弹射的 GGX 天然欠能量，欠多少取决于 $\alpha$。
有一个可以手算的锚点。取 $\alpha=1$、$v=n$、$F\equiv1$：此时 $D\equiv1/\pi$，$G_1(n)=1$，$G_1(l)=2\mu_l/(\mu_l+1)$，而 $h=\text{normalize}(l+n)$ 给出 $n\cdot h=\sqrt{(1+\mu_l)/2}$。代进 $f_s=DGF/\big(4(n\cdot l)(n\cdot v)\big)$ 之后这些项会成片消掉，只剩
$$ f_s=\frac{1}{2\pi,(1+\cos\theta_l)} $$
再积一次（换元 $u=\cos\theta$）：
$$ E=\int_0^{2\pi}!!\int_0^{\pi/2}\frac{\cos\theta}{2\pi(1+\cos\theta)}\sin\theta,d\theta,d\phi=\int_0^1\frac{u}{1+u},du=1-\ln 2=0.306854 $$
所以 $\alpha=1$ 的完全粗糙表面只反射了 30.7% 的能量，剩下 69.3% 全丢在"只碰一次"这个假设里。这是一个闭式解，拿来给整数实现当靶子最合适。
实测
#!/usr/bin/env python3
# file: lantern/sim/cg_shade_06_3_energy.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_shade_06_3_energy.py
"""6.3 能量守恒的验证：半球整数求和、GGX 残差表、单次 vs 多次弹射、kD = 1 - F。
测度走投影圆盘：面积元 dr ds 已含 cos(theta_l) d omega_l，故 ∫ f cosθ dω = Σ f · 格面积。
四个实验：方向反射率 E(mu)、半球平均 E_avg、Kulla-Conty 多次弹射补偿、kD 归一化。
产物：out/cg_shade_06_3_residual.csv、_ms.csv、_kd.csv、_energy.svg、_energy.json。
"""
import json
import os
import sys
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cglib # noqa: E402
from cg_shade_06_pbr import cell_of, disk_dirs, schlick_f, spec_brdf # noqa: E402
from cg_shade_06_util import svg # noqa: E402
from cglib import ONE, fdiv, fmul, fsqrt # noqa: E402
ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
PI, GRID, M, LN2 = 205887, 64, 8, 45426 # ln2 = round(ln 2 * 2^16)
ALPHAS = [("0.05", 164), ("0.10", 655), ("0.20", 2621), ("0.30", 5898),
("0.50", 16384), ("0.75", 36864), ("1.00", 65536)]
MUS = [("1.00", ONE), ("0.50", 32768), ("0.25", 16384), ("0.10", 6554)]
KD_ALPHAS, F0_KD = [("0.25", 4096), ("0.50", 16384), ("1.00", 65536)], 2621
def e_spec(a2, f0, dirs, n, mu):
"""方向反射率 E(mu_v) = ∫ f_s cos(theta_l) d omega_l；v = (sqrt(1-mu^2), 0, mu)。"""
zv, total = fsqrt(ONE - fmul(mu, mu)), 0
for lx, ly, lz in dirs:
hx, hz = lx + zv, lz + mu
hl = fsqrt(fmul(hx, hx) + fmul(ly, ly) + fmul(hz, hz))
if hl:
total += spec_brdf(a2, lz, mu, fdiv(hz, hl),
fdiv(fmul(zv, hx) + fmul(mu, hz), hl), f0)
return fmul(total, cell_of(n))
def lambert_kd(f0, dirs, n, mu, from_f):
"""漫反射项 ∫ (kD/pi) cos(theta_l) d omega_l，albedo = 1。
from_f 为真时 kD = 1 - F(v·h)（方向相关），否则 kD = 1 - f0（常数）。
注意圆盘度量已经含 cos，这里只补 1/pi——多乘一个 cos 会让 0.96 掉到 0.64。
"""
zv, total = fsqrt(ONE - fmul(mu, mu)), 0
for lx, ly, lz in dirs:
hx, hz = lx + zv, lz + mu
hl = fsqrt(fmul(hx, hx) + fmul(ly, ly) + fmul(hz, hz))
vh = fdiv(fmul(zv, hx) + fmul(mu, hz), hl)
total += fdiv(ONE - (schlick_f(f0, vh) if from_f else f0), PI)
return fmul(total, cell_of(n))
def e_avg(a2, dirs, n, m=M, f0=ONE):
"""半球平均反射率 E_avg = 2 ∫_0^1 E(mu) mu d mu（mu 中点法）。"""
total = 0
for j in range(m):
mu = (2 * j + 1) * ONE // (2 * m)
total += fmul(e_spec(a2, f0, dirs, n, mu), mu)
return fmul(2 * ONE // m, total)
def e_ms(a2, dirs, n, mu, eavg):
"""Kulla-Conty 多次弹射项的方向反射率：f_ms = (1-E(mu_l))(1-E(mu_v)) / (pi (1-E_avg))。
恒等式：E_ms(mu) 应当恰好等于 1 - E(mu)。E(mu_l) 按 l 的 z 分量查表，同一 z 只算一次。
"""
ev, total, tab = e_spec(a2, ONE, dirs, n, mu), 0, {}
for _, _, lz in dirs:
if lz not in tab:
tab[lz] = e_spec(a2, ONE, dirs, n, lz)
total += ONE - tab[lz]
return fdiv(fmul(fmul(total, cell_of(n)), ONE - ev), fmul(PI, ONE - eavg))
def main():
out = os.path.join(ROOT, "out")
os.makedirs(out, exist_ok=True)
dirs = disk_dirs(GRID)
resid = [(lab, e_spec(a2, ONE, dirs, GRID, ONE)) for lab, a2 in ALPHAS]
assert all(e <= ONE + ONE // 100 for _, e in resid), "E 超过能量上界 1% 以上"
assert abs(resid[-1][1] - (ONE - LN2)) < 328, "alpha = 1 未落到解析值 1 - ln2"
ms = [(lab, e_avg(a2, dirs, GRID)) for lab, a2 in ALPHAS]
assert all(ms[i][1] > ms[i + 1][1] for i in range(len(ms) - 1)), "E_avg 未随粗糙度下降"
assert all(ea <= ONE for _, ea in ms), "E_avg 超过 1"
# 补上多次弹射项后回到 1：同一个 mu 上比，别拿平均量去比方向量
comp = []
for (lab, ea), (_, a2) in zip(ms, ALPHAS):
e1 = e_spec(a2, ONE, dirs, GRID, ONE)
add = e_ms(a2, dirs, GRID, ONE, ea)
comp.append((lab, ea, e1, add, e1 + add))
assert all(abs(t - ONE) < ONE // 100 for *_, t in comp), "多次弹射未把总能量补回 1"
kd = []
for alab, a2 in KD_ALPHAS:
for mlab, mu in MUS:
es = e_spec(a2, F0_KD, dirs, GRID, mu)
df, dc = lambert_kd(F0_KD, dirs, GRID, mu, True), lambert_kd(F0_KD, dirs, GRID, mu, False)
kd.append((alab, mlab, es, df, dc, es + df, es + dc))
assert all(kd[i][5] <= ONE and kd[i][6] <= ONE for i in (0, 4, 8)), "法线入射就该守住"
assert all(kd[i][5] <= kd[i][6] for i in range(len(kd))), "kD = 1-F 应比 1-F0 更紧"
assert kd[3][5] > ONE, "alpha = 0.25 掠射本该超界——这正是要报的缺陷"
cglib.write_csv(os.path.join(out, "cg_shade_06_3_residual.csv"),
["alpha", "E_q16", "residual_q16"], [(lab, e, e - ONE) for lab, e in resid])
cglib.write_csv(os.path.join(out, "cg_shade_06_3_ms.csv"),
["alpha", "E_avg_q16", "deficit_q16", "E_mu1_q16", "E_ms_q16", "total_q16"],
[(lab, ea, ONE - ea, e1, ad, t) for lab, ea, e1, ad, t in comp])
cglib.write_csv(os.path.join(out, "cg_shade_06_3_kd.csv"),
["alpha", "mu_v", "E_spec_q16", "D_kd_from_F_q16", "D_kd_from_F0_q16",
"total_from_F_q16", "total_from_F0_q16"], kd)
size = svg(os.path.join(out, "cg_shade_06_3_energy.svg"),
[("E(mu=1) 单次弹射", "#ff8a5c", ONE,
[(i * (ONE // (len(resid) - 1)), e) for i, (_, e) in enumerate(resid)]),
("E_avg 半球平均", "#5cd6ff", ONE,
[(i * (ONE // (len(ms) - 1)), e) for i, (_, e) in enumerate(ms)])],
"6.3 GGX 的方向反射率与半球平均（F = 1，能量上界 = 1）")
data = {"grid": GRID, "m": M, "directions": len(dirs), "ln2": LN2,
"residual": [{"alpha": lab, "E_q16": e} for lab, e in resid],
"multiple_scatter": [{"alpha": lab, "E_avg_q16": ea, "E_mu1_q16": e1,
"E_ms_q16": ad, "total_q16": t} for lab, ea, e1, ad, t in comp],
"kd": [{"alpha": al, "mu_v": ml, "E_spec_q16": es, "D_from_F_q16": df,
"D_from_F0_q16": dc, "total_from_F_q16": tf, "total_from_F0_q16": tc}
for al, ml, es, df, dc, tf, tc in kd],
"svg_bytes": size}
with open(os.path.join(out, "cg_shade_06_3_energy.json"), "w", encoding="utf-8") as fh:
json.dump(data, fh, ensure_ascii=False, indent=2)
print(f"投影圆盘 {len(dirs)} 个方向，格面积 {cell_of(GRID)}/65536")
print("残差表（mu=1，F=1，上界 1.0000）：")
for lab, e in resid:
print(f" alpha={lab} E={e / ONE:.4f} 残差={((e - ONE) / ONE):+.4f}")
print(f" 解析锚点 alpha=1：1 - ln2 = {(ONE - LN2) / ONE:.6f}，实算 {resid[-1][1] / ONE:.6f}")
print("单次 vs 多次弹射：")
for lab, ea, e1, ad, t in comp:
print(f" alpha={lab} 单次 E_avg={ea / ONE:.4f}（缺 {((ONE - ea) / ONE):.4f}）"
f" E(mu=1)={e1 / ONE:.4f} + 多次 {ad / ONE:.4f} = {t / ONE:.4f}")
print("kD 归一化（F0 = 0.04，albedo = 1，上界 1.0000）：")
for al, ml, es, _, _, tf, tc in kd:
print(f" alpha={al} mu={ml} E_spec={es / ONE:.4f} kD=1-F {tf / ONE:.4f}"
f" kD=1-F0 {tc / ONE:.4f}" + (" <- 超界" if tf > ONE else ""))
print(f"svg : out/cg_shade_06_3_energy.svg ({size} bytes)")
return 0
if __name__ == "__main__":
sys.exit(main())
投影圆盘 3228 个方向，格面积 64/65536
残差表（mu=1，F=1，上界 1.0000）：
alpha=0.05 E=1.0003 残差=+0.0003
alpha=0.10 E=0.9887 残差=-0.0113
alpha=0.20 E=0.9486 残差=-0.0514
alpha=0.30 E=0.8789 残差=-0.1211
alpha=0.50 E=0.6899 残差=-0.3101
alpha=0.75 E=0.4659 残差=-0.5341
alpha=1.00 E=0.3084 残差=-0.6916
解析锚点 alpha=1：1 - ln2 = 0.306854，实算 0.308426
单次 vs 多次弹射：
alpha=0.05 单次 E_avg=0.9932（缺 0.0068） E(mu=1)=1.0003 + 多次 -0.0007 = 0.9995
alpha=0.10 单次 E_avg=0.9752（缺 0.0248） E(mu=1)=0.9887 + 多次 0.0114 = 1.0000
alpha=0.20 单次 E_avg=0.9176（缺 0.0824） E(mu=1)=0.9486 + 多次 0.0515 = 1.0001
alpha=0.30 单次 E_avg=0.8463（缺 0.1537） E(mu=1)=0.8789 + 多次 0.1214 = 1.0002
alpha=0.50 单次 E_avg=0.6924（缺 0.3076） E(mu=1)=0.6899 + 多次 0.3108 = 1.0007
alpha=0.75 单次 E_avg=0.5172（缺 0.4828） E(mu=1)=0.4659 + 多次 0.5357 = 1.0016
alpha=1.00 单次 E_avg=0.3792（缺 0.6208） E(mu=1)=0.3084 + 多次 0.6938 = 1.0022
kD 归一化（F0 = 0.04，albedo = 1，上界 1.0000）：
alpha=0.25 mu=1.00 E_spec=0.0367 kD=1-F 0.9998 kD=1-F0 0.9999
alpha=0.25 mu=0.50 E_spec=0.0564 kD=1-F 1.0147 kD=1-F0 1.0196 <- 超界
alpha=0.25 mu=0.25 E_spec=0.1096 kD=1-F 1.0602 kD=1-F0 1.0729 <- 超界
alpha=0.25 mu=0.10 E_spec=0.1568 kD=1-F 1.0996 kD=1-F0 1.1201 <- 超界
alpha=0.50 mu=1.00 E_spec=0.0276 kD=1-F 0.9908 kD=1-F0 0.9908
alpha=0.50 mu=0.50 E_spec=0.0371 kD=1-F 0.9955 kD=1-F0 1.0004
alpha=0.50 mu=0.25 E_spec=0.0580 kD=1-F 1.0085 kD=1-F0 1.0212 <- 超界
alpha=0.50 mu=0.10 E_spec=0.0780 kD=1-F 1.0208 kD=1-F0 1.0412 <- 超界
alpha=1.00 mu=1.00 E_spec=0.0123 kD=1-F 0.9755 kD=1-F0 0.9756
alpha=1.00 mu=0.50 E_spec=0.0188 kD=1-F 0.9772 kD=1-F0 0.9821
alpha=1.00 mu=0.25 E_spec=0.0270 kD=1-F 0.9776 kD=1-F0 0.9903
alpha=1.00 mu=0.10 E_spec=0.0356 kD=1-F 0.9783 kD=1-F0 0.9988
svg : out/cg_shade_06_3_energy.svg (1026 bytes)
残差表（$F$ 取 1，剥掉菲涅尔的影响）里有两处要解释。$\alpha=0.05$ 读到 1.0003，超了 0.0003——这不是模型越界，而是方向集的网格分辨率：那一档的 NDF 半高角只有 0.76 度（6.2 量过），$64\times64$ 的圆盘格采不到这个峰，6.1 的 $D$ 归一化在同一个 $\alpha$ 上读到 0.9927 是同一个原因的另一面。$\alpha=1$ 的实算 0.308426 与解析值 0.306854 差 0.0016（0.5%），剩下的差全部来自 3228 个方向对圆盘面积的离散近似。
单次 vs 多次弹射那一组是本节的主结果。$\alpha=0.5$ 时半球平均只剩 0.6924，缺 30.8%；补上 Kulla-Conty 的多次弹射项
$$ f_{ms}=\frac{\big(1-E(\mu_l)\big)\big(1-E(\mu_v)\big)}{\pi\big(1-E_{avg}\big)} $$
之后，方向反射率回到 1.0007。这条恒等式 $E_{ms}(\mu)=1-E(\mu)$ 是解析的，数值验证落在 0.2% 以内——“欠能量"不是模型的缺陷，只是漏了一项。
kD = 1 − F 守得住上界吗
漫反射项一般写成 $k_D\cdot\text{albedo}/\pi$，$k_D$ 取多少有两种流派：常数 $k_D=1-F_0$，或方向相关的 $k_D=1-F(v\cdot h)$。直觉上后者更对，因为它把掠射角变大的菲涅尔反射也算进去了。实测结论有四条：
- 法线入射都守得住：三档 $\alpha$ 的合计分别是 0.9998、0.9908、0.9755。
- 掠射角都会破：$\alpha=0.25$、$\mu_v=0.10$ 时 $k_D=1-F$ 读到 1.0996，$k_D=1-F_0$ 读到 1.1201。
- $k_D=1-F$ 始终更紧：12 行里没有一次例外，最大能省下 2 个百分点。
- 越光滑越危险：$\alpha=1$ 全程不超界，$\alpha=0.25$ 从 $\mu_v=0.5$ 就开始超。
原因在 6.1 那个分母。$1-F$ 只扣掉了"这次没反射出去的能量”，但单次弹射的 GGX 在掠射角还被 $G$ 削掉一大块，$k_D$ 没有为 $G$ 留任何余量。所以严格说：$k_D=1-F$ 是更好的近似，不是能量上界的保证。要保住上界，要么给 $k_D$ 再乘一个遮蔽因子，要么走上面那条多次弹射补偿的路。
P 层
// file: lantern/core/src/lantern/cg_shade_06_3_energy.ts
// tier: P
// verify: node --test 'core/test/cg_shade_06_*.test.ts'
// 6.3 能量守恒验证的 TypeScript 侧：方向反射率、漫反射项、半球平均、Kulla-Conty
// 多次弹射补偿。与 sim/cg_shade_06_3_energy.py 逐字对齐。
//
// 测度全部走投影圆盘：面积元 dr ds 已经含了 cos(theta_l) d omega_l，所以
// ∫ f cos(theta) d omega = Σ f * 格面积——积分式里不再出现余弦因子。
// 这一条是本模块最容易写错的地方：多乘一个 cos 会让漫反射项从 0.96 掉到 0.64。
import { ONE, fdiv, fmul, fsqrt } from "./fixed.ts";
import { type Dir, PI, cellOf, schlickF, specBrdf } from "./cg_shade_06_pbr.ts";
/** 方向反射率 E(mu_v) = ∫ f_s cos(theta_l) d omega_l。v = (sqrt(1-mu^2), 0, mu)。 */
export function eSpec(a2: bigint, f0: bigint, dirs: Dir[], n: number, mu: bigint): bigint {
const zv = fsqrt(ONE - fmul(mu, mu));
let total = 0n;
for (const [lx, ly, lz] of dirs) {
const hx = lx + zv;
const hz = lz + mu;
const hl = fsqrt(fmul(hx, hx) + fmul(ly, ly) + fmul(hz, hz));
if (hl === 0n) continue;
const vh = fdiv(fmul(zv, hx) + fmul(mu, hz), hl);
total += specBrdf(a2, lz, mu, fdiv(hz, hl), vh, f0);
}
return fmul(total, cellOf(n));
}
/** 漫反射项 ∫ (kD / pi) cos(theta_l) d omega_l（albedo = 1）；fromF 为真时 kD = 1 - F。 */
export function lambertKd(
f0: bigint, dirs: Dir[], n: number, mu: bigint, fromF: boolean,
): bigint {
const zv = fsqrt(ONE - fmul(mu, mu));
let total = 0n;
for (const [lx, ly, lz] of dirs) {
const hx = lx + zv;
const hz = lz + mu;
const hl = fsqrt(fmul(hx, hx) + fmul(ly, ly) + fmul(hz, hz));
const vh = fdiv(fmul(zv, hx) + fmul(mu, hz), hl);
total += fdiv(ONE - (fromF ? schlickF(f0, vh) : f0), PI);
}
return fmul(total, cellOf(n));
}
/** 半球平均反射率 E_avg = 2 ∫_0^1 E(mu) mu d mu（mu 中点法）。 */
export function eAvg(
a2: bigint, dirs: Dir[], n: number, m: number = 8, f0: bigint = ONE,
): bigint {
let total = 0n;
for (let j = 0; j < m; j += 1) {
const mu = ((2n * BigInt(j) + 1n) * ONE) / BigInt(2 * m);
total += fmul(eSpec(a2, f0, dirs, n, mu), mu);
}
return fmul((2n * ONE) / BigInt(m), total);
}
/** Kulla-Conty 多次弹射项：f_ms = (1-E(mu_l))(1-E(mu_v)) / (pi (1-E_avg))。 */
export function eMs(
a2: bigint, dirs: Dir[], n: number, mu: bigint, eavg: bigint,
): bigint {
const ev = eSpec(a2, ONE, dirs, n, mu);
const tab = new Map<bigint, bigint>();
let total = 0n;
for (const [, , lz] of dirs) {
let e = tab.get(lz);
if (e === undefined) {
e = eSpec(a2, ONE, dirs, n, lz);
tab.set(lz, e);
}
total += ONE - e;
}
return fdiv(fmul(fmul(total, cellOf(n)), ONE - ev), fmul(PI, ONE - eavg));
}
P 层把四个积分函数逐位复现。测试里有一条专门盯着本节最容易写错的那件事：“圆盘度量已经含余弦”——lambertKd 在法线入射的结果必须落在 0.9 以上；一旦有人补了第二次 $\cos$，它会掉到 0.64 附近，测试立刻红。
阅读导航：上一节：6.2 金属度-粗糙度工作流 · 下一节：7.1 线性空间与 sRGB 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。