5.3 环境贴图与球谐
5.1 和 5.2 的 AO 都只用了可见性，把环境光当成一个常数。但环境光是有方向分布的：天顶比地平亮，太阳是某个方向上的一个亮点。环境贴图把这个分布整个存下来，问题是它分辨率太高——每个方向一个值，没法直接参与逐像素着色。球谐（spherical harmonics, SH）给出一个可解析重建的压缩：用 9 个系数就能重建出光照方向分布的主要形状。实时 GI 里用 SH 存环境光见 实时全局光照 ，本节做的是它的整数投影与重建。
用基函数展开一个球面函数
把环境光写成一个定义在单位球面上的函数 $L(\omega)$，球谐就是它的一组正交基。只取前 9 项：
$$ L(\omega) \approx \sum_{i=0}^{8} L_i,Y_i(\omega) $$
$Y_i$ 是前 9 个实球谐基：1 个常数项（L0 带）、3 个线性项（L1 带，就是 $x,y,z$ 各乘一个常数）、5 个二次项（L2 带）。系数由投影得到：
$$ L_i = \int_{S^2} L(\omega),Y_i(\omega),d\omega $$
离散化时用一组均匀铺在球面上的方向，每个方向代表 $4\pi/N$ 的立体角，于是积分变成一次加权求和。本节用Fibonacci 格点——它把方向螺旋地铺满球面，分布均匀且完全确定（整数 BAM 三角构造），比随机采样更适合逐位复现。
辐照度卷积：为什么 9 个系数够
着色要的不是辐射亮度本身，而是辐照度 $E(n)$——法线为 $n$ 的表面上收到的总入射能量。它等于 $L$ 与一个余弦核的球面卷积。球谐的漂亮之处在于：余弦核的 SH 展开只有常数项，于是卷积在频域里退化成逐带相乘：
$$ E(n) = \sum_i \hat{A}_i,L_i,Y_i(n),\qquad \hat{A}_0=\pi,;\hat{A}_1=\tfrac{2\pi}{3},;\hat{A}_2=\tfrac{\pi}{4} $$
$\hat A_l$ 随 $l$ 迅速衰减，到 L2 已经只剩 $\pi/4$。所以只留 L0/L1/L2 就够——更高频的带对漫反射辐照度的贡献小到可以忽略。这也是为什么环境光的 SH 压缩在实时渲染里如此常见：9 个浮点（这里 9 个 Q16.16 整数）就能驱动一次全屏的环境光。
实测：投影、重建与误差
环境是一张解析可求值的整数天空：天顶到地平的线性梯度（0.25 到 0.75）加一个 15 度角半径、亮度 +3.0 的太阳圆盘。投影用 1024 条 Fibonacci 方向，误差统计用另一组独立的 2048 条。
#!/usr/bin/env python3
# file: lantern/sim/cg_shade_05_3_sh.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_shade_05_3_sh.py
"""5.3 环境贴图的球谐表示：9 个系数（L0/L1/L2）编码环境光，重建方向光照。
环境是一张**解析可求值**的整数天空：天顶到地平的线性梯度（0.25 -> 0.75）加一个
15 度角半径、亮度 +3.0 的太阳圆盘。沿球面的积分用**固定的 Fibonacci 格点方向**
（整数 BAM 构造）离散化，于是投影与重建全程整数，可跨语言逐位复现：
L_i = (4pi/N) * sum_k env(w_k) * Y_i(w_k) # 投影
E(n) = sum_i A_i * L_i * Y_i(n) # 重建（A = 辐照度卷积系数）
关键事实：线性梯度 z 恰好是 L1 带能**精确**表示的分量，故留到 L1 就把背景重建得很好；
太阳圆盘是高频的，9 个系数装不下，L2 带被圆盘拉偏，反而在背景上引入振铃——RMS 单调
下降但逐点 MAE 反升。这就是 SH 环境光的适用边界。
产物：config/cg_shade_05_3_sh.json、out/cg_shade_05_3_sh.json、
out/cg_shade_05_3_sh_recon.png / .rgb、out/cg_shade_05_3_sh_error.csv。
"""
import hashlib
import json
import math
import os
import sys
import numpy as np
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cglib # noqa: E402
from cg_raster_04_trig import cos_bam, sin_bam # noqa: E402
ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
ONE = cglib.ONE
fmul = cglib.fmul
# 实球谐基的 Q16.16 常数（本脚本算出并写进 config，TS 侧从 config 读取以保证同源）。
Y00 = int(round(0.28209479177387814 * ONE))
Y1 = int(round(0.4886025119029199 * ONE))
Y2A = int(round(1.0925484305920792 * ONE))
Y20 = int(round(0.31539156525252005 * ONE))
Y22 = int(round(0.5462742152960396 * ONE))
BASIS = {"y00": Y00, "y1": Y1, "y2a": Y2A, "y20": Y20, "y22": Y22}
# 辐照度卷积系数 A_l（Q16.16）：A0=pi, A1=2pi/3, A2=pi/4。
A0 = int(round(math.pi * ONE))
A1 = int(round(2.0 * math.pi / 3.0 * ONE))
A2 = int(round(math.pi / 4.0 * ONE))
C4PI = int(round(4.0 * math.pi * ONE)) # 4pi，Q16.16
SUN = (0, 39322, 52429) # 太阳方向（Q16.16 单位向量）
COS_DISK = int(round(math.cos(math.radians(15.0)) * ONE))
def fib_dirs(n):
"""Fibonacci 球面格点：整数 BAM 构造，逐位确定。"""
out = []
step_bam = int(round(ONE / ((1.0 + math.sqrt(5.0)) / 2.0))) # 2pi/golden 折算到 BAM
for i in range(n):
z = ONE - ((2 * i + 1) << 16) // n
r = cglib.fsqrt(max(0, ONE - fmul(z, z)))
phi = (i * step_bam) % ONE
out.append((fmul(r, cos_bam(phi)), fmul(r, sin_bam(phi)), z))
return tuple(out)
def env_radiance(d):
"""环境辐射亮度（Q16.16，HDR 上限 3.75）：梯度 0.25->0.75 + 太阳圆盘 +3.0。"""
x, y, z = d
t = (z + ONE) >> 1 # 0（天底）-> ONE（天顶）
base = 16384 + ((32768 * t) >> 16) # 0.25 + 0.5*t
if fmul(x, SUN[0]) + fmul(y, SUN[1]) + fmul(z, SUN[2]) >= COS_DISK:
base += 3 * ONE
return base
def basis(d):
"""9 个实球谐基函数在方向 d 上的取值（Q16.16），顺序 [L0,L1,L2]。"""
x, y, z = d
return (
Y00,
fmul(Y1, y), fmul(Y1, z), fmul(Y1, x),
fmul(Y2A, fmul(x, y)), fmul(Y2A, fmul(y, z)),
fmul(Y20, (3 * fmul(z, z) - ONE)),
fmul(Y2A, fmul(x, z)), fmul(Y22, (fmul(x, x) - fmul(y, y))),
)
def project(dirs):
"""投影：L_i = (4pi/N) * sum_k env_k * Y_i,k（全整数）。"""
n = len(dirs)
acc = [0] * 9
for d in dirs:
e = env_radiance(d)
for i, yv in enumerate(basis(d)):
acc[i] += fmul(e, yv)
return [(C4PI * s) // (n * ONE) for s in acc]
def reconstruct(L, d, irradiance=False):
"""重建：E(n) = sum_i A_i * L_i * Y_i(n)（辐照度）或 sum_i L_i*Y_i(n)（辐射亮度）。"""
coef = (A0, A1, A1, A1, A2, A2, A2, A2, A2) if irradiance else (ONE,) * 9
tot = 0
for i, yv in enumerate(basis(d)):
tot += fmul(fmul(coef[i], L[i]), yv)
return tot
def brute_irradiance(dirs, nvec):
"""直接积分的地面真值辐照度 E(n) = (4pi/N) * sum_k env_k * max(0, n.d_k)。"""
acc = 0
for d in dirs:
dot = fmul(d[0], nvec[0]) + fmul(d[1], nvec[1]) + fmul(d[2], nvec[2])
if dot > 0:
acc += fmul(env_radiance(d), dot)
return (C4PI * acc) // (len(dirs) * ONE)
def is_disk(d):
return fmul(d[0], SUN[0]) + fmul(d[1], SUN[1]) + fmul(d[2], SUN[2]) >= COS_DISK
def main():
dirs = fib_dirs(1024)
L = project(dirs)
ev = fib_dirs(2048)
low = [d for d in ev if not is_disk(d)]
disk = [d for d in ev if is_disk(d)]
bands = {"L0": [L[0]] + [0] * 8, "L0-L1": L[0:4] + [0] * 5, "L0-L2": L}
def stats(c, subset):
e = [reconstruct(c, d) - env_radiance(d) for d in subset]
return math.isqrt(sum(v * v for v in e) // len(e)), sum(abs(v) for v in e) // len(e)
rms, mae, lowmae = {}, {}, {}
for k, c in bands.items():
rms[k], mae[k] = stats(c, ev)
lowmae[k] = stats(c, low)[1]
mae_disk = stats(L, disk)[1]
# 辐照度：SH 重建 vs 直接积分（法线朝天顶），再乘 AO 得到最终环境项
up = (0, 0, ONE)
e_sh, e_true = reconstruct(L, up, irradiance=True), brute_irradiance(dirs, up)
inv_pi, ao, albedo8 = int(round(ONE / math.pi)), 36416, 72 # AO 取自 5.1 的 h=1.5r 构型
albedo_q16 = (albedo8 * ONE) // 255
lit_sh = fmul(fmul(fmul(e_sh, inv_pi), albedo_q16), ao)
lit_gt = fmul(fmul(fmul(e_true, inv_pi), albedo_q16), ao)
# 重建贴图：256 条 Fibonacci 方向按序铺成 16x16 缩略图
thumb_dirs = ev[:256]
rec8 = np.array([min(255, max(0, (reconstruct(L, d) * 255) >> 16)) for d in thumb_dirs],
dtype=np.uint8).reshape(16, 16)
cglib.save_png(os.path.join(ROOT, "out", "cg_shade_05_3_sh_recon.png"),
np.stack([rec8] * 3, axis=-1))
with open(os.path.join(ROOT, "out", "cg_shade_05_3_sh_recon.rgb"), "wb") as fh:
fh.write(rec8.tobytes())
cglib.write_csv(os.path.join(ROOT, "out", "cg_shade_05_3_sh_error.csv"),
["band", "rmsQ16", "maeQ16", "lowMaeQ16"],
[(k, rms[k], mae[k], lowmae[k]) for k in bands])
cfg = {"name": "sh-sky", "dirs": len(dirs), "evalDirs": len(ev), "sunQ16": list(SUN),
"cosDiskQ16": COS_DISK, "basis": BASIS, "A0": A0, "A1": A1, "A2": A2,
"c4pi": C4PI, "coefQ16": L, "thumbDirsQ16": [list(d) for d in thumb_dirs]}
with open(os.path.join(ROOT, "config", "cg_shade_05_3_sh.json"), "w", encoding="utf-8") as fh:
json.dump(cfg, fh, ensure_ascii=False, indent=2)
stats_out = {"coefQ16": L, "rmsQ16": [rms[k] for k in bands],
"maeQ16": [mae[k] for k in bands], "lowMaeQ16": [lowmae[k] for k in bands],
"maeDiskQ16": mae_disk, "lowDirs": len(low), "diskDirs": len(disk),
"irradianceUpShQ16": e_sh, "irradianceUpTrueQ16": e_true,
"irradianceRelErrPermille": abs(e_sh - e_true) * 1000 // max(1, e_true),
"aoQ16": ao, "albedoQ8": albedo8, "albedoQ16": albedo_q16,
"litShQ16": lit_sh, "litGroundTruthQ16": lit_gt,
"litShQ8": min(255, (lit_sh * 255) >> 16),
"reconSha256": hashlib.sha256(rec8.tobytes()).hexdigest()}
with open(os.path.join(ROOT, "out", "cg_shade_05_3_sh.json"), "w", encoding="utf-8") as fh:
json.dump(stats_out, fh, ensure_ascii=False, indent=2)
print(f"方向集 : 投影 {len(dirs)} / 评估 {len(ev)} 条 Fibonacci 格点")
for i, c in enumerate(L):
print(f" L{i} = {c:>8d} ({c / ONE:+.6f})")
print(f"重建 RMS 全向 : L0 {rms['L0'] / ONE:.6f} L1 {rms['L0-L1'] / ONE:.6f} "
f"L2 {rms['L0-L2'] / ONE:.6f}")
print(f"低频背景 MAE ({len(low)} 向): L0 {lowmae['L0'] / ONE:.6f} L1 {lowmae['L0-L1'] / ONE:.6f} "
f"L2 {lowmae['L0-L2'] / ONE:.6f} （L2 振铃反升）")
print(f"太阳圆盘 MAE ({len(disk)} 向) : {mae_disk / ONE:.6f}")
print(f"天顶辐照度 SH / 真值 : {e_sh} / {e_true} = {e_sh / ONE:.6f} / {e_true / ONE:.6f} "
f"({stats_out['irradianceRelErrPermille']}‰)")
print(f"SH 光照 * AO : {lit_sh} / 65536 = {lit_sh / ONE:.6f} "
f"(真值 {lit_gt / ONE:.6f}) 8 位 {stats_out['litShQ8']}")
print(f"重建贴图 sha256 : {stats_out['reconSha256']}")
return 0
if __name__ == "__main__":
sys.exit(main())
方向集 : 投影 1024 / 评估 2048 条 Fibonacci 格点
L0 = 127721 (+1.948868)
L1 = 11998 (+0.183075)
L2 = 49150 (+0.749969)
L3 = 283 (+0.004318)
L4 = 430 (+0.006561)
L5 = 20614 (+0.314545)
L6 = 10859 (+0.165695)
L7 = 453 (+0.006912)
L8 = -7849 (-0.119766)
重建 RMS 全向 : L0 0.438171 L1 0.378906 L2 0.362946
低频背景 MAE (2013 向): L0 0.129684 L1 0.079880 L2 0.098526 （L2 振铃反升）
太阳圆盘 MAE (35 向) : 2.581299
天顶辐照度 SH / 真值 : 168862 / 169240 = 2.576630 / 2.582397 (2‰)
SH 光照 * AO : 8432 / 65536 = 0.128662 (真值 0.128952) 8 位 32
重建贴图 sha256 : 4fae2b1db4606653143cda4dd5c2cf88e8c3d29332ce0940cf8f41b5608bd529
频域边界：梯度装得下，太阳装不下
数字里最值得读的是低频背景的 MAE：L0 是 0.129684，加到 L1 降到 0.079880，再加到 L2 反而升回 0.098526。这不是 bug，是一个干净的教学点。
原因是两件事同时发生。其一，线性梯度 $z$ 恰好是 L1 带能精确表示的分量，所以留到 L1 就把背景重建得很好了。其二，太阳圆盘是高频的，9 个系数装不下它，投影会把 L2 带"拉偏"去凑那个亮点，结果在平滑背景上引入振荡——这正是截断傅里叶级数在跳变处的吉布斯振铃。
但注意 RMS 是单调下降的：0.438171 → 0.378906 → 0.362946。球谐投影是最小二乘意义下的最优拟合，它最小化的是能量（平方误差），不是逐点误差。所以"加一个带让 RMS 变小、却让逐点 MAE 变大"完全不矛盾——它说明 SH 是一个能量最优、逐点可能更差的表示。
太阳圆盘方向上的 MAE 高达 2.581299，是背景误差的 26 倍。要重建一个尖锐的高光，9 个系数远远不够，得上贴图或更高阶的基。
P 层：9 个系数的跨语言复现
// file: lantern/core/src/lantern/cg_shade_05_sh.ts
// tier: P
// verify: node --test 'core/test/cg_shade_05_3_sh.test.ts'
// 5.3 的 P 侧：9 个球谐系数对方向光照的重建，逐字节复现 sim/cg_shade_05_3_sh.py。
// 量级：|L_i| <= 2^18，|Y_i| <= 2^17，累加 9 项 <= 2^39 << 2^53 —— 可用 Number。
// 舍入：一律 Math.floor(a*b/65536)（= Python 的 `>>16`）；禁用 `>>`（转 int32 会截断）。
const ONE = 65536;
const fmul = (a: number, b: number): number => Math.floor((a * b) / ONE);
export type Vec3i = readonly [number, number, number];
export type ShConfig = {
basis: { y00: number; y1: number; y2a: number; y20: number; y22: number };
A0: number;
A1: number;
A2: number;
coefQ16: number[];
};
/** 9 个实球谐基函数在方向 d 上的取值（Q16.16），顺序 [L0,L1,L2]。 */
export function shBasis(d: Vec3i, b: ShConfig["basis"]): number[] {
const [x, y, z] = d;
return [
b.y00,
fmul(b.y1, y), fmul(b.y1, z), fmul(b.y1, x),
fmul(b.y2a, fmul(x, y)), fmul(b.y2a, fmul(y, z)),
fmul(b.y20, 3 * fmul(z, z) - ONE),
fmul(b.y2a, fmul(x, z)), fmul(b.y22, fmul(x, x) - fmul(y, y)),
];
}
/** 重建：辐照度 E(n) = sum A_i L_i Y_i(n)，或辐射亮度 sum L_i Y_i(n)（Q16.16）。 */
export function shReconstruct(cfg: ShConfig, d: Vec3i, irradiance = false): number {
const coef = irradiance
? [cfg.A0, cfg.A1, cfg.A1, cfg.A1, cfg.A2, cfg.A2, cfg.A2, cfg.A2, cfg.A2]
: new Array<number>(9).fill(ONE);
const basis = shBasis(d, cfg.basis);
let tot = 0;
for (let i = 0; i < 9; i += 1) tot += fmul(fmul(coef[i], cfg.coefQ16[i]), basis[i]);
return tot;
}
/** 重建 8 位贴图：与 Python 侧 `min(255, max(0, (L*255)>>16))` 同口径。 */
export function shThumb(cfg: ShConfig, dirs: Vec3i[]): Uint8Array {
const out = new Uint8Array(dirs.length);
for (let i = 0; i < dirs.length; i += 1) {
out[i] = Math.min(255, Math.max(0, Math.floor((shReconstruct(cfg, dirs[i]) * 255) / ONE)));
}
return out;
}
TypeScript 侧读同一份系数与基函数常数（都由 Python 写进 config/），对同一批方向做重建，产出同一张 8 位重建贴图，sha256 必须逐字节相同。断言里还有一条专门盯着频域行为：RMS 必须单调下降、低频 MAE 必须在 L2 反升——把上面那个教学点也钉成了回归测试。
SH 光照与 AO 结合
最后一步是把 SH 环境光和 5.1 的 AO 合起来。SH 给出法线朝上的辐照度 $E(n)$，除以 $\pi$ 得到漫反射的入射亮度，再乘 albedo 和 AO，即 $L_{ambient}=\frac{E(n)}{\pi}\cdot\text{albedo}\cdot AO$。用 5.1 里 h=1.5r 的 AO 值 36416，实测 $L_{ambient}=0.128662$，与用直接积分算的真值 0.128952 只差 0.2%。两个近似——SH 的频域截断、AO 的几何遮蔽——叠加起来仍然落在千分之几，这就是 SH 环境光加 AO 这条组合拳能成为实时管线标配的原因：它把两个高维积分各自压成一个可预算的近似，误差还能互相独立地度量。
阅读导航：上一节：5.2 SSAO 的采样与噪声 · 下一节：6.1 微表面理论 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。