文章目录
使用球谐函数的本质
使用球谐函数的本质就是把一个点上的空间光照信息投影到球谐基向量上,然后再在需要使用的时候,使用投影得到的球谐系数来还原逼近原空间分布。
关于球谐光照,非常非常重要的信息就是,球谐函数描述的是球面分布,不是某个特定的方向,而是 整个球面所有方向所对应的分布,而不同阶数的球谐函数则是各种标准正交基底,即不同的球面分布。我们将空间光照信息投影到球谐的各个正交基上,得到球谐系数。我们使用的时候会带入一个方向 n n n,这个方向是整个球面空间上的某个方向,只要把该方向上对应的所有,球谐系数与球谐函数的乘积加起来,就是该方向上原分布的逼近
所以各个球谐系数其实是整个空间分布在球谐函数上的编码,所以我们在计算球谐系数的时候要对当前顶点的各个方向都进行积分
蒙特卡洛积分
根据大数定理
E
[
f
(
x
)
]
≈
1
N
∑
i
=
1
N
f
(
x
i
)
E[f(x)] \approx \frac{1}{N} \sum_{i=1}^{N} f\left(x_{i}\right)
E[f(x)]≈N1i=1∑Nf(xi)
而计算期望的形式也可以写成,即蒙特卡洛的形式
∫
f
(
x
)
d
x
=
∫
f
(
x
)
p
(
x
)
p
(
x
)
d
x
≈
1
N
∑
i
=
1
N
f
(
x
)
p
(
x
)
\int f(x) d x=\int \frac{f(x)}{p(x)} p(x) d x \approx \frac{1}{N} \sum_{i=1}^{N} \frac{f(x)}{p(x)}
∫f(x)dx=∫p(x)f(x)p(x)dx≈N1i=1∑Np(x)f(x)
我们根据PDF采样得到很多的样本点,最后将它们除以PDF之后的和平均,得到样本期望,如果PDF跟估计的分布相似,那么方差就会比较小。
如果PDF是均匀分布,那么计算就会非常简单,只需要按照均匀分布采样
如果是要均匀的采样一个球面,我们可以使用两个独立的标准随机变量
ξ
x
\xi_x
ξx 和
ξ
y
\xi_y
ξy,把它们这个面积对应的随机点映射到球坐标下,
(
2
arccos
(
1
−
ξ
x
)
,
2
π
ξ
y
)
→
(
θ
,
φ
)
\left(2 \arccos \left(\sqrt{1-\xi_{x}}\right), 2 \pi \xi_{y}\right) \rightarrow(\theta, \varphi)
(2arccos(1−ξx),2πξy)→(θ,φ)
球体的表面积是
4
π
4\pi
4π ,则PDF是
1
4
π
\frac{1}{4\pi}
4π1
采样方法使用 Inverse transform sampling,参考https://www.cnblogs.com/vpegasus/p/sampling.html
假设一个分布的概率密度函数为
f
(
x
)
f(x)
f(x),其累计分布函数为
F
(
x
)
=
P
{
X
<
x
}
F(x)=P\{X<x\}
F(x)=P{X<x},假设有一个样本集来自此分布函数,那么,将每个样本带入累积分布函数,都可以得到对应的
F
(
x
)
F(x)
F(x)的值。如果样本集来自随机采样,与其对应的累积分布函数值也会是随机的,且服从
(
0
,
1
)
(0,1)
(0,1)上的均匀分布。于是,为了得到这个的样本,可以将这个过程颠倒过来,即在
(
0
,
1
)
(0,1)
(0,1)的均匀分布上随机采样,然后带入累积分布函数的反函数中,即可得到服从此分布的样本 。
同时可以参考https://blog.csdn.net/yjr3426619/article/details/102706968
当
y
y
y和
x
x
x的累计分布函数相同,即
P
{
Y
<
=
y
(
x
)
}
=
P
{
X
<
=
x
}
P\{Y<=y(x)\}=P\{X<=x\}
P{Y<=y(x)}=P{X<=x}
对两边求导,
p
(
y
(
x
)
)
d
(
y
)
d
(
x
)
=
p
(
x
)
p(y(x))\frac{d(y)}{d(x)}=p(x)
p(y(x))d(x)d(y)=p(x),
因为微分表面有
d
w
=
s
i
n
(
θ
)
d
θ
d
ϕ
dw = sin(\theta) d \theta d \phi
dw=sin(θ)dθdϕ
所以
p
(
w
)
s
i
n
(
θ
)
=
p
(
θ
,
ϕ
)
p(w)sin(\theta) = p(\theta, \phi)
p(w)sin(θ)=p(θ,ϕ)
因为
p
(
w
)
=
1
4
π
p(w)=\frac{1}{4\pi}
p(w)=4π1,所以
p
(
θ
,
ϕ
)
=
s
i
n
(
θ
)
4
π
p(\theta,\phi)=\frac{sin(\theta)}{4\pi}
p(θ,ϕ)=4πsin(θ),对
θ
\theta
θ 积分得到
p
(
θ
)
=
∫
0
2
π
p
(
θ
,
ϕ
)
d
ϕ
=
s
i
n
(
θ
)
2
p(\theta)=\int_{0}^{2\pi}p(\theta,\phi)d\phi=\frac{sin(\theta)}{2}
p(θ)=∫02πp(θ,ϕ)dϕ=2sin(θ),因为
θ
\theta
θ 和
ϕ
\phi
ϕ 相互独立,那么
p
(
ϕ
)
=
1
2
π
p(\phi)=\frac{1}{2\pi}
p(ϕ)=2π1,取两个均匀分布的变量
ξ
x
\xi_x
ξx 和
ξ
y
\xi_y
ξy,我们应保证累积分布函数是一样的,
所以
ξ
x
=
{
∫
0
θ
p
Θ
(
x
)
d
x
}
−
1
\xi_x=\{\int_{0}^{\theta}p_{\Theta}(x)dx\}^{-1}
ξx={∫0θpΘ(x)dx}−1 则
θ
=
arccos
(
1
−
2
ξ
x
)
\theta=\arccos{(1-2\xi_x)}
θ=arccos(1−2ξx)
,同时
ξ
y
=
{
∫
0
ϕ
p
Φ
(
x
)
d
x
}
−
1
\xi_y=\{\int_{0}^{\phi}p_{\Phi}(x)dx\}^{-1}
ξy={∫0ϕpΦ(x)dx}−1 ,则
ϕ
=
2
π
ξ
y
\phi=2\pi \xi_y
ϕ=2πξy
正交基函数
所有的SH光照paper都假设使用了正交基。使用正交基逼近一个函数,就是把某个函数投影到这些正交基上,然后把投影得到的参数与正交基相乘得到最后还原的函数。
这类正交基有很多很多,但是球谐使用到的则是伴随勒让德多项式(Associated Legendre Polynomials)
伴随勒让德多项式的基地函数有两个变量 l l l 和 m m m,其中 l l l 被称为 band index 只能去正整数,而 m m m 的区间则是 [ 0 , l ] [0,l] [0,l]
伴随勒让德多项式有如下递归关系
P
l
m
=
x
(
2
l
−
1
)
P
l
−
1
m
−
(
l
+
m
−
1
)
P
l
−
2
m
P
m
m
=
(
−
1
)
m
(
2
m
−
1
)
!
!
(
1
−
x
2
)
m
/
2
P
m
+
1
m
=
x
(
2
m
+
1
)
P
m
m
\begin{aligned} P_{l}^{m} &= x(2 l-1) P_{l-1}^{m}-(l+m-1) P_{l-2}^{m}\\ P_{m}^{m} &= (-1)^{m}(2 m-1) ! !\left(1-x^{2}\right)^{m / 2}\\ P_{m+1}^{m} &= x(2 m+1) P_{m}^{m} \end{aligned}
PlmPmmPm+1m=x(2l−1)Pl−1m−(l+m−1)Pl−2m=(−1)m(2m−1)!!(1−x2)m/2=x(2m+1)Pmm
球谐(Spherical Harmonics)
伴随勒让德多项式就是球谐的核心。球谐类似于傅里叶变换,只不过是定义在球面上的。而且SH通常是由虚数定义的,但是我们只关心的它的实数部分。
对于球面上的点
(
sin
θ
cos
φ
,
sin
θ
sin
φ
,
cos
θ
)
→
(
x
,
y
,
z
)
(\sin \theta \cos \varphi, \sin \theta \sin \varphi, \cos \theta) \rightarrow(x, y, z)
(sinθcosφ,sinθsinφ,cosθ)→(x,y,z)
SH函数一般都是用符号
y
y
y 来定义的
y
l
m
(
θ
,
φ
)
=
{
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}(\theta, \varphi)=\left\{\begin{array}{ll} \sqrt{2} K_{l}^{m} \cos (m \varphi) P_{l}^{m}(\cos \theta), & m>0 \\ \sqrt{2} K_{l}^{m} \sin (-m \varphi) P_{l}^{-m}(\cos \theta), & m<0 \\ K_{l}^{0} P_{l}^{0}(\cos \theta), & m=0 \end{array}\right.
ylm(θ,φ)=⎩⎨⎧2Klmcos(mφ)Plm(cosθ),2Klmsin(−mφ)Pl−m(cosθ),Kl0Pl0(cosθ),m>0m<0m=0
其中
P
P
P就是前面说的伴随勒让德多项式,而
K
K
K是标准化系数
K
l
m
=
(
2
l
+
1
)
4
π
(
l
−
∣
m
∣
)
!
(
l
+
∣
m
∣
)
!
K_{l}^{m}=\sqrt{\frac{(2 l+1)}{4 \pi} \frac{(l-|m|) !}{(l+|m|) !}}
Klm=4π(2l+1)(l+∣m∣)!(l−∣m∣)!
为了方便的生成所有的SH函数,我们把
l
l
l和
m
m
m的定义跟勒让德多项相比稍微修改一下,这里的
m
m
m的
−
l
-l
−l任然对应勒让德多项中的0
y
l
m
(
θ
,
φ
)
where
l
∈
R
+
,
−
l
≤
m
≤
l
y_{l}^{m}(\theta, \varphi) \quad \text { where } l \in \mathbf{R}^{+},-l \leq m \leq l
ylm(θ,φ) where l∈R+,−l≤m≤l
y
l
m
(
θ
,
φ
)
=
y
i
(
θ
,
φ
)
where
i
=
l
(
l
+
1
)
+
m
y_{l}^{m}(\theta, \varphi)=y_{i}(\theta, \varphi) \quad \text { where } i=l(l+1)+m
ylm(θ,φ)=yi(θ,φ) where i=l(l+1)+m
球谐投影(SH Projection)
把空间分布函数投影为SH系数则非常简单了,就是函数与球谐基函数在球面的积分。
c
l
m
=
∫
S
f
(
s
)
y
l
m
(
s
)
d
s
c_{l}^{m}=\int_{S} f(s) y_{l}^{m}(s) d s
clm=∫Sf(s)ylm(s)ds
而要根据这些投影到球谐基函数上的系数,重构逼近原函数则只需要把上面的过程反过来,
f
~
(
s
)
=
∑
l
=
0
n
−
1
∑
m
=
−
l
l
c
l
m
y
l
m
(
s
)
=
∑
i
=
0
n
2
c
i
y
i
(
s
)
\tilde{f}(s)=\sum_{l=0}^{n-1} \sum_{m=-l}^{l} c_{l}^{m} y_{l}^{m}(s)=\sum_{i=0}^{n^{2}} c_{i} y_{i}(s)
f~(s)=l=0∑n−1m=−l∑lclmylm(s)=i=0∑n2ciyi(s)
我们可以看到要逼近一个
n
n
n 阶函数需要
n
2
n^2
n2 个系数。当然如果有无数个系数,我们确实可以完美的重构出原函数,但实际中我们只会使用带宽受限的逼近。
我们在进行球面积分的时候,即对每个量化的球面小块,球面微分
d
w
d w
dw 会表示为
sin
(
θ
)
d
ϕ
d
θ
\sin(\theta) d \phi d \theta
sin(θ)dϕdθ
这里的
sin
θ
\sin \theta
sinθ 是因为赤道附近的小块对最终积分的结果影响是最多的。所以我们可以写出单个系数
c
i
=
∫
0
2
π
∫
0
π
light
(
θ
,
φ
)
y
i
(
θ
,
φ
)
sin
θ
d
θ
d
φ
c_{i}=\int_{0}^{2 \pi} \int_{0}^{\pi} \operatorname{light}(\theta, \varphi) y_{i}(\theta, \varphi) \sin \theta d \theta d \varphi
ci=∫02π∫0πlight(θ,φ)yi(θ,φ)sinθdθdφ
把我们之前说到的蒙特卡洛采样拿出来,
∫
S
f
d
s
≈
1
N
∑
j
=
1
N
f
(
x
j
)
w
(
x
j
)
\int_{S} f d s \approx \frac{1}{N} \sum_{j=1}^{N} f\left(x_{j}\right) w\left(x_{j}\right)
∫Sfds≈N1j=1∑Nf(xj)w(xj)
对应于上面的公式,
x
j
x_j
xj是我们之前预采样的点,函数
f
f
f 则是
f
(
x
j
)
=
l
i
g
h
t
(
x
j
)
y
i
(
x
j
)
f\left(x_{j}\right)=light\left(x_{j}\right) y_{i}\left(x_{j}\right)
f(xj)=light(xj)yi(xj)
因为我们在球面上的采样是均匀分布,所以
p
(
x
j
)
=
1
4
π
p(x_j)=\frac{1}{4\pi}
p(xj)=4π1
则采样系数的计算则是如下
c
i
=
1
N
∑
j
=
1
N
light
(
x
j
)
y
i
(
x
j
)
4
π
=
4
π
N
∑
j
=
1
N
light
(
x
j
)
y
i
(
x
j
)
\begin{aligned} c_{i} &=\frac{1}{N} \sum_{j=1}^{N} \operatorname{light}\left(x_{j}\right) y_{i}\left(x_{j}\right) 4 \pi \\ &=\frac{4 \pi}{N} \sum_{j=1}^{N} \operatorname{light}\left(x_{j}\right) y_{i}\left(x_{j}\right) \end{aligned}
ci=N1j=1∑Nlight(xj)yi(xj)4π=N4πj=1∑Nlight(xj)yi(xj)
因为我们取的N不大,所以高频信息都丢失了,最终会导致重构出来的函数会有个负数的鳍,视觉表现上就是在物体的暗测会有环。
SH函数的特性(Properties of SH Functions)
- SH函数是标准正交的(orthonormal),
- SH函数具有旋转不变性(rotationally invariant),
即如果函数 g g g 是函数 f f f 的旋转后的版本,那么在做了SH投影之后
g ~ ( s ) = f ~ ( R ( s ) ) \widetilde{g}(s)=\widetilde{f}(R(s)) g (s)=f (R(s)) - illumination和transfer函数做卷积,都转换到SH系数后,结果就是他们的系数在 SH space 的乘积之和
∫ S L ~ ( s ) t ~ ( s ) d s = ∑ i = 0 n 2 L i t i \int_{S} \widetilde{L}(s) \widetilde{t}(s) d s=\sum_{i=0}^{n^{2}} L_{i} t_{i} ∫SL (s)t (s)ds=i=0∑n2Liti
Rotating Spherical Harmonics
SH 函数具有旋转不变性,
SH光照表面漫反射(SH Lighting Diffuse Surfaces)
Envmap的例子说明
目标是计算一个顶点上的某个方向 n n n 上的 diffuse
We consider the rendering of diffuse objects under distant illumination,
as specified by an environment map.
我们要计算一个顶点处的 Irradiance 即
E
(
n
)
E(\mathbf{n})
E(n),则我们需要对该顶点法线方向对应的半球面上的光照分布进行积分。
E
(
n
)
=
∫
Ω
(
n
)
L
(
ω
)
(
n
⋅
ω
)
d
ω
E(\mathbf{n})=\int_{\Omega(\mathbf{n})} L(\omega)(\mathbf{n} \cdot \omega) d \omega
E(n)=∫Ω(n)L(ω)(n⋅ω)dω
其中
ω
\omega
ω 是半球表面往外的单位向量,
n
\mathbf{n}
n 是该顶点的法线方向,而
L
(
ω
)
L(\omega)
L(ω) 则是光照分布。
这里要分清楚单位。
因为前9个球谐向量方向上的能量已经占据了绝大部分的能量,所以只需要使用前9个球谐向量就可以比较近似的逼近原函数了。
我们只使用球谐函数的整数部分,对于标准球面上的点
(
x
,
y
,
z
)
=
(
sin
θ
cos
ϕ
,
sin
θ
sin
ϕ
,
cos
θ
)
Y
00
(
θ
,
ϕ
)
=
0.282095
(
Y
11
;
Y
10
;
Y
1
−
1
)
(
θ
,
ϕ
)
=
0.488603
(
x
;
z
;
y
)
(
Y
21
;
Y
2
−
1
;
Y
2
−
2
)
(
θ
,
ϕ
)
=
1.092548
(
x
z
;
y
z
;
x
y
)
Y
20
(
θ
,
ϕ
)
=
0.315392
(
3
z
2
−
1
)
Y
22
(
θ
,
ϕ
)
=
0.546274
(
x
2
−
y
2
)
\begin{aligned} (x, y, z) &=(\sin \theta \cos \phi, \sin \theta \sin \phi, \cos \theta) \\ Y_{00}(\theta, \phi) &=0.282095 \\ \left(Y_{11} ; Y_{10} ; Y_{1-1}\right)(\theta, \phi) &=0.488603(x ; z ; y) \\ \left(Y_{21} ; Y_{2-1} ; Y_{2-2}\right)(\theta, \phi) &=1.092548(x z ; y z ; x y) \\ Y_{20}(\theta, \phi) &=0.315392\left(3 z^{2}-1\right) \\ Y_{22}(\theta, \phi) &=0.546274\left(x^{2}-y^{2}\right) \end{aligned}
(x,y,z)Y00(θ,ϕ)(Y11;Y10;Y1−1)(θ,ϕ)(Y21;Y2−1;Y2−2)(θ,ϕ)Y20(θ,ϕ)Y22(θ,ϕ)=(sinθcosϕ,sinθsinϕ,cosθ)=0.282095=0.488603(x;z;y)=1.092548(xz;yz;xy)=0.315392(3z2−1)=0.546274(x2−y2)
根据球谐函数的定义,原函数的逼近等于球谐系数与球谐基向量的乘积之和
L
(
θ
,
ϕ
)
=
∑
l
,
m
L
l
m
Y
l
m
(
θ
,
ϕ
)
E
(
θ
,
ϕ
)
=
∑
l
,
m
E
l
m
Y
l
m
(
θ
,
ϕ
)
\begin{aligned} L(\theta, \phi) &=\sum_{l, m} L_{l m} Y_{l m}(\theta, \phi) \\ E(\theta, \phi) &=\sum_{l, m} E_{l m} Y_{l m}(\theta, \phi) \end{aligned}
L(θ,ϕ)E(θ,ϕ)=l,m∑LlmYlm(θ,ϕ)=l,m∑ElmYlm(θ,ϕ)
我们再定义
A
=
(
n
⋅
ω
)
A=(\mathbf{n} \cdot \omega)
A=(n⋅ω)
A
(
θ
)
=
max
[
cos
θ
,
0
]
=
∑
l
A
l
Y
l
0
(
θ
)
A(\theta)=\max [\cos \theta, 0]=\sum_{l} A_{l} Y_{l 0}(\theta)
A(θ)=max[cosθ,0]=l∑AlYl0(θ)
可以计算得到
E
l
m
=
4
π
2
l
+
1
A
l
L
l
m
E_{l m}=\sqrt{\frac{4 \pi}{2 l+1}} A_{l} L_{l m}
Elm=2l+14πAlLlm
为了方便后面表示,这里定义一个新的变量
A
^
l
=
4
π
2
l
+
1
A
l
\hat{A}_{l}=\sqrt{\frac{4 \pi}{2 l+1}} A_{l}
A^l=2l+14πAl
那么Irradiance就可以通过球谐系数与球谐函数的乘积之和来近似逼近出来
E
~
(
θ
,
ϕ
)
=
∑
l
,
m
A
^
l
L
l
m
Y
l
m
(
θ
,
ϕ
)
\widetilde{E}(\theta, \phi)=\sum_{l, m} \hat{A}_{l} L_{l m} Y_{l m}(\theta, \phi)
E
(θ,ϕ)=l,m∑A^lLlmYlm(θ,ϕ)
其中
l
=
1
A
^
1
=
2
π
3
l
>
1
,
odd
A
^
l
=
0
l
even
A
^
l
=
2
π
(
−
1
)
l
2
−
1
(
l
+
2
)
(
l
−
1
)
[
l
!
2
l
(
l
2
!
)
2
]
\begin{aligned} l=1 & \hat{A}_{1}=\frac{2 \pi}{3} \\ l>1, \text { odd } & \hat{A}_{l}=0 \\ l \text { even } & \hat{A}_{l}=2 \pi \frac{(-1)^{\frac{l}{2}-1}}{(l+2)(l-1)}\left[\frac{l !}{2^{l}\left(\frac{l}{2} !\right)^{2}}\right] \end{aligned}
l=1l>1, odd l even A^1=32πA^l=0A^l=2π(l+2)(l−1)(−1)2l−1[2l(2l!)2l!]
而前几个数值是
A
^
0
=
3.141593
A
^
1
=
2.094395
A
^
2
=
0.785398
A
^
3
=
0
A
^
4
=
−
0.130900
A
^
5
=
0
A
^
6
=
0.049087
\begin{array}{ccc} \hat{A}_{0}=3.141593 & \hat{A}_{1}=2.094395 & \hat{A}_{2}=0.785398 \\ \hat{A}_{3}=0 & \hat{A}_{4}=-0.130900 & \hat{A}_{5}=0 & \hat{A}_{6}=0.049087 \end{array}
A^0=3.141593A^3=0A^1=2.094395A^4=−0.130900A^2=0.785398A^5=0A^6=0.049087
该paper还提供了使用9个参数后的矩阵形式表达,这也是在unity中使用到的形式
对于一个点上的向量n方向上的Irradiance,其中
n
t
=
(
x
,
y
,
z
,
1
)
n^t = (x,y,z,1)
nt=(x,y,z,1)
E
(
n
)
=
n
t
M
n
E(\mathbf{n})=\mathbf{n}^{t} M \mathbf{n}
E(n)=ntMn
其中
M
=
(
c
1
L
22
c
1
L
2
−
2
c
1
L
21
c
2
L
11
c
1
L
2
−
2
−
c
1
L
22
c
1
L
2
−
1
c
2
L
1
−
1
c
1
L
21
c
1
L
2
−
1
c
3
L
20
c
2
L
10
c
2
L
11
c
2
L
1
−
1
c
2
L
10
c
4
L
00
−
c
5
L
20
)
M=\left(\begin{array}{cccc} c_{1} L_{22} & c_{1} L_{2-2} & c_{1} L_{21} & c_{2} L_{11} \\ c_{1} L_{2-2} & -c_{1} L_{22} & c_{1} L_{2-1} & c_{2} L_{1-1} \\ c_{1} L_{21} & c_{1} L_{2-1} & c_{3} L_{20} & c_{2} L_{10} \\ c_{2} L_{11} & c_{2} L_{1-1} & c_{2} L_{10} & c_{4} L_{00}-c_{5} L_{20} \end{array}\right)
M=⎝⎜⎜⎛c1L22c1L2−2c1L21c2L11c1L2−2−c1L22c1L2−1c2L1−1c1L21c1L2−1c3L20c2L10c2L11c2L1−1c2L10c4L00−c5L20⎠⎟⎟⎞
c
1
=
0.429043
c
2
=
0.511664
c
3
=
0.743125
c
4
=
0.886227
c
5
=
0.247708
\begin{array}{c} c_{1}=0.429043 \quad c_{2}=0.511664 \\ c_{3}=0.743125 \quad c_{4}=0.886227 \quad c_{5}=0.247708 \end{array}
c1=0.429043c2=0.511664c3=0.743125c4=0.886227c5=0.247708
拆开后的形式是
E
(
n
)
=
c
1
L
22
(
x
2
−
y
2
)
+
c
3
L
20
z
2
+
c
4
L
00
−
c
5
L
20
+
2
c
1
(
L
2
−
2
x
y
+
L
21
x
z
+
L
2
−
1
y
z
)
+
2
c
2
(
L
11
x
+
L
1
−
1
y
+
L
10
z
)
\begin{aligned} E(\mathbf{n}) &=c_{1} L_{22}\left(x^{2}-y^{2}\right)+c_{3} L_{20} z^{2}+c_{4} L_{00}-c_{5} L_{20} \\ &+2 c_{1}\left(L_{2-2} x y+L_{21} x z+L_{2-1} y z\right) \\ &+2 c_{2}\left(L_{11} x+L_{1-1} y+L_{10} z\right) \end{aligned}
E(n)=c1L22(x2−y2)+c3L20z2+c4L00−c5L20+2c1(L2−2xy+L21xz+L2−1yz)+2c2(L11x+L1−1y+L10z)
消除Ring震荡失真
One way to minimize the ringing artifacts is to multiply in the frequency domain (which is a convolution in the spatial domain) by a kernel with projection coefficients that taper to zero as you approach the cut-off frequency.
这跟傅里叶变换的类似,门信号(矩形波形)做傅里叶变换,在频率空间就是sinc形式,当扩展到图像上,表现的就是幅值出现正负形式的震荡。频域到时域的变换也有相同的关系,这就是傅里叶变换的形式
而球谐变换类似于傅里叶变换,只不过是在球面空间上发生的。所以为了在还原到原空间光照分布的时候,因为我们只取了投影到前9个球谐基的系数,导致球谐信号出现了类似门信号的形式,肯定会出现失真,为了消除Ring震荡失真,我们需要在把这个球谐空间的门信号消解掉,StupidSH36[4]中提到把门形的突变边缘变成一个渐进到0的平滑形状。
这种方法被称为 Windowing。[5]中提到的就是对球谐基乘上一个 windowing函数 s i n c 4 sinc^4 sinc4,来引入额外的衰减,从而在频率空间消除门信号的形式,最终使得恢复的原信号没有负值。
在Deringing Spherical Harmonics[5]提到了解决方案。根据球谐的形式,我们知道球谐系数是可以旋转的,而且具有旋转不变性。[5]中提到,所有的其他球谐球谐表示都是由 zonal harmonics 旋转得到的,zonal harmonics 就是只有在 z 方向上有值,即只有l,m都是0的球谐表示。
[5]给出的最保守的滤波器
(
sin
π
l
w
π
l
w
)
4
\left(\frac{\sin \frac{\pi l}{w}}{\frac{\pi l}{w}}\right)^{4}
(wπlsinwπl)4
而且[5]还给了一个windowing的表,根据需要对频域空间做的windowing的大小,差值得到每个band上的乘数
w | inf | 16.7 | 11.3 | 10 | 9 | 7 | 5.6 |
---|---|---|---|---|---|---|---|
l1 | 1 | 0.9941 | 0.9872 | 0.9634 | 0.9602 | 0.9184 | 0.8915 |
l2 | 1 | 0.9766 | 0.9493 | 0.9355 | 0.9207 | 0.8710 | 0.8030 |
l3 | 1 | 0.9478 | 0.8880 | 0.8584 | 0.8270 | 0.7241 | 0.5904 |
Table 1: Window size and per band scaling coefficients
其中 l l l 是SH band index, w w w 是window函数边界变0的宽度
球谐光照实现例子
具体可以看看我写的软渲染器 SoftRenderer
参考文献
[1] Spherical Harmonic Lighting: The Gritty Details
[2] https://zh.wikipedia.org/wiki/輻照度
[4]StupidSH36