球谐光照与PRT学习笔记（三）：球谐函数

计算机图形学话题下的优秀答主

Spherical Harmonics的"Harmonics"一般译作“调和”，也颇有一些“谐波”的味道在里面，所以SH就叫球谐啦。而[2]调和函数(harmonic function)就是[3]拉普拉斯方程(Laplace's Equation) $\nabla^2 f=0$ 的解，解一般都含有 $\sin$ 和 $\cos$ 项，所以满足拉普拉斯方程的函数就是"Harmonic"的。在傅立叶分析里面，我们有时会把单位圆上的周期函数(periodic function)用调和函数进行展开。把这操作推广到 $n$ 维球面，我们就会得到球面调和(Spherical Harmonics)。

所以类比一维实函数的傅立叶变换的基底 $\sin(nx)$ 和 $\cos(nx)$ ，本文会讲Spherical Harmonics展开的基底，就姑且叫它球谐基函数吧。不过在此之前先引入一下勒让德多项式(Legendre Polynomial)，我们要基于勒让德多项式来给出球谐基函数的解析形式。

（ps：脑壳疼，末流985软件工程学的数学和物理感觉太少了....）

函数的正交

函数也可以基于“基函数”进行展开，大一数学分析接触的一维函数傅立叶变换的“函数基底”就是谐波 $\sin(nx)$ 和 $\cos(nx)$ 。这些基底函数一般都是[4]正交函数(orthogonal function)，这里的正交是在函数空间(function space)&向量空间(vector space)里面的广义的正交，两个正交函数在给定定义域里面的内积为0。如果两个不同的函数 $f$ 和 $g$ 在定义域上正交，那么：

$\langle f, g \rangle= \int f(x)g(x)dx=0$

这跟欧氏几何的向量的垂直与内积很相似。不过其实你要是把函数看成一个无穷维向量，其实这个内积与正交的概念也是可以接受的= =。

所以其实跟欧氏几何也有点类似，欧氏向量可以通过跟标准基点乘得到投影长度（例如三维向量 $(3,4,5)\cdot(1,0,0)=3$ 就得到了在x轴上的投影）。如果给我们一个函数 $f$ ，我们可以通过函数的点积把他们投影到函数基底上。

那么其实正交基函数不只是有 $\sin(nx)$ 和 $\cos(nx)$ 这么一对的，还可以有很多很多种正交函数，只要满足了上面的函数内积性质就好。有梯子的朋友可以上youtube看一下这个关于函数内积与正交函数的一些小例子[5]：

https://www.youtube.com/watch?v=8ZyeHtgMBjkwww.youtube.com

勒让德多项式

有一族正交函数是数学领域比较关注的，就是正交多项式(Orthogonal Polynomial)。正交多项式有一些奇妙的特性（或者说其实这个应该是它的定义吧?），这个特性跟正交基函数族类似：

$\int^1_{-1}F_m(x)F_n(x)dx=\left \{ \begin{aligned} 0 &&(n \ne m)\\ c&&(n=m) \end{aligned} \right.$

特别地，如果 $c=1$ ，那么这组正交多项式 $F(x)$ 就是标准正交(orthonormal)的。有很多种这样的多项式，例如[4]Chebyshev Polynomial，Jacobi Polynomial，Hermite Polynomial等等。而在球谐光照里面，我们最关心的就是勒让德多项式(Legendre Polynomial)，特别是伴随勒让德多项式(Associated Legendre Polynomial，简称ALP吧)。根据wikipedia的说法[6]，勒让德函数是勒让德微分方程(Legendre Differential Equation)的解：

$(1-x^2)\frac{d^2P(x)}{dx^2}-2x\frac{dP(x)}{dx}+n(n+1)P(x)=0$

勒让德方程是物理和工程领域里面常常遇到的一类常微分方程，也是拉普拉斯微分方程的一种变形，当试图在球坐标中求解三维拉普拉斯方程（或者其他偏微分方程的时），问题经常会归结为勒让德方程的求解。当方程满足 $|x|<1$ 时，可以得到有界解（解级数收敛）。并且当 $n\in Z^*$ 时， $x= \pm1$ 处也有有界解。在这种情况下，方程的解随着 $n$ 值变化而变化，构成的一组由正交多项式组成的多项式序列，称为勒让德多项式：

$P_n(x)=\frac{1}{2^n\cdot n!}\frac{d^n}{dx^n}[(x^2-1)^n]$

而[7]伴随勒让德多项式(ALP)则有两个参数 $l,m$ ，可以基于普通勒让德多项式 $P_n(x)$ 来定义：

$P^m_l(x)=(-1)^m(1-x^2)^{m/2}\frac{d^m}{dx^m}(P_l(x))$

Associated Legendre polynomialsen.wikipedia.org

关于伴随勒让德多项式，有几点要注意一下的：

$l,m$ 分别是ALP的"degree/band index"和"order"。伴随勒让德多项式这一族正交多项式是分了很多"层"/"带"的，在每一个band $l$ 里面， $m$ 表示元素在当前band的元素下标。而且注意两个参数的取值范围： $l,m\in Z^*,m\in [0,l]$ 。这意味着前 $n$ 个band总共有 $\frac{n(n+1)}{2}$ 个多项式。大概是这样子：

$P^0_0(x)$

$P^0_1(x),P^1_1(x)$

$P^0_2(x),P^1_2(x),P^2_2(x)$

...........

[7]上面 $P^m_l(x)$ 表达式的 $(-1)^m$ 叫做Condon–Shortley phase，有些人会忽略掉它。（这是因为我们后面的球谐基函数是基于 $P^m_l(x)$ 定义的，这个 $(-1)^m$ 有时会直接定义在球谐基函数那）
[1]普通的勒让德多项式是复数，伴随勒让德多项式则是实数。

一般我们都不会直接求解 $P^m_l(x)$ ,计算复杂度与重复度都很高，在高阶的情况（不能直接hardcoded了）时我们求解会用递推迭代的方法。几条足以求出所有ALP的递推式：

其中eq.2可以算出ALP的斜线的首项，eq.3可以在斜线的基础上往higher order推一项（只有一项），如果这都算不到给定的勒让德多项式，那就用1不断地往高阶继续递推算下去。

还得注意一个坑点，eq.2的"!!"是双阶乘(double factorial)，这东西的意思并不是两个阶乘的复合，而是有其他定义的。[8] n是自然数时，n的双阶乘表示不超过n且与n有相同奇偶性的正整数乘积：

$(2n-1)!!=1\times3\times5\times...\times(2n-1)$

$(2n)!!=2\times 4 \times 6\times ... \times (2n)$

$0!!=1!!=1$

当n是负奇数时，

$[-(2n+1)]!!=\frac{(-1)^n}{|-1|\times|-3|\times|-5|\times...\times|-(2n-1)|}$

n是负偶数时，分母有 $(-2)!!=0!!/0$ ， $n!!$ 无意义。

球谐函数

球谐（基）函数是定义在球面上的。（下面球谐就简写为SH）。SH函数在通用情况下是在复数的基础上定义的，但是我们在图形学里面一般只关心定义在球面的实函数（在球谐光照里面，这个实函数可以是多通道的光强分布，也可以是传输函数，等等），所以这篇文章就只关心实球面调和(Real Spherical Harmonics)。

引入一下球谐函数的解析式吧。先给出单位球上坐标的参数化[1]：

$(\sin\theta \cos\phi,\sin\theta \sin \phi, \cos \theta) \rightarrow(x,y,z)$

那么，一般我们把球谐函数记为 $Y^m_l$ ，表达式是

$\begin{equation} Y^m_l(\theta,\phi)= \left\{ \begin{aligned} \sqrt{2}K^m_l \cos(m\phi)P^m_l(\cos\theta)& & (m>0)\\ \sqrt{2}K^m_l \sin(-m\phi)P^{-m}_l(cos\theta) & & (m<0)\\ K^0_lP^0_l(cos\theta)& & (m=0) \end{aligned} \right. \end{equation}$

其中 $P^m_l$ 就是前文讲的伴随勒让德多项式，突然出现的 $K^m_l$ 是一个缩放系数，用来归一化：

$K^m_l=\sqrt{\frac{2l+1}{4\pi}\cdot \frac{(l-|m|)!}{(l+|m|)!}}$

注意一点，虽然ALP的"degree"参数 $m$ 是非负整数，但是球谐函数的参数 $m$ 则是可以取到负数的，取值区间关于0对称，也就是：

$l\in Z^*,m\in Z,m\in [-l,l]$

有的时候，把SH函数的项展开成1D向量（变成个线性表）是挺有用的，所以我们就可以定义一个单参数的球谐函数项序列：

$y_i(\theta,\phi)=Y^m_l(\theta,\phi)$ ，其中 $i=l(l+1)+m$

明显 $Y^0_0$ 就是第0项， $Y^0_1$ 是第2项， $Y^0_2$ 是第6项， $Y^0_3$ 是第12项....

维基百科[10]有前几阶用 $(x,y,z)$ 参数化的实球谐函数解析式，你可以直接照着他公式来hardcode在代码里面。

https://en.wikipedia.org/wiki/Table_of_spherical_harmonicsen.wikipedia.org

Table of spherical harmonics

https://en.wikipedia.org/wiki/Table_of_spherical_harmonicsen.wikipedia.org

很多论文都会直接把球谐函数的多项式扔给你，没有可视化，非常抽象。上面的图就给出了前五个band的球谐函数 $Y^m_l(\theta,\phi)$ 的可视化结果。绿色是正值，红色是负值，离中心越远的地方绝对值就越大。这些基底你会发现它们都很对称，这个对称主要应该还是来自于谐波函数的对称性( $\sin\theta$ 和 $\cos\theta$ )。[9]提到定义在球面 $S$ 上的函数 $f(\theta,\phi)$ 可以用球谐函数展开成二重广义傅立叶级数。所以笔者认为所以实球谐函数也可以看成是一维傅立叶变换的函数基底 $\sin(nx),\cos(nx)$ 的复杂版本（且函数是定义在单位球面上）。

$f(\theta,\phi)$ 在球面 $S$ 上的展开式为：

$f(\theta,\phi)=\sum^\infty_{l=0}\sum^l_{m=-l}C^m_lY^m_l(\theta,\phi)$

其中 $C^m_l$ 就是频域上的系数了。系数的求解依旧是要走一波积分变换(integral transform)：

$C^m_l=\int _S f(s)Y^m_l(s)ds$

也就是频域上的系数 $C^m_l$ 要求我们求原函数 $f$ 与球谐函数 $Y^m_l$ 的乘积在球面上的积分。一般我们把这个求系数的过程叫做投影(Projection)。在实际操作中，我们写程序时不可能会有对无穷级数进行储存和卷积的操作，一般展开项只能是有限项，也就是：

$\hat f(\textbf s)=\sum ^{n-1}_{l=0}\sum^l_{m=-l}C^m_lY^m_l(\textbf s)=\sum^{n^2}_{i=0}c_iy_i(\textbf s)$

其中n是球谐基band的数量，显然n个band的球谐基数量是 $\frac{(1+2n-1)n}{2}=n^2$ 个（自己翻回去看看前几阶球谐基的图）。这个过程是一个带限的近似(band limited approximation)。在这个语境下，因为频域信号带宽的限制，大于一定阈值的高频信号就被去掉了。所以我们只能用 $n^2$ 个预计算的球谐系数(SH Coefficient)和球谐函数本身近似地重建出原函数：

从图中可以看出，球谐展开阶数越高，能重构出来的信号就越精确。

在实际操作中，球面积分一般也不太能直接求出解析式，毕竟如果场景不同、灯光不同的话，某点上的入射光强的还是会有很多的变化的，我们还是得用数值方法来"投影"。这就是为什么我之前专门铺垫了一篇文章专门讲球面上的均匀采样与Monte-Carlo积分：

回忆一下蒙特卡洛积分估计器：

$\int _S fds\approx \frac{1}{N}\sum ^N_{j=1}f(x_j)w(x_j)$ ，其中权重 $w(x_j)=\frac{1}{p(x_j)}$

那么对于均匀的球面上采样来说，权重就是：

$w(x_j)=\frac{1}{p_(x_j)}=\frac{1}{(\frac{1}{S_{球面}})}=4\pi$

所以球谐系数 (SH Coefficient)的数值估计表达式是：

$\begin{aligned} C_i&=\int _S light(\textbf v)Y_i(\textbf v)ds \\ &\approx \frac{1}{N} \sum^N_{j=1}light(x_j)Y_i(x_j)\cdot 4\pi\\ &=\frac{4\pi}{N} \sum^N_{j=1}light(x_j)Y_i(x_j) \end{aligned}$

在写代码实现的时候，用蒙特卡洛积分求系数的过程大概就是一系列的相乘与求和，伪代码：

void SH_Coefficients()
{
    double weight =4.0 * PI;
    //生成n条光线进行采样
    for(int i=0; i<n_samples; ++i) 
    {
    生成带抖动的无偏采样方向(θ,ϕ)
        for(int n=0; n<n_coeff; ++n)
        {
        //对于某一个light probe，它的每个球谐展开系数c_i就要累加起所有的【某方向上的irradiance * 这个方向上SH函数值】
        result[n] += light(θ,ϕ)* samples[i].SH_basis_coeff[n];
        }
    }
    // 把蒙特卡洛积分的常数项乘上去（恒定的采样权重，总采样数）
    double factor = weight / n_samples;
    for(i=0; i<n_coeff; ++i)
    {
        result[i] = result[i] * factor;
    }
}

上面的result[i]就是某个点上光照分布的球面函数经过”编码”之后得到的球谐系数。我们可以用离线预计算的系数，在runtime通过效率比较高的方式近似重构出原来的光照：

$light(\textbf s)=\sum^{n^2}_{i=0}C_iY_i(\textbf s)$

一般我们可以用SH coefficients来编码低频的环境光（因为边缘是亮度变化剧烈的高频信息，所以用SH来编码看起来糊糊的indirect lighting其实是比较合适的，例如拿来编码天光skylight）。至此，虽然还没涉及到PRT与着色相关的问题，但是你可以用Spherical Harmonics的"投影&重建"这对基本操作做你爱做的事了。

图：阶数越高，一般重建的质量越高。但是很难重构得像原图那么精致了，毕竟只有少量的系数，只能重构出大概的信息

下一篇讲球谐函数的性质以及球谐旋转（SH Rotation）

引用

[1] Green R. Spherical harmonic lighting: The gritty details[C]// Game Developers Conference. 2003

[2] Harmonic function

[3] Laplace's equation

[4] Orthogonal functions

[5] https://www.youtube.com/watch?v=8ZyeHtgMBjk

[6] https://zh.wikipedia.org/wiki/%E5%8B%92%E8%AE%A9%E5%BE%B7%E5%A4%9A%E9%A1%B9%E5%BC%8F

[7] Associated Legendre polynomials

[8] https://baike.baidu.com/item/%E5%8F%8C%E9%98%B6%E4%B9%98

[9] 19. 5 球函数

[10] Table of spherical harmonics