本文主要介绍球谐（Spherical Harmonic，简称SH）函数在光照中的一些计算实现，其内容来自于GDC2003的演讲：Spherical Harmonic Lighting: The Gritty Details

学习总结

球谐函数是一组正交基函数，两两相乘的积分结果是0，而自身相乘的积分结果为1，任意信号都可以通过与球谐函数相乘积分算出其在对应球谐函数上的系数，这个过程可以看成是信号在球谐函数上的投影，通过多个球谐函数按照对应系数累加可以得到原始信号的模拟，参与模拟的球谐函数阶数越高，模拟精度也就越高。

球面坐标系（$\theta, \phi$）下面的球谐函数可以表示任意点到球心的距离，而这个距离也可以解读成强度，从而可以用于实现某点处各个方向上的输入光强。同时，每个点处的输入光强与输出光强的转换关系（BRDF之类）也可以使用球谐函数来表示，实际光照就是上述两个球谐函数相乘的积分输出，而在实际计算中，如果在离线的时候完成两个球谐函数的系数的求取，在运行时只需要一个系数向量点乘即可完成，大大简化了计算量，提升了计算速度。

背景简介

球谐光照（SH Lighting）是Sloan，Kautz以及Snyder大神在2002年的Siggraph上给出的用以实现超真实光照模型渲染的技术。使用SH Lighting可以让我们实现面光源光照模型的计算捕捉存储以及实时重建一个全局光照般的效果。

Siggraph2002的原文实现比较晦涩，对于一些新手来说可能不太友好，而本文则会对SH实现的背景假设以及推导做详细介绍，以达到新手也能快速入门的目的。

光照计算

作为游戏引擎开发者，不能不知道一些光照模型的相关知识。所有光照模型中最简单的一种应该就是diffuse surface reflection model了，因为这个模型中核心的计算就是对光照方向与表面法线方向进行点乘，而得名"dot product lighting"（点乘光照）。

对于场景中存在多个光源，每个光源的光照强度（intensity）都用RGB来表示的话，那么某个像素的最终反射光亮度（Irradiance）就可以表示成如下的公式：

\[ I=surfacecol*\sum_{i=1}^{nlights}lightcol_i*(\mathbf{N}\cdot\mathbf{L}_i) \]

其中$L_i$表示的是第$i$ 个光源的光照方向，$N$ 表示像素的世界法线，$lightcol_i$ 则是第$i$ 个光源的光照强度，surfacecol是当前像素的基色（diffuse/albedo?）。

上述的公式是处于对计算时间的考虑而进行的对真实的物理渲染公式的一个简化处理，真实的物理渲染公式由于过于复杂而无法实现实时计算：

在真实的物理渲染公式中：

