图形程序员入门球谐函数：解锁实时计算机图形学光照模拟新方法！

为啥关注球谐函数在学习计算机图形学过程中你迟早会在论文或代码里碰到球谐函数。它是实用工具用几个系数就能近似表示球面上的给定函数对模拟复杂光照有帮助。虽说球谐函数研究得挺透彻但理解起来有点棘手。不过理解它在实时计算机图形学领域的应用并不难。本文会尽力按期望的方式解释球谐函数助你读完后更好理解相关高级技术资料。阅读需对实时渲染、线性代数和积分有一定了解但无需高深知识不会有严格证明推导也只是基础代数。作为计算机图形从业者任何将量或值与三维空间方向关联的函数都可看作定义在单位球面上的函数之后会交替用“方向函数”“定义在球面上的函数”“定义在单位球面上的函数”表述。实际上球面上的连续函数都能表示为特殊多项式的无限加权和这些多项式就是“球谐函数”。截断无限和为有限和就能近似表示该函数。多项式易计算而球面上的非平凡函数在计算机图形学中常见。比如辐射率 $L_i(p, \vec{\omega_i})$ 是球面上的函数表示从给定方向 $\vec{\omega_i}$ 到达点 $p$ 的光量立方体贴图可看作其“表格化”形式。给定点的辐照度也可看作方向的函数。用多项式求和近似复杂光照环境很有吸引力。而且除光照外很多东西也能表示为球面上的函数。比如球谐函数可近似表面点沿给定方向的网格或体积厚度近似厚度可烘焙到纹理图集用于模拟次表面散射等现象这想法最初是从好朋友、图形程序员 Bagrat Dabaghyan 那听到的。所以球谐函数应用不局限于光照本文主要关注其在光照方面的应用。球谐函数的定义现在该对学习球谐函数充满期待了那就从定义开始。再次强调不会进行严格数学推导目标是提供足够细节助你阅读和理解引用球谐函数的代码和论文。函数空间及其基从线性代数中我们熟悉“线性”或“向量”空间及其“基”的概念。为定义球谐函数需在函数领域引入类似概念Kevin Cassel 的[这个视频](https://www.youtube.com/watch?v177zEpIwI68)对此有很好解释。首先把定义在某个定义域上的一组函数看作“空间”。严格说所有定义要明确函数定义域视频中也提到这点。但在我们的应用场景中定义域始终是单位球面就不再特别提及。在线性代数中一组向量 $\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_{i0}^{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, 2008function 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 0y_11(x, y, z), // l 1y01(x, y, z),y11(x, y, z)]);case 2: return new Float32Array([y00(x, y, z), // l 0y_11(x, y, z), // l 1y01(x, y, z),y11(x, y, z),y_22(x, y, z), // l 2y_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 0y_11(x, y, z), // l 1y01(x, y, z),y11(x, y, z),y_22(x, y, z), // l 2y_12(x, y, z),y02(x, y, z),y12(x, y, z),y22(x, y, z),y_33(x, y, z), // l 3y_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)]);}}编写处理球谐函数的代码时很可能从其他地方复制粘贴基函数这可能出错或原始来源本身有误。实际上准备本文时就遇到过这种情况。幸运的是验证使用的基函数不难。要记住球谐基的基本属性是正交归一任何基函数与自身内积应为 1与其他不同基函数内积应为 0。这些内积是积分可用最简单的[蒙特卡罗](https://en.wikipedia.org/wiki/Monte_Carlo_integration)方法数值计算。以下是验证上述基函数的代码// 辅助函数将球谐系数乘以一个标量值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) {// 这些数组将初始化为 0innerProducts.push(new Float32Array(numBasisFuncs));}for (var s 0; s numSamples;) {// 使用简单的拒绝采样生成球面上的随机点// * 在 [-1, -1, -1] - [1, 1, 1] 立方体中生成一个点