当前位置：首页 > 篮球资讯 > 正文内容

球谐光照——球谐函数

杏彩体育2年前 (2022-12-17)篮球资讯60

早在1877年，Norman Macleod Ferrers就专门写了一本书来介绍球谐函数，后面物理学家把实数球谐函数扩展到复平面上，在复变函数论中作为“特殊函数”来研究，它在物理以及计算化学上有重要的应用，我们主要讨论它在计算机图形渲染上的应用。

球谐函数是拉普拉斯方程的分离 $r$ r变量后，角度部分通解的正交项，那本篇文章就从拉普拉斯方程开始介绍，直至找到我们想要的球谐函数。球谐函数有复数形式和实数形式，我们只关心它的实数形式。

球谐函数有两条重要的性质，正交完备性和旋转不变性。球谐函数构成的函数组，作为正交基，对信号进行投影和重建，例如[7]介绍的辐射度环境贴图，通过9个系数就可以模拟一张环境贴图的漫反射信号，无论是存储还是计算，都有显著的优势。当然，本篇文章希望从数学的角度，来介绍球谐函数，也受限于个人的数学能力，会忽略了一些复杂的推导过程。

文章目录：

球谐始源球谐性质附录参考

球谐始源

球谐函数是拉普拉斯方程的分离 $r$ r变量后，角度部分通解的正交项。这部分将从拉普拉斯方程开始逐项推导，得到最终我们想要知道的球谐函数。内容主要参考“姚端正,梁家宝. 数学物理方法-第4版”和“顾樵. 数学物理方法”两本书。

球面坐标可以表示为

