浪淘沙令·伊吕两衰翁
[宋代][王安石]
伊吕两衰翁，历遍穷通。一为钓叟一耕佣。若使当时身不遇，老了英雄。
汤武偶相逢，风虎云龙。兴王只在谈笑中。直至如今千载后，谁与争功！

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

学习总结

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

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

背景简介

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

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

光照计算

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

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

Diffuse Surface Reflection

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

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

真实物理渲染公式

在真实的物理渲染公式中：
指的是从当前像素反射而出的朝着方向反射的光照强度；
指的是当前像素自发光的强度
指的是当前像素点对应的入射方向为，出射方向为的BRDF（Bilateral Reflectance Distribution Function）。
指的是从其他位置x'沿着方向反射到当前像素上的光照强度
指的是当前像素x与光源x'之间的几何关系函数
指的是x与x'之间的可见性函数。

从这个公式可以看出，当前像素x沿着方向出射的光照强度等于两个部分之和：其一是自发光强度，其二是反射光强。其中反射光强是最为复杂的一项：由各个方向的入射光考虑可见性以及几何关系之后积分而成。

投影角示意图

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

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

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

固体角下的光辐射强度

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

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

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

蒙特卡罗积分

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

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

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

概率积分

对于出现数值小于某个指定数值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]范围内的概率，可以通过积分得到：

积分求概率

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

期望计算

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

期望示例

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

模拟期望计算

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

模拟公式推导

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

蒙特卡罗估计方式

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

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

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

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

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

球面无偏分层随机采样点按照(θ,φ)坐标展示图

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

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

1. 先将(θ,φ)坐标空间先均匀分成M x M个方形区域，每个方形区域中指定其分割的份数：出现概率越高，越大，从而保证在低概率区域采样点数目少于高采样点区域；
2. 在每个方形区域进行进一步的细分，分割成 x 个cell，每个cell生成一个随机点。

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

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

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

struct SHSample 
{
  Vector3d sph;
  Vector3d vec;
  double *coeff;
};

void SH_setup_spherical_samples(SHSample samples[], int sqrt_n_samples)
{
  // fill an N*N*2 array with uniformly distributed
  // samples across the sphere using jittered stratification
  int i=0; // array index
  double oneoverN = 1.0/sqrt_n_samples;
  for(int a=0; a<sqrt_n_samples; a++) 
  {
    for(int b=0; b<sqrt_n_samples; b++) 
    {
      // generate unbiased distribution of spherical coords
      double x = (a + random()) * oneoverN; // do not reuse results
      double y = (b + random()) * oneoverN; // each sample must be random
      double theta = 2.0 * acos(sqrt(1.0 - x));
      double phi = 2.0 * PI * y;
      samples[i].sph = Vector3d(theta,phi,1.0);
      // convert spherical coords to unit vector
      Vector3d vec(sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta));
      samples[i].vec = vec;
      // precompute all SH coefficients for this sample
      for(int l=0; l<n_bands; ++l) 
      {
        for(int m=-l; m<=l; ++m) 
        {
          int index = l*(l+1)+m;
          samples[i].coeff[index] = SH(l,m,theta,phi);
        }
      }
      ++i;
    }
  }
}

正交基函数（Orthogonal Basis Functions）

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

                         f(x) = sum(ci * B_i(x));
                Integral{ B_i(x) * B_j(x) } = 0, for i != j;
                Integral{ B_i(x) * B_j(x) } = 1, for i == j;

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

投影

拟合

拟合结果

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

正交多项式特性

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

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

前六个Associated Legendre Polynomials

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

伴随勒让德多项式命名

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

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

规则1

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

规则2

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

规则3

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

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

double P(int l,int m,double x)
{
    // evaluate an Associated Legendre Polynomial P(l,m,x) at x
    double pmm = 1.0;
    if(m>0) 
    {
        double somx2 = sqrt((1.0-x)*(1.0+x));
        double fact = 1.0;
        for(int i=1; i<=m; i++) 
        {
            pmm *= (-fact) * somx2;
            fact += 2.0;
        }
    }
    // Rule 2
    if(l==m) 
        return pmm;
    
    // Rule 3
    double pmmp1 = x * (2.0*m+1.0) * pmm;
    if(l==m+1) 
        return pmmp1;
    
    // Rule 1
    double pll = 0.0;
    for(int ll=m+2; ll<=l; ++ll) 
    {
        pll = ( (2.0*ll-1.0)*x*pmmp1-(ll+m-1.0)*pmm ) / (ll-m);
        pmm = pmmp1;
        pmmp1 = pll;
    }
    return pll;
}

球谐函数（Spherical Harmonics）

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

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

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

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

坐标转换

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

球谐公式

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

scale因子

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

球谐参数范围

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

球谐展开

第i个分量就是前面勒让德多项式矩阵中第i个元素，其实现代码给出如下：

double K(int l, int m)
{
    // renormalisation constant for SH function
    // 阶乘函数最好提前算好，运行时直接读取
    double temp = ((2.0*l+1.0)*factorial(l-m)) / (4.0*PI*factorial(l+m));
    return sqrt(temp);
}

double SH(int l, int m, double theta, double phi)
{
    // return a point sample of a Spherical Harmonic basis function
    // l is the band, range [0..N]
    // m in the range [-l..l]
    // theta in the range [0..Pi]
    // phi in the range [0..2*Pi]
    const double sqrt2 = sqrt(2.0);
    if(m==0) 
        return K(l,0)*P(l,m,cos(theta));
    else if(m>0) 
        return sqrt2*K(l,m)*cos(m*phi)*P(l,m,cos(theta));
    else 
        return sqrt2*K(l,-m)*sin(-m*phi)*P(l,-m,cos(theta));
}

The first 5 SH bands plotted as unsigned spherical functions by distance from the origin and by colour on a unit sphere. Green (light gray) are positive values and red (dark gray) are negative.

上图给出了前面5个band的球谐函数的图形展示，其中球谐函数的数值用的是到圆心的距离来表示（比如(x) = 1，因此l=0的时候，就是一个等距球面），而由于m的取值范围是[-l, l]，因此，而其中的红色绿色分别表示cos(theta)是正负两种情况下的结果。

