在学习图形渲染的过程中，一直對球谐函数（球谐光照）有一点了解，但没有亲手实现过。

关于球谐函数的数学概念是比较复杂的，笔者也实在无法完全理解，这篇文章只能从做东西的角度来说，公式是如何和代码进行对应的。

笔者的代码是基于DirectX实现的。期间也遇到了一些坑，也顺便记录一下。

在数学中，一个基函数是一个函数空间（Function Space）中的一个基底，就像欧拉空间中的一个坐标轴一样。

在函数空间中，每个连续的函数都可以表示为基函数的线性组合。

例如，对于所有的二次多项式，可以用基函数组{ 1 , t , t 2 1,t,t^2 1,t,t2}进行线性组合得到，即 a + b t + c t 2 a+bt+ct^2 a+bt+ct2。这里{ a , b , c a,b,c a,b,c}即表示各个基函数的系数。

球谐函数是将谐函数限制于球坐标系下的单位球面的一组基函数 y i ( t ) y_i(t) yi​(t)，其各项常数系数为 c i c_i ci​，可以用来近似目标函数 f ( t ) f(t) f(t)。

球谐函数的表达式定义如下：

f ( t ) ≈ f ~ = ∑ l = 0 ∑ m = − l l = + l c l m y l m = ∑ i = 1 N c i y i f(t)\approx \tilde{f} = \sum_{l=0}\sum_{m=-l}^{l=+l}c_l^my_l^m = \sum_{i=1}^{N}c_iy_i f(t)≈f~​=l=0∑​m=−l∑l=+l​clm​ylm​=i=1∑N​ci​yi​

其中，N表示了系数或者说基底函数的数量，一个连续函数，需要无限个基函数的线性组合，才能完全恢复。

而 y l m y_l^m ylm​是一个阶数 ( l , m ) (l,m) (l,m)和角度（法线 n n n）相关的量，称作球谐基。 c l m c_l^m clm​为对应球谐基方向上的系数。

球面空间上的球谐函数 y l m y_l^m ylm​可视化如下：

绿色表示球谐函数的值为正值，而红色表示球谐函数的值为负值；矢径越大球谐函数值的绝对值越大，反之矢径越小球谐函数值的绝对值越小。

而在实际的使用中，一般只取前几阶的球谐函数来近似球面上的函数。

例如，前三阶的球谐基函数如下：

l=0：

Y 0 0 = s = 1 2 1 π Y_0^0 = s = \frac{1}{2} \sqrt{\frac{1}{\pi}} Y00​=s=21​π1​ ​

l=1：

Y 1 − 1 = p y = 3 4 π y r Y_1^-1 = p_y = \sqrt{\frac{3}{4\pi}} \frac{y}{r} Y1−​1=py​=4π3​ ​ry​

Y 1 0 = p z = 3 4 π z r Y_1^0 = p_z = \sqrt{\frac{3}{4\pi}} \frac{z}{r} Y10​=pz​=4π3​ ​rz​

Y 1 1 = p x = 3 4 π x r Y_1^1 = p_x = \sqrt{\frac{3}{4\pi}} \frac{x}{r} Y11​=px​=4π3​ ​rx​

l=2：

Y 2 − 2 = d x y = 1 2 15 π x y r 2 Y_2^-2 = d_{xy} = \frac{1}{2}\sqrt{\frac{15}{\pi}} \frac{xy}{r^2} Y2−​2=dxy​=21​π15​ ​r2xy​

Y 2 − 1 = d y z = 1 2 15 π y z r 2 Y_2^-1 = d_{yz} = \frac{1}{2}\sqrt{\frac{15}{\pi}} \frac{yz}{r^2} Y2−​1=dyz​=21​π15​ ​r2yz​

