【GDC2008】Stupid Spherical Harmon

作者: 离原春草 | 来源:发表于2021-08-15 23:48 被阅读0次

今天介绍的是GDC2008上微软的Peter-Pike Sloan分享的关于SH（球谐）的相关使用技巧，原文PPT与相应的论文链接在参考文献中给出。

整篇文章主要包括上述四个部分：

第一部分是SH相关的概念与背景介绍
第二部分是如何使用SH来对光源进行表达，以及给出SH的数学分析模型
第三部分是关于SH常见的Ringing Artifacts的，会着重介绍产生这个问题的原因以及相应的解决方案
第四个部分则会介绍SH要如何使用

image.png

首先需要知道的是，什么是调和函数？所谓的调和函数，实际上是拉普拉斯方程等于0的解，其物理意义也就是此函数的梯度（gradients，某个标量场的各点切线的方向）的散度（divergence，向量场中各点向量的散度表示的是按照某种公式将向量换算后得到的标量，比如可以用向量强度作为某种散度）为0，在笛卡尔坐标系下用公式可以表示为：
$\Delta f = 0$

散度为0在不同的学科中有不同的意义：

在流体动力学中，散度为0表示的是流体的flow（流动）是irrotational & uncompressible的（不了解流体动力学，不清楚具体指的是什么……）
电磁场与重力场也常常会用到这个概念，具体含义这里没说。

单位圆（unit circle）上散度为0的解是傅里叶基函数（fourier basis function），而在球面上散度为0的解则是球谐（sphere harmonics，简称SH），且解的角度部分是以球面坐标西来表示的，而球谐同时还是量子化学与量子物理中角动量算子（Angular Momentum Operator）的特征函数（eigenfunction）。

上面这些内容看起来可能会有点吃力，不过没关系，这里只是进行简单的背景普及，接下来我们只是会将SH来用于对球面函数（spherical function）的表达，所以这些属性即使不清楚也没关系。

SH在图形学中的应用已经十分广泛了，比如Kajiya在84年的论文中介绍了其使用SH来完成Cloud的散射的计算等。

Ravi在Siggraph 2001上提出了将SH用于实时天光Diffuse计算的方法（以下简称Ravi 01），极大的提升了程序执行的性能。

而本文的作者Sloan在Siggraph 2002上则给出了一个使用SH来对物体表面的“Transfer Function”进行预计算（Precomputed Transfer Function，简称PRT）的方法（以下简称Sloan 02），通过这个方法，可以在实时完成一些复杂的全局光照效果（比如软影，inter-reflections等），在这之后的绝大部分SH应用都是从这两篇论文中衍生出来的。

Ravi 01的算法后面被扩展用于计算从一些特定表面输出的反射效果，而Sloan 02的算法则被扩展用于处理一个比较大范围材质的反射属性，且可以通过添加次表面散射以及对结果进行压缩来提升效果的真实度。

PRT算法还被扩展用于模拟贴图细节以及法线贴图等应用场景。

除了上述扩展应用之外，SH的应用场景包括对SH梯度进行模拟等（左下角那一系列参考文章）。

SH除了用在预计算中，同时还可以用于实时运算，右下角的一系列参考文章就是采用SH来实现对动态场景的阴影以及inter-reflection进行计算与模拟的。

SH是一个很强大的工具，但不是万能的工具，在如下的一些应用场合，SH就不太合适：

高频信号模拟，高频信号需要高阶SH，这会导致性能急剧下降
全局效果支持，全局效果往往意味着高频信号
使用其他基函数具有更好表达效果的情况

对符号进行一下约定。

image.png

上图给出了三阶SH所对应的9个基函数的形状示意图，其中颜色表示正负（红色正数、蓝色负数），到球心的距离表示幅度。

后面SH的基函数用 $y_l^m$ 来表示，其中 $l$ 表示的是band level（阶数），而 $m$ 表示的则是每个band中的序号。

