Global Illumination_Spherical Harmonic Lighting(球谐光照)_球谐系数

首先我们需要知道的是，如何计算环境光shading，一般我们会想到IBL，其实我们也可以使用球谐函数来进行表示，本部分我们就先来了解下如何使用SH来计算环境光照（后续我们也会继续来看一下环境光阴影的计算方式，即PRT方案）。

在计算环境光照之前，我们首先需要了解基函数的概念,如下图简示（闫大神课件内容）:

使用上述基函数的原因主要是因为我们需要定义环境光照，其实环境光照就是一个定义在三维空间上的一个二维函数（极坐标表示）。
其中球谐系数项如下：

下边所有过程都是以上两图为主线思路进行计算求解。

一、球谐光照与IBL区别

IBL相关内容可参照以前文章Vulkan_IBL相关内容。

球谐光照最主要的应用就是用于计算环境光，这里的环境光只天空盒/天空球所发出的光，通常计算环境光的手段有IBL。天空盒/天空球发出的环境光有什么特点？其中最重要的特点天空球无限大，如此一来我们可以忽略掉位置信息而只考虑法线信息。也就是说，不管我们的模型多么大，我们计算模型上某一个具体点的时候我们认为该点处于天空球中心。

不过IBL有些缺点，就是我们需要对整个球面空间进行积分就变成了一次纹理采样即可。不过这个方法最大的问题也就是纹理采样，采用IBL的方法，我们需要存储一张立方体贴图，纹理采样的带宽对于电脑来说可能不算什么，不过对于手机而言就实在是有点奢侈了。

二、球面坐标系

先简单的介绍一下球面坐标系，所谓球面坐标系其实就是由方位角、仰角和距离构成的坐标系（就是把原来的三维空间中的点通过替换成与z轴的夹角和与x轴的夹角以及与坐标原点的距离r 构成）如下：

只需要通过几何变换就可以从直角坐标系与球面坐标系的转化：

而本文讨论的都是位于单位圆上（半径r =1）的坐标，因此后面的讨论中都不考虑r .

二、球谐函数

球谐光照(Spherical Harmonic Lighting)就是基于球面调和(SH, Spherical Harmonics)这个数学工具的一种光照/着色算法。球谐光照也是一种用让我们捕捉(capture)光照、随后进行重新光照(relight)、实时展现全局光照(Global Illumination)风格的区域光源(area light sources)与软阴影的一种技术。

如果你可以理解一维函数的傅立叶变换(Fourier Transform)，那么其实你也很容易理解在图形学领域里面SH的意义，因为在限定在图形学的语境下，SH就是一种定义在球面上的广义二维傅立叶变换，这可以把离散的球面函数参数化。

一般来说，球谐光照可以用有限制带宽的Spherical Harmonics来模拟低频的环境光照明，反射光和高光比较高频，用低阶球谐函数来编码的话精度是不够用的。

如下图可以看出各阶函数的定义：

其中颜色表示变化频率的正负，亮度表示绝对值的大小。
我们也可以从下边动图中直观的来看一下：

具体的渲染方程便不再赘述，可查看以前专栏相关内容。

2.1 球谐函数性质

从球谐的定义可以直观从上图看出来：一个二维函数，阶数越高，频率越高。此外，大家不要被球谐函数吓到了，不要管网上公式有多么复杂，那都是虚的，大家就想象成一个三角函数可。更简单的理解：给它一个向量它就能返回一个均匀且低频的值（一般使用前三/四阶）。
此处我们直接来看一下前四阶基函数的表达：

2.1.1 正交完备性

正交完备性如下，不用管具体Y是什么以及怎么算，只需要知道必须两个球谐函数一模一样才为1否则就为0即可。这个性质是球谐函数用来简化计算的核心。

2.1.2 旋转不变性

这里的旋转不变目的就是环境光照变化之后我们只需要简单的计算就可以得到光源旋转之后的结果。这个性质使得球谐函数适用性更加的强。

2.2 球谐函数定义

