Banner

深入理解球谐光照 (Spherical Harmonics):从数学理论到工业级代码实践

球谐光照 (SH) 是实时渲染里”压缩光照”的经典手段:把一张几 MB 的 HDR 环境图,浓缩为 27 个浮点数 (9 系数 × RGB),运行时只需几个 dot 指令就能还原漫反射环境光。这篇笔记主要讲述了数学推导、工程实现、约定差异「最容易踩坑的部分」。


在进入数学和代码之前,先建立一张全景的 Pipeline 图。后面所有章节都可以对应到这条流水线上的某个节点:

① HDR EnvMap ~MB level / Cubemap or LatLong 投影积分 ② Projection CPU 离线 / Compute Shader cₗₘ = ∫ L(ω) Yₗₘ(ω) dω 蒙特卡洛 / 加权求和 乘 A 系数 ③ 漫反射卷积 A₀=π, A₁=2π/3, A₂=π/4 辐照度系数 9 × Vec3 = 27 floats 运行时 ④ Shader 重建 E(n) ≈ Σ cₗₘ Yₗₘ(n) 仅 3-4 条 dot 指令 → Lₒ = ρ · E(n) / π ━━━━ 离线烘焙 (一次性) ━━━━ ━━ 每帧 / 每像素 ━━

记住这条主线,后面任何一段公式或代码都能找到自己的位置。


一、核心概念与直觉

1.1 一句话定义

SH 光照是用一组定义在球面上的正交基函数(球谐基),把环境光照贴图压缩成极少量系数的技术。运行时不再采样纹理,而是用法线方向 求值这些基函数,加权求和得到辐照度。

1.2 类比:球面上的”傅里叶变换”

一维信号 球面信号
任何周期信号 = 不同频率的正弦/余弦 任何球面函数 = 不同阶的球谐基函数
阶数 = 频率 Band = 频率
系数 = 各频率的能量 = 该球谐方向的能量

漫反射 BRDF 本质上是个低通滤波器——光线打到粗糙表面后被四处散射,高频信息已被磨平。所以漫反射环境光只需要保留低频成分,截断到 2 阶 (L=2) 就足够了

1.3 为什么是 9 个系数?

球谐函数按 阶 (Band ) 组织,每阶有 个基函数:

Band 基函数数量 物理含义
L0 1 整体平均亮度(环境光强度)
L1 3 三个主轴方向的亮度梯度(”光从哪边来”)
L2 5 二阶各向异性(”光线分布的形状”)

累计 1 + 3 + 5 = 9 个基函数 → 9 个系数。每个 RGB 通道独立,共 27 个 float,约 108 字节。

1.4 优缺点速览

维度 SH 光照
存储 极小(108 B vs MB 级 HDR 图)
运行时性能 极高(几条 dot 指令,无纹理采样)
适用范围 漫反射、低频环境光
不适用 镜面高光、锐利阴影、强方向光(→ Ringing)

二、数学基础:投影与重建

2.1 实球谐基函数(L≤2)

渲染里用的不是复数 SH,而是 实球谐 (Real SH)——通过对复数 SH 做 sin/cos 拆分得到,保证基集仍然规范正交。给定单位向量

常见困惑:为什么有些基函数带负号?
这些负号是 复 SH → 实 SH 转换时为保持规范正交性而引入的约定,不代表”负能量”。这是 Sloan Stupid Spherical Harmonics Tricks 中给出的标准实数化形式,Unity / UE / glTF 都遵循它。

2.2 9 个基函数的形状(3D 可视化)

教科书里 SH 函数最经典的画法是球面上鼓起的 lobe——基函数值的绝对值作为半径,正负值用不同颜色着色。L0/L1/L2 三排正好对应 1+3+5 的金字塔结构:

⚡ Three.js · SH 基函数实时渲染 (可交互)
拖拽旋转视图 · 滚轮缩放

可以直观看到:

  • L0:球(常数项),就是环境光的”平均亮度”
  • L1:三个互相垂直的偶极子,每个对应 X / Y / Z 一个主轴方向 ——“光从哪边来”
  • L2:四瓣 + 一根带凹腰的双瓣,编码二阶各向异性 ——“光的分布形状”

如果某个 lobe 红蓝面积大致对称,说明它的积分≈0;这也是为什么 SH 基函数集是正交的。

2.3 投影 (Projection):把环境图分解成 9 个数

投影就是对每个基函数,与原始光照函数在球面上做内积:

工程实现上把球面积分离散成加权求和:

关键点:立体角权重 不能省。
经纬展开图在两极拉伸严重,不加 权重的话,极区像素会被严重过采样。