每阶包含了 $2l+1$ 个基函数，前 $N$ 阶的基函数总和为 $N^2$ 。

如上图所示，球谐函数常常会用球坐标来表示，这是绝大部分教科书以及网站上常用的方式，这种表达方式在进行符号运算（求取解析解）的时候很方便，但是在写代码的时候就不那么好用了。

上面给出的是笛卡尔坐标系下的表示方式，在这个坐标系下，SH就是简单的二次多项式，因此对于代码编写来说会友好很多。

在所有的SH函数中，我们需要关注其中一些比较特殊的，比如上图中用线框标注的，处于每个band中间位置的SH函数，这些函数我们称之为Zonal Harmonics（简称ZH，下同），ZH在散射领域中经常会被用来表示相函数，ZH的笛卡尔坐标表示只受Z轴影响，因此任意投影到SH空间后沿着Z轴圆形对称（即在Z值确定后，其结果也就确定了，不再受XY影响）的函数，将只会在ZH基函数上具有非0系数。

有时候，出于表达简洁的需要，我们会将二维SH转换到一维空间。另外需要注意的是，SH的基函数都是正交的，也就是说任意两个基函数相乘之后的积分，只有在这两个基函数是同一个基函数的情况下才为1，否则为0。

正交基函数的两个重要特性：

任意函数在某个基函数下的投影系数就等于这个函数与基函数相乘之后的积分结果
两个函数相乘之后的积分结果就等于这两个函数的SH系数组成的向量的点积

这里将某个cubemap表示的light probe投影到不同band的SH上，从结果我们可以看到，band越小，投影结果越模糊。

图形学中对光照的计算时通过将输入光与BRDF相乘之后积分来实现的，而BRDF是一个十分平滑的函数，这也就意味着不论是输入光源还是BRDF，我们都只需要一些低阶的SH就能表达，不必担心高阶SH导致的高额计算消耗。

这里展示了一个light probe通过SH Projection对一个球进行照明的demo。

SH照明基函数的一个重要特性就是旋转不变性，也就是输入光源旋转一定角度，其通过SH表达的近似解也是同样旋转一个相同的角度，这也就说明，光照结果在旋转后的SH表达是连续的稳定的，不会因为光照的旋转而导致输出光效的抖动。

这里给出的是SH基函数与SR基函数（sphere radial basis function，这是半条命2中用过的ambient cube方案，相对于SH方案，这个方案具有众多特点，比如内存占用低，计算消耗低等），要对一个cube进行照明的话，我们需要在cube的每个面对应的位置各摆上一个ambient cube（SR基函数方案），每个基函数会随着到中心点的距离成平方衰减（cosine squared falling off，具体可能不是与距离的平方成反比，不过不影响理解，后面有机会再补上相关细节），在这种表示方法下，如果输入光源正好指向的是某个ambient cube的中心点，其结果看起来是正常的，但是当光照方向在多个ambient cube之间转动时，就会出现很明显的能量损耗，而在这种情况下，SH的表现却是一如既往的正常。

从上图我们可以看到，SH旋转矩阵的结构是十分直观的，因为每个band的旋转都是互不相干的，因此整个SH的旋转结果应该是如上图所示的由多个对角线block组成。

右边的一列是笛卡尔坐标系下的SH多项式，可以看到一次函数（X/Y/Z）只需要3x3的矩阵即可，而二次函数则需要5x5的矩阵。

而由于其对角线block结构，整个旋转的计算复杂度为 $O(n^3)$ 其中n表示的是SH的阶数。

相对于SH，对ZH的旋转要简单的多，由于每阶只包含一个非零系数，因此在矩阵中只需要关心一列数据即可，因此其旋转复杂度为 $O(n^2)$ 。

在介绍如何计算这一列数据之前，先来介绍一下上面图片中的一些前置概念。

Delta Function，这个在中文中叫脉冲函数，指的是只在自变量为0的时候具有正无穷的数值，其他时候其数值均为0，这个函数的特征是，跟任何函数相乘之后的基本，都等于这个函数在0位置处的值。