最终在使用的时候，将这些不同参数下的球谐函数作为基函数，通过线性组合的方式就能够进行二维信号的模拟表达（就像前面一维基函数通过线性组合就能够对给出的一维信号进行模拟一样，以作为自变量，到球心的距离作为因变量，看起来像是三维信号）

球谐投影（SH Projection）

在已知球谐函数的前提下，计算某个函数对应于这个函数的系数的过程比较简单：只要将两个函数相乘之后积分即可：

球谐系数求取

在上面这个公式里，积分的自变量s指的是所有采样点，后续会对这个公式进一步具象成球面参数坐标变量作为自变量的形式。

在得到球谐基函数的系数之后，通过系数相乘求和的方式，就能得到原函数的一个近似模拟：

函数重建

从这个公式可以知道，为什么说n阶的模拟实现需要n^2个球谐系数了。另外，当n越大，这个近似函数就越接近于原函数。不过通常使用的是有限n的版本，这种模拟近似也被称作带限模拟（band-limited approximation）。

模拟示例

下面通过一个示例介绍如何通过蒙特卡罗积分实现球谐投影以计算得到各个球谐基函数的系数。第一步需要对球面进行参数化，这里使用的是之前球谐函数使用的的球面坐标系统的参数化方式。第二步则是选择一个低频函数来进行投影（低频则保证了较少的系数就可以得到较为精确的模拟，band越低，频率越低），这里选定了两个相互垂直的大型单色光源作为待模拟函数，下面的公式给出了这个函数的球面坐标表示方式，不过在实际的实现程序中，我们会直接使用光线追踪的方式从几何的角度直接评估此类函数：

待模拟函数

An example lighting function displayed as a color and a spherical plot.

两个max对应的形状

球面积分公式给出如下：

球面积分

大家可能会好奇，这个积分公式中为什么会有一个sin(theta)项。这是因为，所谓的积分其实就是将小块面积进行累加，当小块的面积足够小的时候，就逼近于积分了。而组成球面的一个个小块（按照角度划分）在赤道附近的面积要大于两极位置的面积，而这个sin就是为了来表示这个含义而增加的。

将这个公式跟前面的球谐函数系数求取公式结合起来，可以得到如下的参数化球谐系数投影函数：

SH Projection

这个公式其实已经算是比较明晰了，通过一些常用的数学软件（比如mathematics或者maple）是可以从符号数学的角度推导出结果的，不过我们在实际应用中是没有这些软件的，只能通过蒙特卡罗积分这种数值分析类的方式来求得结果：

蒙特卡罗积分

这个公式中的x_j指的是采样点，而f(xj) = light(xj)yi(xj)，由于我们前面说到，我们的样本都是无偏样本，因此每个样本的概率应该都是相同的，也就是说w(x_j)应该是一个常量（1/(1/4π) = 4π）。而无偏采样的另一个特性就是不论按照什么样的参数化方式，都应该以相同的概率得到相同的采样样本，因此，前面的双重积分可以直接用蒙特卡罗的简化版积分来模拟：

蒙特卡罗积分投影

通过蒙特卡罗积分公式，可以我们需要的球谐基函数进行预计算，其实现代码非常简单：

typedef double (*SH_polar_fn)(double theta, double phi);

void SH_project_polar_function(SH_polar_fn fn, const SHSample samples[],
double result[])
{
    const double weight = 4.0*PI;
    // for each sample
    for(int i=0; i<n_samples; ++i) 
    {
        double theta = samples[i].sph.x;
        double phi = samples[i].sph.y;
        for(int n=0; n<n_coeff; ++n) 
        {
            result[n] += fn(theta,phi) * samples[i].coeff[n];
        }
    }
    // divide the result by weight and number of samples
    double factor = weight / n_samples;
    for(i=0; i<n_coeff; ++i) 
    {
        result[i] = result[i] * factor;
    }
}

通过这段代码，我们可以计算我们之前光源进行投影可以得到以下系数：

球谐系数

重构公式

根据这些系数进行重构得到了如下结果：

重构结果

由于只使用了16个系数的SH模拟，所以这个模拟表现应该算得上非常不错了。不过右边这张图上出现了一个增生，这种增生在渲染中的表现就是背部无光区域会出现一块异常的亮斑，之所以会出现这种异常，是因为我们将SH中高频的部分通过max（clamp）去掉了，导致信号的不连续，当我们使用越来越高阶的SH来进行模拟的时候，这个部分的增生会逐渐消退，而为了能够使用低阶SH模拟的时候也能得到较好效果，就需要在光源的设计上下功夫，比如在进行模拟之前，先将光源通过一个窗口函数（比如高斯函数，低通滤波）干掉高频数据，从而避免上述瑕疵。

SH函数属性

相对于其他可选的基函数，SH函数有很多非常适合用来进行光照计算的属性：

1. 比如SH函数是归一正交的，也就是任意两个球谐函数相乘之后的积分只有当这两个函数完全相同时，结果为1，其他情况为0.
2. SH函数是旋转不变的，也就是说将某个函数f旋转后得到g函数，这两个函数之间存在如下的关系（球谐模拟后的旋转结果等于旋转后进行球谐模拟的结果，这个特性可以保证在光照位置方向等输入发生了平滑旋转变化的时候，输出也是平滑旋转的，不会出现剧烈突变等瑕疵）：

旋转不变性

3.两个球谐模拟函数相乘的积分等于对应系数相乘之后求和：

球谐相乘积分

对于光照计算而言，通常需要将输入光L与物体表面材质的传递函数t（Transfer Function）相乘之后进行球面积分来求取反射输出的光照结果，如果硬算的话基本上是无法做到实时完成的，而使用了球谐函数之后，只需要将L跟t都用球谐模拟的系数表示，根据球谐函数的归一化正交特性，两者相乘之后的积分，就等于对应系数相乘之后求和。

这里给出了球谐函数的一种用法：将两个球谐函数相乘之后的积分，就等于系数相乘之和，其结果是一个标量。除了这种用法之外，球谐函数还可以用来对球谐函数进行转换（transforming）。