球谐函数的性质交代完了我们接下来继续优化光照公式，对于复杂的函数，我们可以通过各种方法比如泰勒展开，广义傅里叶展开。把一个复杂的函数展开成简单的函数如：

fi 为常数， Pi(x) 为其他的函数，大家不用管 P 函数指什么，通过上面的式子可以看出来如果我们采用将函数展开的方式如果展开次数越多也就是 N 越大那么对 f(x) 的还原就越好。

展开的方式多种多样，本文介绍的就是采用球谐函数展开：

其中 cl 为常数也叫做球谐系数，生成球谐系数的过程，也称为投影，其计算方法如下：

这个式子看起来很恐怖，其实就是一个用原函数 f(x) 乘以对应项的球谐函数 Yl 在整个球面空间积分即可，当然这个是提前预计算好的。

现在我们有了球谐函数（此时不必关心 Yl 内部长什么样），也知道了球谐函数的两条性质，接下来就开始简化光照方程。对于光照方程：

我们把积分里面的函数分成两部分即：

这里为了表述方便把 wi 写作 w 。然后我们分别对着两个函数进行球谐函数展开即：

一般来说我们求漫反射环境光只需要L=3即可，也就是N=9（即使用前三阶九个基函数）。接下来我们代回光照积分公式即：

提取常数项 Li，ti：

这个时候只看积分里面就用得上我们的性质一：正交完备性。有且仅有 i == j 的时候才为1，其余都为0。

此时整个光照公式就化简为常数乘积之和。我们接着分析球谐系数 Li，ti：

2.3 球谐系数

2.3.1 Li项

原函数为 L(p,w) ，我们知道在这个漫反射函数里面虽然有着色点p，但是实际与具体着色点无关。因此 Li 这一项我们可以直接对整个天空盒进行积分即可，伪代码如下：

for(pixel &p : Cubemap)
    Li += p.color * Yi(normalise(p.position)) * dw;
1
2

这里注意，我们直接将天空盒的位置信息进行归一化后就作为自变量。因为我们这个积分是球面积分，需要球面上的点，天空盒的位置进行归一化就投影到球面上去了（就好像在球面上去点一样）。

2.3.2 ti项

求 ti 是与具体着色点有关（需要知道法线信息n）。这也就意味着，我们如果需要预计算 ti ，也就需要对每一个方向的法线n 都要算一组 ti 。按照我们i=9计算，也就是我们和IBL方法相比，需要3张（每一张存3个系数）一样大小的纹理才可以！伪代码如下：

for(normal &n: sphere)
{
    for(pixel &p : Cubemap)
        Li[n] += dot(n,normalise(p.position)) * Yi(normalise(p.position)) * dw;
}
1
2
3
4
5

此种求解方式为常规计算方式（Spherical Harmonic Lighting:The Gritty Details的实现原理）。

到此为止我们目前才用了球谐函数的第一个性质，不要忘记球谐函数还是第二个性质：旋转不变性，根据此性质不需要预计算 ti 的方式。

三、球谐函数旋转

首先需要说明的是，所谓旋转不是说直接去改变原函数，而是传入的自变量在三维空间中被旋转之后传入原函数。因此我们的目的是观察传入这个旋转之后的自变量，通过一系列化简之后原函数会怎么变成什么样。

对于球谐函数如上式所示（ u 为三维空间的点），仍然不用关心其中每一项是什么意思。现在大家需要思考一个问题，如果将 u 进行旋转，这个球谐函数会如何变化呢？假设对点 u 进行三个方向的旋转如下所示：

假设进行上述旋转操作为

则球谐函数经过化简之后为：

d为为维格纳D矩阵，具体可百度搜索详细表达式。

而此式子的最大作用就是： 如果想得到 [公式] 旋转之后的值，无需重新预计算 Yl，只需要在以前计算的基础上把已经计算好的值再乘以一个矩阵就好了，这一下就大大简化了计算。

对于函数 f(u) ,将其采用球谐函数展开如下：

同样，也将 u 进行旋转，代入式得：

