Rendering-球谐光照推导及应用

最新推荐文章于 2023-12-21 18:16:32 发布

繁弱

最新推荐文章于 2023-12-21 18:16:32 发布

阅读量2.1k

点赞数 3

分类专栏： UnrealEngine4 Rendering 文章标签： 3d渲染 rendering 球谐 Math unreal

本文链接：https://blog.csdn.net/zzgyy123/article/details/113761598

版权

UnrealEngine4 同时被 2 个专栏收录

27 篇文章 8 订阅

订阅专栏

Rendering

7 篇文章 0 订阅

订阅专栏

这篇草稿在草稿堆里躺了快一年了，去年过年就开始了，一年了中间再没动静，趁今年过年有时间赶紧翻出来捋捋

一直觉得概念或者想法类的东西可以知道一下就可以，但如果背后有其数学表达的，总觉得需要自己手把手推导一遍才觉得踏实，不过说实话，知其然又要知其所以然，还是蛮花时间的，年假到现在，基本每天都在看，饶是如此，还是有些坑没有填到，吾生也有涯，而知也无涯。以有涯随无涯，殆已。

推导

球谐函数是球面坐标系下的拉普拉斯方程(f(x,y,z) = 0 的泊松方程)中角度部分的解，这里参考球谐函数自己手推一遍

先看普通二阶偏微分方程形式的拉普拉斯方程

$\triangledown ^2f = \frac{\partial ^2f}{\partial x^2 } + \frac{\partial ^2f}{\partial y^2 } + \frac{\partial ^2f}{\partial z^2 }$

球面坐标系下的拉普拉斯方程的形式

$\triangledown ^2f = \frac{1}{r^2} \frac{\partial}{\partial r}(r^2\frac{\partial f}{\partial r}) + \frac{1}{r^2sin\theta } \frac{\partial}{\partial \theta }(sin\theta\frac{\partial f}{\partial \theta}) + \frac{1}{r^2sin^2\theta } \frac{\partial ^2f}{\partial\varphi ^2 }$

其中 $\triangledown ^2$ 表示拉普拉斯算子，φ 表示方位角 Azimuth Angle，θ 表示天顶角 Zenith Angle

从笛卡尔坐标系到球坐标系的推导过程，可以参考：https://wenku.baidu.com/view/f6c75dd2da38376baf1faec3.html

然后主要是分离变量法，考虑

$f(r,\theta ,\varphi) = R(r)Y(\theta ,\varphi )= R(r)\Theta (\theta )\Phi (\varphi )$

其中 $Y(\theta ,\varphi )$ 就是我们想要的角度部分的解，球谐函数

然后代入拉普拉斯方程得到

${\Theta \Phi \over r^{2}}{d \over dr}\left(r^{2}{dR \over dr}\right)+{R\Phi \over r^{2}\sin \theta }{d \over d\theta }\left(\sin \theta {d\Theta \over d\theta }\right)+{R\Theta \over r^{2}\sin ^{2}\theta }{d^{2}\Phi \over d\varphi ^{2}}=0$

乘 $r^2$ ，约去得到

${\Theta \Phi}{d \over dr}\left(r^{2}{dR \over dr}\right)+{R\Phi \over \sin \theta }{d \over d\theta }\left(\sin \theta {d\Theta \over d\theta }\right) + {R\Theta \over \sin ^{2}\theta }{d^{2}\Phi \over d\varphi ^{2}} = 0$ .....................(1)

分离变量 r：

${\Theta \Phi}{d \over dr}\left(r^{2}{dR \over dr}\right) = - \frac{1}{R}\left ( {\Phi \over \sin \theta }{d \over d\theta }\left(\sin \theta {d\Theta \over d\theta }\right)+{R\Theta \over \sin ^{2}\theta }{d^{2}\Phi \over d\varphi ^{2}}\right )$

移项得到

