图形学基础|球谐光照（Spherical Harmonics Lighting）

球谐光照（Spherical Harmonics Lighting）

一、前言

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

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

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

二、球谐函数

2.1 基函数

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

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

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

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

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

$$f(t)\approx \tilde{f}=\sum _{l=0}\sum _{m=-l}^{l=+l}{c}_l^m{y}_l^m=\sum _{i=1}^N{c}_i{y}_i$$

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

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

球面空间上的球谐函数$y_l^m$可视化如下：

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

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

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

l=0：

$$Y_0^0=s=\frac{1}{2}\sqrt{\frac{1}{\pi }}$$

l=1：

$$Y_1^{-}1=p_y=\sqrt{\frac{3}{4\pi }}\frac{y}{r}$$

$$Y_1^0=p_z=\sqrt{\frac{3}{4\pi }}\frac{z}{r}$$

$$Y_1^1=p_x=\sqrt{\frac{3}{4\pi }}\frac{x}{r}$$

l=2：

$$Y_2^{-}2=d_{xy}=\frac{1}{2}\sqrt{\frac{15}{\pi }}\frac{xy}{r^2}$$

$$Y_2^{-}1=d_{yz}=\frac{1}{2}\sqrt{\frac{15}{\pi }}\frac{yz}{r^2}$$

$$Y_2^0=d_{z^2}=\frac{1}{4}\sqrt{\frac{5}{\pi }}\frac{-x^2-y^2+2z^2}{r^2}$$

$$Y_2^1=d_{xz}=\frac{1}{2}\sqrt{\frac{15}{\pi }}\frac{zx}{r^2}$$

$$Y_2^2=d_{x^2-y^2}=\frac{1}{4}\sqrt{\frac{15}{\pi }}\frac{x^2-y^2}{r^2}$$

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

2.2 投影与重建

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

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

$$C_l^m={\int }_{S}f(s)Y_l^m(s)ds$$

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

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

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

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

$$w(x_j)=\frac{1}{p(w_j)}=\frac{1}{S_{球面}}=\frac{1}{4\pi }$$

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

$$\begin{array}{rl}{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{array}$$

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

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

$$f(t)\approx \tilde{f}=\sum _{l=0}\sum _{m=-l}^{l=+l}{c}_l^m{y}_l^m=\sum _{i=1}^N{c}_i{y}_i$$

2.3 应用

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

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

可以看出：

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

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

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

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

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

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

三、漫反射环境光

3.1 IrradianceMap

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

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

$$L(p,w_o)={\int }_{\Omega }L(p,w_i)n\cdot w_{i}d{w}_i$$

通常采用预积分的方式来生成一张环境贴图的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={\int }_{\Omega }f(w)Y_l^m(w)dw$$

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

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

$$\{\begin{array}{ll}& light(w)=L(p,w)\\ & t(w)=n\cdot w\end{array}$$

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

$$\{\begin{array}{ll}& light(w)={\sum }_{i=0}{L}_i{Y}_i(w)\\ & t(w)={\sum }_{j=0}{t}_j{Y}_j(w)\end{array}$$

其中，$L_i$和$t_j$都是常数。

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

$$\begin{array}{rl}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{array}$$

由于球谐系数的正交完备性，有且仅有$i==j$的时候，${\int }_{\Omega }{Y}_i(w)Y_j(w)dw=1$。

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

$$L(p,w_o)=\sum _{i=0}{L}_i{t}_i$$

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

我们接着分析球谐系数$L_i$和$t_i$。

分析$L_i$

$$L_i={\int }_{\Omega }L(p,w)Y_i(w)dw$$

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

伪代码如下：

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

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

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

分析$t_i$

$$t_i={\int }_{\Omega }(n\cdot w)Y_i(w)dw$$

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

这意味着如果要预计算$t_i$，也就需要对每一个方向的法线$n$都要算一组$t_i$，若采用三阶球谐（需存储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)=\sum _{l=0}^{\infty }\sum _{m=-l}^l\sqrt{\frac{4\pi }{2l+1}}{L}_l^m{t}_l{Y}_l^m(n)$$

可以分为3项：

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

// 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_i={\int }_{\Omega }L(p,w)Y_i(w)dw$$。

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

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

$$c_i=\sqrt{\frac{4\pi }{2l+1}}{t}_l\cdot {\int }_{\Omega }L(p,w)Y_i(w)dw$$

对这个式子$$L_i={\int }_{\Omega }L(p,w)Y_i(w)dw$$使用蒙特卡洛积分$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)}$近似，有：

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

$$c_i=\sqrt{\frac{4\pi }{2l+1}}{t}_l\cdot \frac{4\pi }{N}\sum L_j{Y}_j$$

则渲染还原时为：

$$L(n)=\sum _{i=0}{c}_i{Y}_i(n)$$

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

3.3 实践

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

而已仅需在预计算时乘以对应的常数项系数，或者在渲染时乘以对应的常数项系数，都可以获得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)

• 球谐光照