博客的公式排版非常不好看,可以移步这里阅读pdf版。

如有错误敬请指正!

图形学中PRT使用球谐函数来预计算光照信息。

网上找到的球谐函数的相关数学原理比较凌乱,有些讲的也不是很清楚,所以自己以数学专业的角度整理了一下。

用到的数学知识:

  • 微积分(数学分析)
  • 线性代数
  • 复变函数
  • 泛函分析(较少)

注意:下文中的单个字母上标,如果不是标在括号上,一律不代表“幂次”。

内积空间

考虑空间L^2 (S^2),即,在三维空间单位球表面上平方可积(模的平方在球面S^2上的积分存在且有限)的函数全体。在该空间上定义一个(复)内积(类比于向量的点乘、内积)
\langle f,\ g \rangle := \int_{S^2} f(w)\overline{g(w)}\,{\rm d} w.
这样的话,\forall f\in L^2(S^2)
\Vert f\Vert^2:=\langle f,f\rangle = \int_{S^2} f(w)\overline{f(w)}\,{\rm d} w\ge 0.
即定义了一个范数\Vert\cdot \Vert(类比于向量的长度)。可以注意到,如上定义的内积总是有这样一些性质:
\begin{aligned}
\langle f,g\rangle &= \overline{\langle g,f\rangle}\\
k\langle f,g\rangle &= \langle kf,g\rangle = \langle f,\overline k g\rangle,\ \forall k\in\mathbb C\\
\langle f_1+f_2,g\rangle &= \langle f_1,g\rangle + \langle f_2,g\rangle
\end{aligned}

类比一般的线性空间,我们想找到一组正交系(正交基),这样就可以引入类似于坐标这样的概念。内积空间中的一组规范正交系,记为(e_i) _ { i \in \mathbb N},需要满足:

\langle e_i, e_j\rangle = \delta _{ij},
其中
\begin{aligned}
\delta_{ij}=\begin{cases}
1,& i=j,\\
0,&i\neq j.
\end{cases}
\end{aligned}

连带勒让德多项式

连带(也称伴随)勒让德多项式是一组在L^2[-1,1]上的正交系。连带勒让德多项式记为P_l^m(x),其中0\le m\le lm,l\in\mathbb N。它没有简单的显式表达,但可以递推地求出。它满足
\int_{-1}^1 P_l^m(x)P_k^m(x) \,{\rm d}x= \frac2{2l+1}\frac{(l+m)!}{(l-m)!}\delta_{kl},
即,在m固定时,P_l^mP_k^m正交,其中l\neq k

对于m<0的情形,推广定义:
P_l^{-m}(x) = (-1)^m \frac{(l-m)!}{(l+m)!} P_l^m(x).

### 球谐函数(复数形式)

定义球谐函数如下:
Y_l^m (\theta,\varphi) = N_l^m P_l^{m}(\cos \theta){\rm e}^{{\rm i}m\varphi}, \ \ \ -l\le m\le l,
其中N_l^m是归一化常数,即,N_l^m是使得Y_l^m与自身的内积为1(范数为1)的常数(总是实数):
(N_l^m)^{-2} = \langle P_l^m(\cos \theta){\rm e}^{{\rm i}m\varphi},\ P_l^m(\cos \theta){\rm e}^{{\rm i}m\varphi}\rangle.
考虑内积:
\langle Y_l^m ,Y_k^n\rangle ,
展开计算(注意到\overline{{\rm e}^{{\rm i}n\varphi}}={\rm e}^{-{\rm i}n\varphi}):
\langle Y_l^m ,Y_k^n\rangle=N_l^m N_k^n\int_{S^2} P_l^m(\cos \theta)P_k^n(\cos\theta){\rm e}^{{\rm i}(m-n)\varphi}\,{\rm d} w.
拆开\theta\varphi相关的项
\langle Y_l^m ,Y_k^n\rangle = N_l^m N_k^n \int_{0}^{2\pi} {\rm e}^{{\rm i}(m-n)\varphi}\,{\rm d}\varphi \int_0^\pi P_l^m(\cos\theta)P_k^n(\cos\theta)\sin\theta\,{\rm d}\theta.

  • 如果m\neq n,那么
    \int_{0}^{2\pi} {\rm e}^{{\rm i}(m-n)\varphi}\,{\rm d}\varphi = \int_0^{2\pi} \cos(m-n)\,{\rm d}\varphi+{\rm i}\int_0^{2\pi} \sin(m-n)\,{\rm d}\varphi = 0.
    m=n时被积项退化为1,积分值为2\pi。于是
    \langle Y_l^m ,Y_k^n\rangle = 2\pi \delta_{mn} N_l^m N_k^n \int_0^\pi P_l^m(\cos\theta)P_k^n(\cos\theta)\sin\theta\,{\rm d}\theta.

  • 现在设m=n,由于连带勒让德函数的性质以及t=\cos\theta换元可以立即得出
    \langle Y_l^m ,Y_k^m\rangle = 2\pi N_l^m N_k^m \delta_{lk}.

