Spherical Harmonic Lighting（球谐光照）

Robin Green

这篇文章只是Spherical Harmonic Lighting这个论文的解释，看的时候请参照原文，原文需要一定的数学与信号处理，对于基础一般的朋友可能稍有难度，所以我想将其理论通俗化的讲出来，降低一些难度，希望对需要的人有所帮助。
    利用球谐光照技术，能让我们实时重现面积光源下3D模型的全局光照效果，在SIG02的一篇论文中曾有介绍（Sloan，Kautz，Snyder）利用这种技术实现模型的超真实光照效果。从这个技术的发展来看，他其实是N多相关技术的合集，游戏开发者可以利用他们来获得VERY GOOD的画面效果。技术本身的代码和结果都是很简单明了的，但这篇文章的目的在于阐述其背后的原理，以飨广大初学者。告诉你WHY，给予你使用SH LIGHTING的所有力量。

 

如果你曾经有开发3D引擎的经验，相信你一定对PHONG光照模型相当的熟悉，镜面高光、漫反射颜色以及环境光照等概念几乎可以成为你的第二感。你使用的最简单的光照模型可能就是“漫反射模型”，通常也被称为“点积光照”。如果场景中有nLights个光源，则模型表面的光照颜色用如下公式计算：N为表面发线，L为光照方向。

这个光照模型只不过是真实物理照光的简化，只有完全的模拟物理才能得到真正照片级的真实。但物理光照公司计算相当复杂，实现实时计算困难。看下面的物理公式，实际上是在半球上对光线方向（V）的亮度（L）函数进行积分。

实现光线跟踪唯一的难点就是这个立体角微分在半球面上的积分，这个教程只是简单的叙述了全局照明的相关概念，我强烈建议读者更深入的学习相关知识。接下来我简单阐述几个热力学方面的概念：
通量：光源发出或表面接受到的辐射功率（W）。
通量密度：单位面积上的通量（W / M2）。
辐射度：表面发射的通量密度。
辐照度：表面接受的通量密度。
BRDF：双向反射分布函数。实际就是描述了空间中表面上的某个点上入射通量密度与出射通量密度的关系函数（微分比）。当然它并不能反映出物质与光的物理交互过程，纯粹只是一个结论的数学抽象。此函数必须具有出射入射互易不变性和归一性等特征。
立体角：在半径为r的球面上面积为A的区域对应的立体角为w ＝ A / r2。由此可知，在球面上极角为Θ，方位角为Φ处的立体角微分为dw = sinΘdΘdΦ。

现在来理解上面物理方程，其实就是出射亮度，等于入射亮度立体角微分在半球上的积分。而半球上的入射亮度由全局光照事先产生，比如利用辐射度或者光线跟踪，这就是所谓的Precomputed Radiance Transform。这个积分是如此难解，基于符号的积分解法在GPU中很难完成，为了实时的完成每个象素的出射亮度积分运算，人们使用了一个近似的积分解法，虽然不能得出精确结果，但却能在效率和效果之间得到一个良好的平衡，它就是Monte Carlo积分法，这个方法又跟概率论有点关系，没办法，我不得不简要的提一下相关知识，但不会扯得太远：）

概率密度描述了一个随机变量的值在区域的概率分布。
  
    期望值则描述了一个以随机变量作为输入的函数值趋势，即平均值。

    而计算期望值的另一种方法就是大数据量采样平均。

综上可得：

      这个就是传说中的Monte Carlo积分估算法。现在回到球上的BRDF与入射亮度积分，要使用Monte Carlo积分法计算这个积分，就必须在球面上进行采样，因为是均匀采样，所以有：

所以P（x）= 1/4π，现在我们需要一对互相独立且取值在[0-1]之间ξx 和ξy的坐标来索引球面上的点，Θ为极角，Φ为方位角。

这个映射比较难理解的是ξx到Θ的关系，为什么他没有使用ξxπ，很显然在上式中ξx = sin2Θ/2，有知道的朋友一定通知我一声。我们现在可以将球面分为NxN份采样区域，这些区域很显然类似地球经纬线将地球表面划分的情况，在区域范围内进行随机采样，这种采样方式被称为分层抽样或类型抽样，关于他们的具体知识请参考概率论书籍，这种方式的优点在于，采样点在球面上分布比较均匀，不会出现簇聚的情况，这样的样本更能反映群体特征。

终于到了SHL的重点难点了，之所以说是重点难点，主要因为这里牵涉到稍多的数学及信号处理方面的知识，需要一定的基础。其实SHL的核心就是球面亮度信号编码和重建。在高数中我们曾学过，信号满足一定条件下，可以分解为一系列正弦谐波的和，谐波频率以倍频增长，这就是所谓的傅立叶级数。经过全局光照计算后，物体表面上每个点会得到一个球面的亮度信号，但我们不可能为每个点都保存一个环境贴图，因此我们需要对这个定义在球面的亮度信号进行编码，而在实时重现时，利用编码快速重建原球面亮度信号，进而计算光照效果。这就是SHL的来龙去脉，现在进行详细阐述。

正交基函数

一个原始信号，可以分解为一系列带缩放谐波之和，如何我们要在以后利用这些谐波（基函数）来重建原始信号，我们必须事先得到每个基函数相对应于原始信号的系数（缩放因子），比如我们要计算信号f(x)中基函数b(x)分量的系数（权重），我们需要在f定义域上对f(x)b(x)进行积分，也就是所谓的卷积，原文章称这一过程为projection。而重建原始信号，则将经相应系数缩放后的所有基函数求和即可。