经过层层变换，最终得到上式，这就是球谐函数的旋转不变性，我觉得旋转不变性描述的不够形象生动，应该叫做旋转可分离性。因为代入旋转后的点并不是最终的函数值不改变，而是我们可以通过层层化简最终将旋转操作拆分为一个函数，从而通过两个函数的乘积就得到了旋转之后的函数值。而根据乘法结合律，本来是应该先计算：

我们转化为先计算：

如此一来就可以不用重新计算球谐系数，而是对原来的球谐系数进行处理就得到旋转之后的球谐系数。

四、漫反射环境光计算

看一下我们简过的漫反射环境光公式：

上式主要在半球空间积分，入射光方向 wi， n 为着色点 p 的法线。这个公式就描述了着色点 p 在整个球面空间中收到的光照总和。

用局部坐标系重新描述漫反射环境光公式（对于不同的坐标系需要明确一点，只是参考系变了），如下所示：

其中，

而局部坐标系是以法线 n 为z轴的坐标系（局部坐标系下的变量都带有 ’ ），此时需要注意的入射光方向 wi 是世界坐标系的，将其转换到以法线 n 为z轴的局部坐标系（这里只有旋转没有平移操作）：

其中,

为将局部坐标系中的 wi’ 旋转到世界坐标系 wi 所需要旋转的角度，也就是法线 n 在世界坐标系的角度（这一点后面会用到，一定要注意角度是怎么来的）。然后分别对着两个函数进行球谐函数展开。

4.1 解析 L(wi’)

根据：

旋转 wi’ 之后得到：

不过还是不能直接用上式，根据上面的局部坐标系图再仔细分析旋转的三个角度 ，对于光照计算公式而言我们只关心光线与法线的夹角也就是 a ，那么对于其他两个角度完全可以不用旋转。接着继续化简上式，根据化简得出结论（注意，这里的推导不是无偏的，所以最后的结果也不是无偏的！）：

• m’!=0 的项都等于0，也就是说只用考虑 m’=0 的情况即可。
• 当 m’=0 时，在只旋转 a 角的情况下，得到式如下：

带回到上边的旋转式得到

最终结合上式就得到旋转 a 之后的 ***L(wi)***：

4.2 解析 t(Q’)

这里不去求世界坐标系下的 t(Q) 是多少而只关心在局部坐标系下的 t(Q’) :

同时，任然只关心与法线的夹角 Q ，其余角度不影响最终结果，那么将其设置为0，如此一来得到 m’！=0 的项都等于0，从而得到展开的公式为：

这样有什么好处呢？好处就在于求 t(Q’) 的球谐系数 ***tl***可以求出具体的值（有人可能会问，为什么求 L(wi’) 不这样做呢，因为旋转的角度是具体的法线与坐标轴的夹角，对于 [公式] 而言是无法提前知道法线究竟是多少的，虽然推导出式中包含法线的夹角，那只是为了后面进一步化简而提前做的准备）。注意 t 的下标是 l 的索引，对于相同的 l不管 m 如何变都相等。

4.3 求漫反射环境光L(n)

将上述计算的 L(wi’) 与 t(Q’) 带入漫反射环境光式中得到：

然后再把与 积分项 wi’ 无关项提出积分中：

再根据推导的正交完备性化简得到只有当 l=k 时积分里面的项才不为0反之为1，如此一来进一步化简得到：漫反射环境光的光照公式

4.4 小结（重中之重）

1. 首先对于上式中的 Ll 而言，仍然需要预计算出来；
2. 对于 tl 可以根据4.2后最后公式所求出的几个常数与 Yl 的乘积；
3. 对于 Yl(n) ，注意他的参数是具体着色点的法线 n。那也就是说完全可以在着色器中根据顶点着色点传过来的法线 n 计算 Yl(n) ,然后再与外部提前预计算好的 Ll 和 tl 进行简单乘积就得到最终的漫反射环境光。

五、高光项环境光计算

接着继续推导镜面反射环境光，在解决了漫反射环境光之后，镜面反射的推导就非常简单了。对于在天空盒影响下的光泽反射环境光的光照公式如下：

像IBL中化简方式一样，我们将上式分解成两个部分（这里会存在误差，不过对于环境光而言没关系）