因此，如果我们想要计算一个沿着Z轴圆周对称的脉冲函数在XY = （0， 0）位置处的SH投影系数，就可以通过乘法积分进行计算，最终得到的结果为 $y_{l0}(z)$ ，如果取z=1，那么根据前面的ZH多项式就可以得到 $y_{l0}(1) = \sqrt{\frac{2l+1}{4\pi}}$ （l=2的时候好像不对？分母应该是16而不是4？不是的，这个结果是正确的，还要考虑 $3z^2-1$ 的赋值结果），而这是一个常数 $d_l$ 。

由于上述脉冲函数是沿着Z轴圆周对称的，因此对其的旋转就只涉及到ZH函数部分，也就是只有ZH项所对应的那一列数据，下面来看一下这一列数据要如何求取。

假设经过某个旋转之后，将Z轴旋转到了 $Z '$ 轴，同样，这个旋转会对脉冲函数造成影响，而从前面关于ZH矩阵下的旋转结构的内容中我们知道，对某个通过ZH表达的函数的旋转就相当于对其每一个基函数乘上一个旋转矩阵的中间一列 $C_{lm}$ ，而脉冲函数乘上这个数值后，在（0， 0， 1）位置处投影的某个band数值就变成了 $C_{lm}d_l$ （函数的SH表达进行一定角度的旋转），而根据旋转不变性，某个函数旋转后的SH表达就等于其SH表达进行同角度的旋转，也就是说，这一项就等于脉冲函数旋转后的SH表达，这个部分可以用 $y_{lm}(z')$ 来表示（原来脉冲函数在Z=(0, 0, 1)方向时的数值等于旋转到Z=z'方向上的数值），据此就可以计算出 $C_{lm}$ 的数值了。

得到ZH旋转矩阵中间列的表达之后，我们就可以将这个旋转矩阵用在任意的沿着Z轴圆周对称的函数上了，这里我们将这个函数用g(s)表示。

通过左侧下方的积分，我们可以计算出其在ZH上的对应系数，这些系数组成了一个列向量。

之后，我们将g(s)从z旋转到z'，得到了新的函数g'(s)，如右上角的图形所示，而这个函数的SH系数则可以通过两步完成：

先计算出SH基函数在z'方向上的数值y(z')（这个数值应该是每个基函数都有一个）
乘上对角矩阵G*