通过以上讨论可以知道,当且仅当l=km=n时上述内积非0。又由N_l^m定义,我们得到
\langle Y_l^m ,Y_k^n\rangle = \delta_{lk}\delta_{mn},
从而全体复球谐函数成为一组规范正交系。

球谐函数直角坐标表示

引理. Y_l^m(x,y,z)是关于x,y,zl次多项式,可以写成关于zl-|m|次多项式与关于x,y|m|次齐次多项式(x\pm{\rm i}y)^{|m|}的乘积。

Proof. 先考虑m\ge 0的情形。
Y_l^m(\theta,\varphi) =N_l^m P_l^{m}(\cos \theta){\rm e}^{{\rm i}m\varphi},
根据递推公式
P_l^m(t) = (-1)^m (1-t^2)^{m/2} \frac{{\rm d}^m}{{\rm d}t^m}P_l(t),
其中P_l(t)为一般的勒让德多项式。于是
P_l^m(\cos \theta) = (-1)^m \sin^m \theta F(\cos \theta),
其中
F(t) :=\frac{{\rm d}^m}{{\rm d}t^m}P_l(t).
因为勒让德多项式
P_l(t) =\frac{1}{2^l l!}\frac{{\rm d}^l}{{\rm d}t^l}((t^2-1)^l)
是一个关于tl次多项式,从而F(t)是一个关于tl-m次多项式。因为在S^2z=\cos\theta,所以F(\cos\theta)是关于zl-m次多项式。于是
\begin{aligned}
Y_l^m(\theta,\varphi) &=N_l^m (-1)^m F(z) \sin^m \theta {\rm e}^{{\rm i}m\varphi}\\
&=N_l^m (-1)^m F(z) (\sin \theta {\rm e}^{{\rm i}\varphi})^m\\
&=N_l^m (-1)^m F(z) (\sin \theta \cos\varphi+{\rm i}\sin\theta \sin \varphi)^m\\
& = N_l^m (-1)^m F(z) (x+{\rm i}y)^m,
\end{aligned}

这就说明了对于m\ge 0的情况成立。对于m<0的情况,与-m同理可得,结果可写为kG(z)(x-{\rm i}y)^{|m|}。证毕。

无穷维坐标与帕塞瓦尔公式

假设我们有一组希尔伯特(范数完备的内积空间)空间X上的规范正交系(e_i) _ {i\ge 0}。对于任意一个点f \in X,我们可以将f写为如下形式(定义与为什么可以这样做的证明稍繁琐,可以参考任意泛函分析教材,不详细展开了):
f = \sum_{j=0}^\infty a_j e_j,
写成无穷维坐标形式即(a_0, a_1, a_2, \dots). 利用规范正交系的性质我们可以得到
\langle f, e_i\rangle = \sum_{j=0}^\infty a_j\langle e_j, e_i\rangle = a_i,
以及内积空间上的勾股定理的推广帕塞瓦尔等式
\Vert f\Vert^2 = \langle f, f\rangle = \sum_{j=0}^\infty \langle a_je_j, a_je_j\rangle = \sum_{j=0}^\infty |a_j|^2.

对于球谐函数系,即X = L^2(S^2),我们可以写成
f = \sum_{l=0}^\infty \sum_{m=-l}^l a_l^m Y_l^m,
以及帕塞瓦尔等式
\Vert f\Vert^2 = \sum_{l=0}^\infty \sum_{m=-l}^l |a_l^m|^2.

球谐系数
a_l^m = \int f(w) \overline{Y_l^m (w)}\,{\rm d}w.

球谐函数(实数形式)

复数形式的球谐函数一般而言并不好用,就比如上面计算系数的积分,如果球谐函数是实值的,能省下很多麻烦。

我们用球谐函数的复数形式的Y_{1}^1, Y_1^0,Y_1^{-1}举例如何构造实值球谐函数。前面我们已经证明了它们满足规范正交系的性质,我们任取Y_l^ml\neq 1。这时候一定有
\begin{aligned}
\langle Y_l^m, Y_1^1\rangle &= 0,\\
\langle Y_l^m, Y_1^0\rangle &= 0,\\
\langle Y_l^m, Y_1^{-1}\rangle &= 0.
\end{aligned}

也就是说,对于函数f,如果满足
f = c_{-1}Y_1^{-1}+c_0Y_1^0 +c_1 Y_1^1,
那么
\langle Y_l^m ,f\rangle = 0.
如果我们从全体可以写为c_{-1}Y_1^{-1}+c_0Y_1^0 +c_1 Y_1^1的函数(即Y_1^m张成的线性子空间)中,重新选出三个相互正交的、实值且范数为1的函数来替代掉原本的三个Y_1^m,这样构造的“球谐函数族”仍然是规范正交基,虽然它可能不再是物理或数学上普遍使用的球谐函数。

实际上这样做的时候,我们总是把Y_l^mY_l^{-m}这一对函数按照上面的方法重新选择。注意到Y_l^mY_l^{-m}形式非常相似,我们总是可以将它们两个作为一对派生出新的另一对,比如相加之后除以\sqrt2再取虚(实)部,相减之后除以\sqrt 2再取实(虚)部。另外,Y_l^0总是实值的,不需要重新选取。以及,-Y_l^m总是可以替换掉Y_l^m