Y 2 0 = d z 2 = 1 4 5 π − x 2 − y 2 + 2 z 2 r 2 Y_2^0 = d_{z^2} = \frac{1}{4}\sqrt{\frac{5}{\pi}} \frac{-x^2-y^2+2z^2}{r^2} Y20​=dz2​=41​π5​ ​r2−x2−y2+2z2​

Y 2 1 = d x z = 1 2 15 π z x r 2 Y_2^1 = d_{xz} = \frac{1}{2}\sqrt{\frac{15}{\pi}} \frac{zx}{r^2} Y21​=dxz​=21​π15​ ​r2zx​

Y 2 2 = d x 2 − y 2 = 1 4 15 π x 2 − y 2 r 2 Y_2^2 = d_{x^2-y^2} = \frac{1}{4}\sqrt{\frac{15}{\pi}} \frac{x^2-y^2}{r^2} Y22​=dx2−y2​=41​π15​ ​r2x2−y2​

其中， r = 1 r=1 r=1，更多公式请参考：wiki-Table_of_spherical_harmonics。

2.1中提到了球谐函数表达式，但还没有介绍其中的系数 c l m c_l^m clm​是如何求解的。

其求解方法如下公式所示：

C l m = ∫ S f ( s ) Y l m ( s ) d s C_l^m = \int_S{f(s)Y_l^m(s)}ds Clm​=∫S​f(s)Ylm​(s)ds

其中， S S S表示积分区间单位球面。这个公式其实就是用原函数 f ( s ) f(s) f(s)和对应的球谐函数 Y l m ( s ) Y_l^m(s) Ylm​(s)在整个球面空间积分即可。这个求解系数（生成球谐系数）的过程，称之为投影。

在实际过程中，投影所得的系数 C l m C_l^m Clm​通过采用蒙特卡洛积分来近似求解。

C l m = ∫ S f ( s ) Y l m ( s ) d s = 1 N ∑ j = 1 N f ( x j ) w ( x j ) C_l^m = \int_S{f(s)Y_l^m(s)}ds = \frac{1}{N} \sum_{j=1}^{N}f(x_j)w(x_j) Clm​=∫S​f(s)Ylm​(s)ds=N1​j=1∑N​f(xj​)w(xj​)

其中，权重 w ( x j ) = 1 p ( w j ) w(x_j)=\frac{1}{p(w_j)} w(xj​)=p(wj​)1​，对于球面均匀采样的而言，有：

w ( x j ) = 1 p ( w j ) = 1 S 球 面 = 1 4 π w(x_j)=\frac{1}{p(w_j)}=\frac{1}{S_{球面}}=\frac{1}{4\pi} w(xj​)=p(wj​)1​=S球面​1​=4π1​

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

C i = ∫ S l i g h t ( v ) Y i ( v ) d s ≈ 1 N ∑ j = 1 N l i g h t ( x j ) Y ( x j ) ⋅ 4 π = 4 π N ∑ j = 1 N l i g h t ( x j ) Y ( x j ) \begin{aligned} C_i = & \int_{S}{light(v)Y_i(v)}ds \\ \approx & \frac{1}{N} \sum_{j=1}^{N}{light(x_j)Y(x_j)\cdot 4\pi}\\ = & \frac{4\pi}{N} \sum_{j=1}^{N}{light(x_j)Y(x_j)} \end{aligned} Ci​=≈=​∫S​light(v)Yi​(v)dsN1​j=1∑N​light(xj​)Y(xj​)⋅4πN4π​j=1∑N​light(xj​)Y(xj​)​

而重建则是通过一系列的系数 C l m C_l^m Clm​近似恢复函数 f f f的过程。

正如前面所描述的，因为我们不会存储无穷阶系数，对 f f f采用 l l l阶的球谐函数进行恢复。

f ( t ) ≈ f ~ = ∑ l = 0 ∑ m = − l l = + l c l m y l m = ∑ i = 1 N c i y i f(t)\approx \tilde{f} = \sum_{l=0}\sum_{m=-l}^{l=+l}c_l^my_l^m = \sum_{i=1}^{N}c_iy_i f(t)≈f~​=l=0∑​m=−l∑l=+l​clm​ylm​=i=1∑N​ci​yi​