对角矩阵可以通过上一页PPT最后的公式计算出来，根据前面的公式，对角矩阵上的缩放因子有
$C_{lm} = y(z')/d_l$
也就是说，旋转矩阵的缩放因子就等于旋转后的z方向向量在SH上的取值，而由于我们目前要求取的是从z'到z的旋转矩阵G*，这个时候的z'就等于(0, 0, 1)，而此时的y(z')就等于 $g_l$ ，因此对角矩阵就等于 $g_l/d_l$ 。

总结一下，一个沿着Z轴圆周对称的函数g，可以通过ZH来表达，假设其在ZH上的投影系数为 $z_l$ ，那么将这个函数从Z轴对称旋转到沿着d方向对称，旋转后的函数叫f，这个f不再是可以通过ZH表达，但是依然可以通过SH来表达，此时其在SH上的投影系数有：

$f_l^m=\sqrt{\frac{4\pi}{2l+1}}z_l y_l^m(d)=\frac{d_l}{z_l} \cdot y_l^m(d)$

这个系数分成两部分，前面一部分 $\frac{d_l}{z_l}$ 含义比较明确，这里不做解释，后面一部分 $y_l^m(d)$ 表示的是什么意思呢？从含义上来说，这一项就是对SH基函数的自变量赋值d的时候的函数值。

关于卷积

$f(t) * h(t) = \int_{-\infty}^{\infty}f(\tau)h(t-\tau)d \tau= \int_{-\infty}^{\infty}f(t-\tau)h(\tau)d \tau$

这里有一个性质，即ZH跟SH的卷积结果是一个卷积（而一个一般的SH跟一个一般的SH的卷积却没有这个特性（得到的是一个manifold in SO(2) Rotation Group(2D Rotation Group，通常用SO(2)表示，这是机械以及几何学中的一个概念，指的是二维欧式空间中关于原点的所有旋转所组成的Group，同理还有对应的SO(3)等)，而这只是一个圆周对称的函数）），而这个卷积就等于在原始SH函数上乘以对应band的ZH函数之后再乘以一个常量（ $\sqrt{\frac{4\pi}{2l+1}}$ ，跟之前旋转ZH所使用的常量 $d_l$ 互为倒数），其实这个公式描述的是一个可以用ZH表达的函数h与一个用SH表达的函数f在卷积之后的结果依然可以用SH来表达，且在各个SH基函数上的投影系数 $(h*f)_l^m$ 可以通过 $h_l^0*f_l^m$ 乘上一个常量来求得。

这个公式的名字叫做Funk-Hecke定理，上述链接来自于论文Lambertian Reflectance and Linear Subspaces，论文中还介绍了不少的SH相关理论，比如对于一个函数的旋转，不会导致某一阶的SH能量损失，只会使得能量在同一阶中的不同基函数上的转移，而不同阶的SH之间不存在能量转移（从对角旋转矩阵也可以感知一二）。

右下角给出了ZH跟SH卷积后的结果示意图，而卷积的物理意义其实就是对于在原始SH函数f上的任意点而言，都相当于是沿着这个点对应的方向（从球心出发）使用一个kernel函数（ZH函数h）对其进行积分，结果相当于是对SH函数做了模糊处理（低通滤波）。

Irradiance Map的计算实际上就可以使用这项特性，其中kernel函数就是一个截断的cosine函数（相当于 $max(0, n \cdot L)$ ），而原始SH函数f则是远距天光，而这里的kernel函数同样是一个低通滤波器，因此只有低频天光信息被保留了。

上图给出的是截断consine函数的示意图，左图是取 $\phi=0$ 绘制的 $\theta$ 作为自变量的函数曲线，右图则是极坐标下的示意图。

其中黑色曲线是原始的截断cosine函数，而红色的则是使用三阶SH近似模拟的截断cosine函数，蓝色使用的是5阶SH近似模拟的效果（如果想要使用HDR光源的话，3阶SH精度就不太够，通常会用5阶SH（第4阶的ZH系数是0，可以直接跳过）），虽然精度有所提升，但是误差依然存在。

下面来介绍一下lighting部分。

上面两图给出了SH系数计算的伪代码，从代码结构来看，计算过程是通过蒙特卡罗积分完成的，积分采用的概率函数也就是权重fWt的倒数，u跟v表示的是cubemap上这个face上的uv坐标（这里的uv坐标不是从-1到1，而是会考虑z值的影响，比如在+X face上，就是y/z），PPT给的解释这里的fWt表示的是差分立体角（也就是单个像素对应的球面积，可以理解为积分公式的ds部分？）。

EvalSHBasis函数用于计算texel位置处的y(s)的值（也就是SH基函数在texel对应方向上的数值，这里的方向需要归一化到单位向量），t(texel)指的应该就是通过对cubemap采样得到的light数据，也就是f(s)的值了。

累加完成后，按照蒙特卡罗积分，还要除以采样数N，不过PPT中说到，在cubemap分辨率较低的时候，采样数较小，导致计算结果精度有较大偏差，而使用归一化后的权重作为被除数得到的结果会更好一点。

这里给出了使用投影后的SH基函数表达的cubemap效果，右下角是原始cubemap，前面对应的是按照顺序排列的基函数乘上对应系数后的效果。

上面介绍的是天光这种没有办法给出解析式的光源的SH投影计算方式，下面来看下具有解析式的光源的投影系数计算方法。对于一些全局光源而言（比如方向光，比如前面的天光），其求取出来的SH投影系数是全场景公用的，而对于一些局部光源而言，求取的SH投影系数则是针对场景中某个特定点的（比如某个模型上的某个顶点），不同位置的SH投影系数是不能共用的，当然，距离较近的点，其SH系数相对而言是比较接近的。

方向光的投影最简单，只需要求取SH基函数在对应方向上的数值即可（基函数在对应方向上的数值只能表示这一个方向上的SH表达结果，其他方向呢？我们需要的是SH基函数的对应系数，又不是SH基函数在某一点的取值啊？）

球形光源的投影系数（这里指的是某个特定点P上的SH投影系数）计算可以通过ZH（因为球形光源肯定是相对于从当前被照明点到球心的连线方向是圆周对称的，不过这个方向不一定是Z轴(0, 0, 1)，按照前面的说法，可以通过乘上一个系数实现转换）来完成，公式给出如上图所示（这个公式是在将Z轴与被照明点到球心连线重叠的情况下给出的，如果不重叠，公式还要做其他处理）。

对于某个点P而言，假设光源中心点在C，半径为r，P到C的距离为d。那么光源对P点的radiance是多少呢？我们可以很容易求出球形光源对于点P的张角 $\theta$ ，那么我们在求取SH系数的时候只需要将积分范围限定在这个角度内部即可，也就是说：

$z_l = \int_{\theta = 0}^a\int_{\phi=0}^{2\pi}y_l^0d\phi\d\theta$

其中a是半个张角（ $sin(a) = \frac{r}{d}$ ，本来应该是要对着P点的整个球形区域进行积分的，之所以收缩到半个张角，是因为超出这个张角的部分数值为0）

这里比较奇怪的是，正常来说，SH系数的求取应该是f跟y相乘之后的积分，这里将f省去了，这里面做了什么处理呢？此外，这个积分除了张角之外，没有任何与球形光源强度相关的输入，难道强度为1跟强度为10得到的投影系数是完全相同的吗，这显然是不合理的，不过PPT跟论文都没有关于这一点的详细解释，因此这里暂时没有办法解决这个疑惑，有其他同学了解背后的原理的还请在评论区帮忙答疑，万分感谢。

一个圆锥光源（cone light source）对于椎顶的ZH系数也可以采用同样的方法计算得到，不过这里就不需要手动求取张角了，这个角度在圆锥定义的时候已经显式给出了。

使用ZH表达的椎体光源与原始光源效果存在较大差距，如上图右上角的图形所示，即使使用6阶的ZH表达与原始的半球光源还是存在较大的误差。

椎体光源通常会在张角之外带上一层衰减层，随着角度增加，光强会逐渐衰减到0，日常中常用的是三次衰减公式（如左侧曲线所示），而这样一种带有软边缘的椎体光源其ZH表达式在模拟原始信号的精度上要比硬边缘的椎体光源ZH表达要高很多。

这里给出了不同尺寸的光源的投影效果。

对于一个沿着某个方向圆周对称的光源而言，计算光源的SH系数可以分成如下两步：

计算出此光源的ZH投影系数，投影系数通常可以用光源的一些参数来表示（比如前面球形光源的张角之类）
将ZH系数按照前面的旋转公式转换成SH投影系数

如果光源的光强是使用[0, 1]之间的浮点数来表达，那么我们可以对反射光强进行归一化，归一化的目的是为了保证能量守恒。

通常来说，当光源到某个点的路径未被其他物体遮挡，那么如果这个点的法线正好指向光源，那么反射的radiance就应该是1。

PPT中说到这里需要使用截断cosine函数的ZH表达（上图中的 $z_l$ ）来计算归一化因子c，而从论文中推断，上图中的 $h_z$ 应该是光源的ZH表达，那么这个归一化式子就相当于：
$\int_{\Omega} L(\omega)cos(\theta)d\omega * c = 1$

也就是光源各个方向上的总体光强输入为1，其照射到某个点后的反射光强应该要等于输入光强1，所以c就等于上式中积分的倒数，而如果将光源跟cosine项都用ZH表达，那么就可以将积分转换为点乘了。

前面介绍了如何用SH来表示光源，下面来看下如何从SH表达的光源中解析出原始光源的强度信息，而之所以会有这种需求，有如下几个原因：

一些老式硬件不支持irradiance environment map（SH表达的光源已经不再需要irradiance environment map了吧，这个从何说起？），因此需要从SH中解析出方向光+环境光来（这两种光源已经能够应对绝大部分的光照效果）
对于一些光滑材质而言，我们需要计算高频的specular反射，而SH表达的光源都只有低频信息
我们有时候需要通过SH来移除一些高频信息，从而得到一些更为柔和的光源效果。

上图中的E(c)是方向光强度c求取时的误差函数，其中积分是针对所有法线进行的，积分函数是一个差值平方，差值的前一项指的是方向光打在法线N上的反射光部分，而后一项是一个积分，表示的是外部环境光源入射到法线为N的面片上的反射光强。

由于外部反射光强都来自于方向光，因此理论上这两者应该完全一致，而这也就构成了误差函数的理论基础。

而如果将前面的误差公式用SH表达式来翻译一下就会更清楚一点，如上图右上角的公式所示，其中 $L_e$ 表示的是irradiance environment map的SH表达，而 $L_d$ 则是方向光的SH表达。

而这个误差公式实际上就是两个函数相乘之后的积分，而这种情况实际上可以用这个函数的SH表达的系数相乘来简化计算。

这是一个二次公式，对这个简化后的公式求导数，导数为0的c就是误差函数最小时候的c了。

球谐系数中的第一个系数 $f_0^0$ 表示的是函数f在球面投影的平均值，这个值有一个专业的术语叫做DC项。

DC项是专门供环境光使用的，也就是说环境光相关的SH系数只有DC项，而DC项也只影响环境光。

因此如果我们希望解算在环境光照明下的方向光，只需要先将DC项移除之后就可以按照之前的方式计算方向光数据了，在计算完方向光之后，就可以进行环境光的计算了，得到结果如上图所示。

拿到方向光强度之后，下一步就是计算方向光的方向。而最匹配的方向的常用的一种计算方式就是使用Optimal Linear方法，这个方法就是使用线性SH系数（一阶SH系数）添加上对应的符号，从中挑选出球面上亮度最高的一个线性SH系数向量作为方向光的方向。

文章还介绍了如何求取多个方向光叠加情况下的光照方向，不过也只是一带而过，这里就不介绍了，下面介绍一下SH表达方法中的一些问题。

Ringing也叫Gibbs现象，是信号处理中十分常见的问题，当一个在空间域不连续的函数投影到有限阶数的傅里叶基函数上时，最终组合出来的近似模拟结果就会出现Ringing问题（之所以会有这个表现，是因为空间上不连续的信号实际上在频域上对应的就是高频信号，而通过傅里叶投影或者SH投影得到系数，通常只会记录其低频的几阶数据，而高频信息则是丢弃的，因此导致重建回来的信号跟原始信号相比，是缺失了高频信号的，从而导致震荡，当然，随着投影保留的数据阶数的增加，信号也越来越接近原始信号，但是由于不连续导致的信号频率是无穷的，因此还是无法完全模拟原始信号），具体示例下面会说。

一些本身不连续的函数（比如在取值上存在突变的阶跃函数）是无法使用一些有限的平滑的函数按照线性组合来模拟的。

如果将一个阶跃信号投影到傅里叶基函数中，会产生两个问题：1. 在不连续的位置产生一个超出原信号幅度的波峰，2. 则是在那之后的连续的震荡。

而这两个问题就是Ringing问题的表征。

上图展示的5阶傅里叶基函数模拟的结果。

随着阶数的增加，模拟的精度也随之上升，前面的两个问题还依然存在，而这种不连续位置的幅度增幅比例在傅里叶投影下是一个常量，数值接近18%。

在信号处理中，解决Ringing问题最常用的方式就是在频域中添加一个窗口函数（windowing，也叫sigma-factors），这个窗口函数的作用是过滤掉原始信号中的高频信息，从而可以使用较低的阶数完成对窗口函数截断后的信号的重建，而最常用的窗口函数就是sinc（这是因为频域中的sinc函数对应到空间域中就是一个box函数，而box函数正好可以用于对重建后的信号进行截断处理，从而避免震荡的发生），如上图中公式所示，其图形展示如上图右下角所示。

通常情况下，我们会选取sinc函数最主要的一个lobe，也就是从原点开始的第一个lobe，第一个lobe下降到0之后的震荡区域在频域中被抹去不用。

由于空间域中box函数的傅里叶变换得到的频域函数就是sinc函数，如上图所示，红色表示的是空间域的box函数，绿色曲线对应的则是频域的sinc函数，且box越宽，sinc lobe就越窄，而我们这里的sinc lobe相对来说比较宽，因此其对应的box就比较窄。

而在频域对一个信号乘上sinc函数，就相当于在空间域对信号的空间域函数用box函数进行卷积（而这个过程就相当于对重建信号进行截断处理），而得到的信号再用于傅里叶展开，就不会出现前面的不连续的问题了，卷积结果如上图右侧所示。

频率越高，box越窄，卷积之后的结果也就越接近此前的阶跃信号。

上图是跟box函数卷积之前的投影结果。

如上图所示，作为对比，跟box函数卷积之后的投影结果就消除了此前的两个问题，得到的平滑曲线就是卷积后的结果，不再存在震荡问题，不过这里需要注意的是，跟box函数卷积之前的投影信号得到的累计平方误差要小一些。

不断增加阶数，可以看到信号精度越来越高，不过不管怎样，误差总还是存在的。

在这种情况继续添加一个窗口函数，也就是在频域再次乘以一个sinc函数，在空间域再次与一个box函数进行卷积，实际上两个box函数的卷积是一个tent函数，这也是生成b样条曲线的一种方式，而与tent函数的卷积则可以消除前面的误差。

除了sync函数之外，还有其他的窗口函数，比如这里提到的Hanning函数（绿色），这个函数在降低高频噪音上面比sync函数（红色）见效更快，但是又比sync函数的平方（黄色）要慢。

sinc函数也叫sigma factor，而使用截断的sinc函数（只保留第一个lobe）常被叫做Lancosz sigma factors。

在实际使用过程中，选用哪个窗口函数并不是最重要的，更为重要的是如何做到在不至于对原始信号造成过于强烈的模糊的前提下消除Ringing问题，也就是做到Ringing问题与模糊问题之间的平衡。

上图跟下图给出的是脉冲函数与窗口函数在频域相乘后被投影到6阶SH上的表现，上图使用的是Lancosz函数，下图使用的则是Hanning函数，分别展示了不同频率的窗口函数的效果比对。

可以看到，在不同的频率下，窗口函数对Ringing问题的表现是不同的，而调整窗口函数的频率也是平衡Ringing问题与模糊问题常用的方法。

这是Hanning函数的投影结果，在同样的宽度下，Hanning函数的模糊程度更严重一点（高频信号消除速度当然也更快一点）

这里将多种窗口函数的效果放在一起对比，红色的是没有添加窗口函数的情况，蓝绿色则是对应前面不同频率窗口函数的表现。

前面展示的都是一些数字信号下的表现，上图给出的是采用脉冲光源的渲染效果在不同窗口函数下的表现，左边是没有添加窗口函数的情况，可以看到在明暗交界处存在光线的异常明亮，而添加了不同频率的Hanning窗口函数之后，这个问题得到了有效的扼制。

中间使用的是频率为6（自变量上的振幅，不是真的频率）的Hanning函数，其对高频信号的衰减强度较高，虽然在移除Ringing上的表现较好，但是同样也造成了阴影信号的不必要模糊。

这里给出的是Ringing问题在光照颜色上的异常表现，在这个示例中，光源为环境光+一盏黄色的方向光，可以看到在没有添加窗口函数的情况下，背光面（方向光）存在着黄色的漏光现象，且在最暗的地方将红绿颜色抽出，只留下了蓝色，使得阴影发蓝，而添加了窗口函数之后，这个问题就得到了有效的缓解。

另外一种添加窗口函数的方法则是改变大家对于投影方式的认知，我们知道，一个函数在SH上的投影，实际上是为了最小化上图中最上方的积分误差。

但是实际上，我们可以换一种误差计算方式，比如上图中间给出了另一种误差计算公式，是在原始的误差函数基础上添加了一个补偿项（penalty），而这个补偿项是拉普拉斯函数平方的积分乘上一个 $\lambda$ （这是球面上常用的曲度测量（curvature measure）算子），这个补偿项的作用是抑制拟合过程中的震荡（为啥这一项可以抑制震荡？），而拉普拉斯函数平方的积分通过SH投影可以很方便的计算出来。

在给定 $\lambda$ 的情况下，我们可以求出这个公式的闭合解（解析符号解），其中 $f_l^m$ 是直接投影到不做任何处理下的SH基函数上的系数，而 $c_l^m$ 则是对误差函数进行修正后的投影系数。

可以看到，最终的系数将取决于 $\lambda$ ，当 $\lambda$ 为0的时候，就得到了原始的投影系数。当 $\lambda$ 接近无穷大时，此时只有DC项（ $f_0^0$ ）是非0的，此时的曲度（curvature）是0。

而要想找到合适的 $\lambda$ ，可以考虑对平方拉普拉斯函数项按照一个固定的比例进行衰减（比如每次减半，即二分查找）来确定哪个 $\lambda$ 表现最好（有点类似于牛顿寻根法）。

上图给出了不同 $\lambda$ 下的投影表现，最右侧将几个投影表现放在一起进行对比。

前面给定窗口函数是作用在全信号上的，但是在图形学中，我们的输入信号本身是具有有限作用域的，也就是说，我们需要计算SH投影的图形学函数通常是一些本身就比较平滑，比较紧凑（Compact）的。

而这些函数如果有一个主要的波峰，那么我们只需要在震荡会对渲染结果造成影响的情况下添加窗口函数即可。

这里给出的一个光滑的物体被脉冲光源照明的案例，可以看到那一圈黑色阴影中心出现了白色的Ringing问题，而通过全局窗口函数，我们可以消除这个问题，但是在阴影边缘上则不可避免的出现了阴影的模糊，而最右侧给出的Content Sensitive的窗口函数则不存在上述两个问题。

在介绍SH Product的意义以及具体做法之前，先来介绍一下应用场景：

这个东西可以用于给环境光（天光cubemap）开一个口子，比如天空飞过的带有天光效果的飞行物，或者近处的一栋大楼挡住了天光，或者消除地表下的天光影响等
有几篇论文也在动态场景中使用这个东西来进行visibility累加
还可以用来模拟天光环境下云层的移动

如上图最上方的公式所示，所谓的SH Product实际上就是将两个SH函数（图中的f跟g，从上标推测应该是使用SH基函数表达的函数f跟函数g）投影到SH基函数上，这个过程涉及到了三重乘法。

对三个SH基函数相乘的积分得到的是学术上称之为clebsch gordan coupling coefficient 的系数，这个系数可以用三重积(triple product)tensor与两个函数对应的系数的乘积来表示，而如果两个函数之中有一个是固定值（比如f是固定的，而g可能会是其他经常变化的函数），那么我们可以为f构建一个矩阵M，而构建了矩阵之后，计算复杂度就从原来的2577个乘法减少到1296个（依然很复杂）。