2.4 重建 (Reconstruction):用 9 个数还原方向亮度

给定单位法线 ,重建公式:

注意:上面的 只是”压缩后还原的环境亮度”,而不是漫反射所需的辐照度 。两者差一个余弦核卷积,这就是下一节要补的关键一步。


三、漫反射卷积:神奇的 A 系数从哪来?

3.1 为什么直接重建不是辐照度?

漫反射表面接收到的辐照度是入射辐亮度对法线半球做余弦加权积分:

如果直接拿 重建,得到的是模糊的环境图本身,而不是已经被余弦核卷积过的辐照度。差别就在那个 项。

3.2 Ramamoorthi 的关键洞察

Ramamoorthi & Hanrahan 在 SIGGRAPH 2001 论文 An Efficient Representation for Irradiance Environment Maps 中指出:

漫反射余弦核 是 zonal(轴对称)函数,与 SH 卷积的结果可以闭式求解,且能量在 之后衰减极快。

也就是说,漫反射卷积可以写成按阶逐项缩放的形式:

是只跟阶数有关的常数。

3.3 推导(一遍过,跟得上的话很短)

利用两个工具:

(1)Legendre 展开:余弦核 展开为
,其中

(2)SH 加法定理

代入卷积积分后整理出:

对漫反射核 算 Legendre 展开系数:

0
1
2
3 (奇数 都为 0) 0 0
4 极小 极小 极小

物理意义 这阶整个为零( 是奇函数,对半球积分自动消掉), 起的能量已经几乎可以忽略。这就是为什么截断到 L=2(9 个系数)就够用的根本原因。

最终物理版的 A 系数是:

工程版(预除 ,下文章节四会解释为什么):


四、工程踩坑指南 (Gotchas) ★

这是整篇笔记的核心。SH 数学不难,但工程约定不统一——同样一段代码,A 系数有时是 ,有时是 ,有时干脆消失。下面把这些坑集中盘清楚。

坑 1:A 系数到底乘在哪一步?两套主流约定

约定 A:原始投影 + 重建时乘 A

1
2
投影:  c_lm = ∫ L(ω) · Y_lm(ω) dω         // 原始环境投影系数
重建: E(n) = Σ A_l · c_lm · Y_lm(n) // 在 Shader 里乘 A
  • 优点:系数定义和数学公式一一对应,干净。
  • 缺点:Shader 里多 9 次乘法。

约定 B:投影时就预卷积 + 重建时不乘 A

1
2
投影:  c'_lm = A_l · ∫ L(ω) · Y_lm(ω) dω   // 预卷积
重建: E(n) = Σ c'_lm · Y_lm(n) // Shader 干净
  • 优点:运行时省指令,Shader 只做 dot
  • 缺点:存储的不是”纯投影系数”,调试不方便。

Unity / UE / 大多数游戏引擎用约定 B,所以你看到的 ShadeSH9 里通常没有 A 系数。

坑 2:那段神秘的 color / PI 又是什么?

很多 GLSL/HLSL 重建代码末尾会有一句:

1
color = color.rgb / PI;

这是因为 Lambert BRDF 是 。最终着色:

工程上把 提前合到环境光里——重建时就直接输出 ,那么后面只要 Lo = albedo * envLight,省一次除法。这就是工程版 A 系数 (即 )的由来。

一个快速判断技巧:拿一张 常量亮度 的环境图烘焙。

  • 如果重建结果 ≈ :你在用”E/π” 约定(工程版)。
  • 如果重建结果 ≈ :你在用真实辐照度约定(物理版)。

坑 3:基函数符号要一致!

computeCoefficients 用的基函数和 ShadeSH9 用的基函数必须是同一套约定(包括所有负号)。

  • 投影:coeff += color * Y_lm(dir) * weight —— 这里 带的负号会”吃进” coeff。
  • 重建:color += coeff * Y_lm(normal) —— 同样的负号再次出现。
  • 两个负号相乘抵消,结果正确。

但如果投影端用的是 Y_lm 而 Shader 端不带负号,结果方向反了,环境光会出现”上下颠倒”或”左右翻转”的视觉异常。

坑 4:坐标系陷阱(Y-Up vs Z-Up)

Unity / UE 用 Y-Up(场景物体 transform.up = Y);
但很多论文和 OpenCV 代码用 Z-Up(球坐标 ,Z 向上)。

如果你直接把论文里的 SH 系数表搬到 Unity 中,天空盒会”躺下”。修正方法:在生成方向向量时调换 y/z,或者在 Shader 里 normal.xzy

