球谐光照（Spherical Harmonics Lighting）

文章目录

• 球谐光照（Spherical Harmonics Lighting）
• • 一、前言
  • 二、球谐函数
  • • 2.1 基函数
    • 2.2 投影与重建
    • 2.3 应用
  • 三、漫反射环境光
  • • 3.1 IrradianceMap
    • 3.2 数学推导
    • 3.3 实践
  • 参考博文

一、前言

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

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

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

二、球谐函数

2.1 基函数

在数学中，一个基函数是一个函数空间（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.2 投影与重建

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)的数值估计表达式是：

KaTeX parse error: No such environment: align* at position 7: \begin{̲a̲l̲i̲g̲n̲*̲}̲ C_i = & \int_…

而重建则是通过一系列的系数 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​

2.3 应用

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

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

可以看出：

• 当球谐函数的阶数越高，还原的效果越好。
• 球谐函数的结果是比较粗糙的，只能模拟低频信号。

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

• 效果不好，根本原理是球谐函数只能模拟低频信号，而HDR图像的数值范围较大，无法模拟好。

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

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

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

三、漫反射环境光

3.1 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。

3.2 数学推导

这里要对球谐函数近似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​都是常数。

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

KaTeX parse error: No such environment: align* at position 7: \begin{̲a̲l̲i̲g̲n̲*̲}̲ L(p,w_o) = & \…

由于球谐系数的正交完备性，有且仅有 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，但实际与具体的着色点无关。因此这一项我们可以直接对整个天空盒进行积分。

伪代码如下：

for(pixel p : cubeMap)
    Li += p.color * Yi(normalise(p.position))*dw;

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

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

分析 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个系数）。

• 还要和IBL的CubeMap一样大小的纹理。

伪代码如下：

// t与法线相关，只能把所有的法线都计算了
for(normal &n: sphere)
{
    for(pixel &p : Cubemap)
        ti[n] += dot(n,normalise(p.position)) * Yi(normalise(p.position)) * dw;
}

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

球谐函数旋转

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

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

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

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​ 为一系列常数。下面列出了前三阶的系数。

// Convolve with SH-projected cosinus lobe
float ConvolveCosineLobeBandFactor[] =
{
    PI,
    2.0f * PI/3.0f, 2.0f * PI/3.0f, 2.0f * PI/3.0f,
    PI/4.0f, PI/4.0f, PI/4.0f, PI/4.0f, PI/4.0f
}

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。

3.3 实践

相较于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呢？

原因很简单：

• 对于一个大的场景而言，每个位置的环境光都可能不同，如果每个点的环境光都采用一张Cubemap来存储，并且每次在Shader进行采样，那么这种方法无疑是非常昂贵的。
• 而使用球谐函数的话，可以将Cubemap减少到只有9个SH系数（如果是三阶的话）。

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

原图：

预积分的IrradianceMap：

SH恢复的IrradianceMap：

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

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

DirectX::ScratchImage image;
FCommandQueue& Queue = g_CommandListManager.GetQueue(D3D12_COMMAND_LIST_TYPE_DIRECT);
hr = DirectX::CaptureTexture(Queue.GetD3D12CommandQueue(), m_Resource.Get(), true/*isCubeMap*/, image, m_CurrentState, m_CurrentState);

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

const Image img = dstSImg.GetImages()[i];
{
    size_t rowPitch = 0;
    size_t slicePitch = 0;
    ComputePitch(img.format, img.width, img.height, rowPitch, slicePitch);
    for (int rowIndex = 0; rowIndex < img.height; rowIndex++)
    {
        uint16* dst = (uint16*)(img.pixels + rowPitch * rowIndex);
        for (int colIndex = 0; colIndex < img.width; colIndex++)
        {
            // ... 这样并不能读取到正确的数据
            uint16 R = dst[4*colIndex+0];
        }
    }
}

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

/*
* InputImage 输入的R16G16B16A16_FLOAT格式Cubemap
* Width、Height 图像分辨率
* SampleNum 随机采样的总数量
* OutSamples 输出的所有样本点（法线方向、颜色）
*/
class Vertex
{
public:
	Vector3f pos, color;
};
void RandomSample(const DirectX::ScratchImage& InputImage, int Width, int Height,int SampleNum, std::vector<Vertex>& OutSamples)
{
    // 通过 DirectX::_ConvertFromR16G16B16A16 将其转成DXGI_FORMAT_R32G32B32A32_FLOAT格式
	DirectX::ScratchImage dstSImg;
	dstSImg.Initialize2D(DXGI_FORMAT_R32G32B32A32_FLOAT, InputImage.GetMetadata().width, InputImage.GetMetadata().height, InputImage.GetMetadata().arraySize, InputImage.GetMetadata().mipLevels);
	for (int i = 0; i < InputImage.GetImageCount(); i++)
	{
		DirectX::_ConvertFromR16G16B16A16(InputImage.GetImages()[i], dstSImg.GetImages()[i]);
    }

    OutSamples.clear();
	OutSamples.resize(SampleNum);

	for (int i = 0; i < SampleNum; ++i)
	{
		float x, y, z;
		do {
			x = NormalRandom();
			y = NormalRandom();
			z = NormalRandom();
		} while (x == 0 && y == 0 && z == 0);

		Vertex vex;

		Vector3f pos(x, y, z);
		vex.pos = pos.Normalize();

		CubeUV cubeuv = XYZ2CubeUV(pos);

		int colIndex = (int)(cubeuv.u * (Width - 1));
		int rowIndex = (int)((1.f - cubeuv.v) * (Height - 1));
		
		const DirectX::Image* images = dstSImg.GetImages();
		int Lod0FaceIndex = cubeuv.index * dstSImg.GetMetadata().mipLevels;
		const DirectX::Image image = images[Lod0FaceIndex];

		size_t rowPitch = 0;
		size_t slicePitch = 0;
		ComputePitch(image.format, image.width, image.height, rowPitch, slicePitch);

        // 注意强转为float类型
		float* dst = (float*)(image.pixels + rowPitch * rowIndex);
		float R = dst[colIndex * 4 + 0];
		float G = dst[colIndex * 4 + 1];
		float B = dst[colIndex * 4 + 2];

		//printf("[%f,%f,%f]\n", R, G, B);
		vex.color = { R,G,B };
		OutSamples[i] = vex;
	}
}

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

//Face0: Mip0, Mip1, Mip2, ...
//Face1: Mip0, Mip1, Mip2, ...
//...
//Face5: Mip0, Mip1, Mip2, ...

参考博文

• 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)

• 球谐光照