$\frac{1}{R}{d \over dr}\left(r^{2}{dR \over dr}\right) = - \frac{1}{{\Theta \Phi}}\left ({\Phi \over \sin \theta }{d \over d\theta }\left(\sin \theta {d\Theta \over d\theta }\right)+{R\Theta \over \sin ^{2}\theta }{d^{2}\Phi \over d\varphi ^{2}}\right )= \lambda$ .....................(2)

注意前面的 ΘΦ 就是我们要求的球谐部分 Y。

分离变量 φ：

${R\Theta \over \sin ^{2}\theta }{d^{2}\Phi \over d\varphi ^{2}} = -\Phi\left ( {\Theta }{d \over dr}\left(r^{2}{dR \over dr}\right)+{R \over \sin \theta }{d \over d\theta }\left(\sin \theta {d\Theta \over d\theta }\right)\right )$

同样移项得到

${1 \over \Phi }{d^{2}\Phi \over d\varphi ^{2}} = -\left ( { \sin ^{2} \theta \over R }{d \over dr}\left(r^{2}{dR \over dr}\right)+{\sin \theta \over \Theta}{d \over d\theta }\left(\sin \theta {d\Theta \over d\theta }\right)\right ) = -m^2$ ....................(3)

分离变量 θ：

代入(2),(3)到(1)式，得到只剩 θ 的表达式

${\Theta \Phi R}\lambda +{R\Phi \over \sin \theta }{d \over d\theta }\left(\sin \theta {d\Theta \over d\theta }\right) - { {\Theta \Phi R} m^2 \over \sin ^{2}\theta } = 0$

约去 $f = {\Theta \Phi R}$ ，化简得到

$\lambda +{1 \over \Theta \sin \theta }{d \over d\theta }\left(\sin \theta {d\Theta \over d\theta }\right) - { m^2 \over \sin ^{2}\theta } = 0$ ....................(4)

各自转换形式得到

(2)式