在渲染中，球谐函数可以用来重建辐射亮度（Radiance），一般为环境贴图Cubemap。

以下是1阶到6阶SH重建辐射亮度的效果图：

可以看出：

笔者输入了一张HDR的Cubemap，重建辐射亮度如下：

如果采用一张较为低频的图像，则是以下的结果：

当然，重建辐射亮度（Radiance）并非球谐函数最主要的作用。

其最重要的作用是第三节中要介绍的恢復IrradianceMap！

这里所说的环境光指的是天空盒/天空球（无限远）所发出的光，如此一来我们可以忽略掉位置信息而只考虑法线信息。

漫反射光照计算公式如下：

L ( p , w o ) = ∫ Ω L ( p , w i ) n ⋅ w i d w i L(p,w_o)=\int_{\Omega} L(p,w_i)n\cdot w_idw_i L(p,wo​)=∫Ω​L(p,wi​)n⋅wi​dwi​

通常采用预积分的方式来生成一张环境贴图的IrradianceMap（下图右侧）。

在实时渲染时，仅需使用一个法线向量对立方体贴图进行采样，就可以获取该方向上的场景辐照度。

由于辐照度图（IrradianceMap）变化较为缓慢，看起来有点像环境的平均颜色或光照图，可以以较低的分辨率进行存储，例如64x64x6。

这里要对球谐函数近似IrradianceMap的数学理论进行简单的介绍。

这里的数学推导摘錄自以下兩篇博客，寫的實在太棒了！建议观看原文！

筆者在這僅僅摘录了其中的结果！

【论文复现】Spherical Harmonic Lighting:The Gritty Details

【论文复现】An Efficient Representation for Irradiance Environment Maps

首先，对漫反射光照计算公式按照球谐函数进行展开：

KaTeX parse error: No such environment: align* at position 7: \begin{̲a̲l̲i̲g̲n̲*̲}̲ f(x) = & \sum…

其中，系数为 c l m = ∫ Ω f ( w ) Y l m ( w ) d w c_l^m = \int_{\Omega } f(w)Y_l^m(w)dw clm​=∫Ω​f(w)Ylm​(w)dw

接下来，对光照方程进行简化。

将积分里的函数分为两个部分即：