当然有很多种基函数可以用来重建信号，比如正弦信号（傅立叶级数）等，但我们最感兴趣的是一种在数学上称为正交多项式的基函数。他们有一个相当有意思的特性，如果你对任何两个多项式进行卷积运算，如两多项式相同，你会得到一个常数，如果不同你会得到0值，直观上讲，就是虽然这些多项式定义在相同作用域中，但他们缺可以互不影响。我们要讨论的是一种叫做伴随勒让德多项式的正交函数，通常用P表示这类多项式，伴随勒让德多项式具有两个参数l和m，定义域为[-1，1]，返回值为实数（一般勒让德多项式返回的是复数，不要搞混淆了），l和m两个参数将这些多项式归为若干组，l便是组的编号，取值为大于0的整数，而m取值为[0，l]之间的整数，同一组中的任意2个多项式之间卷积为一个固定的常数项，而不同组中的任意2个多项式之间卷积又是另外一个固定的常数。很显然一个n组勒让德多项式一共有n(n+1)项，如果用这些多项式进行信号projection可以得到n(n+1)个参数。因为勒让德多项式定义本身就非常麻烦，所以很少用它来重建一维信号。勒让德级数原始数学定义很是复杂，又是导数又是符号变换又是虚数，不适合浮点计算，因此我们用递归的方式来定义它，只需要3条规则即可。

三条规则见原文……

以上三条规则定义了一维勒让德多项式，但如何将它扩展在二维球面呢？伴随勒让德多项式是球谐的核心所在，其意义与傅立叶变换一样，只不过它的定义域为球面坐标。本来SH（球谐）函数是用来定义虚数的，但现在我们只关心其实数域，我们称这些球谐函数为实域球谐函数，以下我们讨论的都是这类球谐函数。实域球谐函数定义见原文，需要注意的是这里的勒让德多项式定义需要稍加改变，m的取值范围要从[0 - l]改为[-l - l]，以下有一个表示实域球谐函数的形象图形，红色是负数，绿色是正数，到原点距离为返回值的绝对值大小，单位小球表示了球面上正负值分布情况。

很显然，当l = 0时，球谐函数可以用来近似环境光照，可以用在AO计算中。而l = 1时，因为球面上值呈COS分布，故可以用来近似漫反射计算中的COS项，欧了！现在我们要利用这些球谐函数对球面上的原始入射亮度信号进行编码了，只要在球面上做卷积运算就可以了。

    而重建原始信号更简单，事先计算的卷积值与对应谐函数相乘的和即为近似的原始信号。

显然，一个n阶球谐近似，需要n的2次方个系数，因为2 * （0+1+2+…+（n - 1））/ 2 + n = n*n。理论上我们需要无穷项的球谐函数，才能完美重建原始信号，所以我们无法做到完美，只能近似，这里我们只取有限的低频谐函数，而将高频函数忽略掉，这就是所谓的带宽限制，这注定我们重建的信号会丢失掉很多高频信号，即亮度信号的细节变化。

     从上图可见，即使用了100个参数（10阶球谐），也丢失了很多的细节，细节越多（高频），丢失越严重，如第二个信号（因为细节多的信号其能量在高频部分比例很大）。
    篇外：这里有篇小波重建信号的文章，据说效果比球谐效果好很多，基于小波的SHL改进 (http://graphics.stanford.edu/papers/allfreq/)，有兴趣的朋友可以研究一下。
      基础理论到此基本就讲完了，可能有点晕乎，现在我们理论结合实践，来看看这些理论是如何发挥作用的。利用MONTE CARLO积分法将一信号投影成球谐函数的系数。我们利用之前确定的球面坐标系在定位点（极角与方位角），同时使用比较低阶的球谐函数，这样就不会有太多的系数需要计算，现在假设一个简单光照环境，有两个略为垂直的单色光源，我们可以用如下公式来表示这个光照环境：

    计算球谐系数的积分公式如下计算球谐系数的积分公式如下:

    这个积分公式可以利用Maple或者Mathematica符号积分运算工具计算，但我们还是需要用MONTE CARLO数字估算法，计算如下：
  
    为什么Sin项没有了？因为是对球面积分，后面 r²sinΘdΘdΦ即是球面微分，换成一个二维微分单位就是dΩ，推导可参看立体角的定义。重建亮度信号方式如下:

    看上图重建结果，对于16个系谐波重建还是比较满意的，虽然丢失了很多高频信息，但毕竟我们将计算量缩减得如此厉害。细心的朋友可能发现在右图中原本应该黑暗的区域出现了一个亮度突起，这回产生错误的光照效果，物体本来黑暗的背面可能突兀的出现一小块亮斑，在重建信号中出现这个突起的原因主要是原始信号里能量比较高的高频对低频产生了干扰，这个高频其实就是MAX函数将负数CLAMP成0导致的，虽然用更高阶数了球谐函数来重建原始信号可惜消除这个干扰，但显然会增加不少的计算量，其实有个等效的做法，我们可以事先消除信号中的一部分高频，然后再进行卷积编码，这个消除过程就是利用低通滤波器过滤原始亮度信号，高斯卷积核其实就是一个不错的低通滤波器，对原始信号采样前做个高斯模糊就OK了，虽然能降低甚至消除高频干扰，但带来的负面影响就是细节和突变更少了，也就是说信号更模糊了。

上部基本上把SHL的基础知识都过了一遍，中部主要会讲SH的一些属性，而下部将介绍SHL程序实现。