[笔记]球谐函数
前言
此笔记主要用来记录阅读Stupid Spherical Harmonics (SH) Tricks所总结的要点,主要用来理清自己胡乱学来的球谐函数知识;
看不懂很正常,建议多看几遍,文章中的蓝色链接里面有大佬对相应知识的详细介绍,点进去学习可以解除图形学习人员的很多疑惑;
我自己也不知道看了多少遍,每次看都会有新的理解,并能解除之前几遍所留下来的疑惑;
背景
定义
球谐函数定了一组在球面上的正交基,球面可用S表示,即:
S
=
(
x
,
y
,
z
)
=
(
sin
θ
cos
ϕ
,
sin
θ
sin
ϕ
,
cos
θ
)
S = (x,y,z) = (\sin \theta \cos \phi, \sin \theta \sin \phi, \cos \theta)
S=(x,y,z)=(sinθcosϕ,sinθsinϕ,cosθ)
那么这组正交基则可表示为:
Y
l
m
(
θ
,
ϕ
)
=
K
l
m
e
i
m
ϕ
P
l
∣
m
∣
(
cos
θ
)
,
l
∈
N
,
−
l
<
=
m
<
=
l
Y_l^m(\theta, \phi) = K_l^m \space e^{im\phi} \space P_l^{|m|}(\cos \theta), l \in N,-l<=m<=l
Ylm(θ,ϕ)=Klm eimϕ Pl∣m∣(cosθ),l∈N,−l<=m<=l
这里
P
l
m
P_l^m
Plm表示伴随勒让多德多项式(associated Legendre polynomials ),
K
l
m
K_l^m
Klm表示归一化常数,即:
K
l
m
=
(
2
l
+
1
)
(
l
−
∣
m
∣
)
!
4
π
(
l
+
∣
m
∣
)
!
K_l^m = \sqrt{\frac{(2l+1)(l-|m|)!}{4\pi(l+|m|)!} }
Klm=4π(l+∣m∣)!(2l+1)(l−∣m∣)!
此组正交基实际上是针对复数来进行定义的,在图形学中通常只使用其实数表示形式,即:
y
l
m
=
{
2
R
e
(
Y
l
m
)
,
m
>
0
2
I
m
(
Y
l
m
)
,
m
<
0
Y
l
0
,
m
=
0
=
{
2
K
l
m
cos
m
ϕ
P
l
m
(
cos
θ
)
,
m
>
0
2
K
l
m
sin
∣
m
∣
ϕ
P
l
∣
m
∣
(
cos
θ
)
,
m
<
0
K
l
0
P
l
0
(
cos
ϕ
)
,
m
=
0
y_{l}^{m} = \begin{cases} \sqrt{2} Re(Y_l^m), m>0 \\\sqrt{2} Im(Y_l^m), m<0 \\Y_l^0, m = 0 \end{cases} = \begin{cases} \sqrt{2} K_l^m \space \cos m\phi \space P_l^m (\cos \theta ), m>0 \\\sqrt{2} K_l^m \space \sin |m|\phi \space P_l^{|m|} (\cos \theta ), m<0 \\K_l^0P_l^0(\cos \phi ), m = 0 \end{cases}
ylm=⎩⎪⎨⎪⎧2Re(Ylm),m>02Im(Ylm),m<0Yl0,m=0=⎩⎪⎨⎪⎧2Klm cosmϕ Plm(cosθ),m>02Klm sin∣m∣ϕ Pl∣m∣(cosθ),m<0Kl0Pl0(cosϕ),m=0
其中的L通常称之为“band”,每一个band相当于多项式中相应次数的项,即band0表示常数项,band1表示线性项;
每一个相应的band会有
2
l
+
1
2l+1
2l+1个基函数;
一个n阶的SH展开需要用到band(n-1)内的所有基函数;
在球谐基函数下,使用球谐坐标来进行积分是非常方便的,当然也可以使用xyz的多项式来进行表示;
实际上,在工程使用过程中也确实大多都是用xyz多项式的表达形式,毕竟方向使用xyz来进行表示最为方便;参考StupidSH的附录Appendix A2 Polynomial Forms of SH Basis;
可视化
球谐函数的可视化有多种方式
最常见的就是扭曲一个单位球来进行显示,通过轴向缩放球面上的点来构建模型,球谐函数的值来表示颜色,通常正负用不同颜色来表示;如下图所示;
另外一种可以使用cubemap展开的形式来进行显示,cubamap展开后只需在相应点上显示相应颜色即可,相当于去除了几何建模这一步;如下图所示;
投影与重建
由于球谐函数是正交的,一个定义在S表面下的标量函数f,在球谐函数的最小二次误差投影,可以直接用该函数与相应的球谐基函数的积分来计算;投影结果为球谐基函数下的球谐系数,即:
f
l
m
=
∫
f
(
s
)
y
l
m
(
s
)
d
s
f_l^m = \int f(s)y_l^m(s)ds
flm=∫f(s)ylm(s)ds
该函数的重建为:
f
~
(
s
)
=
∑
l
=
0
n
∑
m
=
−
l
l
f
l
m
y
l
m
(
s
)
\widetilde{f}(s) = \sum_{l=0}^{n}\sum_{m=-l}^{l}f_l^m y_l^m(s)
f
(s)=l=0∑nm=−l∑lflmylm(s)
由于n阶球谐函数需要
n
2
n^2
n2个球谐系数来表示,因此可以转化为一个数组形式的表示,即:
f
~
(
s
)
=
∑
i
=
0
n
2
f
i
y
i
(
s
)
\widetilde{f}(s) = \sum_{i=0}^{n^2}f_i y_i(s)
f
(s)=i=0∑n2fiyi(s)
其中的
f
0
f_0
f0通常代表原函数在整个球域的平均值,常称之为DC term;
基础特性
球谐函数各基函数之间具有正交性;
m=0下的球谐函数,及中间那一列的球谐函数具有绕z轴的旋转对称性;此列基函数常称之为zonal harmonics (ZH) ;
旋转不变性,假如一个球面函数F,通过一个旋转矩阵操作,得到了另外一个球面函数G;则G的球谐系数(G在球谐基函数下的投影)等于F的球谐系数进行同样旋转操作后的数值;
旋转不变性意味着光照围绕着物体旋转将不会产生抖动及撕裂现象;假设原光照为F,对应的球谐系数为f,旋转后的光照为G,旋转操作为R,则直接对f进行R操作即可得到G对应的球谐系数;
卷积
假设一频域核函数h(z)具有圆对称性(绕z轴),那个任一函数的与此函数的卷积结果的球谐系数可以由一下公式进行计算:
(
h
∗
f
)
l
m
=
(
4
π
)
2
l
+
1
h
l
0
f
l
m
(h*f)_l^m = \sqrt{\frac{(4\pi)}{2l+1} }h_l^0 f_l^m
(h∗f)lm=2l+1(4π)hl0flm
个人觉得此性质主要归功于下面内容中的Zonal Harmonics的性质;
后面内容的Irridiance Enviroment Maps就是采用此性质进行推导,详细可参考论文An Efficient Representation for Irradiance Environment Maps;
旋转
球谐函数各band之间具有旋转独立性,意味着将球谐系数转换为旋转后的球谐系数,这个操作过程中,同一band下的球谐系数结果仍需要统一band下的球谐系数来进行变换计算;
通常同一band下的旋转操作使用(2l+1)x(2l+1)的矩阵来进行操作,对于低阶下的计算;高阶下的计算仍需要进行分解计算才更高效;
Zonal Harmonics
假如某一SH投影具有轴对称性质,则称这种SH为ZH,实际上为m=0下的SH;
ZH的n阶表示只需要n个ZH系数即可;
在散射介质的相函数表示中,由于相函数的轴对称性,经常使用ZH来表示;参考PBRT的章节Phase Functions
ZH的旋转比普通SH要简单很多,给定一组ZH系数,旋转至D方向后的ZH系数为:
f
l
m
=
4
π
2
l
+
1
z
l
y
l
m
(
d
)
f_l^m = \sqrt{\frac{4\pi}{2l+1} }z_l y_l^m(d)
flm=2l+14πzlylm(d)
SH Products
n阶球谐函数f与n阶球谐函数g的乘积为p,求p的球谐系数;最终的结果为2n-1阶的球谐函数,每个球谐系数为
p
k
=
∑
i
j
Γ
i
j
k
f
j
g
j
p_{k} = \sum_{ij}^{} \Gamma _{ijk}f_jg_j
pk=ij∑Γijkfjgj
这里
Γ
i
j
k
\Gamma _{ijk}
Γijk是一个三阶张量;
Γ
i
j
k
=
∫
y
i
(
s
)
y
j
(
s
)
y
k
(
s
)
d
s
\Gamma _{ijk} = \int y_i(s)y_j(s)y_k(s)ds
Γijk=∫yi(s)yj(s)yk(s)ds
最终结果的阶数2n-1实际使用中略微偏高,一般会取较低阶数的结果;
如果函数f是固定的的,可以进行Product Matrix的构建,来减少计算的成本;即
(
M
f
)
i
j
=
∑
k
Γ
i
j
k
f
k
(M_f)_{ij} = \sum_{k}^{} \Gamma _{ijk}f_k
(Mf)ij=k∑Γijkfk
Irridiance Enviroment Maps
Irridiance Enviroment Maps的创建过程实际上就是环境贴图与一个clamp cos函数的卷积过程,最后再除以PI进行归一化;实际上就是环境漫反射的预处理贴图,参考LearnOpenGL的章节Diffuse irradiance;
这个过程可以使用SH进行处理及rendering,使用3阶SH就可以足够的精确与高效;
具体思路需要先将Cubemap进行投影(参考下面Projection from Cubemaps),然后再进行一些变换处理(这一部分会考虑cos的影响),转换为光照使用的Light coefficient,才能进行高效的rendering,具体做法参考StupidSH的附录Appendix A10 Shader/CPU code for Irradiance Environment Maps ;
关于这一部分的理论介绍参考论文An Efficient Representation for Irradiance Environment Maps;
Sébastien Lagarde的文章PI or not to PI in game lighting equation也对这一块进行了介绍;
模拟的结果为,使用3阶SH模拟,在光照方向会多出1/16,在反方向多一个-1/16;使用5阶模拟,同样也有一个反向的lobe,正向与反向lobe的比例为31:1;如下图所示,红色表示3阶模拟,蓝色表示5阶模拟;
之所以叫Irridiance Enviroment Maps,是因为对于环境漫反射,需要对球面上的每个点积分整个半球上的Radiance光照作用;相应的概念在皮肤的次表面实时渲染中也有提到Irridiance Maps,参看GPU Gems3中的章节Advanced Techniques for Realistic Real-Time Skin Rendering;
Lighting Models
使用SH可以模拟很多光源模型,最常见的是模拟Cubemap;也可以以解析的形式去模拟一些物理光源,比如一些解析的大气光照模型;
Projection from Cube Maps
透射cubemap至sh伪代码如下:
float f[],s[];
float fWtSum=0;
Foreach(cube map face)
Foreach(texel)
float fTmp = 1 + u^2+v^2;
float fWt = 4/(sqrt(fTmp)*fTmp); // 这里不太懂
EvalSHBasis(texel,s);
f += t(texel)*fWt*s; // vector
fWtSum += fWt;
f *= 4*Pi/fWtSum; // area of sphere // 这里也不太懂
代码中u与v代表面中当前像素的uv坐标,t(texel)表示贴图颜色,texel表示归一化方向;最后的归一化有时可以省略(但是要替换为除以采样的数量),因为低分辨率下,体积角的积分与4pi会有些偏离;
这里对 float fWt = 4/(sqrt(fTmp)fTmp);以及 f = 4Pi/fWtSum;进行解释
fWt按照上述理论应该代表微分体积角,而单个像素对应微分体积角为:
4pi/(6widthwidth) * cos theta / radius^2;
4pi代表整个球面体积角,6widthwidth代表像素个数,cos theta代表当前像素朝向与轴向的夹角,即1/(sqrt(fTmp),radius^2代表当前像素距球心距离的平方,即fTmp;
4pi/(6width*width)在计算的过程中会抵消掉,所以实际上fWt应该为1/(sqrt(fTmp)*fTmp);
实际上微分角的和不一定为4pi,所以f = 4Pi/fWtSum;先除以微分角的和,在乘以4pi,表示归一化;
Analytic Models
除了数值上去积分cubemap,还可以将一些解析光照模型,以解析的方式去计算其SH;
对于Directional Light,其计算会简单一下,但是StupidSH36里面并没有给出具体的结果;平行光的理论模型是一个Delta函数,只在平行光的光照方向光强为1,其它方向光强为0;
对于Spherical Light sources,即球形光源,也是能比较容易的去计算其SH,将Z轴对齐与光源中心及原点连线,就能发现这其实也是一个径向对称模型;此模型下,只有球形光源对应的体积角下存在光照,积分下来的SH结果为附录Appendix A3 ZH Coefficients for Spherical Light Source;
对于Spherical Light sources,可以改进成为Smooth Cone光源,此光源在中心轴处强度为1,在锥形角的边缘强度为0;最终的结果由于精度问题,只能适用于锥夹角大于8度的情况,积分下来的SH结果为附录Appendix A4 ZH Coefficients for Smooth Cone;
从6阶SH的测试结果来看,Smooth Cone的SH拟合结果更加接近假设模型;Spherical Light sources的SH拟合模型会有较大下次,毕竟其光照模型本身就存在高频信息,么得办法;结果如下图,上一行是球形光源,下一层是Smooth Cone:
Normalization
解析光源的SH可以乘以一个归一化系数,来使得最终的结果与clamp cos的积分值为1;
clamp cos函数对应的SH相对好求一些,由于clamp cos是径向对称的,一次可以使用ZH去表示,前6个band为:
1
2
π
,
3
3
π
,
5
8
π
,
0
,
−
1
16
π
,
0
\frac{1}{2\sqrt{\pi } } ,\frac{\sqrt{3} }{3\sqrt{\pi } } ,\frac{\sqrt{5} }{8\sqrt{\pi } } ,0,\frac{-1 }{16\sqrt{\pi } } ,0
2π1,3π3,8π5,0,16π−1,0
对于Smooth Cone,其归一化系数为$\frac{1}{sin^2a}
;
对
于
平
行
光
,
其
归
一
化
系
数
为
; 对于平行光,其归一化系数为
;对于平行光,其归一化系数为\frac{16\pi}{17} $;
Extracting Lights from SH
对于一个SH,是可以抽离出一个平行光成分加一个环境光成分的,也可以抽离出多个光照成分;
这一节看不太懂,后面再看
Ringing
当使用SH去拟合不连续函数时,在不连续的地方,就会出现overshoot与undershoot现象;当使用截断的SH去拟合连续函数时,也会出现同样的问题;
有两种方法来解决这个问题:
- 使用windowing技术;
- 最小化变分函数的其它形式,而不仅仅是平方误差;
Windowing
所谓Windowing技术,即对截断的SH内,针对不同的结束进行平滑缩放,并在截断频率内缩放为0;
常用的是Lanczos window(
s
i
n
π
x
w
π
x
w
\frac {sin \frac {\pi x} w} {\frac {\pi x} w}
wπxsinwπx)与Hanning window(
1
+
c
o
s
π
x
w
2
\frac {1 + cos \frac {\pi x} w} {2}
21+coswπx),相比Lanczos,Hanning会具有更好的模糊效果;两者的衰减图像如下所示,红色为Lanczos window,蓝色为Hanning window:
最小化误差函数
在游戏以及图形领域,一般都使用低阶的SH,所以很少去使用这种方法来减小ringing问题;
该方法的一种简单形式是寻找新的SH系数表示形式,来尽可能接近平方误差(满足此条件的就是原生的SH系数)的同时,还要最小化weighted squared Laplacian;常用的逼近形式为:
c
l
m
=
f
l
m
(
1
+
λ
l
2
(
l
+
1
)
2
c_l^m = \frac {f_l^m} {(1+\lambda l^2(l+1)2}
clm=(1+λl2(l+1)2flm
这里新产生的windowing函数依赖于
λ
\lambda
λ的数值,当
λ
\lambda
λ为0时,即为原生的SH系数;
一般来说计算irradiance environment maps时,windowing一般是不需要的;但是对于真实的场景中,法线的变化都是不平滑的,这是ringring现象就比较严重;另外,对于使用HDR时的irradiance environment maps,由于亮光处的问题,使用windowing有时也是必要的;
另一个经常伴随着ringring现象的问题是color artifacts,如下图所示,为一黄色平行光与一白色环境光产生的Irradiance enviroment maps,可以明显的看到未使用windowing产生的color artifact;
Content Sensitive Windowing
由于windowing技术是全局的,这样就会导致消除ringing问题的同时,也会模糊没有ringing问题的区域;一个比较好的方法是在整个球形区域的部分区域使用原生的SH,有ring问题的区域使用windowing SH;
对于一个平行光,可以采用下面的公式来决定原生SH区域所占的权重:
w
a
=
(
m
a
x
(
0
,
(
n
⋅
r
)
−
c
t
1
−
c
t
)
)
p
w_a = (max(0, \frac {(n·r)-c_t} {1-c_t}))^p
wa=(max(0,1−ct(n⋅r)−ct))p
式中的
c
t
c_t
ct与
p
p
p分别决定了blend的位置与blend的范围;
对于其他复杂的light,可以采用前面所说的光源抽离,来算出大约的主平行光,然后再使用上面的公式进行blend;
使用主平行光抽离的方法对于静态的光源没有较大问题,但对于动态光源,则可能会产生temporal artifact;
SH Products
使用SH来去计算两个函数的乘积是非常的有用的,使用场景有:
- 在一个skylight model上添加移动的物体;
- 光照时乘以可见性函数;
- 缩放或修改一个SH light probe,乘以0-1的某常数用来模拟云;
SH product的计算需要计算三阶张量,这里只谈论一些特殊的情况;
Products with varying Orders
在product的阶数比较低得情况下,可以极大程度的降低计算的复杂度;