$L(x, \vec{w_o})$ 指的是从当前像素反射而出的朝着 $w_o$ 方向反射的光照强度；
$L_e(x, \vec{w_o})$ 指的是当前像素自发光的强度;
$f_r(x, w_i\rightarrow\vec{w_o})$指的是当前像素点对应的入射方向为$w_i$，出射方向为 $w_o$ 的BRDF（Bilateral Reflectance Distribution Function）;
$L(x^{'}, \vec{w_i})$指的是从其他位置$x^{'}$ 沿着 $w_i$ 方向反射到当前像素上的光照强度;
$G(x, x^{'})$指的是当前像素 $x$ 与光源 $x^{'}$ 之间的几何关系函数;
$V(x, x^{'})$指的是 $x$ 与 $x^{'}$ 之间的可见性函数。
从这个公式可以看出，当前像素 $x$ 沿着 $w_o$ 方向出射的光照强度等于两个部分之和：其一是自发光强度，其二是反射光强。其中反射光强是最为复杂的一项：由各个方向的入射光考虑可见性以及几何关系之后积分而成。

如上图所示，假设不考虑时间的流动，在一个装满光子的空间volume中，可以认为每个位置的光子密度是均匀的常量；而如果假设这个volume中的光子都沿着某个固定的方向匀速流动，那么在一定的时间段内，就会有一定数量的光子穿过流动方向成一定夹角的截面，为了进行定量的分析，我们需要检测单位时间内穿过截面的光子数量，这个检测的结果称为flux，通常用joules/second(J/s，也就是watt)为单位。

而从空间几何关系计算，可以得知，flux的数值与截面法线方向跟流动方向夹角的余弦成正比，夹角越大（比如前面两个方向越接近垂直的话），余弦越小，flux数值越小。

为了表征光子密度流速以及截面之间的关系，我们将flux除以截面的面积，就得到了单位时间流经单位截面面积的光子数目，这个数值通常用光辐射照度irradiance来表示，其单位为w/m2。

irradiance考虑的是单位面积下的光照强度，不过按照主观意识，可以知道，光强是随着距离而递减的，如果将光源看成是球面输出模型的话，那么就可以理解成，距离越远，球面表面积就越大，此时均分得到的光强自然越弱，但实际上一束光线沿着某个固定的方向射出，不考虑途中的损耗，其光强应该是不随距离的变化而变化的，为了表征这个特性，引入单位固体角下的光辐射照度这个表达项，这个数值通常称为光照亮度radiance，其单位为w/(m2*sr)，sr是固体角的单位。

如果对上图中的volume的边长进行缩小逼近，我们可以得到一个微分 ($dA, dt, d\Omega$) 的volume模型,，这个模型可以用来近似单束光的数值表现，之后根据需要进行积分运算就能得到实际需要的各种结果。

这里简要介绍了物理真实的光照渲染模型的一些基本背景以及实现公式，之后要考虑的就是怎样将这个实现公式用在实时渲染中，下面将会介绍实际渲染中使用的各种模拟逼近方法。

蒙特卡罗积分

我们现在已经有一个输入光强作为变量的积分函数光照模型，不过由于目前这个模型还存在很多不明确的部分，使得我们无法通过符号运算的方法直接给出积分结果。对于这种情况来说，蒙特卡罗积分方法是一个比较合适的计算方法，蒙特卡罗算法主要的计算过程都是通过概率论的方法完成的，所以不需要了解那些不明确部分的具体数值或规律。

先来介绍一下概率论的一些基础知识。一个随机变量指的是在一个指定的定义域中按照某种分布规律输出对应数值的变量。比如我们要抛一个六面的骰子，那么出现1~6中任意一个数值的概率都是1/6，这个我们抛出的数值x就是一个随机变量。

随机变量x出现的规律我们称之为概率分布函数(probability distribution function, PDF)，简称$P(x)$，$P(x)$ 是一个输出数值永远都为正数的函数，且对相应的随机变量按照全定义域进行积分或者求和，得到的结果一定为1。

\[ \int^{+\infty}_{-\infty}P(x)=1, \ \ \ where \ P(x)\ge0 \]

对于出现数值小于某个指定数值 $x$ 的概率称之为累计概率函数（cumulative distribution function， CDF），简称$C(x)$。$P(x)$ 与 $C(x)$ 之间的关系是导数与积分的关系：对 $C(x)$ 进行微分（连续）或者差分（离散）就得到了 $P(x)$，对 $P(x)$ 进行积分（连续）或者求和（离散）就得到了 $C(x)$。在刚才的骰子实验中，$P(x) = 1/6, C(x) = x/6$（此处假设分布式连续的，而非离散变量）。此外对于出现任一数值的概率都相同的分布函数，我们称这个随机变量符合均匀分布（uniform distribution），这个随机变量成为均匀随机变量。其中定义域为[0,1)的连续均匀随机变量由于在采样理论中会经常出现，我们这里给其设定一个专有的名字：归一化均匀随机变量(canonical random variable)，且为其设置一个专有符号 ξ（/ksi/）。

随机变量x出现在[a,b]范围内的概率，可以通过积分得到：

\[ P(x\in[a, b])=\int^b_ap(x)dx \]

每个由随机变量 $x$ 作为自变量的函数 $f(x)$ 都能够计算得到一个平均值，这个平均值我们称之为 $f(x) $的期望值 $E[f(x)]$。

\[ E[f(x)]=\int f(x)p(x)dx \]

假设某个均匀随机变量的定义域为[0,2]（可以计算出P(x) = 1/2），而以x作为自变量的函数$f(x) = 2 - x$；可以计算得到其期望值：

\[ E[2-x]=\int_0^2 (2-x)/2dx=1 \]

除了上述直接通过符号化的公式计算函数的期望值之外，我们还可以通过采样足够多的的数据（大数定理）来进行数值的模拟（蒙特卡罗估计方式）：

\[ E[f(x)]\approx\frac{1}{N}\sum_{i=1}^Nf(x_i) \]

根据这个期望的计算公式，我们可以推导出一个积分公式如何转换为一个求和计算公式，这个公式就称之为蒙特卡罗估计公式：

\[ \int f(x)dx=\int \frac{f(x)}{p(x)}p(x)dx \approx\frac{1}{N}\sum_{i=1}^N\frac{f(x_i)}{p(x_i)} \]

这个公式表达的意思是，一个积分公式，可以转换为一个采样点除以采样值出现概率之后的平均期望。而从这个公式可以看到，当积分函数f(x)与概率函数p(x)的一致性越高，那么$f(x)/p(x)$的方差就越低，期望计算时所需要的采样点数目也就越少（保证计算结果精度的前提下），在极端情况，比如 $f(x) = p(x)$ 的时候（这里说的应该是基本上等于，否则根本无需采样，直接计算就完了），可能只需要一个采样点就足以计算出f(x)的结果。蒙特卡罗估计方式的另一种写法是将$1/p(x)$ 用 $w(x)$来表示（权重相乘求和）：

\[ \int f(x)dx\approx\frac{1}{N}\sum_{i=1}^Nf(x_i)w(x_i) \]

这个公式的直观意义是，出现概率越大的采样数值，在积分中贡献的权重越低（因为出现概率越大，说明在采样结果列表中出现的频次越高，而积分的意义却是每个值只出现一次，所以从这个角度来看，这种处理是符合实际需要的）。

按照这个公式，假如我们能够保证概率分布函数$p(x)$ 是一个均匀概率函数的话，此时$p(x) = 1/N$，前面的求和公式就可以近似为$\sum f(x_i)$，那么我们在计算图形积分的时候，只需要通过点采样（point sample）后求和，就能够得到积分的近似结果了。

现在回到前面的物理真实的渲染公式，这个渲染公式是一个对半球的积分，这个积分可以通过在半球上均匀散布的采样点（更专业一点的说法叫做无偏随机采样点：unbiased random samples）进行采样求和计算来近似。这里需要两个相互独立的均匀随机变量 $ξ_x$ 以及 $ξ_y$，这两个变量分别代表了横竖垂直的两个采样方向，其表征的范围是一个方形区域，通过下面的公式，可以转换成球面坐标：

\[ (2arccos(\sqrt{1-\xi_x}), 2\pi \xi_y)\rightarrow(\theta, \varphi) \]

平面坐标到球面坐标的转换

由于对于整个球面的采样点是均匀分布的，也就意味着各个采样点之间的概率对于整个球面分布而言是相同的常量，而球面的表面积公式为 $4\pi R^2$ ，其中R = 1，可以简化为4π，球面上任意一点出现的概率就应该是1/4π，对应的公式计算中的权重函数$w(x) = 4π$. 根据这种分布方式得到的采样点分布大致如下图所示：

蒙特卡罗积分公式是对期望积分的模拟与简化，而采样点相对于期望值越分散（方差越大），表明对应函数曲线突变越明显，此时模拟的精度就越差，而在实际使用中，通常可以使用一种分层采样（stratified sampling）的方法来降低方差。这种方法是如何实现的呢？那就是将输入的方块区域分割成N x N个cell单元组成的细小网格，之后在每个cell中生成随机采样位置（cell中的采样点称为抖动采样点，jittered sample）进行采样。

这里有个问题，如果cell都是等大的，那么怎样才能保证不同区域的采样点密度不一样呢？实际上，之所以称为分层采样，是因为采样点并不是统一在一起的，而是分为多层进行，比如这里的这个例子，就可以分成两层进行采样：

所谓的分层采样，是根据一些先验条件对采样位置进行人为干预以降低方差的方法，为什么这种方式可以降低方差呢，其实就跟曲线拟合是差不多的，在变化陡峭的地方增加一些采样点，在变化平缓的地方减少采样点，从而保证最终的曲线拟合精度足够高。

经过证明，使用分层采样方式得到的采样结果中，每个cell的采样方差之和不会大于未分层之前的采样方差，实际情况中通常要远远低于这个值。当然，在实际应用中还有许多小trick可以增加蒙特卡罗积分的精确性，不过对于基本的SH光照计算来说，分层采样就已经够用了。

下面给出为某个方块生成抖动采样点（每行cell数目为sqrt_n_samples）并计算每个采样点的SH系数的示例代码：

正交基函数（Orthogonal Basis Functions）

正交基函数是一组相互正交的基础函数（类似于相互垂直的坐标轴，这里的正交指的是任意两个基函数相乘后的结果的积分都是0，而基函数本身相乘后的积分结果为1），通过这些基函数的线性组合，可以实现一些复杂函数（信号）的模拟：

\[ f(x)=\sum_ic_i*B_i(x)\\ \int B_i(x)*B_j(x)dx=0,\ for\ i\neq j\\ \int B_i(x)*B_j(x)dx=1,\ for\ i= j \]

$B_i(x)$ 为基函数，$c_i$ 是此基函数的加权系数。在给定基函数的前提下，求取 $f(x)$ 对应的各个基函数的系数 $c_i$ 的过程，称为投影。根据基函数的特性，我们可以通过对$f(x) * B_i(x)$进行积分来求取$c_i$。

候选基函数有许多，而其中一些非常有意思的函数被数学家单拎出来组成了一个基函数组，并被冠名正交多项式(orthogonal polynomial，这个基函数组内的基函数存在以下的积分特性：

\[ \int_{-1}^1F_m(x)F_n(x)dx=\begin{cases}0&\text{for} n\neq m\\c&\text{for} n=m\end{cases} \]

而如果我们对上述公式中返回的结果做进一步限制，比如将返回值设定为0，或者1（而非c）的话，那么这个基函数组就缩小为归一化正交基函数(othonormal basis function)，这类多项式正交基函数通常会以它们的研究者命名，其中最让人感兴趣的一组函数被称作"勒让德多项式(Legendre Polynomials，此处我们指的是更为细分的伴随勒让德多项式，Associated Legendre Polynomials)"。

勒让德多项式通常用符号 $P$ 表示，而伴随勒让德多项式一共包含两个参数 $l$ 与 $m$，这些多项式的定义域为$[-1, 1]$，其中普通勒让德多项式返回的结果是复数域的，而伴随勒让德多项式返回的结果是实数域的。

伴随勒让德多项式的两个参数 $l$ 跟$m$ 会将这个基函数组分割成更为细分的单位，其中处于同一个 $l$ 中的多个基函数称为这些函数处于一个函数带中（band of functions），而 $l$ 被称为带系数(band index)，其取值范围为任意非负整数，确定了 $l$ 之后，$m$ 的取值范围就变成了 $[0, l]$，按照这个命名规则，伴随勒让德多项式矩阵应该是一个三角形状：

处于同一个带内的勒让德多项式其正交的结果（也就是相乘后的积分结果）为某一个常数c1，而不同带内的勒让德多项式的正交结果为另一个常数c2.

勒让德多项式的处理过程通常比较耗时，因此很少用于对一维函数的模拟，且这个系列的基函数通常是面向导数定义的，且处理过程非常繁琐，这个过程对于浮点数计算非常不友好。因此，出于简化处理的考虑，我们此处会利用前后多项式的关系（比如按照递归方式）来根据前面的生成结果生成后续的多项式，这种生成方式只需要遵循三条规则：

\[ ( l- m) P_l^m= x( 2l- 1) P_{l- 1}^m- ( l+ m- 1) P_{l- 2}^m \]

这是我们给出的多项式生成规则中最主要的一条，利用前面两个band的多项式直接计算出当前band对应的多项式，不过注意，这个公式中由于存在 $P^m_{l-2}$ 部分，也就隐式包含了其使用限制：l-2 >= m。也就是说无法处理$m = l$ 以及 $m + 1 = l$ 这两种情况，因此给出了下面的两条补充规则：

\[ P_m^m=\begin{pmatrix}-1\end{pmatrix}^m\begin{pmatrix}2m-1\end{pmatrix}!!\begin{pmatrix}1-x^2\end{pmatrix}^{m/2} \]

这个规则主要用于处理规则1中的m无法等于l的漏洞，且这个公式计算得到的多项式不需要前面band的多项式结果。注意，公式中的两个感叹号指的是对这个值进行双阶乘（注意，负奇数的双阶乘是其绝对值的双阶乘的倒数，负偶数的双阶乘是无意义的），其结果是所有不大于2m -1的奇数相乘，将m = 0代入公式得到 $P_0^0(x=1)$。

\[ P_{m+ 1}^m= x( 2m+ 1) P_m^m \]

规则3给出了$l = m + 1$ 的补充规则。

有了这三条规则，我们就能够生成整个勒让德多项式矩阵了，其使用方式是先根据规则2生成 $P_m^m$

（对应于之前的勒让德多项式矩阵的斜边），之后根据规则3生成 $p_{m+1}^m$（对应于勒让德多项式矩阵的次斜边，也就是斜边各元素正下方的第一个元素），之后根据规则1，填充剩余的各个元素。注意，使用规则1计算得到的多项式结果其四舍五入误差要小于按照规则3计算得到的结果（怎么看出来的？另外，这两个公式的应用范围没有重叠，也没有办法比较误差吧？），计算代码给出如下：

def P(l, m, x):
    # evaluate an Associated Legendre Polynomial P(l,m,x) at x
    pmm = 1.0
    if m > 0:
        somx2 = np.sqrt((1.0 - x) * (1.0 + x))
        fact = 1.0
        for i in range(1, m + 1):
            pmm *= (-fact) * somx2
            fact += 2.0

    # Rule 2
    if l == m:
        return pmm

    # Rule 3
    pmmp1 = x * (2.0 * m + 1.0) * pmm
    if l == m + 1:
        return pmmp1

    # Rule 1
    pll = 0.0
    for ll in range(m + 2, l + 1):
        pll = ((2.0 * ll - 1.0) * x * pmmp1 - (ll + m - 1.0) * pmm) / (ll - m)
        pmm = pmmp1
        pmmp1 = pll

    return pll


# 测试例子
l = 3
m = 2
x = 0.5
print(P(l, m, x))

球谐函数（Spherical Harmonics）

前面介绍的都是一些基础的数学知识，那么这些知识跟球谐函数有什么关系呢？从前面的公式中可以看出，勒让德多项式表达的是对一维信号的模拟实现，而对于二维三维以及更高维的信号，比如二维球面上的数据信号，又该如何表达呢？

实际上，球面二维信号有一种经典的表达方式，这种方式叫做Spherical Harmonic，简称SH（球谐函数），这种数学系统其实跟傅里叶变换比较类似，不过是直接定义在球面上的，其核心的算法就跟勒让德多项式有关。

球谐函数原本是定义在虚数空间，不过我们此处只关心实数域的模拟（比如球面上的光照强度），因此本文后续介绍的所有内容都是跟实数域SH相关的，如果不做解释，默认指的实数域的SH。

单位半径球面上的点的笛卡尔坐标与标准参数化方式得到的球面坐标存在着以下关系：

\[ \begin{pmatrix}\sin\theta\cos\varphi,\sin\theta\sin\varphi,\cos\theta\end{pmatrix}\to\begin{pmatrix}x,y,z\end{pmatrix} \]

经典的球谐函数公式形式给出如下，可以看到其核心的部分是以$\cos(\Theta)$ 为自变量的勒让德多项式，前面提到的一维多项式怎么表达二维信号的问题，从这个公式中可以得到解决，即将一维勒让德多项式作为核心组成部分，分别揉入二维信号的自变量 (两个角度$θ,φ$), 就得到了两维信号的表达：

$$
y_l^m(\theta,\varphi)=\begin{cases}\sqrt{2}K_l^m\cos\bigl(m\varphi\bigr)P_l^m\bigl(\cos\theta\bigr),&m>0\[2ex]\sqrt{2}K_l^m\sin\bigl(-m\varphi\bigr)P_l^{-m}\bigl(\cos\theta\bigr),&m<0\[2ex]K_l^0P_l^0\bigl(\cos\theta\bigr),&m=0\end{cases}

其中$P_l^m$ 项是伴随勒让德多项式，$K_l^m$ 是用于进行归一化处理的scale因子：

\[ K_l^m=\sqrt{\frac{\left(2l+1\right)}{4\pi}\frac{\left(l-\left|m\right|\right)!}{\left(l+\left|m\right|\right)!}} \]

为了能够表达所有的球谐函数，此处的I跟m的定义跟前面伴随勒让德函数中的定义有所区别：
m不再是$[0,l]$,而是$[-l,l]$。

\[ y_{l}^{m}( \theta , \varphi )\ \ \ where\ l\in\mathbf{R}^{+},-l\leq m\leq l \]

$y_l^m$ 就是对应于伴随勒让德多项式 $P_l^m$ 的一个分量，将 $y_0^0$~$y_n^n$ 的多个分量加起来，得到的就是球面信号的一个近似模拟，而由于各个分量是相互独立的，因此可以将球谐信号看成是各个球谐分量按照顺序逐一发生的信号，从而可以展开成多个一维的信号（其实就是将之前的三角形排布的信号按照顺序放到同一行中）：