比如说，我们有一个光源的球谐函数表示a(s)，以及一个对表面上的点产生遮挡阴影的球谐函数（用于表征当前点周边的遮挡情况的）表示b(s)，如果直接将a(s)用作输入光源的SH函数显然是不合适的，因为中间隔着一个b(s)，那么是否可以先将b(s)的作用计算出来，在最终进行光照计算的时候利用这个作用先对a(s)进行转换，得到遮挡处理后的球谐函数表示c(s)呢？其实是可以做到的，通过光线追踪的方式，我们能够将b(s)的作用表示成一个传递矩阵（transfer matrix）：

传递矩阵

这里有两个疑问尚未解决：1.这个公式是怎么推导出来的？2.因为这里相当于三个基函数相乘之后的积分，其结果不能简单通过系数相乘得到，那么这个公式要怎么计算出结果？难道是代入伴随勒让德多项式硬算？

之后将这个矩阵与a(s)的系数进行线性运算，就能得到c(s)的系数：

传递矩阵使用方法

这里的积分是对j进行的，是不是意味着可将ai提取出来，放到求和之外，所以整个公式就变成了:?

这里用一个实际的例子来阐明这种做法的效果：简单考虑，假如我们的阴影函数就是某个SH基函数，那么矩阵元素求取公式就变成了：

矩阵元素

这里计算得到的结果是非常稀疏的（0多1少），按照25*25来进行计算，得到的矩阵结果用1表示黑点，0则什么也不绘制，那么其结果图如下所示（难道真的是硬算的？感觉不能做到实时完成？不过因为整个计算过程不需要输入光源的数据，倒是可以通过离线完成）：

传递矩阵示例

将这个矩阵应用在光源的SH函数上，得到的新的SH函数的系数，其结果就跟将光源的球谐表示与阴影的球谐表示相乘之后的结果差不多，如下图所示：

示例结果

有了这种方法之后，我们可以在不进行最终的光照与投影物之间的交互计算的前提下完成阴影数据的记录。

旋转球谐函数（Rotating SH）

球谐函数的旋转不变特性是最难以理解的，首先需要搞明白的球谐投影函数的旋转是什么意思：按照欧拉角进行旋转？以什么轴序进行？通过3x3旋转矩阵完成？还是通过quaternion完成？

从正交的角度，我们能够理解的旋转过程应该是一个线性操作，且不同band之间的系数是无交互的。按照这种说法，从实际计算的角度出发，我们可以表示为：我们可以通过n^2 x n^2的块状对角矩阵（如下图所示的格式，块状是因为band之间的系数是存在交互的）将某个SH向量（系数组成的向量）转换成另一个SH向量：

转换矩阵

在旋转矩阵中元素已知的情况下实现旋转是一个简单地过程，不过矩阵元素的计算却不简单。为了解释这个过程，我们来看下球谐基函数表示的另一种形式——前面用的是球面坐标形式，现在转换到笛卡尔坐标形式：

笛卡尔坐标SH基函数形式

这种格式的SH当然也可以如球面坐标形式一样用来进行投影以求取对应的系数，不过更常用的做法是用来进行符号运算。通过对旋转过的SH样本与未旋转过的样本相乘之后进行球面符号积分，我们可以构建出前面提到的旋转矩阵：

旋转矩阵元素

示例

什么原理？

通过这个n^2 x n^2矩阵，我们可以将未旋转的样本转换为旋转后的样本。如果取n = 3,我们将得到由前三个band组成的SH系统的绕着z轴旋转的9x9旋转矩阵：

z轴旋转alpha

按照上面矩阵的规律，我们可以猜测将这个矩阵扩充到N band的话，得到的最边缘的非零元素中的角度应该是Nalpha（准确来说，应该是（N-1)alpha）。

如果只使用低阶的SH函数的话，这种技术就显得非常有用了，我们可以将旋转拆分成多个简单的旋转步骤，之后在得到计算结果后再重新组合出来，在实际使用中，如果SH函数高过两阶，这种技术使用起来就很成问题（什么问题，消耗过高？）。

在使用这个技术之前，我们需要考虑使用哪种旋转拆分方式，可以用最少的步骤达到想要的目的。比如说，如果我们采用ZYZ的旋转方式（欧拉角的旋转方式一共有十二种，采用三个轴组合旋转的有A(3,3) = 6种，采用两个轴组合旋转的有C(3,2) * 2 = 6种，注意，两轴旋转需要满足相同两轴并不连续，比如不能出现XXY，只能是XYX），最终可能只需要进行两次旋转就能够达到目标（其中一次已经是我们想要的格式的旋转矩阵？每太看懂这段话，后面再回过头来分析一下），那么要怎么完成对Y轴的旋转呢？这个旋转可以拆分成对X轴的90度旋转后面接一个普通的z轴旋转，之后再接一个x轴的-90度的旋转：

x轴的旋转矩阵

ZYZ形式的旋转，按照上述拆分可以变成如下格式：

矩阵拆分