${d \over dr}\left(r^{2}{dR \over dr}\right) = \lambda R \Rightarrow \mathbf{r^{2}R''+2rR'-\lambda R=0}$ ....................(5)

(3)式

${1 \over \Phi }{d^{2}\Phi \over d\varphi ^{2}} = -m^2 \Rightarrow \mathbf{\Phi ''+m^{2}\Phi =0}$ ....................(6)

(4)式，左右同乘 $\Theta \sin^2 \theta$

$\lambda \Theta \sin ^2 \theta +{\sin \theta }{d \over d\theta }\left(\sin \theta {d\Theta \over d\theta }\right) - { m^2 \Theta } = 0 \Rightarrow \mathbf{\sin \theta {\dfrac {d}{d\theta }}(\sin \theta \Theta ')+(\lambda \sin ^{2}\theta -m^{2})\Theta =0}$ ....................(7)

我们关心的球谐函数是方程的角度部分，因此我们需要求的是 $\Theta (\theta )\Phi (\varphi )$ 这两个函数的解，也就是(6)式和(7)式

求解 Θ，已知(6)式，且 φ 是球形坐标里的方位角，是一个2π 的周期函数，因此 m 必须是整数

$\Phi ''= -m^{2}\Phi \quad and \quad \Phi (\phi) = \Phi(\phi + 2\Pi )$ ，可解得虚数域上的解

$\Phi =e^{{im\phi }}, m \in \mathbb {Z}$ ....................(8)

其中i是虚数单位，之后计算球谐的话会换算到实数域上的解(转换的推导见下方注释)

求解 Φ ，对(7)式变量替换， $t= cos \theta, dt = -sin\theta d \theta, |t|\leq 1$ .

我们关注的 θ 是球形坐标系里的天顶角，取值范围[0,π]，所以cosθ 范围是 [-1,1], t 的绝对值小于等于1，也就是我们的目标函数是有界的

$-sin^2 \theta {\dfrac {d}{dt}}(-\sin^2 \theta {d\Theta \over dt})+(\lambda \sin ^{2}\theta -m^{2})\Theta =0$

$(t^2 -1 ) {\dfrac {d}{dt}}((t^2 -1) {d\Theta \over dt})+(\lambda( t ^{2} - 1) -m^{2})\Theta =0$

除 $1-t^2$

$-{\dfrac {d}{dt}}((t^2 -1) {d\Theta \over dt})+(\lambda - {m^{2} \over ( t ^{2} - 1)})\Theta =0$

整理得到

$(1-t^2) {d^2\Theta \over dt^2}-2t {d\Theta \over dt}+(\lambda - {m^{2} \over ( t ^{2} - 1)})\Theta =0$

再根据上面的(2)式，令

$\frac{1}{R}{d \over dr}\left(r^{2}{dR \over dr}\right) = \lambda = \ell(\ell +1)$

参考一下伴随勒让德多项式的标准形式

$(1-x^2)\,\frac{d^2\,y}{dx^2} -2x\frac{dy}{dx} + \left(\ell[\ell+1] - \frac{m^2}{1-x^2}\right)\,y = 0$

正好 Φ 项就是伴随勒让德多项式的形式，这也是为什么很多地方一提到这个球谐函数就会提到这个多项式的原因

而且从上面链接能看到，这个方程在方程仅当 $\ell$ 和 m 均为整数且满足 $0\leq m\leq \ell$ 时，才在区间 [−1, 1] 上有非奇异解，参考wiki，此时我们能得到多项式的解

$\Theta =P_{\ell }^{m}(t) =P_{\ell }^{m}(cos \theta),l\in {\mathbb {N}},l\geqslant |m| \geqslant 0$ ，l 表示“degree”，m表示“order”

而求解伴随勒让德多项式，一种方法是可以根据它的形式，对相应的勒让德多项式求m 次导得来（下面的自变量 x，即上面的t，也即 cos θ）

$P_l^m(x)=(1-x^2)^{m/2}P_l^{(m)}(x)$

对应的勒让德多项式的形式是这样

$P_{\ell }(x)={1 \over 2^{\ell }\ell !}{d^{\ell } \over dx^{\ell }}(x^{2}-1)^{l}$

但可以看到，直接这样求解是比较复杂的，因此我们一般会基于递推公式来求解伴随勒让德多项式，它有如下几个性质

$\\P^m_m=(-1)^m(2m-1)!!(1-x^2)^{m/2} \\ P^m_{m+1}=x(2m+1)P^m_m \\ (l-m)P^m_l=x(2l-1)P^m_{l-1}-(l+m-1)P^m_{l-2} \\$

其中 !! 不是普通的阶乘，而是表示双阶乘，正整数的双阶乘表示不超过这个正整数且与它有相同奇偶性的所有正整数的乘积。

比如

由第一个公式可以求首项，得到 $P^0_0 = 1$ , $P^1_1 = -(1-x^2)^{1/2}$ , $P^2_2 = 3(1-x^2)$ ...

由第二个公式可以在每个degree “l” 推高一个 order "m"，得到 $P^0_1 = x$ , $P^1_2 =x*3*P^1_1 = -3x (1-x^2)^{1/2}$ , $P^2_3 =x*5*P^2_2 = 15x (1-x^2)$ ...

由第三个公式，可以根据前两项，计算得到第三项，开始递归补足前面两步中间留下的，

$(2-0)P^0_2 = x(2*2-1)P^0_1 - (2+0-1)*P^0_0 = 3x^2 -1$ ，得到 $P^0_2 = \frac{3x^2 -1}{2}$

$(3-0)P^0_3 = x(2*3-1)P^0_2 - (3+0-1)*P^0_1 = 5x\frac{3x^2 -1}{2}-2x$ ，得到 $P^0_3 = 5x\frac{3x^2 -1}{6}-\frac{2x}{3} = \frac{5x^3-3x}{2}$ ...

这里简单的表达一下递归流程，更多的计算结果可以参考 https://en.wikipedia.org/wiki/Associated_Legendre_polynomials

现在两项都有了，整合起来就是最终的球谐函数的表达式

$Y_{\ell }^{m}(\theta ,\varphi )=A_{l}^m\Phi (\varphi )\Theta (\theta )=A_{l}^m\,e^{{im\varphi }}\,P_{\ell }^{m}(\cos {\theta })$ ，其中 $A_{l}^m$ 是归一化因子，用于让 Y 的整个球面积分等于 1，它的推导可以参考下面的 Ref 8，经归一化后表示为

$Y_{\ell }^{m}(\theta ,\ \varphi )=(-1)^{m}{\sqrt {{(2\ell +1) \over 4\pi }{(\ell -|m|)! \over (\ell +|m|)!}}}\,P_{\ell }^{m}(\cos {\theta })\,e^{im\varphi }$

其中， $(-1)^m$ 是 Condon–Shortley Phase，用于物理学里的计算，防止一个相位计算两次，这里计算球谐可以忽略

另外后面的Θ项，可以从虚数域转到实数域，参考Spherical_harmonics，表示为(其中sin,cos 的推导见注释)

$Y^m_{ \ell}(\theta,\varphi)= \left\{\begin{matrix} {\sqrt {{2\ell +1 \over 4\pi }{(\ell -|m|)! \over (\ell +|m|)!}}}\ P_{\ell }^{|m|}(\cos \theta )\ \sin(|m|\varphi ) & {\mbox{if }}m<0 \\ {\sqrt {2\ell +1 \over 4\pi }}\ P_{\ell }^{0}(\cos \theta )&{\mbox{if }}m=0 \\{\sqrt {{2\ell +1 \over 4\pi }{(\ell -m)! \over (\ell +m)!}}}\ P_{\ell }^{m}(\cos \theta )\ \cos(m\varphi )&{\mbox{if }}m>0 \end{matrix}\right.$

注释：

有些地方能看到下面这种实数域的球谐公式表示

$Y^m_{ \ell}(\theta,\varphi)= \left\{\begin{matrix} (-1)^{m}{\sqrt {2}}{\sqrt {{2\ell +1 \over 4\pi }{(\ell -|m|)! \over (\ell +|m|)!}}}\ P_{\ell }^{|m|}(\cos \theta )\ \sin(|m|\varphi ) & {\mbox{if }}m<0 \\ {\sqrt {2\ell +1 \over 4\pi }}\ P_{\ell }^{0}(\cos \theta )&{\mbox{if }}m=0 \\(-1)^{m}{\sqrt {2}}{\sqrt {{2\ell +1 \over 4\pi }{(\ell -m)! \over (\ell +m)!}}}\ P_{\ell }^{m}(\cos \theta )\ \cos(m\varphi )&{\mbox{if }}m>0 \end{matrix}\right.$

$(-1)^m$ 解释过了，可以不要， $\sqrt2$ 的来源其实是球谐公式从虚数域到实数域转化时进行归一化引入的

虚数域的公式有如下推导，参考 https://math.stackexchange.com/questions/1095011/in-the-real-spherical-harmonics-where-does-the-sqrt2-factor-come-from

$Y_{l}^m = \left\{\begin{matrix} \frac{1}{i \sqrt{2}}( Y_l^{-m} - (-1)^mY_l^m) & \text{if } m < 0 \\ Y_l^m & \text{if } m = 0 \\\frac{1}{\sqrt{2}} ( Y_l^m + (-1)^mY_l^{-m} ) & \text{if } m > 0 \end{matrix}\right.$

以m>0的情况为例，将虚数域转到实数域，需要用到欧拉公式，由于

$\\e^{ix} = cosx+isinx \Rightarrow \\\\cosmx = \frac{e^{imx}+e^{-imx}}{2} \\sinmx = \frac{e^{imx}-e^{-imx}}{2}$

所以有

$\\Y_{l}^m = \frac{1}{\sqrt{2}} ( Y_l^m + (-1)^mY_l^{-m} ) \\ = A_l^m P_l^m\frac{1}{\sqrt{2}} (e^{imx}+e^{-imx})\Rightarrow (Euler Formula) \\ = A_l^m P_l^m\frac{1}{\sqrt{2}} (cosmx+isin(mx)+cos(-mx)+isin(-mx)) \\= A_l^m P_l^m\frac{1}{\sqrt{2}} *2*cosmx \\= A_l^m P_l^m{\sqrt{2}} *cosmx$

这就得到了上面实数域公式里的 cos(mx) 以及 $\sqrt 2$ ，对于的，也能得到 m<0 时的 sin(mx)

到此，得到了球谐函数的基本形式，

代入上面的伴随勒让德多项式的各阶公式，就可以得到最终的球谐函数的各阶公式(另外代入 $x=rsin\theta cos\varphi \quad y=rsin\theta sin\varphi \quad z=r cos\theta$ ，可转到笛卡尔坐标系，用于实际应用，因为一般传入的也是一个 normal vector）

$Y_{0 }^{0}(\theta ,\ \varphi )={\sqrt {{1\over 4\pi }}}P^{0}_{0}(\cos {\theta })cos(0) = \sqrt{\frac{1}{4\pi}}$

$Y_{1}^{0}(\theta ,\ \varphi )={\sqrt {{3\over 4\pi }}}P^{0}_{1}(cos\theta)= \sqrt{\frac{3}{4\pi}}cos\theta = \frac{1}{2}\sqrt{\frac{3}{\pi}}\frac{z}{r}$

$Y_{1}^{1}(\theta ,\ \varphi )={\sqrt {{3\over4 \pi }}}P^{1}_{1}(cos\theta) cos\varphi= \sqrt{\frac{3}{4\pi}}sin\theta cos\varphi= \frac{1}{2}\sqrt{\frac{3}{\pi}}\frac{x}{r}$ ...

更多的可以在Table_of_spherical_harmonics 看到

把公式里前面的数字单拎出来

$\\\sqrt{\frac{1}{4\pi}} = 0.28209,\quad \frac{1}{2}\sqrt{\frac{3}{\pi}} = 0.488602,\quad \frac{1}{4}\sqrt{\frac{5}{\pi}} = 0.315391, \\\quad \frac{1}{2}\sqrt{\frac{15}{\pi}} = 1.09254,\quad \frac{1}{4}\sqrt{\frac{15}{\pi}} = 0.54627$

这些数字，经常会被引擎作为魔法数字硬编码到代码里，第一眼看到，可能会不明所以，其实它们就是来自这里

比如 UE 里的 SHBasisFunction3，能看到它基本和上面直角坐标系中的表达式一一对应，传进去normal，然后返回出来三阶球谐

FThreeBandSHVector SHBasisFunction3(half3 InputVector)
{
	FThreeBandSHVector Result;
	// These are derived from simplifying SHBasisFunction in C++
	Result.V0.x = 0.282095f; 
	Result.V0.y = -0.488603f * InputVector.y;
	Result.V0.z = 0.488603f * InputVector.z;
	Result.V0.w = -0.488603f * InputVector.x;

	half3 VectorSquared = InputVector * InputVector;
	Result.V1.x = 1.092548f * InputVector.x * InputVector.y;
	Result.V1.y = -1.092548f * InputVector.y * InputVector.z;
	Result.V1.z = 0.315392f * (3.0f * VectorSquared.z - 1.0f);
	Result.V1.w = -1.092548f * InputVector.x * InputVector.z;
	Result.V2 = 0.546274f * (VectorSquared.x - VectorSquared.y);

	return Result;
}

性质

开始应用之前，首先需要了解一下球谐函数的几个很重要的性质

1.正交完备性

正交是指，球谐函数的所有基任意两个都是正交的，它满足

$\int _{\theta =0}^{\pi }\int _{\varphi =0}^{2\pi }Y_{\ell }^{m}\,Y_{\ell '}^{m'}{}^{*}\,d\Omega =\delta _{\ell \ell '}\,\delta _{mm'}$

$\delta _{ij}$ 是克罗内克δ函数，它表示如果 i=j，函数就 =1，否则 =0，

dΩ 是球面积分， dΩ = sin θ dφ dθ

完备是说由于基是无限的，理论上它可以表达所有的球面函数 F(x)

2.旋转不变性

哇，这又是一个大话题，不过这个已经有很多的文章都解释了，所以具体的推导这里先不提，要是再捋这个，我怕这篇草稿又得搁这躺一年...咱有缘再会.....

这里先就简单的直接使用一下它的结论，它的旋转不变不是说怎么转都不变，而是说假定它的旋转变换是 r(x)，我可以先计算旋转再计算球谐 f(r(x))，也可以先计算球谐之后再旋转，也就是 f(r(x)) = r(f(x))

这个的应用一个在于本身采样球是可以旋转的，旋转之后，不用重复计算球谐系数，另外一个在于可以进行空间变换，比如local space 到 world space 的旋转变换，比如有些应用会把球谐系数记录在模型的顶点上，由于是预处理的，一般是local space的，运行时需要转换到 world space，而由于球谐是旋转不变的，因此光照结果不会出错

3.低频高频的问题

经常看到是说球谐适合模拟低频部分的光照信息，比如间接光照，不适合表达高光或者反射等高频的光照信息，对于低频高频，我个人的理解是如下图所示(图像来源 Associated_Legendre_polynomials)

图像来源 wiki

像图中，所示 m=0 的伴随勒得让多项式的不同的 l 阶的曲线，阶数越高，频率越高，变化越高，就越能表达更高频的信息(频率高的基底，也只能表达高频信息)，而由于一般实际应用，最多应用3到6阶的球谐，并不足以表达高频信息，所以一般会说球谐适合表达低频信息，高频的可以考虑用小波变换来表示

这里也可以搬出这张经典的图(图像来源 Spherical_harmonics)，浅色表示 $Y_l^m$ 的正值，深色表示负值，表面到原点的距离表示 $Y_l^m$ 在角度 (θ ,φ)的绝对值大小，整个球面函数就是它们乘各自系数叠加起来，从中也能看到越高阶对越高频信息的贡献区域

4.光照函数和传输函数的积分可以转变成它们之间球谐系数的点积，将积分变点乘，大幅简化了计算

传输函数 Tranfer Function 这个概念经常在 PRT(Precompute Radiance Transfer) 提到，我个人的理解是它是对入射照度 Light(x) 进行变换操作的函数，不管是 Lambert 里面的一个 cos，还是标准 BRDF公式里的 FDG，它们都可以变成传输函数，只是说有些传输函数本身需要的是高频信息，而一般的球谐函数只用有限阶数，适合来表现一些低频信息，不太适合表现高频信息，所以有些函数不会去用球谐表示，而这个特性要做的是说，每种传输函数本来计算都是要组合起来然后积分的，现在都可以各自预计算各自的球谐系数，然后运行时可以把积分变成系数之间的点积操作，大幅简化了运行时的计算量。

总之有了传输函数，我们就可以表现更加复杂的光照情况，比如遮挡关系等等，而且不会带来很大计算量的开销

应用

了解了基本的性质，就可以开始应用的两步走了，先是投影，然后是重建

先是投影，一般是预计算的，

我们从最后到底想要的是什么出发，其实我们想要的是一个球面函数 $F(\theta,\phi)$ ，运行时只要传进去角度或者法线向量，就可以得到球面这点的光照情况，而为了得到这个函数，需要先得到它的函数表达

为了得到函数表达，就涉及到了投影以及傅里叶变换的概念

傅里叶变换的概念是说对于任意函数 F(x)，都可以把它展开成一组正交基底乘上各自系数的和，也就是使用这组正交基底，我可以表示这个函数域内的所有函数

而投影是说，因为正交完备，只有这个基底自己和自己相乘，结果才为1，能够保留，因此相当于相乘后就能得到这个基底对于原函数的贡献量(系数)，这个过程就是投影，所以用原函数F(x)往每个基底上投影，就可以得到每个基底对于这个函数的系数 c1,c2,c3...

之后，还原的时候，只要用得到的系数乘以基底，然后加起来就可以还原回来原函数 F(x)

所以，经投影，对于傅里叶变换，原函数可以表示为 $F(x) = \Sigma c*f(x)$ ，使用的基底越多，阶数越高，能够还原的越精确，理论上，只要取的够多，就能完全还原

球谐函数同理，只不过它是拓展之后的二重广义傅里叶展开，

$\\F(\theta,\phi)=\sum^\infty_{l=0}\sum^l_{m=-l}C^m_lY^m_l(\theta,\phi)\\=C^0_0Y^0_0(\theta,\phi) + C^0_1Y^0_1(\theta,\phi) + C^{-1}_1Y^{-1}_1(\theta,\phi)+C^1_1Y^1_1(\theta,\phi)+......$

它的所有基 $Y_l^m$ 也都是正交的，也可以投影

$\\F(\theta,\phi)*Y^0_0(\theta,\phi)=C^0_0 \\F(\theta,\phi)*Y^0_1(\theta,\phi)=C^0_1 \\F(\theta,\phi)*Y^{-1}_1(\theta,\phi)=C^{-1}_1 \\F(\theta,\phi)*Y^1_1(\theta,\phi)=C^1_1$

因此我们可以用这些球谐系数去还原球面函数F(θ,φ)

现在有 $Y_l^m$ 这组正交基底，只要得到不同的系数 $C_l^m$ ，就可以表示出来所有的球面函数

换句话说，只要有每个位置预先算好的系数 $C_l^m$ ，运行时只要传进去角度，就能得到这点的光照了

于是现在问题变成了怎么在预计算的阶段，算出来这个球面函数的球谐系数 $C_l^m$ ，此时我们要求解的目标函数变为上述傅里叶展开后求解系数的公式

$C^m_l=\int _S F(\theta,\phi)Y^m_l(s)d\Omega = \int _{\theta =0}^{\pi }\int _{\varphi =0}^{2\pi }Y_{\ell }^{m}\, F(\theta,\phi)Y^m_l(s)sin\theta d\varphi d\theta$ （dΩ 是球面积分，dΩ = sin θ dφ dθ，这个公式的得来从公式表象就能看出来，它表示用 S 球面每个点的光照结果 F(θ,φ) 都往 $Y_l^m$ 投影一次，然后把球面无限个点的贡献量积分积起来就是 $C_l^m$ ）

要求解它，就需要用到蒙特卡洛积分了，把上述公式里求积分的过程转变成用概率密度函数求期望的过程

蒙特卡洛积分是说

假定 x 服从概率密度函数 p(x), 那函数 f(x) 的数学期望是

$E\left [f(x)\right ] = \int f(x)p(x)dx$

如果我们需要求

$I = \int f(x)dx$

就可以转换成求数学期望的形式，求 N 次 f(x)/p(x) 的均值结果作为 I 值结果

$I = \int f(x)dx = \int \frac{f(x)}{p(x)}p(x)dx = E\left [ \frac{f(x)}{p(x)}\right ] \approx \frac{1}{N}\sum^N_{i=1}\frac{f(x_i)}{p(x_i)}$

令 $w(x_i) = \frac{1}{p(x_i)}$ ，表示 $x_i$ 的权重，则最终可以表示为

$I = \int f(x)dx \approx \frac{1}{N}\sum^N_{i=1}\frac{f(x_i)}{p(x_i)} = \frac{1}{N}\sum^N_{i=1}f(x_i)w(x_i)$

对于球谐同理，不过我们需要的是球面的均匀采样，假定(0,1)的随机值 $\xi _x,\xi_y$

球面均匀采样公式是，推导可以参考下面 Ref6

$\left\{\begin{matrix} \theta = 2arccos\sqrt{1-\xi _x} \\\varphi = 2\pi\xi_y \end{matrix}\right.$

于是应用蒙特卡洛积分后，上述球谐系数的求解函数变成了

$C^m_l=\int _S F(\theta,\phi)Y^m_l(\theta,\phi)d\Omega =\frac{1}{N}\sum^N_{i=1}F(\theta,\phi)Y^m_l(\theta,\phi)w(\theta,\phi)$

其中，N 是采样个数，F(θ,φ) 是这点采样到的光照结果 light(θ,φ), $Y^m_l(\theta,\phi)$ 是球谐函数， $w(\theta,\phi)$ 是这点的概率权重

由于是球面采样，球面积 = 4ⲡ，所以概率是 1/4ⲡ，公式变为

$C^m_l=\frac{4\Pi }{N}\sum^N_{i=1}F(\theta,\phi)Y^m_l(\theta,\phi)$

因此离线时，根据质量标准和算力决定采样个数后，采样周围光照结果乘以球谐函数，最后累加起来，除以采样个数，就得到了这点的球谐系数。

此外，还可以考虑传输函数，假定传输函数 t(θ,φ)，现在light function F(θ,φ) 有了自己的球谐系数 $C_l^m$ ，传输函数 t(θ,φ) 也可以单独计算自己的球谐系数 $T_l^m$ ，运行时点积相乘就可以

这里也有一些细节待补充，可以结合 UE4的ILC或者 VLM 一起看看。

上述是离线预计算阶段的，下面就可以到实时阶段了

可以开始运行时重建，根据预计算好的球谐系数，只需要简单的计算，就可以还原原来的球面函数

$F(\theta,\phi)=\sum^\infty_{l=0}\sum^l_{m=-l}C^m_lY^m_l(\theta,\phi)$

运行时传进去θ,φ，或者法线向量，就可以得到这点的光照结果F(θ,φ)= light(θ,φ)。

文中的应用部分有一些并没亲手实现一下，如有错误，欢迎指正。

嗯......

那大概这样。

Reference：

https://zh.wikipedia.org/wiki/%E7%90%83%E8%B0%90%E5%87%BD%E6%95%B0
https://en.wikipedia.org/wiki/Spherical_harmonics
https://zh.wikipedia.org/wiki/%E6%8B%89%E6%99%AE%E6%8B%89%E6%96%AF%E6%96%B9%E7%A8%8B
https://zh.wikipedia.org/wiki/%E6%B3%8A%E6%9D%BE%E6%96%B9%E7%A8%8B
https://wenku.baidu.com/view/f6c75dd2da38376baf1faec3.html
https://www.zhihu.com/column/c_1046337272613527552
https://zhuanlan.zhihu.com/p/66989673
https://wuli.wiki/changed/SphHar.html
https://en.wikipedia.org/wiki/Kronecker_delta
https://math.stackexchange.com/questions/1095011/in-the-real-spherical-harmonics-where-does-the-sqrt2-factor-come-from
https://zh.wikipedia.org/zh-hans/%E6%AC%A7%E6%8B%89%E5%85%AC%E5%BC%8F
https://zh.wikipedia.org/wiki/%E5%85%8B%E7%BD%97%E5%86%85%E5%85%8B%CE%B4%E5%87%BD%E6%95%B0
https://zhuanlan.zhihu.com/p/153352797
https://en.wikipedia.org/wiki/Spherical_harmonics
Spherical Harmonic Lighting - the gritty details
Stupid Spherical Harmonics (SH) Tricks
An Efficient Representation for Irradiance Environment Maps

繁弱

关注

3
点赞
踩
8

收藏

觉得还不错? 一键收藏
0
评论
Rendering-球谐光照推导及应用

球谐函数是球面坐标系下的拉普拉斯方程(f(x,y,z) = 0 的泊松方程)中角度部分的解，这里参考球谐函数自己手推一遍先看球面坐标系下的拉普拉斯方程的形式其中表示拉普拉斯算子，表示方位角 Azimuth Angle，表示天顶角 Zenith AngleRef： https://zh.wikipedia.org/wiki/%E7%90%83%E8%B0%90%E5%87%BD%E6%95%B0 ...
复制链接

扫一扫