注意这里与IBL 文章中的分解方式不一样，n.wi 被分到前面一部分了（IBL是分到后面一部分的），而前面一部分就是4.1中推导的公式。

对于另一部分针对brdf项的积分在split sum approximation中我们知道该如何分解化简。这里唯一区别就是没有 n.wi ，导致预计算出来的BRDF LUT有些许不同，如下：

但是球谐函数处理低频信息更好（意思是更适合处理漫反射或者是粗糙度比较大的光泽反射），因为对于纯镜面反射而言，我们需要非常高阶的球谐函数才能进行逼近，那计算量实在是太大了(如果需要，你可以直接使用反射光采样天空盒来模拟)。

六、代码实现

我们首先来看一下预计算好的前四阶数据（此处为预计算的天空盒数据，你也可以使用别的天空盒来）

1.12892	1.33143	1.51915
-0.230839	-0.25334	-0.247597
0.5335	0.650908	0.775002
-0.132986	-0.139497	-0.134058
0.127024	0.106568	0.0791303
-0.0833286	-0.0890604	-0.0880843
-0.29704	-0.285707	-0.217588
-0.058678	-0.0539401	-0.0443014
-0.0854982	-0.0481909	-0.0101895
0.025283	0.0135456	0.0011381
0.0629861	0.0414612	0.0173981
0.0223243	0.0187693	0.00901232
-0.228068	-0.236081	-0.226528
-0.00804389	-0.00198417	0.00426891
-0.0345589	-0.020069	-0.0065323
0.0822039	0.0685155	0.0529989
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

可视化成贴图数据后如下：

我们主要看一下渲染物体使用的shader。

顶点着色器：

#version 430 core
layout(location = 0) in vec3 _Position;
layout(location = 1) in vec3 _Normal;
layout(location = 2) in vec2 _TexCoord;

layout(std140, binding = 0) uniform u_Matrices4ProjectionWorld
{
	mat4 u_ProjectionMatrix;
	mat4 u_ViewMatrix;
};

uniform mat4 u_ModelMatrix;

out vec2 v2f_TexCoords;
out vec3 v2f_Normal;
out vec3 v2f_FragPosInViewSpace;