坑 5:Ringing 振铃伪影

强光源(如 HDRI 中的太阳)属于高频信号,硬截断到 L=2 会引发 Gibbs 现象——某些方向的系数变成负值,画面上出现:

  • 太阳对面出现一圈”暗环”;
  • 法线指向太阳时辐照度爆掉(>1 倍正常值);
  • 极少数情况下出现负辐照度(被 max(color, 0) 截掉)。

缓解方案

  • Windowing:对系数加一个低通窗(如 Hanning),代价是整体变暗、对比度降低;
  • 预卷积:先把 HDR 图低通滤波再投影;
  • Spherical Gaussians:高频用 SG,低频用 SH(混合方案,参考 MJP 博客)。

坑 6:sRGB / Linear 空间

SH 必须在 Linear 空间计算。如果输入贴图是 sRGB(普通天空盒),导入设置必须关闭 sRGB 标记,或在采样时手动 pow(color, 2.2)。HDR .hdr / .exr 默认是 Linear,不用管。

坑 7:归一化常数因子 千万不要漏

经纬展开图遍历时,赤道附近像素和极点附近像素代表的球面面积不同。如果直接像素均匀加权求和,极点会被过采样(画面会有”上下亮中间暗”的偏差)。必须乘 权重。


五、工程实战

5.1 CPU 离线版(C++):标准参考实现

最经典的 SH 投影,用于离线烘焙或工具链。这个版本是约定 B(预卷积)

SphericalHarmonics.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
class SphericalHarmonics {
public:
using SHCoefficients = std::array<Vec3, 9>;
static constexpr float PI = 3.14159265359f;

// 工程版 A 系数 (已预除 π)
static constexpr std::array<float, 9> A = {
1.0f, // l=0
2.0f/3.0f, 2.0f/3.0f, 2.0f/3.0f, // l=1
1.0f/4.0f, 1.0f/4.0f, 1.0f/4.0f, // l=2
1.0f/4.0f, 1.0f/4.0f
};

// 直接由方向向量求 9 个基函数值(避免 sin/cos)
static std::array<float, 9> evaluateBasis(const Vec3& d) {
float x = d.x, y = d.y, z = d.z;
return {
0.282095f,
-0.488603f * y,
0.488603f * z,
-0.488603f * x,
1.092548f * x * y,
-1.092548f * y * z,
0.315392f * (3.0f * z * z - 1.0f),
-1.092548f * x * z,
0.546274f * (x * x - y * y)
};
}

// 投影:环境图 → 9 个 SH 系数
static SHCoefficients project(const std::vector<std::vector<Vec3>>& env,
int width, int height) {
SHCoefficients coeffs{};
float totalWeight = 0.0f;

for (int y = 0; y < height; ++y) {
float v = (y + 0.5f) / height;
float theta = v * PI;
float sinTheta = std::sin(theta);

for (int x = 0; x < width; ++x) {
float u = (x + 0.5f) / width;
float phi = u * 2.0f * PI;

Vec3 dir(sinTheta * std::cos(phi),
std::cos(theta), // Y-up
sinTheta * std::sin(phi));

Vec3 radiance = env[y][x];
float weight = sinTheta; // 立体角权重 (坑 7)
auto basis = evaluateBasis(dir);

for (int i = 0; i < 9; ++i)
coeffs[i] += radiance * (basis[i] * weight);

totalWeight += weight;
}
}

// 归一化:让加权和等价于 ∫ dω = 4π
float norm = 4.0f * PI / totalWeight;
for (auto& c : coeffs) c = c * norm;
return coeffs;
}

// 重建:根据法线方向求辐照度(已乘 A,输出为 E/π)
static Vec3 reconstruct(const SHCoefficients& coeffs, const Vec3& n) {
auto basis = evaluateBasis(n);
Vec3 irradiance{0,0,0};
for (int i = 0; i < 9; ++i)
irradiance += coeffs[i] * (A[i] * basis[i]);
return irradiance; // 注:实际是 E/π,配合 BRDF ρ/π 直接用
}
};

5.2 GPU 加速版(Unity Compute Shader):分块归约

CPU 版对一张 2048×1024 的 HDRI 至少要 200 万次循环,离线烘焙能接受,但工具链需要”按一下立即更新”时就太慢了。Compute Shader 利用 GPU 大规模并行可以快上百倍

核心难点是归约求和 (Reduction)——几十万像素的贡献最终要累加成 9 个数。我们采用分块归约:每个 ThreadGroup(8×8 = 64 线程)算出该块的局部和,CPU 端再做最后的全局累加。

