球谐光照(Spherical Harmonics Lighting)

球谐光照(Spherical Harmonics Lighting)

文章目录

一、前言

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

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

笔者的代码是基于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​可视化如下:

球谐光照(Spherical Harmonics Lighting)

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

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

例如,前三阶的球谐基函数如下:

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重建辐射亮度的效果图:

球谐光照(Spherical Harmonics Lighting)

可以看出:

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

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

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

球谐光照(Spherical Harmonics Lighting)

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

球谐光照(Spherical Harmonics Lighting)
球谐光照(Spherical Harmonics Lighting)
当然,重建辐射亮度(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(下图右侧)。

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

球谐光照(Spherical Harmonics Lighting)

由于辐照度图(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。

原图:

球谐光照(Spherical Harmonics Lighting)

预积分的IrradianceMap:

球谐光照(Spherical Harmonics Lighting)

SH恢复的IrradianceMap:

球谐光照(Spherical Harmonics Lighting)

笔者在实践中使用的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, ...

参考博文

上一篇:ARM-Linux(内核)修改开机启动logo


下一篇:LED And Incandescent, Who Is Suitable For Holiday Lighting?