$(1){x=rsinθcosφy=rsinθsinφz=rcosθ$ \left\{ \begin{matrix} x=r\sin \theta \cos \varphi \\ y=r\sin \theta \sin \varphi \\ z=r\cos \theta \\ \end{matrix} \right. \tag{1}\\

三维空间下的拉普拉斯（Laplace）方程可以表示为

$(2)\nabla2=\partial2\partialx2+\partial2\partialy2+\partial2\partialz2=0$ {{\nabla }^{\text{2}}}\text{=}\frac{{{\partial }^{2}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}}{\partial {{z}^{2}}}=0 \tag{2}\\

把球面坐标代入拉普拉斯方程（参见“Laplacian in Spherical Coordinates”），可以得到

$(3)1r2\partial\partialr(r2\partialf\partialr)+1r2sinθ\partial\partialθ(sinθ\partialf\partialθ)+1r2sin2θ\partial2f\partialφ2=0$ \frac{1}{{{r}^{2}}}\frac{\partial }{\partial r}\left( {{r}^{2}}\frac{\partial f}{\partial r} \right)+\frac{1}{{{r}^{2}}\sin \theta }\frac{\partial }{\partial \theta }\left( \sin \theta \frac{\partial f}{\partial \theta } \right)+\frac{1}{{{r}^{2}}{{\sin }^{2}}\theta }\frac{{{\partial }^{2}}f}{\partial {{\varphi }^{2}}}=0 \tag{3}\\

我们的目标就是求解拉普拉斯方程的解。首先，把表示距离的变数 $r$ r跟表示方向的变数 $θ$ \theta和 $φ$ \varphi分离，即

$(4)f(r,θ,φ)=R(r)Y(θ,φ)$ f\left( r,\theta ,\varphi \right)=R\left( r \right)Y\left( \theta ,\varphi \right) \tag{4}\\

$Y(θ,φ)$ Y\left( \theta ,\varphi \right)表示角度部分，把等式（4）代入等式（3）可以得到

$Yr2\partial\partialr(r2\partialR\partialr)+Rr2sinθ\partial\partialθ(sinθ\partialY\partialθ)+Rr2sin2θ\partial2Y\partialφ2=0$ \frac{Y}{{{r}^{2}}}\frac{\partial }{\partial r}\left( {{r}^{2}}\frac{\partial R}{\partial r} \right)+\frac{R}{{{r}^{2}}\sin \theta }\frac{\partial }{\partial \theta }\left( \sin \theta \frac{\partial Y}{\partial \theta } \right)+\frac{R}{{{r}^{2}}{{\sin }^{2}}\theta }\frac{{{\partial }^{2}}Y}{\partial {{\varphi }^{2}}}=0\\

等式两边乘以 $r2/RY$ {{r}^{2}}/RY，并移项，可得

$1R\partial\partialr(r2\partialR\partialr)=-1Ysinθ\partial\partialθ(sinθ\partialY\partialθ)-1Ysin2θ\partial2Y\partialφ2$ \frac{1}{R}\frac{\partial }{\partial r}\left( {{r}^{2}}\frac{\partial R}{\partial r} \right)=-\frac{1}{Y\sin \theta }\frac{\partial }{\partial \theta }\left( \sin \theta \frac{\partial Y}{\partial \theta } \right)-\frac{1}{Y{{\sin }^{2}}\theta }\frac{{{\partial }^{2}}Y}{\partial {{\varphi }^{2}}}\\

左边是 $r$ r的函数，跟 $θ$ \theta和 $φ$ \varphi无关；右边是 $θ$ \theta和 $φ$ \varphi的函数，跟 $r$ r无关。两边相等，显然是不可能的。除非两边实际上是同一个常数，通常把这个参数记为 $l(l+1)$ l\left( l+1 \right)。即

$1R\partial\partialr(r2\partialR\partialr)=-1Ysinθ\partial\partialθ(sinθ\partialY\partialθ)-1Ysin2θ\partial2Y\partialφ2=l(l+1)$ \frac{1}{R}\frac{\partial }{\partial r}\left( {{r}^{2}}\frac{\partial R}{\partial r} \right)=-\frac{1}{Y\sin \theta }\frac{\partial }{\partial \theta }\left( \sin \theta \frac{\partial Y}{\partial \theta } \right)-\frac{1}{Y{{\sin }^{2}}\theta }\frac{{{\partial }^{2}}Y}{\partial {{\varphi }^{2}}}=l\left( l+1 \right)\\

这就分解为两个方程

$\partial\partialr(r2\partialR\partialr)-l(l+1)R=0$ \frac{\partial }{\partial r}\left( {{r}^{2}}\frac{\partial R}{\partial r} \right)-l\left( l+1 \right)R=0\\

$(5)1sinθ\partial\partialθ(sinθ\partialY\partialθ)+1sin2θ\partial2Y\partialφ2+l(l+1)Y=0$ \frac{1}{\sin \theta }\frac{\partial }{\partial \theta }\left( \sin \theta \frac{\partial Y}{\partial \theta } \right)+\frac{1}{{{\sin }^{2}}\theta }\frac{{{\partial }^{2}}Y}{\partial {{\varphi }^{2}}}+l\left( l+1 \right)Y=0 \tag{5}\\

当然，我们只关心角度部分的解，即方程（5）的解，这个方程也称为球函数方程。进一步采用分离变数法，以

$Y(θ,φ)=Θ(θ)Φ(φ)$ Y\left( \theta ,\varphi \right)=\Theta \left( \theta \right)\Phi \left( \varphi \right)\\

代入球函数方程，得

$Φsinθ\partial\partialθ(sinθ\partialΘ\partialθ)+Θsin2θ\partial2Φ\partialφ2+l(l+1)ΘΦ=0$ \frac{\Phi }{\sin \theta }\frac{\partial }{\partial \theta }\left( \sin \theta \frac{\partial \Theta }{\partial \theta } \right)+\frac{\Theta }{{{\sin }^{2}}\theta }\frac{{{\partial }^{2}}\Phi }{\partial {{\varphi }^{2}}}+l\left( l+1 \right)\Theta \Phi =0\\

在方程两边乘以 $sin2θ/ΦΘ$ {{\sin }^{2}}\theta /\Phi \Theta并移项，即可得

$sinθΘ\partial\partialθ(sinθ\partialΘ\partialθ)+l(l+1)sin2θ=-1Φ\partial2Φ\partialφ2$ \frac{\sin \theta }{\Theta }\frac{\partial }{\partial \theta }\left( \sin \theta \frac{\partial \Theta }{\partial \theta } \right)+l\left( l+1 \right){{\sin }^{2}}\theta =-\frac{1}{\Phi }\frac{{{\partial }^{2}}\Phi }{\partial {{\varphi }^{2}}}\\

左边是 $θ$ \theta的函数，跟 $φ$ \varphi无关；右边是 $φ$ \varphi 的函数，跟 $θ$ \theta 无关。两边相等显然是不可能的，除非两边实际上是同一个常数，这个常数记作 $λ$ \lambda，

$sinθΘ\partial\partialθ(sinθ\partialΘ\partialθ)+l(l+1)sin2θ=-1Φ\partial2Φ\partialφ2=λ$ \frac{\sin \theta }{\Theta }\frac{\partial }{\partial \theta }\left( \sin \theta \frac{\partial \Theta }{\partial \theta } \right)+l\left( l+1 \right){{\sin }^{2}}\theta =-\frac{1}{\Phi }\frac{{{\partial }^{2}}\Phi }{\partial {{\varphi }^{2}}}=\lambda\\

这就又分解为两个常微分方程

$(7)\partial2Φ\partialφ2+λΦ=0$ \frac{{{\partial }^{2}}\Phi }{\partial {{\varphi }^{2}}}+\lambda \Phi =0 \tag{7}\\

$(8)sinθ\partial\partialθ(sinθ\partialΘ\partialθ)+[l(l+1)sin2θ-λ]Θ=0$ \sin \theta \frac{\partial }{\partial \theta }\left( \sin \theta \frac{\partial \Theta }{\partial \theta } \right)+\left[ l\left( l+1 \right){{\sin }^{2}}\theta -\lambda \right]\Theta =0 \tag{8}\\

对于常微分方程（7），它有一个隐含的“自然的周期条件”（ $Φ(φ+2π)=Φ(φ)$ \Phi \left( \varphi +2\pi \right)=\Phi \left( \varphi \right)），两者构成本征值问题。即

$(9)λ=m2,(m=0,\pm1,\pm2,\dots)$ \lambda ={{m}^{2}},\left( m=0,\pm 1,\pm 2,\cdots \right) \tag{9}\\

它的周期解用复数形式，可以表示为

$(10)Φ(φ)=eimφ,m=0,\pm1,\pm2,\dots$ \Phi \left( \varphi \right)={{e}^{im\varphi }},m=0,\pm 1,\pm 2,\cdots \tag{10}\\

严格来说，这里忽略了常数项，且是两个cos和sin函数的混合，或者是两个正负虚数的混合，但不影响最终的通解。

复数形式在复变函数论里面，有一个如下所示的转换关系，它们都是复平面上坐标的不同表现形式。由于本篇文章并不是为了严格的数学理论推导，而是为了梳理球谐函数在计算机上从理论到应用这条线，最终我们采用的是实数形式的球谐函数，所以看看就可以了。

$cosmφ+isinmφ=eimφ$ \cos m\varphi +i\sin m\varphi ={{e}^{im\varphi }}\\

再看常微分方程（8），它可以改写为

$1sinθ\partial\partialθ(sinθ\partialΘ\partialθ)+[l(l+1)-m2sin2θ]Θ=0$ \frac{1}{\sin \theta }\frac{\partial }{\partial \theta }\left( \sin \theta \frac{\partial \Theta }{\partial \theta } \right)+\left[ l\left( l+1 \right)-\frac{{{m}^{2}}}{{{\sin }^{2}}\theta } \right]\Theta =0\\

设 $x=cosθ$ x=\cos \theta,则

$\partialΘ\partialθ=\partialΘ\partialx\partialx\partialθ=-sinθ\partialΘ\partialx$ \frac{\partial \Theta }{\partial \theta }=\frac{\partial \Theta }{\partial x}\frac{\partial x}{\partial \theta }=-\sin \theta \frac{\partial \Theta }{\partial x}\\

代入上式，化简可得

$(11)(1-x2)\partial2Θ\partialx2-2x\partialΘ\partialx+[l(l+1)-m21-x2]Θ=0$ \left( 1-{{x}^{2}} \right)\frac{{{\partial }^{2}}\Theta }{\partial {{x}^{2}}}-2x\frac{\partial \Theta }{\partial x}+\left[ l\left( l+1 \right)-\frac{{{m}^{2}}}{1-{{x}^{2}}} \right]\Theta =0 \tag{11}\\

这个方程就是l次连带勒让德方程，也称为缔合勒让德方程，其中， $m=0$ m=0的特例，即

$(12)(1-x2)\partial2Θ\partialx2-2x\partialΘ\partialx+l(l+1)Θ=0$ \left( 1-{{x}^{2}} \right)\frac{{{\partial }^{2}}\Theta }{\partial {{x}^{2}}}-2x\frac{\partial \Theta }{\partial x}+l\left( l+1 \right)\Theta =0 \tag{12}\\

叫作l次勒让德方程。

后面只考虑连带勒让德方程，它的解就称为连带勒让德函数，只有当 $λ=l(l+1)$ \lambda =l\left( l+1 \right)， $l=0,1,\dots$ l=0,1,\cdots 时才有有界周期解，用 $Plm(x)$ P_{l}^{m}\left( x \right)表示，即

$(13)Θ(θ)=Plm(cosθ){m=0,\pm1,\dots,\pml}$ \Theta \left( \theta \right)=P_{l}^{m}\left( \cos \theta \right)\left\{ m=0,\pm 1,\cdots ,\pm l \right\} \tag{13}\\

经过数学家论证，连带勒让德函数表示为

Pml(x)=(−1)m(1−x2)m22ll!dl+mdxl+m(x2−1)l

P_{l}^{m}\left( x \right)=\frac{{{\left( -1 \right)}^{m}}{{\left( 1-{{x}^{2}} \right)}^{\frac{m}{2}}}}{{{2}^{l}}l!}\frac{{{d}^{l+m}}}{d{{x}^{l+m}}}{{\left( {{x}^{2}}-1 \right)}^{l}} \tag{14}\\

这叫称为l次m阶连带勒让德函数，“次”的英文是“degree”，“阶”的英文是“order”。当l">m>l。连带勒让德函数里面有一个m+l次导数计算，在计算机上这个很难处理，但是有递归关系[3]，即

\left\{ \begin{aligned} & \left( l-m \right)P_{l}^{m}\left( x \right)=x\left( 2l-1 \right)P_{l-1}^{m}\left( x \right)-\left( l+m-1 \right)P_{l-2}^{m}\left( x \right) \\ & P_{m}^{m}\left( x \right)={{\left( -1 \right)}^{m}}\left( 2m-1 \right)!!{{\left( 1-{{x}^{2}} \right)}^{m/2}} \\ & P_{m+1}^{m}\left( x \right)=x\left( 2m+1 \right)P_{m}^{m}\left( x \right) \\ \end{aligned} \right. \tag{15}\\

两个!!表示双阶乘，即\left( 2m-1 \right)!!=1\cdot 3\cdot 5\cdots \left( 2m-1 \right)。

连带勒让德函数的递归关系，保证了计算机实现的基础。

此外，再给一个l次m阶连带勒让德函数的关系等式

P_{l}^{m}\left( x \right)={{\left( -1 \right)}^{m}}\frac{\left( l+m \right)!}{\left( l-m \right)!}P_{l}^{-m}\left( x \right) \tag{16}\\

回到球函数方程（5）的求解，它的Y\left( \theta ,\varphi \right)通解的复数形式表示为

Y\left( \theta ,\varphi \right)=\sum\limits_{l=0}^{\infty }{\sum\limits_{k=-l}^{l}{P_{l}^{k}\left( \cos \theta \right){{e}^{im\varphi }}}},m=0,\pm 1,\pm 2,\cdots \tag{17}\\

严格来说，由于\Phi \left( \varphi \right)忽略了常数项，这里也是忽略常数项的情况。

一般的l次m阶球谐函数{{Y}_{lm}}\left( \theta ,\varphi \right)的复数形式可以表示为

{{Y}_{lm}}\left( \theta ,\varphi \right)={{P}_{lm}}\left( \cos \theta \right){{e}^{im\varphi }},m=0,\pm 1,\pm 2,\cdots \tag{18}\\

l表示球谐函数的次数，m表示球谐函数的阶数

球谐函数的模长可以表示为

{{\left( N_{l}^{m} \right)}^{2}}={{\iint\limits_{S}{{{Y}_{lm}}\left( x \right)\left[ {{Y}_{lm}}\left( x \right) \right]}}^{*}}\sin \theta d\theta d\varphi =\frac{2}{2l+1}\frac{\left( l+\left| m \right| \right)!}{\left( l-\left| m \right| \right)!}2\pi\\

归一化的球谐函数Y_{l}^{m}\left( \theta ,\varphi \right)的复数形式可以表示为

Y_{l}^{m}\left( \theta ,\varphi \right)=K_{l}^{m}{{Y}_{lm}}\left( \theta ,\varphi \right) \tag{19}\\

其中

K_{l}^{m}=\frac{1}{N_{l}^{m}}=\sqrt{\frac{2l+1}{4\pi }\frac{\left( l-\left| m \right| \right)!}{\left( l+\left| m \right| \right)!}} \tag{20}\\

注意区分两种数学表示的含义，{{Y}_{lm}}\left( \theta ,\varphi \right)表示一般形式的球谐函数，Y_{l}^{m}\left( \theta ,\varphi \right)表示归一化的球谐函数。

当m>0时采用实数cos部分，当m<0时采用虚数sin部分，则归一化的球谐函数的实数形式可以表示为

0 \\ \sqrt{2}K_{l}^{m}\sin \left( -m\varphi \right)P_{l}^{-m}\left( \cos \theta \right) & m<0>Y_{l}^{m}\left( \theta ,\varphi \right)=\left\{ \begin{matrix} \sqrt{2}K_{l}^{m}\cos \left( m\varphi \right)P_{l}^{m}\left( \cos \theta \right) & m>0 \\ \sqrt{2}K_{l}^{m}\sin \left( -m\varphi \right)P_{l}^{-m}\left( \cos \theta \right) & m<0 \\ K_{l}^{0}P_{l}^{m}\left( \cos \theta \right) & m=0 \\ \end{matrix} \right. \tag{21}

根据上述计算公式，可以得到前面4次的球谐函数，为

前4次球谐函数

参见[6]，已经推导出前6次的球谐函数。

球谐函数可视化，前面几次的三维图像如图1所示

图1. 前4次的球谐函数三维图像

至此，你应该理解两条重要结论：

球谐函数是拉普拉斯方程分离r变量后，角度部分通解的正交项（后面介绍正交性）。如何计算球谐函数，可以参见等式（15）（20）（21）。

球谐性质

接着，讨论归一化的球谐函数的性质，它具备两条重要的性质构成了它应用的基石：

正交完备性旋转不变性

正交完备性

对于任意两个归一化的球谐函数在球面上的积分有

\iint\limits_{S}{Y_{l}^{m}\left( \theta ,\varphi \right)Y_{k}^{n}\left( \theta ,\varphi \right)\sin \theta d\theta d\varphi }=\left\{ \begin{matrix} 0 & m\ne n,or,l\ne k \\ 1 & m=n,l=k \\ \end{matrix} \right. \tag{21}\\

这就表示由球谐函数构成的函数组\left\{ Y_{l}^{m}\left( \theta ,\varphi \right) \right\}是正交归一化的。

以某一正交归一函数组为基，把一个给定的函数用这些函数的线性组合来表示，这就是一种重要的展开，这种用正交函数组展开为级数的一个显著的例子就是傅里叶变换。

任意一个球面函数f\left( \theta ,\varphi \right)可以用正交归一的球函数Y_{l}^{m}\left( \theta ,\varphi \right)进行展开，这种展开类似于傅里叶展开，称为广义傅里叶展开

f\left( \theta ,\varphi \right)=\sum\limits_{l=0}^{\infty }{\sum\limits_{m=-l}^{l}{C_{l}^{m}Y_{l}^{m}\left( \theta ,\varphi \right)}} \tag{22}\\

其中，广义傅里叶系数C_{l}^{m}为

C_{l}^{m}=\int_{0}^{2\pi }{\int_{0}^{\pi }{f\left( \theta ,\varphi \right)Y_{l}^{m}\left( \theta ,\varphi \right)\sin \theta d\theta d\varphi }} \tag{23}\\

当次数l\to \infty 的时候，展开的级数和会平均收敛于f\left(\theta ,\varphi \right)。换句话说，当次数l越大，那么级数和就会越趋近于被展开的函数f\left(\theta ,\varphi \right)，就称\left\{ Y_{l}^{m}\left( \theta ,\varphi \right) \right\}为完备函数组。平均收敛，并不代表收敛，只是表示趋近于的含义。

从计算机的角度来说，如等式（22）所示的函数展开，n的取值不可能是无穷大，往往取一个给定系数，则可以确定球谐函数组，例如n=2，那么球谐函数组就是

\left\{ Y_{l}^{m}\left( \theta ,\varphi \right) \right\}=\left\{ Y_{0}^{0},Y_{1}^{-1},Y_{1}^{0},Y_{1}^{1} \right\}\\

任意给定n，得到的球谐函数组的个数为

S=1+\text{3}+5+\cdots 2n-1={{n}^{2}}\\

那么，广义傅里叶系数相当于这样一个排列

C_{0}^{0},C_{1}^{-1},C_{1}^{0},C_{1}^{1},C_{2}^{-2},C_{2}^{-1},\cdots\\

类似的球谐函数也可以构成这样一个类似的排列，若我们用一个普通的系数{{c}_{k}}来表示上面的广义傅里叶系数，用一个函数{{y}_{k}}\left( \theta ,\varphi \right)来表示球谐函数，那么等式（22）可以换成另外一种形式

f\left( \theta ,\varphi \right)=\sum\limits_{k=0}^{{{n}^{2}}-1}{{{c}_{k}}{{y}_{k}}\left( \theta ,\varphi \right)} \tag{24}\\

这种形式的展开与等式（22）是完全一样的，它只是把球谐函数的次数展开，用一个系数来表示，但是它隐藏了一个条件：取的系数个数必需是{{n}^{2}}。

回过来，球谐函数组相当于一组正交基，将函数f\left( \theta ,\varphi \right)表示为这组正交基的线性组合，生成线性组合系数的过程就称为投影（Projection），例如一个函数可以表示为

f\left( \theta ,\varphi \right)\approx aY_{0}^{0}+bY_{1}^{-1}+cY_{1}^{0}+dY_{1}^{1}\\

生成系数\left\{ a,b,c,d \right\}的过程，就是投影，等式（23）就确定了投影的方法。相反，利用这组系数和正交基组合，得到原函数的过程，就称为重建（Reconstruction）。

投影的过程就是计算函数积分，计算消耗较大，可以采用离线处理来生成广义傅里叶系数；在实时渲染时，就要简单的线性组合，就可以重建原始函数。当然，由于是有限个系数，就必然存在误差。

我们再来看下连带勒让德函数，如图2所示，随着次数的增加，函数的振动频率会越快。对于函数的展开而言，振动频率越大的基底，它就只能表示越高频的信息，往往一个函数里面的高频信息量是较少的。

图2. 连带勒让德函数曲线图，次数越大，振频越快

类似的，球谐函数也具备这种随着次数增加，振动频率增加的的特性。它就使得n的取值不需要很大时，就可以得到很好的重建效果，当然只能还原出低频信息。根据Robin[3]得到的数据，如图3所示，当n > 6时，就能还原出整体效果，但是边缘棱角这些高频信息是无法还原出来的。

图3. 不同n取值下球谐函数的重建效果

由球谐函数构成的函数组构成正交归一的基底，对球面上的函数进行投影和重建，也就是广义傅里叶展开，数学上的完备性，保证了展开结果会趋近于被展开的函数。

旋转不变性

第一个问题是：什么叫旋转不变性。

任意一个球面上的函数f\left( \theta ,\varphi \right)可以用球谐函数组作为基底展开，需要根据等式（23）计算广义傅里叶系数C_{l}^{m}。如果我们对原函数进行旋转操作的话，设旋转变换表示为R\left( \theta ,\varphi \right)，我们就得到了一个新的函数f\left( R\left( \theta ,\varphi \right) \right)。对新函数进行展开的话，我们需要重新计算广义傅里叶系数，设为B_{l}^{m}，这就有点为难了。在图形渲染中，广义傅里叶系数的生成是离线实现的，它的消耗很大，这就表示，一旦光源发生了旋转后，由于原函数的改变导致提前生成的系数失效。旋转不变性，表示原函数发生了旋转，只需要对生成的广义傅里叶系数进行变换，就能保证变换后的系数能等价还原出新函数。在图形渲染上的表现就是，当光源发生旋转后，我们只要同步的计算出变换后的广义傅里叶系数，就能保证画面的光照效果不会抖动跳变。旋转不变性，并不是表示源函数发生旋转后，对重建结果没有影响，而是表示通过对系数与匹配的旋转进行变换后，能等价的还原出旋转后的函数。

举个[6]实验的例子，如图4所示，球谐函数表示的光照发生旋转后，仍然能等价重建新的变换函数，但是采用Ambient Cube的方法效果就出现了异常。

图4. SH表示球谐基底，HL2表示Ambient Cube基底

第二个问题是：怎么对生成的系数进行变换。

针对这个问题，这里写些自己的理解，不做深入的研究。

对于l次球谐函数，就会有2l + 1个系数，表示为

{{C}_{l}}=\left\{ C_{l}^{-l},C_{l}^{-l+1},\cdots ,C_{l}^{l-1},C_{l}^{l} \right\} \tag{25}\\

设变换矩阵为R_{SH}^{l}，它是一个(2l + 1)*(2l + 1)的矩阵，那么系数的变换就可以表示为

B_{l}^{m}=\sum\limits_{k=-l}^{k=l}{M_{l}^{m,k}C_{l}^{k}} \tag{26}\\

或者用向量与矩阵的乘积形式，表示为

{{B}_{l}}={{C}_{l}}\cdot R_{SH}^{l} \tag{27}\\

那么，经过旋转变换后的函数的展开就可以表示为

f\left( R\left( \theta ,\varphi \right) \right)=\sum\limits_{l=0}^{\infty }{\sum\limits_{m=-l}^{l}{B_{l}^{m}Y_{l}^{m}\left( \theta ,\varphi \right)}} \tag{28}\\

唯一的区别就是，系数由{{C}_{l}}变成了{{B}_{l}}。

想表达的一点是，系数的变换是基于球谐函数的次数，即第3次球谐函数的系数{{B}_{3}}，只能由第3次球谐函数的系数{{C}_{3}}变换而来。

若取前3次的球谐函数构成正交基，函数组共有0，1，2次三类球谐函数，若采用等式（24）的形式，则3个子矩阵需要整合成一个完整的变换矩阵，对于前3次球谐函数的例子，就组成一个9x9的变换矩阵，它的形状如下所示。

9x9的球谐系数变换矩阵

考虑低维的情况[3]。旋转可以用旋转矩阵、欧拉角、四元数等方式表示，任意一个旋转矩阵R可以用{{Z}_{\alpha }}{{Y}_{\beta }}{{Z}_{\gamma }}型的欧拉角表示，它们间的变换关系[5]表示为

\left( \begin{matrix} {{R}_{0,0}} & {{R}_{0,1}} & {{R}_{0,2}} \\ {{R}_{1,0}} & {{R}_{1,1}} & {{R}_{1,2}} \\ {{R}_{2,0}} & {{R}_{2,1}} & {{R}_{2,2}} \\ \end{matrix} \right)=\left( \begin{matrix} {{c}_{\alpha }}{{c}_{\beta }}{{c}_{\gamma }}-{{s}_{\alpha }}{{s}_{\gamma }} & {{c}_{\alpha }}{{s}_{\gamma }}+{{s}_{\alpha }}{{c}_{\beta }}{{c}_{\gamma }} & -{{s}_{\beta }}{{c}_{\gamma }} \\ -{{s}_{\alpha }}{{c}_{\gamma }}-{{c}_{\alpha }}{{c}_{\beta }}{{s}_{\gamma }} & {{c}_{\alpha }}{{c}_{\gamma }}-{{s}_{\alpha }}{{c}_{\beta }}{{s}_{\gamma }} & {{s}_{\beta }}{{s}_{\gamma }} \\ {{c}_{\alpha }}{{s}_{\beta }} & {{s}_{\alpha }}{{s}_{\beta }} & {{c}_{\beta }} \\ \end{matrix} \right) \tag{29}\\

其中，c表示cos，s表示sin。

有了这个变换关系后，就很容易计算出欧拉角\alpha ,\beta ,\gamma，表示为

\begin{aligned} & \sin \beta =\sqrt{1-R_{2,2}^{2}} \\ & \left\{ \begin{aligned} & \alpha \text{=atan2f}\left( {{R}_{2,1}}/\sin \beta ,{{R}_{2,0}}/\sin \beta \right) \\ & \beta =\text{atan2f}\left( \sin \beta ,{{R}_{2,2}} \right) \\ & \gamma =\text{atan2f}\left( {{R}_{1,2}}/\sin \beta ,-{{R}_{0,2}}/\sin \beta \right) \\ \end{aligned} \right. \\ \end{aligned}\\

对于，{{R}_{2,2}}=1的退化情况，欧拉角表示为

\left\{ \begin{aligned} & \alpha \text{=atan2f}\left( {{R}_{0,1}},{{R}_{0,0}} \right) \\ & \beta =0 \\ & \gamma =0 \\ \end{aligned} \right.\\

那么，相应l次的球谐系数的变换矩阵可以表示为

R_{SH}^{l}\left( \alpha ,\beta ,\gamma \right)={{Z}_{\gamma }}{{Y}_{-90}}{{Z}_{\beta }}{{Y}_{+90}}{{Z}_{\alpha }}\\

对于第0次的球谐变换矩阵为

R_{SH}^{0}\left( \alpha ,\beta ,\gamma \right)=\left( 1 \right) \tag{30}\\

其它维度的矩阵推导比较麻烦，就推导了第1次的球谐变换矩阵，它可以表示为

R_{SH}^{1}\left( \alpha ,\beta ,\gamma \right)=\left( \begin{matrix} {{c}_{\alpha }}{{c}_{\gamma }}-{{s}_{\alpha }}{{c}_{\beta }}{{s}_{\gamma }} & -{{s}_{\beta }}{{s}_{\gamma }} & -{{s}_{\alpha }}{{c}_{\gamma }}-{{c}_{\alpha }}{{c}_{\beta }}{{s}_{\gamma }} \\ -{{s}_{\alpha }}{{s}_{\beta }} & {{c}_{\beta }} & -{{c}_{\alpha }}{{s}_{\beta }} \\ {{c}_{\alpha }}{{s}_{\gamma }}+{{s}_{\alpha }}{{c}_{\beta }}{{c}_{\gamma }} & {{s}_{\beta }}{{c}_{\gamma }} & {{c}_{\alpha }}{{c}_{\beta }}{{c}_{\gamma }}-{{s}_{\alpha }}{{s}_{\gamma }} \\ \end{matrix} \right)=\left( \begin{matrix} {{R}_{1,1}} & -{{R}_{1,2}} & {{R}_{1,0}} \\ -{{R}_{2,1}} & {{R}_{2,2}} & {{R}_{2,0}} \\ {{R}_{0,1}} & -{{R}_{0,2}} & {{R}_{0,0}} \\ \end{matrix} \right)\\

变换矩阵R_{SH}^{l}的计算可以参见附录D3D的实现，实现了前6次的球谐系数的旋转，对于图形渲染来说，已经够用了。

对于高维矩阵的构造方法非常的复杂，采用的是魏格纳d矩阵（Wigner d-matrices），可以参见文献[4]的讨论，网上也有这个算法的高效实现，有兴趣可以研究研究，参见SHTns。

附录

根据等式（15）的递归关系，就可以很容易计算出连带勒让德函数[3]。

double P(int l,int m,double x) { // evaluate an Associated Legendre Polynomial P(l,m,x) at x double pmm = 1.0; if(m>0) { double somx2 = sqrt((1.0-x)*(1.0+x)); double fact = 1.0; for(int i=1; i<=m; i++) { pmm *= (-fact) * somx2; fact += 2.0; } } if(l==m) return pmm; double pmmp1 = x * (2.0*m+1.0) * pmm; if(l==m+1) return pmmp1; double pll = 0.0; for(int ll=m+2; ll<=l; ++ll) { pll = ( (2.0*ll-1.0)*x*pmmp1-(ll+m-1.0)*pmm ) / (ll-m); pmm = pmmp1; pmmp1 = pll; } return pll; }

根据等式（20）（21），可以计算出球谐函数[3]。

double K(int l, int m) { // renormalisation constant for SH function double temp = ((2.0*l+1.0)*factorial(l-m)) / (4.0*PI*factorial(l+m)); return sqrt(temp); } double SH(int l, int m, double theta, double phi) { // return a point sample of a Spherical Harmonic basis function // l is the band, range [0..N] // m in the range [-l..l] // theta in the range [0..Pi] // phi in the range [0..2*Pi] const double sqrt2 = sqrt(2.0); if(m==0) return K(l,0)*P(l,m,cos(theta)); else if(m>0) return sqrt2*K(l,m)*cos(m*phi)*P(l,m,cos(theta)); else return sqrt2*K(l,-m)*sin(-m*phi)*P(l,-m,cos(theta)); }

在D3D中实现的球谐系数的旋转D3DXSHRotate的实现为：

FLOAT* WINAPI D3DXSHRotate(FLOAT *out, UINT order, const D3DXMATRIX *matrix, const FLOAT *in) { FLOAT alpha, beta, gamma, sinb, temp[36], temp1[36]; TRACE("out %p, order %u, matrix %p, in %p\n", out, order, matrix, in); out[0] = in[0]; if ((order > D3DXSH_MAXORDER) || (order < D3DXSH_MINORDER)) return out; if (order <= 3) { out[1] = matrix->u.m[1][1] * in[1] - matrix->u.m[2][1] * in[2] + matrix->u.m[0][1] * in[3]; out[2] = -matrix->u.m[1][2] * in[1] + matrix->u.m[2][2] * in[2] - matrix->u.m[0][2] * in[3]; out[3] = matrix->u.m[1][0] * in[1] - matrix->u.m[2][0] * in[2] + matrix->u.m[0][0] * in[3]; if (order == 3) { FLOAT coeff[]={ matrix->u.m[1][0] * matrix->u.m[0][0], matrix->u.m[1][1] * matrix->u.m[0][1], matrix->u.m[1][1] * matrix->u.m[2][1], matrix->u.m[1][0] * matrix->u.m[2][0], matrix->u.m[2][0] * matrix->u.m[2][0], matrix->u.m[2][1] * matrix->u.m[2][1], matrix->u.m[0][0] * matrix->u.m[2][0], matrix->u.m[0][1] * matrix->u.m[2][1], matrix->u.m[0][1] * matrix->u.m[0][1], matrix->u.m[1][0] * matrix->u.m[1][0], matrix->u.m[1][1] * matrix->u.m[1][1], matrix->u.m[0][0] * matrix->u.m[0][0], }; out[4] = (matrix->u.m[1][1] * matrix->u.m[0][0] + matrix->u.m[0][1] * matrix->u.m[1][0]) * in[4]; out[4] -= (matrix->u.m[1][0] * matrix->u.m[2][1] + matrix->u.m[1][1] * matrix->u.m[2][0]) * in[5]; out[4] += 1.7320508076f * matrix->u.m[2][0] * matrix->u.m[2][1] * in[6]; out[4] -= (matrix->u.m[0][1] * matrix->u.m[2][0] + matrix->u.m[0][0] * matrix->u.m[2][1]) * in[7]; out[4] += (matrix->u.m[0][0] * matrix->u.m[0][1] - matrix->u.m[1][0] * matrix->u.m[1][1]) * in[8]; out[5] = (matrix->u.m[1][1] * matrix->u.m[2][2] + matrix->u.m[1][2] * matrix->u.m[2][1]) * in[5]; out[5] -= (matrix->u.m[1][1] * matrix->u.m[0][2] + matrix->u.m[1][2] * matrix->u.m[0][1]) * in[4]; out[5] -= 1.7320508076f * matrix->u.m[2][2] * matrix->u.m[2][1] * in[6]; out[5] += (matrix->u.m[0][2] * matrix->u.m[2][1] + matrix->u.m[0][1] * matrix->u.m[2][2]) * in[7]; out[5] -= (matrix->u.m[0][1] * matrix->u.m[0][2] - matrix->u.m[1][1] * matrix->u.m[1][2]) * in[8]; out[6] = (matrix->u.m[2][2] * matrix->u.m[2][2] - 0.5f * (coeff[4] + coeff[5])) * in[6]; out[6] -= (0.5773502692f * (coeff[0] + coeff[1]) - 1.1547005384f * matrix->u.m[1][2] * matrix->u.m[0][2]) * in[4]; out[6] += (0.5773502692f * (coeff[2] + coeff[3]) - 1.1547005384f * matrix->u.m[1][2] * matrix->u.m[2][2]) * in[5]; out[6] += (0.5773502692f * (coeff[6] + coeff[7]) - 1.1547005384f * matrix->u.m[0][2] * matrix->u.m[2][2]) * in[7]; out[6] += (0.2886751347f * (coeff[9] - coeff[8] + coeff[10] - coeff[11]) - 0.5773502692f * (matrix->u.m[1][2] * matrix->u.m[1][2] - matrix->u.m[0][2] * matrix->u.m[0][2])) * in[8]; out[7] = (matrix->u.m[0][0] * matrix->u.m[2][2] + matrix->u.m[0][2] * matrix->u.m[2][0]) * in[7]; out[7] -= (matrix->u.m[1][0] * matrix->u.m[0][2] + matrix->u.m[1][2] * matrix->u.m[0][0]) * in[4]; out[7] += (matrix->u.m[1][0] * matrix->u.m[2][2] + matrix->u.m[1][2] * matrix->u.m[2][0]) * in[5]; out[7] -= 1.7320508076f * matrix->u.m[2][2] * matrix->u.m[2][0] * in[6]; out[7] -= (matrix->u.m[0][0] * matrix->u.m[0][2] - matrix->u.m[1][0] * matrix->u.m[1][2]) * in[8]; out[8] = 0.5f * (coeff[11] - coeff[8] - coeff[9] + coeff[10]) * in[8]; out[8] += (coeff[0] - coeff[1]) * in[4]; out[8] += (coeff[2] - coeff[3]) * in[5]; out[8] += 0.86602540f * (coeff[4] - coeff[5]) * in[6]; out[8] += (coeff[7] - coeff[6]) * in[7]; } return out; } if (fabsf(matrix->u.m[2][2]) != 1.0f) { sinb = sqrtf(1.0f - matrix->u.m[2][2] * matrix->u.m[2][2]); alpha = atan2f(matrix->u.m[2][1] / sinb, matrix->u.m[2][0] / sinb); beta = atan2f(sinb, matrix->u.m[2][2]); gamma = atan2f(matrix->u.m[1][2] / sinb, -matrix->u.m[0][2] / sinb); } else { alpha = atan2f(matrix->u.m[0][1], matrix->u.m[0][0]); beta = 0.0f; gamma = 0.0f; } D3DXSHRotateZ(temp, order, gamma, in); rotate_X(temp1, order, 1.0f, temp); D3DXSHRotateZ(temp, order, beta, temp1); rotate_X(temp1, order, -1.0f, temp); D3DXSHRotateZ(out, order, alpha, temp1); return out; } static void rotate_X(FLOAT *out, UINT order, FLOAT a, FLOAT *in) { out[0] = in[0]; out[1] = a * in[2]; out[2] = -a * in[1]; out[3] = in[3]; out[4] = a * in[7]; out[5] = -in[5]; out[6] = -0.5f * in[6] - 0.8660253882f * in[8]; out[7] = -a * in[4]; out[8] = -0.8660253882f * in[6] + 0.5f * in[8]; out[9] = -a * 0.7905694842f * in[12] + a * 0.6123724580f * in[14]; out[10] = -in[10]; out[11] = -a * 0.6123724580f * in[12] - a * 0.7905694842f * in[14]; out[12] = a * 0.7905694842f * in[9] + a * 0.6123724580f * in[11]; out[13] = -0.25f * in[13] - 0.9682458639f * in[15]; out[14] = -a * 0.6123724580f * in[9] + a * 0.7905694842f * in[11]; out[15] = -0.9682458639f * in[13] + 0.25f * in[15]; if (order == 4) return; out[16] = -a * 0.9354143739f * in[21] + a * 0.3535533845f * in[23]; out[17] = -0.75f * in[17] + 0.6614378095f * in[19]; out[18] = -a * 0.3535533845f * in[21] - a * 0.9354143739f * in[23]; out[19] = 0.6614378095f * in[17] + 0.75f * in[19]; out[20] = 0.375f * in[20] + 0.5590170026f * in[22] + 0.7395099998f * in[24]; out[21] = a * 0.9354143739f * in[16] + a * 0.3535533845f * in[18]; out[22] = 0.5590170026f * in[20] + 0.5f * in[22] - 0.6614378691f * in[24]; out[23] = -a * 0.3535533845f * in[16] + a * 0.9354143739f * in[18]; out[24] = 0.7395099998f * in[20] - 0.6614378691f * in[22] + 0.125f * in[24]; if (order == 5) return; out[25] = a * 0.7015607357f * in[30] - a * 0.6846531630f * in[32] + a * 0.1976423711f * in[34]; out[26] = -0.5f * in[26] + 0.8660253882f * in[28]; out[27] = a * 0.5229125023f * in[30] + a * 0.3061861992f * in[32] - a * 0.7954951525f * in[34]; out[28] = 0.8660253882f * in[26] + 0.5f * in[28]; out[29] = a * 0.4841229022f * in[30] + a * 0.6614378691f * in[32] + a * 0.5728219748f * in[34]; out[30] = -a * 0.7015607357f * in[25] - a * 0.5229125023f * in[27] - a * 0.4841229022f * in[29]; out[31] = 0.125f * in[31] + 0.4050463140f * in[33] + 0.9057110548f * in[35]; out[32] = a * 0.6846531630f * in[25] - a * 0.3061861992f * in[27] - a * 0.6614378691f * in[29]; out[33] = 0.4050463140f * in[31] + 0.8125f * in[33] - 0.4192627370f * in[35]; out[34] = -a * 0.1976423711f * in[25] + a * 0.7954951525f * in[27] - a * 0.5728219748f * in[29]; out[35] = 0.9057110548f * in[31] - 0.4192627370f * in[33] + 0.0624999329f * in[35]; } FLOAT * WINAPI D3DXSHRotateZ(FLOAT *out, UINT order, FLOAT angle, const FLOAT *in) { UINT i, sum = 0; FLOAT c[5], s[5]; TRACE("out %p, order %u, angle %f, in %p\n", out, order, angle, in); order = min(max(order, D3DXSH_MINORDER), D3DXSH_MAXORDER); out[0] = in[0]; for (i = 1; i < order; i++) { UINT j; c[i - 1] = cosf(i * angle); s[i - 1] = sinf(i * angle); sum += i * 2; out[sum - i] = c[i - 1] * in[sum - i]; out[sum - i] += s[i - 1] * in[sum + i]; for (j = i - 1; j > 0; j--) { out[sum - j] = 0.0f; out[sum - j] = c[j - 1] * in[sum - j]; out[sum - j] += s[j - 1] * in[sum + j]; } if (in == out) out[sum] = 0.0f; else out[sum] = in[sum]; for (j = 1; j < i; j++) { out[sum + j] = 0.0f; out[sum + j] = -s[j - 1] * in[sum - j]; out[sum + j] += c[j - 1] * in[sum + j]; } out[sum + i] = -s[i - 1] * in[sum - i]; out[sum + i] += c[i - 1] * in[sum + i]; } return out; }

参考

[1] 姚端正, 梁家宝. 数学物理方法-第4版..

[2] 顾樵. 数学物理方法.

[3] Robin Green. "Spherical harmonic lighting: The gritty details." Archives of the Game Developers Conference. Vol. 56. 2003.

[4] Joseph Ivanic, and Klaus Ruedenberg. "Rotation matrices for real spherical harmonics. Direct determination by recursion." The Journal of Physical Chemistry 100.15, 6342-6347, 1996.

[5] Wikipea. Euler angles

[6] Peter-Pike Sloan. "Stupid spherical harmonics (sh) tricks." Game developers conference. Vol. 9. 2008.

[7] Ravi Ramamoorthi, and Pat Hanrahan. "An efficient representation for irradiance environment maps." Proceedings of the 28th annual conference on Computer graphics and interactive techniques. 2001.