5.1 环境光遮蔽
环境光遮蔽（ambient occlusion, AO）回答一个非常具体的问题：一个着色点周围，有多大比例的半球方向能"看见"天空。它把环境光这个本来要沿整个半球积分的高维量，压成一个标量乘子——这正是它便宜的原因，也是它近似的边界。实时 GI 的工程做法见 实时全局光照 ，本节只做一件事：把 AO 的定义写成一个可解析求解、可逐像素数值估计、且两边能互相验证的定点算式。
定义：半球可见性积分
设着色点为 $P$、法线为 $n$，$\Omega$ 是 $n$ 周围的所有方向。AO 定义为
$$ AO(P) = \frac{1}{\pi}\int_{\Omega} V(P,\omega),(n\cdot\omega),d\omega $$
其中 $V$ 是可见性函数：方向 $\omega$ 没被挡住取 1，被挡住取 0。分母 $\pi$ 是余弦加权的归一化因子，于是全可见时 $AO=1$，全遮挡时 $AO=0$。注意它是一个纯几何量，与光源无关——这一点决定了它可以被预计算、被复用。
一个能解析求解的最小场景
任意场景里这个积分都没有闭式解，但有一个构型例外：无限大的平面（法线 $+z$），正上方 $h$ 处悬一个半径 $r$ 的球，球心落在平面的垂足上。问平面正下方那一点 $P$ 的 AO。
球从 $P$ 看去张出一个半角 $\alpha$ 的圆锥，且 $\sin\alpha = r/h$。圆锥内的方向被球挡住，圆锥外全可见。余弦加权的可见比例恰好是
$$ AO = 1 - \sin^2\alpha = 1 - (r/h)^2 $$
这个式子只在球心正上方那一点严格成立（离开垂足后，被挡住的区域不再是完整圆锥），但它给了我们一个精确的锚点：任何数值估计器都必须在这个点上逼近它。脚本把它算成一条扫描线，h/r 取 1.2 到 4.0。
#!/usr/bin/env python3
# file: lantern/sim/cg_shade_05_1_ao.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_shade_05_1_ao.py
"""5.1 环境光遮蔽：半球可见性积分的解析解与定点数值估计。
最小构型：地面 z=0（法线 +z），半径 r 的球心悬在正上方 h 处。对地面正下方那一点，
球在法线周围张出半角 alpha 的锥（sin alpha = r/h），余弦加权的遮蔽占比恰为
sin^2(alpha)，于是解析真值 **AO = 1 - (r/h)^2**。
数值估计用一组**固定的整数方向集**（64 条余弦加权方向，Q16.16，见 DIRS_Q16）逐方向
做射线-球求交，统计未被遮蔽的比例。误差来源：方向集的离散量化（估计值只能是 1/N 的
整数倍）与方向分布偏差。
关键化简：地面点在球外（c = |P-C|^2 - r^2 > 0），故「射线命中球」等价于
`b < 0 且 b*b >= c`（b = d·(P-C)），**无需开方**，全整数、可跨语言逐位复现。
产物：config/cg_shade_05_1_ao.json、out/cg_shade_05_1_ao.json、
out/cg_shade_05_1_ao_map.png / .rgb、out/cg_shade_05_1_ao_shaded.png、
out/cg_shade_05_1_ao_scan.csv。
"""
import hashlib
import json
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
# 固定整数方向集：64 条余弦加权半球方向（极轴 = +z），Q16.16，单位化后取整。
# 由 Hammersley 序列 + 余弦半球映射确定性生成，**逐字内联**故跨语言同源。
DIRS_Q16 = (
(5793, 0, 65279), (-10033, 0, 64763), (0, 12953, 64243), (0, -15326, 63719),
(12288, 12288, 63190), (-13585, -13585, 62657), (-14768, 14768, 62119),
(15864, -15864, 61576), (22066, 9140, 61029), (-23327, -9663, 60477),
(-10158, 24524, 59919), (10631, -25666, 59357), (11084, 26758, 58789),
(-11519, -27808, 58215), (-28820, 11938, 57636), (29797, -12342, 57051),
(32637, 6492, 56459), (-33611, -6686, 55862), (-6874, 34558, 55258),
(7057, -35480, 54647), (20607, 30840, 54030), (-21103, -31583, 53405),
(-32309, 21588, 52773), (33019, -22063, 52134), (33715, 22527, 51486),
(-34396, -22983, 50830), (-23429, 35064, 50166), (23867, -35719, 49492),
(8532, 42893, 48809), (-8680, -43639, 48117), (-44372, 8826, 47415),
(45094, -8970, 46702), (46477, 4578, 45977), (-47186, -4647, 45242),
(-4716, 47885, 44494), (4784, -48574, 43733), (31397, 38258, 42959),
(-31825, -38778, 42171), (-39292, 32246, 41368), (39799, -32662, 40548),
(45978, 24576, 39712), (-46542, -24877, 38858), (-25175, 47099, 37985),
(25470, -47650, 37091), (15863, 52294, 36175), (-16041, -52879, 35235),
(-53457, 16216, 34270), (54028, -16389, 33276), (54594, 16561, 32252),
(-55154, -16731, 31194), (-16899, 55708, 30099), (17065, -56257, 28963),
(27981, 52348, 27780), (-28246, -52844, 26545), (-53336, 28509, 25249),
(53823, -28769, 23884), (47599, 39064, 22435), (-48019, -39408, 20886),
(-39749, 48434, 19212), (40087, -48847, 17378), (6246, 63412, 15326),
(-6297, -63934, 12953), (-64452, 6348, 10033), (64965, -6399, 5793),
)
def reference_dirs(side=32):
"""参考方向集：side×side 分层余弦半球方向，全整数构造（用整数 BAM 三角函数）。
与 DIRS_Q16 的构造**不同源**（分层网格 vs Hammersley 序列），故可作为独立参考。
"""
out = []
for i in range(side):
u1 = ((2 * i + 1) << 16) // (2 * side)
r = cglib.fsqrt(u1)
z = cglib.fsqrt(ONE - u1)
for j in range(side):
phi = ((2 * j + 1) << 16) // (2 * side)
out.append((cglib.fmul(r, cos_bam(phi)), cglib.fmul(r, sin_bam(phi)), z))
return tuple(out)
def occlusion_mask(xs, ys, cz, r, dirs):
"""对 (xs, ys, 0) 上一批地面点逐方向返回遮蔽布尔矩阵，形状 (n_dirs, G, G)。"""
o0, o1, o2 = xs, ys, -cz
c = (((o0 * o0) >> 16) + ((o1 * o1) >> 16) + ((o2 * o2) >> 16)) - ((r * r) >> 16)
occ = np.zeros((len(dirs),) + xs.shape, dtype=bool)
for k, (dx, dy, dz) in enumerate(dirs):
b = ((dx * o0) >> 16) + ((dy * o1) >> 16) + ((dz * o2) >> 16)
occ[k] = (b < 0) & (((b * b) >> 16) >= c)
return occ
def ao_from_mask(occ, n_dirs):
"""未被遮蔽比例（Q16.16）：((N - 遮蔽数) << 16) // N，分子非负故 floor 与截断一致。"""
vis = n_dirs - occ[:n_dirs].sum(axis=0)
return (vis.astype(np.int64) << 16) // n_dirs
def main():
r = ONE
grid, extent = 64, 4 * ONE
step = (2 * extent) // grid
coord = np.array([-extent + step * i for i in range(grid)], dtype=np.int64)
xs, ys = np.meshgrid(coord, coord)
ref_dirs = reference_dirs(64) # 4096 条独立参考方向
# --- 解析扫描：正下方点的 AO = 1 - (r/h)^2 ---
scan = []
origin = np.zeros((1, 1), dtype=np.int64)
for hr_num, hr_den in [(12, 10), (15, 10), (20, 10), (30, 10), (40, 10)]:
h = r * hr_num // hr_den
q = (r << 16) // h # r/h in Q16.16
analytic = ONE - ((q * q) >> 16) # 1 - (r/h)^2
est64 = int(ao_from_mask(occlusion_mask(origin, origin, h, r, DIRS_Q16), 64)[0, 0])
est1k = int(ao_from_mask(occlusion_mask(origin, origin, h, r, ref_dirs), len(ref_dirs))[0, 0])
scan.append((hr_num * 10 // hr_den, analytic, est64, est1k))
# --- 二维 AO 图（h = 1.6 r）---
h = r * 16 // 10
ao64 = ao_from_mask(occlusion_mask(xs, ys, h, r, DIRS_Q16), 64)
ao_ref = ao_from_mask(occlusion_mask(xs, ys, h, r, ref_dirs), len(ref_dirs))
abs_err = int(np.abs(ao64 - ao_ref).sum())
mae_q16 = abs_err // ao64.size
mae_permille = abs_err * 1000 // (ao64.size * ONE)
# --- AO 作为乘性项接入着色：ambient * AO + direct * max(0, n·L) ---
albedo, ambient, direct = (72, 82, 66), 13107, 58982
lz = 37837 # 方向光 z 分量，n = +z 故 n·L = lz
k = np.clip((ambient * ao64) >> 16, 0, ONE)
k = np.clip(k + ((direct * lz) >> 16), 0, ONE)
shade = np.stack([(albedo[i] * k) >> 16 for i in range(3)], axis=-1).astype(np.int64)
out = os.path.join(ROOT, "out")
cfg = {"name": "ao-plane-sphere", "radiusQ16": r, "heightQ16": h, "grid": grid,
"extentQ16": extent, "dirs": len(DIRS_Q16), "referenceDirs": len(ref_dirs),
"referenceSide": 64, "ambientQ16": ambient, "directQ16": direct,
"albedoQ8": list(albedo), "dirsQ16": [list(d) for d in DIRS_Q16]}
with open(os.path.join(ROOT, "config", "cg_shade_05_1_ao.json"), "w", encoding="utf-8") as fh:
json.dump(cfg, fh, ensure_ascii=False, indent=2)
ao8 = np.clip(ao64 * 255 >> 16, 0, 255).astype(np.uint8)
cglib.save_png(os.path.join(out, "cg_shade_05_1_ao_map.png"), np.stack([ao8] * 3, axis=-1))
with open(os.path.join(out, "cg_shade_05_1_ao_map.rgb"), "wb") as fh:
fh.write(ao8.tobytes())
cglib.save_png(os.path.join(out, "cg_shade_05_1_ao_shaded.png"), shade)
cglib.write_csv(os.path.join(out, "cg_shade_05_1_ao_scan.csv"),
["h_over_r_x10", "analyticQ16", "est64Q16", "est1024Q16"], scan)
stats = {"analyticCenterQ16": scan[1][1], "est64CenterQ16": scan[1][2],
"est1024CenterQ16": scan[1][3], "scan": scan,
"aoMapMaeQ16": mae_q16, "aoMapMaePermille": mae_permille,
"aoMapMinQ16": int(ao64.min()), "aoMapMaxQ16": int(ao64.max()),
"aoMapCenterQ16": int(ao64[grid // 2, grid // 2]),
"aoMapCenterQ8": int(ao8[grid // 2, grid // 2]),
"aoMapSha256": hashlib.sha256(ao8.tobytes()).hexdigest(),
"dirs": len(DIRS_Q16), "referenceDirs": len(ref_dirs)}
with open(os.path.join(out, "cg_shade_05_1_ao.json"), "w", encoding="utf-8") as fh:
json.dump(stats, fh, ensure_ascii=False, indent=2)
print(f"方向集 : 固定 {len(DIRS_Q16)} 条 / 参考 {len(ref_dirs)} 条")
for hr, an, e64, e1k in scan:
print(f" h/r={hr / 10:.1f} 解析 {an / ONE:.6f} 64 向 {e64 / ONE:.6f} "
f"4096 向 {e1k / ONE:.6f} 64 向误差 {(e64 - an) / ONE:+.6f}")
print(f"AO 图 64 向 vs 4096 向: MAE {mae_q16} / 65536 = {mae_q16 / ONE:.6f} "
f"({mae_permille}‰)")
print(f"AO 值域 : [{ao64.min() / ONE:.4f}, {ao64.max() / ONE:.4f}]")
print(f"AO 图 sha256 : {stats['aoMapSha256']}")
return 0
if __name__ == "__main__":
sys.exit(main())
方向集 : 固定 64 条 / 参考 4096 条
h/r=1.2 解析 0.305573 64 向 0.312500 4096 向 0.312500 64 向误差 +0.006927
h/r=1.5 解析 0.555573 64 向 0.562500 4096 向 0.562500 64 向误差 +0.006927
h/r=2.0 解析 0.750000 64 向 0.750000 4096 向 0.750000 64 向误差 +0.000000
h/r=3.0 解析 0.888901 64 向 0.890625 4096 向 0.890625 64 向误差 +0.001724
h/r=4.0 解析 0.937500 64 向 0.937500 4096 向 0.937500 64 向误差 +0.000000
AO 图 64 向 vs 4096 向: MAE 543 / 65536 = 0.008286 (8‰)
AO 值域 : [0.5938, 1.0000]
AO 图 sha256 : b23b26e8821ec5b62ae27ebb23f62c0292071158843f50f9137579cf5984e887
数值估计：固定方向集
数值做法是蒙特卡洛式的：撒 $N$ 条方向，数有多少条被挡住。本节刻意不用随机数，而是用一组固定的整数方向——64 条余弦加权半球方向，Q16.16，逐字内联在脚本的 DIRS_Q16 里。固定方向集有三个好处：跨语言逐位可复现、误差里不含随机成分、收敛曲线干净。
每条方向做一次射线-球求交。这里有一个关键化简：地面点在球外，记 $o = P-C$、$c = |o|^2 - r^2 > 0$、$b = d\cdot o$，则"射线命中球"等价于
$$ b < 0 \quad\text{且}\quad b^2 \ge c $$
不需要开方，全程整数，跨语言可逐位复现。
误差有两个来源。其一是离散量化：$N$ 条方向的估计值只能是 $1/N$ 的整数倍，所以 h/r=1.5 时真值 0.555573 只能落到 0.562500（36/64），误差 +0.006927；h/r=2.0 恰好整除，误差为 0。其二是方向分布的偏差：固定方向集在极点附近密、赤道附近疏，于是整幅 AO 图对 4096 向参考的 MAE 是 543/65536 ≈ 8‰。
P 层：同一张 AO 图的逐字节复现
// file: lantern/core/src/lantern/cg_shade_05_ao.ts
// tier: P
// verify: node --test 'core/test/cg_shade_05_1_ao.test.ts'
// 5.1 的 P 侧：环境光遮蔽的定点估计，逐字节复现 Python 侧 sim/cg_shade_05_1_ao.py 的 AO 图。
//
// 逐像素层量级：地面坐标 |o| <= 4*2^16 = 2^18，其平方 <= 2^36，三项累加 <= 2^37.6 < 2^53，
// 故本文件全程用 Number 即可（不需要 BigInt）。但**不能用 `>> 16`**——JS 的位移先做
// ToInt32，2^36 量级会被静默截断；一律用 `Math.floor(x / 65536)`，对负数也是 floor，
// 与 Python 的算术右移 `>> 16` 同义。方向集由 config JSON 提供，与 Python 侧同源。
export type Vec3i = readonly [number, number, number];
const ONE = 65536;
/** 地面点 o 沿方向 d 是否命中球心 (0,0,cz)、半径 r 的球（b<0 且 b^2>=c，无需开方）。 */
export function occluded(ox: number, oy: number, oz: number, d: Vec3i, c: number): boolean {
const b = Math.floor((d[0] * ox) / ONE) + Math.floor((d[1] * oy) / ONE) + Math.floor((d[2] * oz) / ONE);
return b < 0 && Math.floor((b * b) / ONE) >= c;
}
/** grid×grid 的地面 AO 图（8 位），与 Python 的 ao_from_mask + 量化对齐。 */
export function aoMap(
radiusQ16: number, heightQ16: number, grid: number, extentQ16: number, dirs: Vec3i[],
): Uint8Array {
const step = Math.floor((2 * extentQ16) / grid);
const coord = Array.from({ length: grid }, (_, i) => -extentQ16 + step * i);
const out = new Uint8Array(grid * grid);
const r2 = Math.floor((radiusQ16 * radiusQ16) / ONE);
for (let row = 0; row < grid; row += 1) {
const oy = coord[row];
for (let col = 0; col < grid; col += 1) {
const ox = coord[col];
const oz = -heightQ16;
const c = Math.floor((ox * ox) / ONE) + Math.floor((oy * oy) / ONE) + Math.floor((oz * oz) / ONE) - r2;
let blocked = 0;
for (const d of dirs) if (occluded(ox, oy, oz, d, c)) blocked += 1;
const ao = Math.floor(((dirs.length - blocked) * ONE) / dirs.length); // 分子非负
out[row * grid + col] = Math.min(255, Math.floor((ao * 255) / ONE));
}
}
return out;
}
TypeScript 侧读同一份方向表（由 Python 写进 config/），用同一套整数判定，产出同一张 AO 图。四条断言全绿，其中最关键的一条是 sha256 比对——一旦运算顺序或舍入方式不同，哈希立刻不同，这是可证伪的。负例把遮蔽判定写成 b > 0，哈希随即改变，证明这个判据确实有区分力。
把 AO 接进着色
AO 的用法是乘性的：把环境项乘以 AO，直射项不动。
$$ L = k_a\cdot AO\cdot \text{albedo} + k_d\cdot \text{albedo}\cdot \max(0, n\cdot L) $$
脚本里以 h=1.6r 的 AO 图为例算了一遍着色（环境强度 0.2、方向光强度 0.9），输出 out/cg_shade_05_1_ao_shaded.png。乘性接入的代价是：AO 只用了可见性的一阶信息——它知道"挡住了多少"，不知道"挡住的东西有多亮"。
要更准，就得把入射光的方向分布也带上，那正是 5.3 的球谐环境光；而如果场景里没有解析解，可见性本身也得靠屏幕空间的深度缓冲去估，那就是 5.2 的 SSAO。
阅读导航：上一节：4.3 软阴影与 PCF · 下一节：5.2 SSAO 的采样与噪声 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。