整个计算过程需要涉及到4次9x9矩阵的乘法（还没算三角函数的计算消耗），如果矩阵运算的时间复杂度为O(n3），那么整个运算过程需要消耗2916个MAD指令，虽然我们可以利用矩阵的稀疏性将5阶旋转矩阵（这里应该说的不是5x5，也不是5次相乘，而是band5（也就是25个系数））的计算消耗降低到340个指令，但是这个消耗依然不低。

那么是否可以将上述的复杂的旋转与乘法操作整合到一个矩阵中呢？实际的结果矩阵比较吓人，比如上述计算结果的二阶表示：

二阶结果

这个结果只能用来进行调试，如果想在实际游戏运行中使用，就得将引擎中的旋转都转换成ZYZ的形式，这就无法避免因此而带来的gimbal lock的问题。那么是否因此就要采用quaternion（四元数）来进行旋转计算呢？其实不需要，3D旋转矩阵中拥有众多的对称性（有吗？），我们可以利用这些对称性来完成旋转计算的加速。下面是一个普通的3x3旋转矩阵：

旋转矩阵

ZYZ-(α,β,γ)

我们可以通过ZYZ欧拉角(α,β,γ)的三角旋转矩阵的形式，得到上述矩阵各个元素与角度之间的关系（这个地方给出的结果有误，看起来像是将ZYZ矩阵进行了转置后的结果，不过不妨碍后续的推导与结论）：

旋转矩阵与角度

上述这个公式在β=0的时候是不存在的，此时需要按照之前的格式给出结果如下：

β=0时的SH矩阵

而如果γ=0，我们可以得到如下的简化公式（依然是之前的ZYZ的转置矩阵推导的结果）：

简化公式

使用上述公式能够轻易计算出前面几阶的SH函数的系数，不过对于更加高阶的SH系数，可能就难以为继了，而这正是当前研究中的一个空白，一个痛点。为了能够达成这个目的，就需要一种能够阐明旋转时高阶函数与低阶函数之间的迭代关系的方法，从而实现高阶对低阶的递归调用。关于这种递归关系，目前有三篇文章给出了叙述：

1. Ivanic J and Ruedenberg K, “Rotation Matrices for Real Spherical Harmonics, Direct Determination by Recursion”, J. Phys Chem. A, Vol. 100, 1996, pp 6342-6347. See also “Additions and Corrections: Rotation Matrices for Real Spherical Harmonics”, J. Phys Chem. A, Vol. 102, No.45, 1998, pp 9099-9100
2. Choi, Cheol Ho et al, “Rapid and stable determination of rotation
   matrices between spherical harmonics by direct recursion”, J. Chem.Phys. Vol 111, No. 19, 1999, pp 8825-8831
3. Blanco, Miguel A et al, “Evaluation of the rotation matrices in the
   basis of real spherical harmonics”, J. Molecular Structure (Theochem), 419, 1997, pp19-27

这个GDC的作者实现了Blanco的方法，不过据他介绍Choi的方法实现效率更高（尤其是是在复数领域，不过在复数领域可能会需要进行一次复数到实数的转换）。

使用SH计算diffuse表面光照(SH Lighting Diffuse Surface)

下面介绍使用SH函数来实现光照，这里会给出多种不同的光照技术的实现方案。

现在给定一个场景，场景由多个顶点构成，每个顶点都有自己的法线，通过对顶点-法线进行处理，生成所谓的“光照点”，这些光照点的光照结果会通过raytracer来实现。之后对于每个光照点，会计算其对应的传递函数（transfer function），传递函数的作用是什么呢？从后面的内容来看，传递函数是光在传播过程中经过与场景的交互后对光照的一种影响，比如直接光照中的BRDF就是传递函数的重要组成部分。将传递函数跟跟输入光强函数相乘后积分得到的结果就是这个光照点的输出光照结果。

对于漫反射表面（diffuse surface）而言，总共有三种传递函数可供选择，这三种函数的复杂度是层层递进的。

无阴影漫反射传递（Diffuse Unshadowed Transfer）

如果我们将前面的物理真实的渲染计算公式精简到极致，就变成了如下的一种实现：某个朝向的表面上的某个点，通过直接光照的作用，得到一个无阴影的渲染结果。

无阴影渲染示意图

其中：
L(x,ωo)表示的是点x沿着ωo方向出射的光强；
fr(x, ωo, ωi)表示的是两个方向的BRDF；
Li(x, ωi)表示的是点x沿着方向ωi的入射光强
H(x, ωi)表示的是几何关系（跟之前的描述一致，其实就是光照方向与表面法线的余弦关系）

由于diffuse表面的特性是沿着任意方向的反射光的光强都是相同的，因此可以对上述公式进行简化处理：
1.任意方向输出光强相同，意味着反射光跟观察方向无关，因此可以去除ωo
2.BRDF可以用一个常量来代替

简化结果

其中：
ρx表示的是点x处的表面反射率（albedo）
Nx表示的是x点处的表面法线
DU表示的是Diffuse Unshadowed

这里需要注意的是，此处使用反射率ρx来表示Lambertian Diffuse表面的反射幅度，通常来说，反射率指的是出射radiance跟irradiance的比值，对于漫反射表面而言，这个值通常取值范围为[0, π]，为了简化处理，这里给出的ρx取值范围直接就限定在[0,1]之间，也正好跟我们通常的颜色范围相一致。

将上述公式中的光源部分剔除掉，就得到了传递函数部分M_DU，这个传递函数就是我们准备通过球谐函数来进行模拟近似的部分了：

传递函数部分

蒙特卡罗积分计算SH投影系数

按照蒙特卡罗积分计算SH投影系数公式，计算传递函数的球谐系数的实现代码给出如下：

int red_offset = 0;
int green_offset = 1;
int blue_offset = 2;
// for each sample
for(int i=0; i<n_samples; ++i) 
{
    // calculate cosine term for this sample
    double Hs = DotProduct(sample[i].vec, normal);
    if(Hs > 0.0) 
    {
        // ray inside upper hemisphere so...
        // SH project over all bands into the sum vector
        for(int j=0; j<n_coeff; ++j) 
        {
            //sample[i].coeff[j] represents for the result of j-th Spherical Harmonics
            //Base function takes the i-th sample as input 
            value = Hs * sample[i].coeff[j];
            //ρx for albedo_r/g/b
            result[j + red_offset ] += albedo_red * value;
            result[j + green_offset] += albedo_green * value;
            result[j + blue_offset ] += albedo_blue * value;
        }
    } 
    else 
    {
        // ray not in upper hemisphere
    }
}
// divide the result by probability / number of samples
// area is usually 4π
double factor = area / n_samples;
for(i=0; i<3*n_coeff; ++i) 
{
    coeff[i] = result[i] * factor;
}

应用球谐无阴影漫反射传递函数计算得到的光照效果给出如下图所示，看起来就跟只计算了光照方向与法线方向直接点乘得到的效果类似，不同的是，由于实现方式使用的是球谐传递，因此任意通过球谐函数表达的面光源都能通过简单的计算得到其光照效果：

5阶（25个系数）无阴影漫反射传递实现效果

这里有一个问题，球谐表达的传递函数其系数是固定的，球谐表达的光照也对应于一套常数系数，那么前面说过，两个球谐函数相乘的积分，其结果就相当于系数相乘之后求和，那么按照刚才的说法，传递函数跟光照的积分得到的结果就应该是一个固定值了，怎么表现成光照明暗变化呢？难道前面的猜测都是错误的，其实这些系数是带有变量的？但是如果是这样，又不符合正交基函数的定义了。

之所以会有这个问题，是因为对于传递函数的理解错误，前面给出的传递函数M_DU实际上是单点的传递函数，而非整个场景的传递函数，也就是说上面一段代码计算得到的只是一个采样点（一个顶点）的传递函数的系数，对于一个场景而言，会存在众多的采样点，因此需要对很多采样点（通常是模型上的每个顶点用作一个采样点）分别计算一套系数，在实际使用的时候，计算得到的光照只是此点的光照，不同点的光照效果是不同的，从而才会使得整个场景呈现出明暗变化。

不需要担心对于每个顶点计算一套SH系数会有很大的消耗，因为这个工作不需要光源的参与，可以离线完成。

另一个问题，离线计算的模型的顶点球谐系数，在模型发生了刚体运动（平移或者旋转）之后，是否还依然有效？

从代码给出的计算过程看到，发生平移与旋转后参与计算的各项数值可以看成是不变的（法线虽然变化，但是各个采样点的位置向量也跟随变化，相对而言可以理解成不变的），因此其计算结果应该也是不变的，从而使得在发生刚体运动之后，球谐函数的各个系数依然是有效的。

带阴影漫反射传递(Shadowed Diffuse Transfer)

在前面的无阴影漫反射光照积分公式中添加一项遮挡项（投影物的作用项）之后，就可以得到带阴影的光照公式：

带阴影光照积分公式

其中V(ω_i)是一个开关项，如果方向ω_i上的入射光被当前像素所在物体的其他部分所遮挡，则返回0，否则返回1。

由于光照方向来自四面八方，因此这里有个假设，只有当遮挡物足够近，这个方向上的入射光作用才会接近于0，更进一步假设，只有几何体自身才被认为是比较近的遮挡物（实际情况中当然不是如此，不过这里只是做简化处理罢了），不知这个理解是否存在问题？

前面无阴影漫反射其实是默认几何体是一个无穷大的平面，不会受到任何遮挡物的投影干扰。而此处带阴影的漫反射则开始考虑自身导致的光照遮挡作用，其对应的传递函数（Diffused Shadow Transfer）M_DS可以表示成以下方式：

带阴影漫反射传递

这个遮挡项其实就是平时说的AO（Ambient Occlusion）或者OA（Occluded Ambient），这一项其实就是需要计算当前点在不同方向的遮挡情况，不过这里只需要知道是否遮挡，而不需要知道其具体的遮挡物远近的详细信息，其实现代码给出如下：

int red_offset = 0;
int green_offset = 1;
int blue_offset = 2;
// for each sample
for(int i=0; i<n_samples; ++i) 
{
    // calculate cosine term for this sample
    double Hs = DotProduct(sample[i].vec, normal);
    if(Hs > 0.0) 
    {
        // ray inside upper hemisphere...
        if(self_shadow(pos,sample[i].vec))
          continue;
        // ray inside upper hemisphere so...
        // SH project over all bands into the sum vector
        for(int j=0; j<n_coeff; ++j) 
        {
            //sample[i].coeff[j] represents for the result of j-th Spherical Harmonics
            //Base function takes the i-th sample as input 
                        // Why should here multiply Samples'Coeff?Means The SH of sample is taken as a point light
                        // how to compute sample[i].coeff[j]?Calculate Offline
            value = Hs * sample[i].coeff[j];
            //ρx for albedo_r/g/b
            result[j + red_offset ] += albedo_red * value;
            result[j + green_offset] += albedo_green * value;
            result[j + blue_offset ] += albedo_blue * value;
        }
    } 
    else 
    {
        // ray not in upper hemisphere
    }
}
// divide the result by probability / number of samples
// area is usually 4π
double factor = area / n_samples;
for(i=0; i<3*n_coeff; ++i) 
{
    coeff[i] = result[i] * factor;
}

可以看到，相对于无阴影漫反射，此处代码中增加了一个self_shadow()函数用以表征AO项，这个函数的具体实现取决于使用哪种光线追踪技术，不过在计算这一项的时候有一些tips可以帮助加速：

1.在使用光线追踪技术时候，通常会使用一些空间划分数据结构（比如说hierarchical bounding box或者BSP等）来进行加速，而使用这类结构来计算self_shadow的时候，生成的光线的起点其实是位于几何体bounding box内部的（这个点应该就是被投影点吧），因此不需要额外的计算来搜寻光线的起点（难道在其他的光线追踪中，还需要这样的一个过程？不是很能理解）。

2.self_shadow的计算需要从几何体上的顶点沿着各个方向投射光线，并计算这些光线是否与模型相交从而判断是否被遮挡。不过在这里会存在一些问题，由于顶点是被多个面所共享的，那么就会存在很多光线与这些面所平行，这些光线与当前几何体的交点要么是光线的出射原点或者沿着光线的方向上的某一点，或者存在很多光线是穿过这些面片的，这些光线的交点就是原点，而这样计算得到的碰撞点实际上是会导致阴影结果异常的，因此，为了得到更好的结果，可以在实际实现中保存面片与顶点的邻接信息，在计算self_shadow的时候，将当前顶点相邻的面片从碰撞检测中剔除出去，如下图所示，这样前面提到的那些光线就会与下一个面（背面）碰撞（前提是这个几何体是一个闭合的形状，而不是单一面片）：

self_shadow出射光线在剔除相邻面片后的碰撞结果示意图

3.这里的shadow计算需要在闭合模型上进行，而不能在单面模型上进行，否则会得到错误的结果；此外，模型本身的面片也不能交叉，否则在两个面的交叉线上是没有顶点的，从而导致以顶点作为射线出射点的做法会得到不正确的结果。且self_shadow的碰撞点距离射线起点的远近是没有关系的（比如上图中的B点），其最终结果都是正确的。

带阴影漫反射模型

间接光漫反射（Diffuse Interreflected Transfer）

前面两种模型，入射光都是直接取自光源，下面的间接光漫反射模型，则在前面二者的基础上增加了从其他可见位置反射而来的光线，其渲染公式如下：

间接光漫反射公式

其中L_DS(x)就是前面一节中带阴影漫反射公式的输出结果，V(ω_i)则是前面给出的可见性检测结果，L(x',ω_i)则是新增部分的核心，指的是从当前模型上（注意这个修饰）的某点x'朝着点x反射的光强。

间接光漫反射的传递函数用数学公式就不是很好描述了，且其意义也不是特别直观，不过其计算方式倒是非常简单，总共可以分成四步：

1.对于模型上的任一点x，先利用前面一节的漫反射模型（比如带阴影漫反射）计算出直接光照射的漫反射传递函数。

2.从当前点x发出射线，此射线与模型上另一个面片的交点x'，从x'反射向x的传递函数可以通过对x'所在面片的三个顶点的直接光球谐函数按照质心公式（barycentric）插值得到。

3.将第二步得到的反射数值乘上光线方向与x点法线的点乘结果，并累加到一个初始值为空的球谐函数中，当从x点出发的所有射线的结果都累加完成之后，再将这个部分除以累加点数并乘以对应的权重，得到蒙特卡罗积分形式

4.当模型上的所有特征点都完成上述步骤之后，我们就得到一次反弹的间接光漫反射结果（球谐函数），之后将刚才得到的各个特征点（也就是面片顶点）的球谐表示作为初始的光照强度，开启新一轮的积分模拟，从而得到二次反弹光照模型，循环往复以至于得到N次反弹的光照模型。将直接光（可以看成0次反弹光照模型）以及从1~N次反弹的间接光照模型都加起来，就得到了一个新的球谐函数列表。

从几何角度来看，间接光漫反射模型是很简单的，将当前点可见的其他点的反射光作用通过SH相加的方式计算出来，之所以可以通过球谐相加，是因为球谐表示的坐标系是相同的。这里有个假设：光源模型的球谐表示结果在同一个模型上的各个点可以认为是近似不变的（不用考虑光损等各种复杂因素）。在这个假设下，当前点的反射光强，就不需要考虑不同反射点的光源不同的因素，从而可以直接相加。

这个模型的代码给出如下：

void self_transfer_sh()
{
    const double area = 4.0*PI;
    double *sh_buffer[n_bounces+1]; // list of light bounce buffers.
    // allocate and clear buffers for self transferred light
    sh_buffer[0] = sh_coeff; // already calculated from direct lighting
    for(int i=1; i<=n_bounces; ++i) 
    {
        sh_buffer[i] = new double[n_lighting * 3 * n_coeff];
        memset(sh_buffer[i], 0, n_lighting*3*n_coeff*sizeof(double));
    }
    // for each bounce of light
    for(int bounce=1; bounce<=n_bounces; ++bounce) 
    {
        // loop through all lighting points redistributing self light
        for(int i=0; i<n_lighting; ++i) 
        {
            // find rays that hit self
            bitvector::iterator j;
            int n = 0;
            double u = 0.0, v = 0.0, w = 0.0;
            Face *fptr = 0;
            double sh[3*n_coeff];
            // get the surface albedo of the lighting point.
            double albedo_red = mlist[plist[i].material].kd.x / PI;
            double albedo_green = mlist[plist[i].material].kd.y / PI;
            double albedo_blue = mlist[plist[i].material].kd.z / PI;
            // loop through boolean vector looking for a ray that hits self…
            for(j=hit_self[i].begin(); j!=hit_self[i].end(); ++n,++j) 
            {
                if(*j) 
                {
                    // calc H cosine term about surface normal
                    float Hs = DotProduct(sample[n].vec, plist[i].norm);
                    // if ray inside hemisphere, continue processing.
                    if(Hs > 0.0) 
                    {
                        // trace ray to find tri and (u,v,w) barycentric coords of hit
                        u = v = w = 0.0;
                        fptr = 0;
                        bool ret = raytrace_closest_triangle(plist[i].pos,
                        sample[n].vec,
                        face_ptr, u, v);
                        // if (surprise, surprise) the ray hits something...
                        if(ret) 
                        {
                            // lerp vertex SH vector to get SH at hit point
                            w = 1.0 - (u+v);
                            double *ptr0 = sh_buffer[bounce-1] + face_ptr->vert[0]*3*n_coeff;
                            double *ptr1 = sh_buffer[bounce-1] + face_ptr->vert[1]*3*n_coeff;
                            double *ptr2 = sh_buffer[bounce-1] + face_ptr->vert[2]*3*n_coeff;
                            for(int k=0; k<3*n_coeff; ++k) 
                            {
                                sh[k] = u*(*ptr0++) + v*(*ptr1++) + w*(*ptr2++);
                            }
                            // sum reflected SH light for this vertex
                            for(k=0; k<n_coeff; ++k) 
                            {
                                sh_buffer[bounce][i*3*n_coeff + k+0*n_coeff] += albedo_red * Hs * sh[k+0*n_coeff];
                                sh_buffer[bounce][i*3*n_coeff + k+1*n_coeff] += albedo_green * Hs * sh[k+1*n_coeff];
                                sh_buffer[bounce][i*3*n_coeff + k+2*n_coeff] += albedo_blue * Hs * sh[k+2*n_coeff];
                            }
                        } // ray test
                    } // hemisphere test
                } // hit self bit is true
            } // loop for bool vector
        } // each lighting point
        // divide through by n_samples
        const double factor = area / n_samples;
        double *ptr = sh_buffer[bounce];
        for(int j=0; j<n_lighting * 3 * n_coeff; ++j)
            *ptr++ *= factor;
    } // loop over all bounces
    // sum all bounces of self transferred light back into sh_coeff
    for(i=1; i<=n_bounces; ++i) 
    {
        double *ptra = sh_buffer[0];
        double *ptrb = sh_buffer[i];
        for(int j=0; j<n_lighting * 3 * n_coeff; ++j)
            *ptra++ += *ptrb++;
    }
    // deallocate SH buffers
    for(i=1; i<=n_bounces; ++i) 
    {
        delete[] sh_buffer[i];
    }
    return;
}

上述的这种算法就是所谓的概率累加算法（probabilistic gathering solution）：对于模型上的每个顶点，将此位置可见的所有顶点的反射效果累加起来。如果每次迭代只考虑一个顶点的反射结果的话，那么通过这种算法，我们可以看到当前的点是如何被一点点的点亮的。除了这种算法之外，还有所谓的non-probabilistic solution以及shooting solution，都有各自的优缺点，此处就不再赘述，有兴趣请自行查阅相关资料。

相对于self_shadow对于模型闭合性的要求而言，间接光漫反射模型对于模型的要求更高。如果在计算某个处于孔洞中的顶点得到的漫反射结果存在比较大的问题，那么可以通过将此顶点沿着其邻面的法线的平均值往上提一个epsilon的范围来快速修复。

间接光漫反射模型效果

漫反射表面SH渲染

经过上面的处理之后，场景中的每个顶点都有了一套SH系数，那么其光照输出按照前面的表述，只需要通过乘法运算就能得到：

相乘的积分等于SH系数的乘法

上述公式得到的是顶点的光照结果，整个几何体面片上的其他像素光照可以通过插值得到。此外需要注意的是，上述计算过程是在模型空间中完成的，如果需要重新调整模型的朝向（旋转模型）或者对光照传递函数进行旋转操作，如果想要保证光照效果的正确性，就需要先将光照转换到模型空间中。

为什么需要先将光照转换到模型空间中？因为模型上各顶点的SH是在模型空间计算得到的，如果不将光照转换到模型空间中计算，那么模型转换之后得到的结果就跟旋转之前得到的结果无区别，这显然是不符合实际情况的。

如果考虑到光照的颜色，那么每个顶点的输出光照应该分成三个通道计算其光照亮度：

for(int j=0; j<n_coeff; ++j) 
{
  vertex[i].red += light_red[j] * vertex[i].sh_red[j];
  vertex[i].green += light_green[j]* vertex[i].sh_green[j];
  vertex[i].blue += light_blue[j] * vertex[i].sh_blue[j];
}

上面的公式以及代码片段中都没有用到法线数据，是因为法线已经集成在各个顶点的SH中了，同样，AO数据也一并集成在SH中，所以上述计算得到的结果会看到AO的效果。

在实际应用中，如果对于阴影的柔和度要求不高的话（比如角色面部），那么可以将上面的3阶（9个系数）的SH替换成2阶（4个系数）的版本。这种做法可以节省内存与计算资源的消耗，同时还可以通过一个SIMD指令完成相关计算。

另一个优化就是将彩色的光照进行简化处理，比如上面的代码就可以改写成如下的样子,此时，对于光照与顶点的SH数据，都只需要存储一份（而不需要存储RGB三份）就可以了（只是不知道效果怎么样）：

for(int j=0; j<n_coeff; ++j) 
{
  vertex[i].red += k_red * light[j] * vertex[i].sh[j];
  vertex[i].green += k_green * light[j] * vertex[i].sh[j];
  vertex[i].blue += k_blue * light[j] * vertex[i].sh[j];
}

由于光照之间的叠加运算是线性的，因此我们可以将不同的算法得到的不同的光照部分结果累加起来，对于不同的光照部分采用最优的算法进行计算，比如我们可以球谐来计算Ambient与Diffuse（这部分的SH计算非常简单，但是效果又非常好），并通过其他算法（比如预生成环境贴图的方式）计算出Specular部分，之后将二者通过Alpha Blend（能量守恒原则，Specular与Diffuse等部分的能量之和是有限的，不超出输入光亮度，因此会呈现一个此消彼长的趋势）累加起来，就能得到兼顾各种光照效果的结果，比如开篇之初给出的汽车效果，就是通过这种方式计算得到的：

SH-Diffuse + Prefiltered Specular

到目前为止，我们上面的所有计算都是建立在光源距离几何表面无穷远的假设上的，对于局部光源（比较近的如点光聚光灯等）而言，由于其不满足模型上的所有点接收到的光照区别较小的假设，因此不能使用间接光漫反射传递模型，不过我们可以一种只考虑self-shadow而不考虑self-transfer（内部光线反射）的方式来进行模拟。

局部光源的模拟是通过从模型表面上一些空间分布位置提前设计过得的点为局部光照环境构建出少数（比如6~8个）采样点，之后就对这些采样点进行光照计算，并完全忽略模型的影响。

模型上各个顶点的局部光照可以通过对前面的采样点的光照进行加权累加来求得，权重与顶点到采样点的距离有关。Sloan， Kautz以及Snyder的球谐公式使用的权重公式如下所示：

权重公式

i指的是第i个顶点，j指的是第j个采样点，n在Sloan的文章中取值为10.

除了上面的模拟方法之外，还有人给出了其他实现更快消耗更低的方法。比如我们可以通过如下的方式计算各个顶点的光照：
1.跟前面一样，选择出N个采样点，并提前计算出这N个点的光照
2.将N个采样点组合成一个由三角面片组成的多面体
3.根据到上述多面体的哪个面更近，将当前模型的各个顶点分割到不同的Group中
4.将每个Group中的每个顶点投影到对应的三角面片平面上，根据质心公式计算三角面片的权重，并根据权重对三角面片三个顶点的光照进行加成累加，将结果作为模型上的当前顶点的光照输出。
这种方法还有一个优点是可以通过硬件的光栅化来完成计算。

关于根据SH重构光源的作用，本文只给出了一些浅显的表述，以此寄希望于读者能够迸现出更多的灵感与思想火花。当前遭遇的最大最有价值的挑战在于实现带动作的模型的SH光照，对于前面的三种传递模型中，只有最简单的无阴影漫反射传递模型可以通过对多个SH进行线性插值实现，而带阴影的传递模型由于存在太多的不确定性，使用线性插值的表现并不好，后续还需要更多的研究来突破这个限制。

Creating Lighting Source

这里先给出一个阴天模型，其光强给出如下：

阴天模型

其中：
Lz指的是天空正上方的天空照明，其单位用kcd/m2表示（千坎德拉每平方米）
β指的是Z跟P之间的夹角（与中心点连线的夹角）,弧度。

阴天模型光照强度分布

CIE给出的晴朗天气时的天空模型就相对来说要复杂一些，因为需要考虑太阳所在的位置，以及大气对于阳光的散射：

晴天模型

其中：
L(θ ,ϕ)指的是球面上在(θ ,ϕ)角度的点P的光照强度（kcd/m2）.
Lz跟前面的阴天模型一样，指的是中心正上方的天光强度(kcd/m2)
θ指的是Z跟P的夹角，弧度
ϕ 指的是P点跟太阳所在位置投影到地平面的夹角，弧度
S指的是太阳跟Z点的夹角，弧度
γ指的是太阳跟P点的平面夹角，弧度

使用上述公式模拟得到的光源实施效果非常逼真：

晴天模型光强分布

CIE给出的有云天气模型公式给出如下：

有云天气模型

相关变量的含义跟之前晴天模型的一致，其光强分布模拟图如下：

有云天气光强分布

另一种可以用来表示光照的模型是由Paul Debevec带火的HDR Lighting Probe模型。HDR贴图是使用专门的相机经过标定与多次曝光而采集到的现实世界真实光强能量数值的浮点数据位图。通过自由软件HDRShop，可以为HDR场景生成一组HDR角度投影贴图（angle map projection，这个是啥，看起来像是cubemap？），如下图所示：

angle map projection示例

angle map projection不管是通过球面坐标还是笛卡尔坐标，访问都非常方便。可以将其直接集成到SH投影（求取各个基的系数）函数中来求取光照函数的SH系数，其代码给出如下：

typedef Vector3d (*SH_vector_fn_rgb)(float dx, float dy, float dz);
Vector3d hdr_lightsource(Vector3d *hdr_image, int image_size, float dx,float dy, float dz)
{
    // assume angle map projection
    const float one_over_pi = 1.0f / PI;
    float invl = 1.0f / sqrtf(dx*dx+dy*dy);
    float r = one_over_pi * acosf(dz) * invl;
    float u = dx * r; // -1..1
    float v = dy * r; // -1..1
    // map to pixel coordinates
    int x = int(u * image_size + image_size) >> 1;
    int y = int(v * image_size + image_size) >> 1;
    // return the float RGB value at (x,y)
    return hdr_image[y*image_size + x];
}

void SH_project_vector_function_rgb( SH_vector_fn_rgb fn,
int n_bands,
int sqrt_n_samples,
const SHSample sh_samples[],
Vector3d result[] )
{
    int n_coeff = n_bands*n_bands;
    int n_samples = sqrt_n_samples*sqrt_n_samples;
    const double area = 4.0*PI;
    // for each sample
    for(int i=0; i<n_samples; ++i) 
    {
        Vector3d color = fn(sh_samples[i].vec.x,
        sh_samples[i].vec.y,
        sh_samples[i].vec.z);
        
        for(int n=0; n<n_coeff; ++n) 
        {
            result[n] += color * float(sh_samples[i].coeff[n]);
        }
    }
    // divide the result by number of samples
    double factor = area / n_samples;
    for(i=0; i<n_coeff; ++i) 
    {
        result[i] = result[i] * factor; // NOTE: vector-scalar multiply
    }
}

使用angular map projection的一个不足之处在于，对于HDR的曝光来说，没有一种标准的归一化方案，因此在实现的时候，必须给出一套调整曝光度的设置方案，方便随时根据需要对曝光度进行调整。

通过上述方法产生的是完全的合成光源，其表达的方式为SH系数。通过对球面函数f(theta, phi)进行符号积分，可以得到如下的结果（在theta小于t的时候，其结果为1）：

光照结果

从而得到上图所示的一个Z轴对称的圆曲面光源，而对这个光源，可以通过符号积分，得到此光源SH投影的一个分析解：

光源SH系数公式

根据这个公式就可以直接生成对应的光源了，唯一需要做的就是将光源旋转到我们需要的方向（可以根据本文给出的光源旋转代码实现），并在最终输出的时候与其他光源的作用一起累加进行输出：

inline double is_positive(double x) { return x>0 ? 1.0 : 0.0; }
void synth_light(double cutoff, double sh[])
{
    // clear all values to 0.0
    memset(sh, 0, 16*sizeof(double));
    // symbolic integration automatically generated by Maple.
    double t3 = is_positive(-cutoff + PI);
    double t5 = cos(cutoff);
    double t6 = t3*t5;
    double t9 = is_positive(-cutoff);
    double t11 = t9*t5;
    double t14 = sin(cutoff);
    double t15 = t14*t14;
    double t16 = t3*t15;
    double t18 = t9*t15;
    double t21 = t5*t5;
    double t22 = t21*t5;
    double t31 = t21*t21;
    sh[0] = 3.544907702 - 1.772453851*t3 - 1.772453851*t6 –
    1.772453851*t9 + 1.772453851*t11;
    sh[2] = 1.534990062*t16 - 1.534990062*t18;
    sh[6] = -1.98166365*t3*t22 + 1.98166365*t6 + 1.98166365*t9*t22 –
    1.98166365*t11;
    sh[12] = 2.930920062*t3 - 2.930920062*t3*t31 - 3.517104075*t16 -
    2.930920062*t9 + 2.930920062*t9*t31 + 3.517104075*t18;
}

LionHead公司的Alex Evans提出可以使用相似的方法为SH光源生成假影（proxy shadow）：
1.为每个模型生成对应的光照函数方程
2.根据光照方程，计算出模型上的阴影区域（跟前面计算光照区域刚好相反）
3.将其他模型转换到当前模型空间，并计算出从两者的连线转换到Z轴（为什么是Z轴？）的旋转矩阵
4.将模型的光照函数应用上述的旋转矩阵（这样做的意义何在呢）
5.将旋转后的光照函数方程与光源SH相乘，就得到了当前模型在这个方向下的光强（这里有两种方法可以实现这个相乘计算，一种方式是通过循环卷积circular convolution，另一种是通过符号积分生成一个分析式的传递矩阵，不管使用哪种方式，都需要对光照函数进行旋转，这样才能保证阴影的方向是沿着Z轴的，这种乘法完成之后，再将光照函数旋转回去，此时的阴影就正好处于两者的连线方向上，换句话说，当前物体的投影正好落在其他物体上，此时再继续进行物体的光照计算，总的来说，每个物体每进行一次投影，就需要进行两次光照旋转操作）
6.将上述结果减去此方向上的光源光强，就得到了此方向上的阴影数据。

这个过程比较抽象，还不是能完全理解，看看后面能不能逐步明白

总结

本文给出了球谐光照是如何实现静态与动态物件的真实的光影效果的技术阐述与说明。介绍了如何实现光照在静态模型上是如何完成模型内部反射的预计算的，以及如何为漫反射阴影生成对应的传递函数。

球谐光照可以用于替代静态物体的光照实时计算，其中有许多可选的方案可以实现实时的全局光与局部光照明。

附录

附录给出了如何实现SH旋转的快速实现方案，有兴趣的同学可以自行翻阅。