上面的讨论说明了实值球谐函数的选取方式非常多,主要看实际应用中你想如何定义它。实值球谐函数的例子,可以参考en wiki的实值球谐函数表。

投影

对于任意一个定义在球面S^2上的函数f,我们有分解
f = \sum_{l=0}^\infty \sum_{m=-l}^l a_l^m Y_l^m,
如果我们只取其中的前几层,比如l\le 2,就得到了f对子空间的投影
f_p = \sum_{l=0}^2 \sum_{m=-l}^l a_l^m Y_l^m.
f_p是在l\le 2的所有球谐函数张成的空间里,最接近 f的那个。因为根据帕塞瓦尔等式
\Vert f-f_p\Vert^2 = \sum_{l=3}^\infty \sum_{m=-l}^l |a_l^m|^2,
对于l\le 2球谐函数张成的子空间,这一差值的范数至少是\Vert f-f_p\Vert^2.可以注意到,如果层数l取值越大,f_pf的逼近就越好。

在图形学中我们一般用投影f_p来代替原本的f,这样做一是因为我们不可能算出无穷个球谐系数a_l^m,二是取前几层得到的投影在大多数情况下,已经非常近似于原本的f

球谐函数的旋转

根据施密特正交化,我们总是能轻松构造出一组规范正交基,除了球谐函数,当然也可以选择别的规范正交基来做分解。使用球谐函数是因为它们的另一个优秀性质:快速旋转。

考虑一个旋转Rf与旋转R的复合,仍然是一个定义在S^2的函数。

对于Y_l^mR的复合,必然也有
Y_l^m\circ R(\theta,\varphi) = \sum_{k=0}^\infty \sum_{n=-k}^k b_{l,m,k}^n(R) Y_k^n(\theta, \varphi).
但是,对于球谐函数而言,实际上可以写成
Y_l^m\circ R(\theta,\varphi) = \sum_{n=-l}^l b_{mn} Y_l^n(\theta, \varphi),
其中b_{mn}R有关。也就是说,第l层的球谐函数与任意旋转的复合,仍然在第l层球谐函数张成的子空间内。这一性质与三维旋转群SO(3)的不可约表示密不可分,严格的证明可以参考物理学中的群论(第2版)第122页。

至于如何解出上述的a_l^n,实际上我们可以把它看成解一个2l+1元线性方程组,那么我们自然需要2l+1个方程。取2l+1个不同的球面上的采样点(\theta_i, \varphi_i),\ i=-l,\dots, l,得到2l+1个等式:(有的说法是将采样点投影为球谐系数,这一说法是非常不严谨的。球谐函数分解的是L^2(S^2)上的函数,而不是S^2上的点)
Y_l^m\circ R(\theta_i, \varphi_i) = \sum_{n=-l}^l b_{mn} Y_l^n(\theta_i, \varphi_i),\ i=-l,\dots, l.
记矩阵A = \left(Y_l^n(\theta_i, \varphi_i)\right) _ {i,n},向量b_{m} = (b_{m,-l},\ \dots,\ b_{ml})^{\rm T},向量c_m = (Y_l^m\circ R(\theta_{-l}, \varphi_{-l}), \dots, Y_l^m\circ R(\theta_{l}, \varphi_{l}))^{\rm T}。则得到方程组
Ab_m = c_m.
于是(若A可逆)
b_m = A^{-1} c_m.
这个等式只解出了Y_l^m的旋转。现在假设我们待旋转的f在第l层上的球谐系数是a_l^m,这一层的投影f_p可以分解为
f_p= \sum_{m=-l}^l a_l^m Y_l^m,

f_p\circ R = \sum_{m=-l}^l a_l^m \sum_{n=-l}^l b_{mn} Y_l^n = \sum_{n=-l}^l Y_l^n \sum_{m=-l}^l a_l^m b_{mn}.
把所有列向量b_m合并,写成矩阵形式B = (b_m) _ {m=-l}^l,同理C=(c_m) _ {m=-l}^l,得到
B = A^{-1}C
Y_l^n项新的系数是\sum_{m=-l}^l a_l^m b_{mn},写成一列,即为
Ba_l = A^{-1}C a_l.
实际应用中,矩阵A的值是可以预先算好的,因为A与输入的旋转R没有关联,也就是说,只要选取的2l+1个采样点,能够让算出来的A是可逆矩阵,那么这些采样点就是良好的。我们可以直接计算出A^{-1},需要进行旋转时,只需要用R计算出C,即可得到旋转后的球谐系数。

如果你在别处见到的结果是写成SA^{-1}的,而不是A^{-1}C,那是因为,在上面的讨论中采用的A,相对于那种写法是转置过的,这样的目的在于,最终计算时用的是矩阵左乘列向量(数学上比较通用的写法),而不是矩阵右乘行向量。

其他分解方式

球面上当然也可以选取其他的方式分解,比如Zonal Harmonics以及Wavelets. Zonal Harmonics是一些特殊的球谐函数旋转后的结果。Wavelets不适合快速旋转,但它应用于高频信息的保留时效果比球谐基好。


引き籠り 絶対 ジャスティス!