void main()
{

	vec4 FragPosInViewSpace = u_ViewMatrix * u_ModelMatrix * vec4(_Position, 1.0f);
	gl_Position = u_ProjectionMatrix * FragPosInViewSpace;
	v2f_TexCoords = _TexCoord;
	v2f_Normal =  normalize(mat3(transpose(inverse(u_ViewMatrix * u_ModelMatrix))) * _Normal);
	v2f_FragPosInViewSpace = vec3(FragPosInViewSpace);
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

片元着色器：

#version 430 core

in  vec3 v2f_FragPosInViewSpace;
in  vec2 v2f_TexCoords;
in  vec3 v2f_Normal;

layout (location = 0) out vec4 Albedo_;

const float PI = 3.1415926535897932384626433832795;
uniform vec3 u_Coef[16];
uniform vec3 u_DiffuseColor;
uniform sampler2D u_BRDFLut;


vec3 FresnelSchlickRoughness(float cosTheta, vec3 F0, float roughness)
{
    return F0 + (max(vec3(1.0 - roughness), F0) - F0) * pow(max(1.0 - cosTheta, 0.0), 5.0);
}  

void main()
{	
	if((abs(v2f_Normal.x) < 0.0001f) && (abs(v2f_Normal.y) < 0.0001f) && (abs(v2f_Normal.z) < 0.0001f))
	{
		Albedo_ = vec4(0, 0, 0, 1);
		return;
	}

	float Basis[9];
	float x = v2f_Normal.x;
	float y = v2f_Normal.y;
	float z = v2f_Normal.z;
	float x2 = x * x;
	float y2 = y * y;
	float z2 = z * z;
    
    Basis[0] = 1.f / 2.f * sqrt(1.f / PI);
    Basis[1] = 2.0 / 3.0 * sqrt(3.f / (4.f * PI)) * z;
    Basis[2] = 2.0 / 3.0 * sqrt(3.f / (4.f * PI)) * y;
    Basis[3] = 2.0 / 3.0 * sqrt(3.f / (4.f * PI)) * x;
    Basis[4] = 1.0 / 4.0 * 1.f / 2.f * sqrt(15.f / PI) * x * z;
    Basis[5] = 1.0 / 4.0 * 1.f / 2.f * sqrt(15.f / PI) * z * y;
    Basis[6] = 1.0 / 4.0 * 1.f / 4.f * sqrt(5.f / PI) * (-x2 - z2 + 2 * y2);
    Basis[7] = 1.0 / 4.0 * 1.f / 2.f * sqrt(15.f / PI) * y * x;
    Basis[8] = 1.0 / 4.0 * 1.f / 4.f * sqrt(15.f / PI) * (x2 - z2);

	vec3 Diffuse = vec3(0,0,0);
	for (int i = 0; i < 9; i++)
		Diffuse += u_Coef[i] * Basis[i];


	vec3 F0 = vec3(0.2,0.2,0.2);
	float Roughness = 0.2;
	vec3 N = normalize(v2f_Normal);
	vec3 V = -normalize(v2f_FragPosInViewSpace);
	vec3 R = reflect(-V, N); 
	vec3 F        = FresnelSchlickRoughness(max(dot(N, V), 0.0), F0, Roughness);
	vec2 EnvBRDF  = texture(u_BRDFLut, vec2(max(dot(N, V), 0.0), Roughness)).rg;
	vec3 LUT = (F * EnvBRDF.x + EnvBRDF.y);

	x = R.x;
	y = R.y;
	z = R.z;
	x2 = x * x;
	y2 = y * y;
	z2 = z * z;
	Basis[0] = 1.f / 2.f * sqrt(1.f / PI);
    Basis[1] = 2.0 / 3.0 * sqrt(3.f / (4.f * PI)) * z;
    Basis[2] = 2.0 / 3.0 * sqrt(3.f / (4.f * PI)) * y;
    Basis[3] = 2.0 / 3.0 * sqrt(3.f / (4.f * PI)) * x;
    Basis[4] = 1.0 / 4.0 * 1.f / 2.f * sqrt(15.f / PI) * x * z;
    Basis[5] = 1.0 / 4.0 * 1.f / 2.f * sqrt(15.f / PI) * z * y;
    Basis[6] = 1.0 / 4.0 * 1.f / 4.f * sqrt(5.f / PI) * (-x2 - z2 + 2 * y2);
    Basis[7] = 1.0 / 4.0 * 1.f / 2.f * sqrt(15.f / PI) * y * x;
    Basis[8] = 1.0 / 4.0 * 1.f / 4.f * sqrt(15.f / PI) * (x2 - z2);

	vec3 Specular = vec3(0,0,0);
	for (int i = 0; i < 9; i++)
		Specular += u_Coef[i] * Basis[i];


	Albedo_ = vec4(1.0 * Diffuse + 0.5 * Specular * LUT,1.0);
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82

在片元着色器中我们使用了前三阶的数据，你也可以使用更高阶的来模拟高频信息。
更直观的你可以看到普通渲染器与使用SH进行环境光模拟的模型如下所示。

你也可以尝试使用别的贴图数据：

1.50002	1.30591	0.978739
0.147293	0.100121	0.0407086
-0.293406	-0.243371	-0.162667
-0.0326943	-0.0538553	-0.072974
0.0171384	0.00874806	0.00898071
-0.00190327	-0.0483498	-0.0816752
0.109316	0.0956668	0.0549429
-0.0168429	-0.0114185	-0.0111622
-0.0330506	-0.0336814	-0.0358358
-0.0410804	-0.0277784	-0.00791818
0.00992956	0.0083802	0.00807565
0.0160037	0.0151785	0.0067624
-0.00266439	-0.0188293	-0.0268928
-0.000100907	-0.00546308	-0.006031
-0.03704	-0.0292321	-0.017578
0.00265965	-0.000746732	-0.00970712

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

可视化成贴图数据后如下：

运行如下所示：