图形程序员入门球谐函数：解锁实时计算机图形学光照模拟新方法！
为啥关注球谐函数？
在学习计算机图形学过程中，你迟早会在论文或代码里碰到球谐函数。它是实用工具，用几个系数就能近似表示球面上的给定函数，对模拟复杂光照有帮助。虽说球谐函数研究得挺透彻，但理解起来有点棘手。不过，理解它在实时计算机图形学领域的应用并不难。本文会尽力按期望的方式解释球谐函数，助你读完后更好理解相关高级技术资料。阅读需对实时渲染、线性代数和积分有一定了解，但无需高深知识，不会有严格证明，推导也只是基础代数。作为计算机图形从业者，任何将量或值与三维空间方向关联的函数，都可看作定义在单位球面上的函数，之后会交替用“方向函数”“定义在球面上的函数”“定义在单位球面上的函数”表述。实际上，球面上的连续函数都能表示为特殊多项式的无限加权和，这些多项式就是“球谐函数”。截断无限和为有限和，就能近似表示该函数。多项式易计算，而球面上的非平凡函数在计算机图形学中常见。比如，辐射率 $L_i(p, \vec{\omega_i})$ 是球面上的函数，表示从给定方向 $\vec{\omega_i}$ 到达点 $p$ 的光量，立方体贴图可看作其“表格化”形式。给定点的辐照度也可看作方向的函数。用多项式求和近似复杂光照环境很有吸引力。而且，除光照外，很多东西也能表示为球面上的函数。比如，球谐函数可近似表面点沿给定方向的网格或体积厚度，近似厚度可烘焙到纹理图集，用于模拟次表面散射等现象（这想法最初是从好朋友、图形程序员 Bagrat Dabaghyan 那听到的）。所以，球谐函数应用不局限于光照，本文主要关注其在光照方面的应用。
球谐函数的定义
现在该对学习球谐函数充满期待了，那就从定义开始。再次强调，不会进行严格数学推导，目标是提供足够细节，助你阅读和理解引用球谐函数的代码和论文。
函数空间及其基
从线性代数中，我们熟悉“线性”或“向量”空间及其“基”的概念。为定义球谐函数，需在函数领域引入类似概念，Kevin Cassel 的[这个视频](https://www.youtube.com/watch?v=177zEpIwI68)对此有很好解释。首先，把定义在某个定义域上的一组函数看作“空间”。严格说，所有定义要明确函数定义域（视频中也提到这点）。但在我们的应用场景中，定义域始终是单位球面，就不再特别提及。在线性代数中，一组向量 $\vec{v_0}, \vec{v_1}, ..., \vec{v_n}$ 中若没有向量能表示为其他向量的线性组合，就称这组向量“线性无关”。类似地，一组函数 $f_0(x), f_1(x), ..., f_n(x)$ 中若没有函数能表示为其他函数的线性组合，这组函数就线性无关。在线性代数中，有内积（或点积）的概念。对于两个向量 $\vec{p} = (p_0, p_1, ..., p_n)$ 和 $\vec{q} = (q_0, q_1, ..., q_n)$，内积 $\vec{p} \cdot \vec{q}$ 定义为 $\vec{p} \cdot \vec{q} = \sum_{i=0}^{n}p_iq_i$。对于函数 $f$ 和 $g$，内积 $\langle f,g \rangle$ 定义为它们在定义域上乘积的积分：$\langle f, g \rangle = \int f(\omega)g(\omega) \,d\omega$。直观地说，把函数看作无限维向量，积分看作求和的类似操作，两种内积概念的相似性就很明显。同样，函数的范数是该函数与其自身内积的平方根：$\left\| f \right\| = \sqrt{\langle f, f \rangle }$。由此，能理解“正交归一函数集”的概念：一组函数中，任何函数与自身内积为 1，与集合中其他函数内积为 0。对于线性空间，“正交归一基”是一组正交、线性无关的向量，空间中任何向量都能唯一表示为基向量的线性组合。类似地，对于函数空间，“正交归一基”是一组正交、线性无关的函数，空间中任何函数都能唯一表示为基函数的线性组合。不过要注意，函数空间的基集总是无限的！把函数看作无限维向量就好理解了：N 维向量空间的基含 N 个元素，“无限维向量空间”有无限基是合理的。
球谐函数
球谐函数（简称“球谐”）是定义在球面上的特殊函数，构成定义在球面上的所有连续函数空间的正交归一基。回顾定义可知，球面上的任何连续函数，无论多复杂难算，都能表示为球谐函数的无限加权和。若球谐函数易计算，且有有效方法找权重/系数，这将很有实用价值。如前文所说，球谐函数是多项式，易计算，也有方法找系数，后面会讨论。但在计算机上无法处理无限个系数，只能存有限数量的系数。没关系，可接受近似值。但该存多少系数，选哪些系数？丢弃系数会损失哪些信息？下节解答。
球谐函数的阶和次
介绍球谐函数具体形式前，需了解函数的组织方式。球谐函数分为编号组，通常叫“频带”。每个频带有关联数字 $\ell \in {0,1,2,...}$，叫该频带内函数的“阶”。阶为 $\ell$ 的频带含 $2\ell + 1$ 个函数。按惯例，阶为 $\ell$ 的频带内函数从 $-\ell$ 到 $\ell$ 索引，索引叫函数的“次”。要注意，有些资料用“次”表示这里的“阶”，用“相位”表示“次”。我们不用这种术语，但要注意避免混淆！阶为 $\ell$、次为 $m$ 的球谐函数通常表示为 $Y_{\ell}^m$。按频带和次对球谐函数可视化展示很有启发性：图片来自《实时渲染》（Realtime Rendering）一书，原始可视化由 Robin Green 完成。原始版本用红色和绿色可视化，此版本为提高可读性做了色调调整。
从图可见，阶数小的球谐函数捕捉大（“低频”）细节，阶数大的球谐函数捕捉小（“高频”）细节。回顾之前的问题：截断无限和会损失哪些信息？图很清楚表明，会损失试图近似的函数中的小尺度变化信息，即“细节”。若试图近似的函数是“低频”的（变化缓慢），只需存前几个频带的系数。在实时计算机图形学实际应用中，通常不用阶数超个位数的频带，很多应用中，阶数 $\ell = 2$ 就够了。
球谐多项式的形式
到目前为止，只是讨论了球谐函数，还没看到具体形式。现在看看前三个频带内球谐函数的多项式形式，接下来的示例会用到。以下是计算这些多项式的代码，用 JavaScript 编写，因为示例将在浏览器中运行。
// 定义常量，方便编写基函数
const RECIP_PI = 1/Math.PI;
const C = [
Math.sqrt(RECIP_PI) * 0.5,
Math.sqrt(3 * RECIP_PI) * 0.5,
Math.sqrt(15 * RECIP_PI) * 0.5,
Math.sqrt(5 * RECIP_PI) * 0.25,
Math.sqrt(15 * RECIP_PI) * 0.25,
Math.sqrt(70 * RECIP_PI) * 0.125,
Math.sqrt(105 * RECIP_PI) * 0.5,
Math.sqrt(42 * RECIP_PI) * 0.125,
Math.sqrt(7 * RECIP_PI) * 0.25,
Math.sqrt(105 * RECIP_PI) * 0.25
];
// 阶数 l 最大为 3 的球谐基函数
// 球谐基函数定义的来源：
// "Stupid Spherical Harmonics Tricks", Peter - Pike Sloan, 2008
function y00(x, y, z) { return C[0]; }
function y_11(x, y, z) { return C[1] * y; }
function y01(x, y, z) { return C[1] * z; }
function y11(x, y, z) { return C[1] * x; }
function y_22(x, y, z) { return C[2] * y * x; }
function y_12(x, y, z) { return C[2] * y * z; }
function y02(x, y, z) { return C[3] * (3 * z * z - 1.0); }
function y12(x, y, z) { return C[2] * x * z; }
function y22(x, y, z) { return C[4] * (x * x - y * y); }
function y_33(x, y, z) { return C[5] * y * (3 * x * x - y * y); }
function y_23(x, y, z) { return C[6] * z * (y * x); }
function y_13(x, y, z) { return C[7] * y * (5 * z * z - 1); }
function y03(x, y, z) { return C[8] * z * (5 * z * z - 3); }
function y13(x, y, z) { return C[7] * x * (5 * z * z - 1); }
function y23(x, y, z) { return C[9] * z * (x * x - y * y); }
function y33(x, y, z) { return C[5] * x * (x * x - 3 * y * y); }
// 计算阶数小于等于 l 的球谐基函数在给定方向 d 上的值，结果以 Float32 数组形式返回
// 仅支持 l <= 3 的值
function evalSHBasis(d, l) {
const x = d[0];
const y = d[1];
const z = d[2];
switch (l) {
case 0: return new Float32Array([y00(x, y, z)]);
case 1: return new Float32Array([
y00(x, y, z), // l = 0
y_11(x, y, z), // l = 1
y01(x, y, z),
y11(x, y, z)
]);
case 2: return new Float32Array([
y00(x, y, z), // l = 0
y_11(x, y, z), // l = 1
y01(x, y, z),
y11(x, y, z),
y_22(x, y, z), // l = 2
y_12(x, y, z),
y02(x, y, z),
y12(x, y, z),
y22(x, y, z)
]);
default: return new Float32Array([
y00(x, y, z), // l = 0
y_11(x, y, z), // l = 1
y01(x, y, z),
y11(x, y, z),
y_22(x, y, z), // l = 2
y_12(x, y, z),
y02(x, y, z),
y12(x, y, z),
y22(x, y, z),
y_33(x, y, z), // l = 3
y_23(x, y, z),
y_13(x, y, z),
y03(x, y, z),
y13(x, y, z),
y23(x, y, z),
y33(x, y, z)
]);
}
}// 辅助函数：将球谐系数乘以一个标量值
function mulScalarBySHCoeffs(scalar, coeffs) {
return coeffs.map(function(c) { return scalar * c; });
}
// 辅助函数：将两组球谐系数相加
function addSHCoeffs(coeffs0, coeffs1) {
return coeffs1.map(function(c, i) {
return c + (coeffs0.length == 0 ? 0 : coeffs0[i]);
});
}
// 使用简单的蒙特卡罗积分验证基函数是否正交归一
function testBasisFunctions() {
const numSamples = 100000;
const numBasisFuncs = 16; // 支持 l 最大为 3，共 16 个函数
// innerProducts[i] 包含第 i 个球谐基函数与所有球谐基函数（包括自身）的内积
var innerProducts = [];
for (var b = 0; b < numBasisFuncs; ++b) {
// 这些数组将初始化为 0
innerProducts.push(new Float32Array(numBasisFuncs));
}
for (var s = 0; s < numSamples;) {
// 使用简单的拒绝采样生成球面上的随机点：
// * 在 [-1, -1, -1] - [1, 1, 1] 立方体中生成一个点；