SHProjection.compute

SHProjection.compute
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
#pragma kernel CSMain
#define PI 3.14159265359

Texture2D<float4> _InputTexture;
RWStructuredBuffer<float3> _OutputBuffer; // 每组写 9 个 float3
uint _TexWidth;
uint _TexHeight;

void GetSHBasis(float3 d, out float basis[9]) {
basis[0] = 0.282095;
basis[1] = -0.488603 * d.y;
basis[2] = 0.488603 * d.z;
basis[3] = -0.488603 * d.x;
basis[4] = 1.092548 * d.x * d.y;
basis[5] = -1.092548 * d.y * d.z;
basis[6] = 0.315392 * (3.0 * d.z * d.z - 1.0);
basis[7] = -1.092548 * d.x * d.z;
basis[8] = 0.546274 * (d.x * d.x - d.y * d.y);
}

groupshared float3 gs_cache[64][9];

[numthreads(8, 8, 1)]
void CSMain(uint3 id : SV_DispatchThreadID,
uint3 gtid : SV_GroupThreadID,
uint3 gid : SV_GroupID) {
uint local = gtid.y * 8 + gtid.x;

// 初始化共享内存
for (int k = 0; k < 9; k++) gs_cache[local][k] = float3(0,0,0);

if (id.x < _TexWidth && id.y < _TexHeight) {
float u = (id.x + 0.5) / _TexWidth;
float v = (id.y + 0.5) / _TexHeight;
float theta = (1.0 - v) * PI;
float phi = u * 2.0 * PI;
float st = sin(theta);

// Y-up 球面方向
float3 dir = float3(st * cos(phi), cos(theta), st * sin(phi));
float3 color = _InputTexture[id.xy].rgb;
float weight = st; // sin(theta) 权重,dφ·dθ 留到 CPU 乘

float basis[9];
GetSHBasis(dir, basis);
for (int i = 0; i < 9; i++)
gs_cache[local][i] = color * basis[i] * weight;
}

GroupMemoryBarrierWithGroupSync();

// 组内归约:让 0 号线程汇总 64 个结果
if (local == 0) {
float3 sum[9];
for (int i = 0; i < 9; i++) sum[i] = float3(0,0,0);
for (int t = 0; t < 64; t++)
for (int i = 0; i < 9; i++) sum[i] += gs_cache[t][i];

uint groupsPerRow = (_TexWidth + 7) / 8;
uint outIdx = (gid.y * groupsPerRow + gid.x) * 9;
for (int i = 0; i < 9; i++) _OutputBuffer[outIdx + i] = sum[i];
}
}

SHSolver.cs 调度端

SHSolver.cs
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
public static Vector3[] ComputeSH(Texture2D envMap, ComputeShader shader) {
int width = envMap.width, height = envMap.height;
int kernel = shader.FindKernel("CSMain");

int groupsX = Mathf.CeilToInt(width / 8.0f);
int groupsY = Mathf.CeilToInt(height / 8.0f);
int totalGroups = groupsX * groupsY;

var buffer = new ComputeBuffer(totalGroups * 9, sizeof(float) * 3);
shader.SetTexture(kernel, "_InputTexture", envMap);
shader.SetBuffer(kernel, "_OutputBuffer", buffer);
shader.SetInt("_TexWidth", width);
shader.SetInt("_TexHeight", height);
shader.Dispatch(kernel, groupsX, groupsY, 1);

// 读回 → CPU 归约
var partial = new Vector3[totalGroups * 9];
buffer.GetData(partial);
buffer.Release();

var coeffs = new Vector3[9];
for (int g = 0; g < totalGroups; g++)
for (int k = 0; k < 9; k++)
coeffs[k] += partial[g * 9 + k];

// 补回 dφ·dθ
float dOmega = (2f * Mathf.PI / width) * (Mathf.PI / height);
for (int k = 0; k < 9; k++) coeffs[k] *= dOmega;

// 漫反射卷积:乘 A 系数(物理版含 π;如果 Shader 端最后还要除 π,请用工程版)
float[] A = {
Mathf.PI,
2f*Mathf.PI/3f, 2f*Mathf.PI/3f, 2f*Mathf.PI/3f,
Mathf.PI/4f, Mathf.PI/4f, Mathf.PI/4f, Mathf.PI/4f, Mathf.PI/4f
};
for (int k = 0; k < 9; k++) coeffs[k] *= A[k];
return coeffs;
}

性能数据:2048×1024 的 HDRI,CPU 单线程 ~600 ms,GPU 版 ~3 ms。差距 200×。

