2.1 BRDF 与能量守恒
站内的 实时 PBR 与阴影贴图 讲的是把 BRDF 塞进 GPU 着色器、再用阴影贴图出帧的工程做法；本节反过来，把 BRDF 当成一个可以在 CPU 上逐项求和的整数函数，用一组固定的半球采样方向，把「半球反射率 ≤ 1」这条物理约束做成可跑的断言。第一章的 Lambert、Phong、逐像素插值都是具体的光照公式，本章要给它们一个共同的度量——BRDF。
从辐亮度到 BRDF
BRDF 描述「一束入射光有多少被反射到某个出射方向」，定义是出射辐亮度与入射辐照度之比：
$$ f_r(\omega_i, \omega_o) = \frac{\mathrm{d}L_o(\omega_o)}{\mathrm{d}E_i(\omega_i)},\quad \mathrm{d}E_i(\omega_i) = L_i(\omega_i)\cos\theta_i,\mathrm{d}\omega_i $$
它的单位是球面度的倒数（sr⁻¹），因为辐照度里已经含了一个立体角。两个性质决定它能不能用：
- 互易性：$f_r(\omega_i, \omega_o) = f_r(\omega_o, \omega_i)$。交换光源与观察者，结果不变。
- 能量守恒：对任意入射方向 $\omega_i$，半球反射率不超过 1：
$$ \int_{\mathcal{H}} f_r(\omega_i, \omega_o),\cos\theta_o,\mathrm{d}\omega_o \le 1 $$
本节验证的就是右边这个积分——它是一个可以离散成整数求和的量。
固定的半球采样集
要在整数上做这个积分，采样方向必须确定、可复现、且立体角已知。方案是 64 个极角带 × 8 个方位：
- 极角在 $\cos\theta$ 上等分成 64 段，取中点 $\cos\theta_k = (2k+1)/128$。这是精确有理数，在 Q16.16 里就是 $(2k+1)\cdot 512$，没有任何舍入。
- $\sin\theta_k = \sqrt{1-\cos^2\theta_k}$ 用
fsqrt的整数平方根求出。 - 方位只用 ${0, \pm 1, \pm \sqrt{2}/2}$ 八个方向的整数常量拼出（45° 间隔），完全不调用浮点三角函数。
- 每个方向带相等立体角 $\Delta\omega = 2\pi/512$。
各向同性归一化：为什么中点法则恰好精确
入射方向取法线时，Lambert 与 Phong（反射向量 = 法线）的被积函数都各向同性——只随 $\theta_o$ 变，8 个方位各给同一个值，方位求和因此精确。而极角方向用的是中点法则，对 $\cos\theta$ 这种线性函数也精确：
$$ \sum_k \cos\theta_k \cdot \Delta\omega = \frac{2\pi}{N_Z}\cdot\frac{N_Z}{2} = \pi = \int_{\mathcal{H}}\cos\theta,\mathrm{d}\omega $$
也就是说，这个 512 方向集对「常数 × cosθ」的积分是逐位精确的。实跑里 $\int\cos\theta,\mathrm{d}\omega$ 的整数估计是 205887，与 $\pi$ 的 Q16.16 表示完全相同，残差 0。
S 层：整数半球积分
下面的脚本搭采样集、跑 Lambert 能量、并把结果落进 out/cg_shade_02_1_energy.json：
#!/usr/bin/env python3
# file: lantern/sim/cg_shade_02_1_brdf_energy.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_shade_02_1_brdf_energy.py
"""2.1 BRDF 与能量守恒（S 层）。
把 BRDF 的**半球反射率** ∫ f_r(ω_i, ω_o)·cosθ_o dω_o 写成**整数估计**：
1. 一个**固定的半球采样方向集**：64 个极角带 × 8 个方位 = 512 个方向；
cosθ_k = (2k+1)/128 是**精确有理数**（Q16.16 里就是 (2k+1)·512），
sinθ_k = fsqrt(1-cos²θ_k)，方位用 {0, ±1, ±√2/2} 四个整数常量拼出，
全程不调用任何浮点三角函数。
2. 每个方向带**相等立体角** Δω = 2π/512（在 (cosθ, φ) 上分层，故 Δω 精确相等）。
3. 估计值 = Σ f_r(ω_o)·cosθ_o·Δω，断言 ≤ ONE。
当入射方向取**法线**、BRDF 各向同性（Lambert、Phong 的反射向量 = 法线）时，
被积函数只随 θ_o 变，8 个方位各给同一个值，方位求和**精确**——这就是
「各向同性归一化」在本采样集上成立的原因：midpoint 法则对 cosθ 是线性的，
Σ_k cosθ_k·Δω = π 逐位精确（见输出里的 cos_integral）。
产物：out/cg_shade_02_1_energy.json。
"""
import json
import os
import sys
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cglib # noqa: E402
ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
ONE = cglib.ONE
TWO_PI = 411775 # round(2π · ONE) = round(6.283185307 · 65536)
INV_PI = 20861 # round(ONE / π) = round(65536 / 3.14159265)
SQRT_HALF = 46341 # round(√2/2 · ONE) = round(0.707106781 · 65536)
N_Z = 64 # 极角带数（cosθ 上的等分）
N_PHI = 8 # 方位数（45° 间隔）
# 8 个方位单位向量，只用 {0, ±ONE, ±SQRT_HALF} 三个整数常量
AZIMUTH = [
(ONE, 0), (SQRT_HALF, SQRT_HALF), (0, ONE), (-SQRT_HALF, SQRT_HALF),
(-ONE, 0), (-SQRT_HALF, -SQRT_HALF), (0, -ONE), (SQRT_HALF, -SQRT_HALF),
]
def build_hemisphere():
"""固定半球采样集：N_Z 个极角带 × N_PHI 个方位，返回 (dx, dy, cosθ) 列表。
cosθ_k = (2k+1)/(2·N_Z) 是第 k 个 cosθ 区间的**中点**——正是让 midpoint
法则对线性被积函数精确的那一步。z 分量直接取 cosθ_k（不用归一化），
故积分里的 cosθ 因子逐位等于该整数。
"""
samples = []
step = ONE // (2 * N_Z) # 1/(2·N_Z) 的 Q16.16（N_Z 为 2 的幂时精确）
for k in range(N_Z):
c = (2 * k + 1) * step # (2k+1)/(2·N_Z) 的 Q16.16
s = cglib.fsqrt(ONE - cglib.fmul(c, c)) # sinθ_k = √(1-cos²θ_k)
for cx, cy in AZIMUTH:
samples.append((cglib.fmul(s, cx), cglib.fmul(s, cy), c))
return samples
def hemisphere_integral(samples, brdf):
"""Σ f_r(ω_o)·cosθ_o·Δω 的整数估计，Δω = 2π/N（N = 采样数）。
Δω 不逐项取整：先求和 Σ f·cosθ，再整体乘 2π 后整除 N，
避免 N 次舍入累积。分子 total ≥ 0、TWO_PI > 0，故 `//` 与截断一致。
"""
total = 0
for dx, dy, cz in samples:
total += cglib.fmul(brdf(dx, dy, cz), cz)
return cglib.fmul(total, TWO_PI) // len(samples)
def lambert(rho):
"""Lambert：f_r = ρ/π，与方向无关。返回 Q16.16 常量。"""
return cglib.fmul(rho, INV_PI)
def main():
samples = build_hemisphere()
# ---- 各向同性归一化的推导：Σ_k cosθ_k 与 Σ cosθ·Δω ----
step = ONE // (2 * N_Z)
cos_sum = sum(cz for _, _, cz in samples) # 8 方位 × Σ_k cosθ_k
band_sum = sum((2 * k + 1) * step for k in range(N_Z)) # Σ_k cosθ_k
cos_integral = hemisphere_integral(samples, lambda dx, dy, cz: ONE) # ∫ cosθ dω 的估计
# ---- Lambert 的能量估计（ρ = 1，理想漫反射）----
f = lambert(ONE)
energy = hemisphere_integral(samples, lambda dx, dy, cz: f)
residual = ONE - energy
record = {
"hemisphere": {
"bands": N_Z, "azimuths": N_PHI, "samples": len(samples),
"cos_theta_k": [(2 * k + 1) * (ONE // (2 * N_Z)) for k in range(N_Z)],
"solid_angle_q16": TWO_PI // (N_Z * N_PHI),
"two_pi_q16": TWO_PI,
},
"isotropic_normalization": {
"sum_cos_theta_bands": band_sum,
"sum_cos_theta_samples": cos_sum,
"cos_integral_q16": cos_integral, # 应 ≈ π = 205887
"pi_q16": 205887,
"identity": "∫_hemi cosθ dω = π",
},
"lambert": {
"inv_pi_q16": INV_PI,
"rho_q16": ONE,
"f_r_q16": f,
"energy_q16": energy,
"energy_over_one_permille": round(energy * 1000.0 / ONE, 4),
"residual_q16": residual,
"assert_le_one": energy <= ONE,
},
}
out = os.path.join(ROOT, "out")
os.makedirs(out, exist_ok=True)
with open(os.path.join(out, "cg_shade_02_1_energy.json"), "w", encoding="utf-8") as fh:
json.dump(record, fh, ensure_ascii=False, indent=2)
print(f"半球采样 : {N_Z} 带 × {N_PHI} 方位 = {len(samples)} 方向")
print(f"cosθ 带中点 (前 4) : {record['hemisphere']['cos_theta_k'][:4]} … "
f"{record['hemisphere']['cos_theta_k'][-1]}")
print(f"Σ_k cosθ_k : {band_sum} (= (N_Z/2)·ONE = {N_Z // 2 * ONE})")
print(f"∫ cosθ dω 的估计 : {cos_integral} (π = 205887, 残差 {cos_integral - 205887})")
print(f"1/π (Q16.16) : {INV_PI}")
print(f"Lambert f_r = ρ/π : {f} (ρ = ONE)")
print(f"能量估计 ∫f·cosθ dω : {energy} ({record['lambert']['energy_over_one_permille']}‰ of ONE)")
print(f"能量守恒残差 : {residual} (Q16.16)")
print(f"断言 estimate ≤ ONE : {record['lambert']['assert_le_one']}")
return 0
if __name__ == "__main__":
sys.exit(main())
实跑输出（/opt/homebrew/bin/python3 sim/cg_shade_02_1_brdf_energy.py）：
半球采样 : 64 带 × 8 方位 = 512 方向
cosθ 带中点 (前 4) : [512, 1536, 2560, 3584] … 65024
Σ_k cosθ_k : 2097152 (= (N_Z/2)·ONE = 2097152)
∫ cosθ dω 的估计 : 205887 (π = 205887, 残差 0)
1/π (Q16.16) : 20861
Lambert f_r = ρ/π : 20861 (ρ = ONE)
能量估计 ∫f·cosθ dω : 65533 (999.9542‰ of ONE)
能量守恒残差 : 3 (Q16.16)
断言 estimate ≤ ONE : True
理想漫反射（$\rho = 1$）的 $f_r = 1/\pi = 20861$（Q16.16），能量估计 65533，是 $\text{ONE}=65536$ 的 999.9542‰，残差只有 3（约 $4.6\times10^{-5}$）。这个残差全部来自 $f_r$ 与 $\cos\theta$ 相乘时的逐项右移取整，不是模型误差。
P 层：跨宿主的同一组整数
TypeScript 侧用同一套整数算术重算采样集与能量，node --test 断言与 S 侧 JSON 逐位相同，并带一个负例（把 $f_r$ 放大 2 倍后能量必须越界）：
// file: lantern/core/test/cg_shade_02_1_brdf_energy.test.ts
// tier: P
// verify: node --test 'core/test/cg_shade_02_1_brdf_energy.test.ts'
// 2.1 的 P 层：TypeScript 侧用同一套整数算术重算半球采样集与 Lambert 能量，
// 断言与 S 侧 sim/cg_shade_02_1_brdf_energy.py 的 out/cg_shade_02_1_energy.json 逐位相同。
//
// 金标准由 S 侧独立产出：
// /opt/homebrew/bin/python3 sim/cg_shade_02_1_brdf_energy.py
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 { buildHemisphere, hemisphereIntegral, lambert } from "../src/lantern/cg_shade_02_material.ts";
const root = join(dirname(fileURLToPath(import.meta.url)), "..", "..");
const g = JSON.parse(
readFileSync(join(root, "out", "cg_shade_02_1_energy.json"), "utf8"),
) as {
hemisphere: { samples: number; cos_theta_k: number[]; two_pi_q16: number };
isotropic_normalization: { sum_cos_theta_bands: number; cos_integral_q16: number; pi_q16: number };
lambert: { inv_pi_q16: number; f_r_q16: number; energy_q16: number; residual_q16: number; assert_le_one: boolean };
};
const samples = buildHemisphere();
test("半球采样集规模与 S 侧一致（64 带 × 8 方位 = 512）", () => {
assert.equal(samples.length, g.hemisphere.samples);
assert.equal(samples.length, 512);
});
test("Σ_k cosθ_k 与 S 侧逐位相同（midpoint 对 cosθ 精确）", () => {
let band = 0n;
for (let k = 0; k < 64; k += 1) band += BigInt(g.hemisphere.cos_theta_k[k]);
assert.equal(band, BigInt(g.isotropic_normalization.sum_cos_theta_bands));
});
test("【核心】∫ cosθ dω 的整数估计 = π（残差 0）", () => {
const cosInt = hemisphereIntegral(samples, () => ONE);
assert.equal(cosInt, BigInt(g.isotropic_normalization.cos_integral_q16));
assert.equal(cosInt, BigInt(g.isotropic_normalization.pi_q16));
assert.equal(cosInt, 205887n);
});
test("【核心】Lambert 能量估计与 S 侧逐位相同，且 ≤ ONE", () => {
const f = lambert(ONE);
assert.equal(f, BigInt(g.lambert.f_r_q16));
assert.equal(f, 20861n);
const energy = hemisphereIntegral(samples, () => f);
assert.equal(energy, BigInt(g.lambert.energy_q16));
assert.equal(energy, 65533n);
assert.ok(energy <= ONE, "能量必须 ≤ ONE");
assert.equal(ONE - energy, BigInt(g.lambert.residual_q16));
});
test("负例：f_r 放大 2 倍后能量必须越界（判据可被拒绝）", () => {
const over = hemisphereIntegral(samples, () => 2n * lambert(ONE));
assert.ok(over > ONE, "2 倍 BRDF 的能量必须 > ONE");
assert.notEqual(over, BigInt(g.lambert.energy_q16));
});
node --test 'core/test/cg_shade_02_1_brdf_energy.test.ts' 的实测输出：
✔ 半球采样集规模与 S 侧一致（64 带 × 8 方位 = 512）
✔ Σ_k cosθ_k 与 S 侧逐位相同（midpoint 对 cosθ 精确）
✔ 【核心】∫ cosθ dω 的整数估计 = π（残差 0）
✔ 【核心】Lambert 能量估计与 S 侧逐位相同，且 ≤ ONE
✔ 负例：f_r 放大 2 倍后能量必须越界（判据可被拒绝）
ℹ tests 5 ℹ pass 5 ℹ fail 0
5 条全绿。最后一条负例是关键：一个对任何改动都无所谓的断言证明不了守恒——把 BRDF 放大 2 倍后能量必须立刻越界，说明这条判据真的对量级敏感。有了这个可断言的积分器，下一节就能把三种常见 BRDF 分别塞进去，逐个量它们的能量残差。
阅读导航：上一节：1.3 着色频率与插值 · 下一节：2.2 常见 BRDF 的定点化 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。