{   l i g h t ( w ) = L ( p , w )   t ( w ) = n ⋅ w \begin{cases} &\text{ } light(w) = L(p,w) \\ &\text{ } t(w) = n \cdot w \end{cases} {​ light(w)=L(p,w) t(w)=n⋅w​

然后我们分别对着两个函数进行球谐函数展开即：

{   l i g h t ( w ) = ∑ i = 0 L i Y i ( w )   t ( w ) = ∑ j = 0 t j Y j ( w ) \begin{cases} &\text{ } light(w) = \sum _{i=0} L_i Y_i(w) \\ &\text{ } t(w) = \sum _{j=0} t_j Y_j(w) \end{cases} {​ light(w)=∑i=0​Li​Yi​(w) t(w)=∑j=0​tj​Yj​(w)​

其中， L i L_i Li​和 t j t_j tj​都是常数。

将上述带回到光照积分公式中：

L ( p , w o ) = ∫ Ω l i g h t ( w ) t ( w ) d w = ∫ Ω ( ∑ i = 0 L i Y i ( w ) ) ⋅ ( ∑ j = 0 t j Y j ( w ) ) d w = ∫ Ω ∑ i = 0 ∑ j = 0 ( L i Y i ( w ) t j Y j ( w ) ) d w = ∑ i = 0 ∑ j = 0 L i t j ∫ Ω Y i ( w ) Y j ( w ) d w \begin{aligned} L(p,w_o) = & \int_{\Omega} light(w) t(w)dw \\ = & \int_{\Omega} (\sum _{i=0} L_i Y_i(w)) \cdot (\sum _{j=0} t_j Y_j(w)) dw \\ = & \int_{\Omega} \sum _{i=0}\sum _{j=0}(L_i Y_i(w) t_j Y_j(w)) dw \\ = & \sum _{i=0}\sum _{j=0} L_i t_j \int_{\Omega} Y_i(w)Y_j(w) dw \end{aligned} L(p,wo​)====​∫Ω​light(w)t(w)dw∫Ω​(i=0∑​Li​Yi​(w))⋅(j=0∑​tj​Yj​(w))dw∫Ω​i=0∑​j=0∑​(Li​Yi​(w)tj​Yj​(w))dwi=0∑​j=0∑​Li​tj​∫Ω​Yi​(w)Yj​(w)dw​

由于球谐系数的正交完备性，有且仅有 i = = j i == j i==j的时候， ∫ Ω Y i ( w ) Y j ( w ) d w = 1 \int_{\Omega} Y_i(w)Y_j(w) dw=1 ∫Ω​Yi​(w)Yj​(w)dw=1。

则光照积分公式可以简化为：

L ( p , w o ) = ∑ i = 0 L i t i L(p,w_o) = \sum_{i=0} L_i t_i L(p,wo​)=i=0∑​Li​ti​

此时整个光照函数就简化为了常数积分之和。

我们接着分析球谐系数 L i L_i Li​和 t i t_i ti​。

分析 L i L_i Li​

L i = ∫ Ω L ( p , w ) Y i ( w ) d w L_i= \int_{\Omega } L(p,w)Y_i(w)dw Li​=∫Ω​L(p,w)Yi​(w)dw

函数里面虽然有着色亲P，但实际与具体的着色点无关。因此这一项我们可以直接对整个天空盒进行积分。

伪代码如下：

这里注意，我们直接将天空盒的位置信息进行归一化后就作为自变量。

因为我们这个积分是球面积分，需要球面上的点，天空盒的位置进行归一化就投影到球面上去了（就好像在球面上去点一样）。

分析 t i t_i ti​

t i = ∫ Ω ( n ⋅ w ) Y i ( w ) d w t_i= \int_{\Omega } (n \cdot w)Y_i(w)dw ti​=∫Ω​(n⋅w)Yi​(w)dw

可惜，求 t i t_i ti​时域具体的着色器是有关的。因为我们需要知道法线的信息 n n n。

这意味着如果要预计算 t i t_i ti​，也就需要对每一个方向的法线 n n n都要算一组 t i t_i ti​，若采用三阶球谐（需存储9个系数），那就需要三张CubeMap（每一张存3个系数）。

伪代码如下：

直接采用球谐系数替代原有光照的问题，即：针对每一个方向的法线都需要预计算出一组球谐系数。

球谐函数旋转

所谓的旋转，不是说直接去改变原函数，而是将传入自变量在三维空间中被旋转之后传入原函数。

// … 这里有一系列公式的推导，笔者也看不完全懂就不摘录了。

经过一系列推导，光照函数表示为：

L ( n ) = ∑ l = 0 ∞ ∑ m = − l l 4 π 2 l + 1 L l m t l Y l m ( n ) L(n)=\sum_{l=0}^{\infty } \sum_{m=-l}^{l} \sqrt{\frac{4\pi }{2l+1}} L_l^m t_l Y_l^m(n) L(n)=l=0∑∞​m=−l∑l​2l+14π​ ​Llm​tl​Ylm​(n)

可以分为3项：

1） 4 π 2 l + 1 t l \sqrt{\frac{4\pi }{2l+1}} t_l 2l+14π​ ​tl​ 为一系列常数。下面列出了前三阶的系数。

2） L l m L_l^m Llm​ 为预计算的， L i = ∫ Ω L ( p , w ) Y i ( w ) d w L_i= \int_{\Omega } L(p,w)Y_i(w)dw Li​=∫Ω​L(p,w)Yi​(w)dw。

3） Y l m ( n ) Y_l^m(n) Ylm​(n) 为球谐基函数，参数是具体着色点的法线 n n n。可以根据着色器中的法线进行实时计算。

在预计算 L l m L_l^m Llm​ 时，其实可以将第一项也包含进来，一起进行存储 c i c_i ci​。如下公式所示：

c i = 4 π 2 l + 1 t l ⋅ ∫ Ω L ( p , w ) Y i ( w ) d w c_i = \sqrt{\frac{4\pi }{2l+1}} t_l \cdot \int_{\Omega } L(p,w)Y_i(w)dw ci​=2l+14π​ ​tl​⋅∫Ω​L(p,w)Yi​(w)dw

对这个式子 L i = ∫ Ω L ( p , w ) Y i ( w ) d w L_i= \int_{\Omega } L(p,w)Y_i(w)dw Li​=∫Ω​L(p,w)Yi​(w)dw使用蒙特卡洛积分 I = ∫ f ( x ) d x = ∫ f ( x ) p ( x ) p ( x ) d x ≈ 1 N ∑ i = 1 N f ( X ) p ( 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)}{p(x)} I=∫f(x)dx=∫p(x)f(x)​p(x)dx≈N1​∑i=1N​p(x)f(X)​近似，有：

对于均匀的球面上采样，概率 p = 1 4 π p=\frac{1}{4\pi} p=4π1​，则有：

c i = 4 π 2 l + 1 t l ⋅ 4 π N ∑ L j Y j c_i = \sqrt{\frac{4\pi }{2l+1}} t_l \cdot \frac{4\pi}{N}\sum L_jY_j ci​=2l+14π​ ​tl​⋅N4π​∑Lj​Yj​

则渲染还原时为：

L ( n ) = ∑ i = 0 c i Y i ( n ) L(n)=\sum_{i=0} c_i Y_i(n) L(n)=i=0∑​ci​Yi​(n)

此时求出的为SH近似出来的IrradianceMap。

相较于2.2对Radiance的重建，3.2重建IrradianceMap仅多出来一项，其他完全相同。即 4 π 2 l + 1 t l \sqrt{\frac{4\pi }{2l+1}} t_l 2l+14π​ ​tl​ 常数项。

而已仅需在预计算时乘以对应的常数项系数，或者在渲染时乘以对应的常数项系数，都可以获得IrradianceMap。

但是为什么要使用SH替代CubeMap呢？

原因很简单：

下面给出笔者实现的SH重建的Irradiancemap。

原图：

预积分的IrradianceMap：

SH恢复的IrradianceMap：

笔者在实践中使用的DirectX，期间遇到一个读取数据的问题，也记录一下。

笔者通过DirectXTex库的CaptureTexture函数获取了CubeMap在CPU中的数据存储于DirectX::ScratchImage类型中。

但由于格式为DXGI_FORMAT_R16G16B16A16_FLOAT，以下是笔者尝试的代码，并不能正确地读取到图像数据。

最终，笔者通过将DXGI_FORMAT_R16G16B16A16_FLOAT转为DXGI_FORMAT_R32G32B32A32_FLOAT格式，才可以顺利读取到正确的像素值。代码如下：

同时需要注意的是对于Cubemap，其Mipmap的存储方式为：

SphericalHarmonicsLighting

【论文复现】Spherical Harmonic Lighting:The Gritty Details

【论文复现】An Efficient Representation for Irradiance Environment Maps

球谐光照技术 —— 原理+程序（全网最通俗易懂的版本）

Unity渲染编程(灯光篇)【第一卷：Spherical Harmonics Lighting】

PI or not to PI in game lighting equation

球谐光照-实验篇

球谐光照-应用篇

球谐(SH)和预计算辐射度(PRT)

球谐光照