5.3 Shader 运行时重建(HLSL)

得到 9 个系数后,运行时只需要:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
// SH[0..8] 已经是预卷积过的 (含 A 系数)
float3 ShadeSH9(float3 n, float3 SH[9]) {
float x = n.x, y = n.y, z = n.z;
float Y[9] = {
0.282095,
-0.488603 * y,
0.488603 * z,
-0.488603 * x,
1.092548 * x * y,
-1.092548 * y * z,
0.315392 * (3.0 * z * z - 1.0),
-1.092548 * x * z,
0.546274 * (x * x - y * y)
};
float3 c = 0;
[unroll] for (int i = 0; i < 9; i++) c += SH[i] * Y[i];
return max(c, 0.0); // 防止 ringing 产生负值
}

// PBR 漫反射最终着色
float3 diffuse = albedo * ShadeSH9(N, _SH); // 注意:约定 B 已含 1/π

Unity URP 的 ShadeSHPerPixel 实际上更紧凑——把基函数项和系数合并进 4 个 float4 里,用 3-4 条 dot4 完成全部计算。


六、交互式演示:SH 重建可视化

下面这个小工具让你直观感受 9 个系数如何”塑造”环境光。点击预设场景,或拖动滑条调节系数,左侧球面展示重建出的辐照度,右侧是经纬展开图

球面 (拖动旋转)
经纬展开图

试着观察:

  • 只开 (按”均匀”预设):球完全均匀,等同环境光常量项;
  • 加大 系数(z 方向):光源出现在”上方”;
  • L2 项 会引入”鞍形” / “条带”形态——这正是为什么 L2 能比 L1 更精细地拟合环境分布。

七、扩展应用与延伸阅读

7.1 现代渲染管线中的 SH

应用场景 角色
Unity Light Probes / UE Indirect Lighting Cache 每个探针存 27 个 float,运行时三线性插值后 dot
Sky / Ambient 光 烘焙天空盒为 SH,作为漫反射兜底
DDGI (Dynamic Diffuse GI) 探针网格 + SH 编码方向辐照度
Lumen (UE5) 局部用 SH,远场用屏幕追踪
VRS / Probe Volumes SH 是基础数据格式之一

7.2 SH 的局限与替代

  • Spherical Gaussians (SG):MJP 系列博客详细对比了 SH/SG 在表达 specular 和高频光时的优劣;
  • Wavelet 球面基:能表达更高频,但运行时开销大;
  • PRT (Precomputed Radiance Transfer):SH 的”重型版”,把 visibility / inter-reflection 也烘焙进去;
  • MixLight (2024):SH(低频)+ SG(高频)混合方案。

7.3 推荐阅读

  • Sloan, Stupid Spherical Harmonics Tricks, GDC 2008 —— 工业界最常被引用的 SH 工程指南。
  • Ramamoorthi & Hanrahan, An Efficient Representation for Irradiance Environment Maps, SIGGRAPH 2001 —— 9 系数足够漫反射的根本论证。
  • Robin Green, Spherical Harmonic Lighting: The Gritty Details, 2003 —— 每一步推导都给出完整数学。
  • MJP’s Blog, SG Series —— SH 视角下的 Spherical Gaussians 对比。

八、一图总结

SH Lighting · 全流程速查 HDR EnvMap 几 MB 投影 (CPU 或 GPU) cₗₘ = ∫ L(ω) Yₗₘ(ω) sin θ dθ dφ ★ 立体角权重 sin θ 不可省 乘 A 系数 (漫反射卷积) A = {π, 2π/3, π/4} 物理版 {1, 2/3, 1/4} 工程版 9×Vec3 108 字节 Light Probe ━━━ 运行时 (Per Pixel) ━━━ Shader 重建 E(n)/π = Σ cₗₘ Yₗₘ(n) ~ 3-4 dot 指令 Lambert BRDF Lₒ = ρ · E(n)/π 直接乘 albedo 即可 最终颜色 输出 sceneColor ⚠ 坑点提醒:约定差异 (A 是否预乘 / π 是否预除) · 坐标系 (Y-up vs Z-up) · 振铃 · sRGB

结语

球谐光照的精彩之处不在于数学有多深,而在于它用 27 个浮点数解决了实时漫反射环境光这件事——这种”压缩到极致还能用”的优雅,是计算机图形学里最经典的近似艺术之一。

下次看到 Shader 里那一串看似神秘的 0.2820950.488603color / PI,希望你能笑着说一句:”哦,这就是把环境图傅里叶变换之后的低频系数嘛。”