球谐函数的简单介绍和PlenOctree中的代码实现

• 关于球谐函数的具体内容有很多很棒的介绍：

  • https://zhuanlan.zhihu.com/p/34343836?utm_id=0

  • https://zhuanlan.zhihu.com/p/351289217

    这里主要介绍实际的计算过程以及介绍代码为啥这样写，Plenoxel中的球谐函数是用CUDA写的，不过PlenOctree: https://github.com/sxyu/plenoctree/ 用的Python，还是比较容易理解的。

• 球谐函数在PlenOctree和Plenoxel以及大火的3D Gaussian中都用了，用来表示和观察方向相关的颜色。代码实现部分是将球谐系数转换为具体的RGB值。
  球谐函数就是用一组球谐系数对于球谐基加权求和，球谐基就是关于球坐标系下$(\theta ,\varphi )$的函数，下面的$y_l^m$代表第$l$阶的的SH Basis的第$m$个分量：
  （下面这个公式摘抄自PlenOctree：https://arxiv.org/pdf/2103.14024.pdf 的补充材料中)

$$c(d;k)=Sigmoid(\sum _{l=0}^{l_{max}}\sum _{m=-l}^l{K}_l^m{y}_l^m(d))$$

---

• STEP 1 $d\to (\theta ,\phi )$:

  $d$是$(x,y,z)$的单位向量，在单位球坐标系下：

• STEP 2 球谐函数的公式：( https://zhuanlan.zhihu.com/p/351289217 )

y l m ( θ , φ ) = { 2 K l m c o s ( m φ ) P l m ( c o s θ ) m > 0 2 K l m s i n ( − m φ ) P l − m ( c o s θ ) m < 0 K l 0 P l 0 ( c o s θ ) m = 0 y_l^m(\theta, \varphi) = \left\{

$$\begin{matrix} \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^0P_l^0(cos\theta) & m=0 \end{matrix}$$ \right. ylm​(θ,φ)=⎩ ⎨ ⎧​2 ​Klm​cos(mφ)Plm​(cosθ)2 ​Klm​sin(−mφ)Pl−m​(cosθ)Kl0​Pl0​(cosθ)​m>0m<0m=0​

$$\begin{gathered} P_n(x)=\frac{1}{2^n\cdot n!}\frac{d^n}{d{x}^n}[(x^2-1{)}^n] \\ {P}_l^m=(-1{)}^m(1-x^2{)}^{\frac{m}{2}}\frac{d^m}{d{x}^m}(P_l(x)) \\ {K}_l^m=\sqrt{\frac{(2l+1)(l-|m|)!}{4\pi (l+|m|)!}} \end{gathered}$$

• STEP 3 开始代入不同的$l,m$：

P 0 ( x ) = 1 { P 0 0 ( x ) = 1 P_0(x) = 1 \\ \left\{

$$\begin{matrix} P_0^0(x) = 1 \end{matrix}$$ \right. P0​(x)=1{P00​(x)=1​

P 1 ( x ) = x { P 1 0 ( x ) = x P 1 1 ( x ) = − ( 1 − x 2 ) 1 2 P_1(x) = x \\ \left\{

$$\begin{matrix} P_1^0(x) = x \\ P_1^1(x) = -(1-x^2)^{\frac{1}{2}} \end{matrix}$$ \right. P1​(x)=x{P10​(x)=xP11​(x)=−(1−x2)21​​

P 2 ( x ) = 3 x 2 − 1 2 { P 2 0 ( x ) = 3 x 2 − 1 2 P 2 1 ( x ) = − 3 x ( 1 − x 2 ) 1 2 P 2 2 ( x ) = 3 ( 1 − x 2 ) P_2(x) = \frac{3x^2-1}{2} \\ \left\{

$$\begin{matrix} P_2^0(x) = \frac{3x^2-1}{2} \\ P_2^1(x) = -3x(1-x^2)^{\frac{1}{2}} \\ P_2^2(x) = 3(1-x^2) \end{matrix}$$ \right. P2​(x)=23x2−1​⎩ ⎨ ⎧​P20​(x)=23x2−1​P21​(x)=−3x(1−x2)21​P22​(x)=3(1−x2)​

• STEP 4 求SH basis，即$y_l^m$:

  { − y 0 0 = K 0 0 \left\{

  $$\begin{matrix}- y_0^0 = K_0^0 \end{matrix}$$ \right. {−y00​=K00​​

  { y 1 − 1 = − 2 K 1 − 1 s i n ( φ ) ∣ s i n θ ∣ = − 2 K 1 − 1 y y 1 0 = K 1 0 c o s θ = K 1 0 z y 1 1 = − 2 K 1 − 1 c o s ( φ ) ∣ s i n θ ∣ = − 2 K 1 − 1 x \left\{

  $$\begin{matrix} y_1^{-1}= -\sqrt{2} K_1^{-1} sin(\varphi)|sin\theta| = -\sqrt{2} K_1^{-1}y \\ y_1^{0} = K_1^0 cos\theta = K_1^0 z \\ y_1^{1} = -\sqrt{2} K_1^{-1} cos(\varphi)|sin\theta| = -\sqrt{2} K_1^{-1}x \end{matrix}$$ \right. ⎩ ⎨ ⎧​y1−1​=−2 ​K1−1​sin(φ)∣sinθ∣=−2 ​K1−1​yy10​=K10​cosθ=K10​zy11​=−2 ​K1−1​cos(φ)∣sinθ∣=−2 ​K1−1​x​

  { y 2 − 2 = 3 2 K 2 − 2 s i n ( 2 φ ) s i n 2 θ = 6 2 K 2 − 2 x y y 2 − 1 = − 3 2 K 2 − 1 s i n ( φ ) ∣ s i n θ ∣ c o s θ = − 3 2 y z y 2 0 = 3 c o s 2 θ − 1 2 = K 2 0 2 z 2 − x 2 − y 2 2 y 2 1 = − 3 2 K 2 1 c o s ( φ ) ∣ s i n θ ∣ c o s θ = − 3 2 K 2 1 x z y 2 2 = 3 2 K 2 2 c o s ( 2 φ ) s i n 2 θ = 3 2 K 2 2 ( x 2 − y 2 ) \left\{

  $$\begin{matrix} y_2^{-2} = 3\sqrt{2} K_2^{-2} sin(2\varphi){sin^2\theta} = 6\sqrt{2} K_2^{-2}xy \\ y_2^{-1} = -3\sqrt{2} K_2^{-1} sin(\varphi) |sin\theta| cos\theta = -3\sqrt{2} yz \\ y_2^{0} = \frac{3 cos^2 \theta - 1}{2} = K_2^0 \frac{2z^2-x^2-y^2}{2} \\ y_2^{1} = -3\sqrt{2} K_2^1 cos(\varphi) |sin\theta| cos\theta = -3\sqrt{2} K_2^1 xz \\ y_2^{2} = 3\sqrt{2} K_2^2 cos(2 \varphi) sin^2\theta = 3\sqrt{2}K_2^2 (x^2-y^2) \end{matrix}$$ \right. ⎩ ⎨ ⎧​y2−2​=32 ​K2−2​sin(2φ)sin2θ=62 ​K2−2​xyy2−1​=−32 ​K2−1​sin(φ)∣sinθ∣cosθ=−32 ​yzy20​=23cos2θ−1​=K20​22z2−x2−y2​y21​=−32 ​K21​cos(φ)∣sinθ∣cosθ=−32 ​K21​xzy22​=32 ​K22​cos(2φ)sin2θ=32 ​K22​(x2−y2)​

• STEP 5 求$K_l^m$:

  $$K_l^m=\sqrt{\frac{(2l+1)(l-|m|)!}{4\pi (l+|m|)!}}$$

  K 0 0 = 1 4 π = 0.282094 2 K 1 − 1 = 2 K 1 1 = K 1 0 = 3 4 π = 0.4886 { 6 2 K 2 − 2 = 1.0925 − 3 2 K 2 − 1 = − 1.0925 K 2 0 2 = 0.31539 − 3 2 K 2 1 = − 1.0925 3 2 K 2 2 = 0.54627 K_0^0 = \frac{1}{\sqrt{4\pi}} = 0.282094\\ \sqrt{2}K_1^{-1} = \sqrt{2}K_1^{1} = K_1^0 = \sqrt{\frac{3}{4\pi}} = 0.4886 \\ \\ \left\{

  $$\begin{matrix} 6\sqrt{2} K_2^{-2} = 1.0925 \\ -3\sqrt{2}K_2^{-1} = -1.0925 \\ \frac{K_2^0}{2} = 0.31539 \\ -3\sqrt{2}K_2^{1} = -1.0925 \\ 3\sqrt{2}K_2^{2} = 0.54627 \end{matrix}$$ \right. K00​=4π ​1​=0.2820942 ​K1−1​=2 ​K11​=K10​=4π3​ ​=0.4886⎩ ⎨ ⎧​62 ​K2−2​=1.0925−32 ​K2−1​=−1.09252K20​​=0.31539−32 ​K21​=−1.092532 ​K22​=0.54627​

• 因此PlenOctree的代码

  ( https://github.com/sxyu/plenoctree/blob/master/nerf_sh/nerf/sh.py )

  就解释的通啦：

C0 = 0.28209479177387814
C1 = 0.4886025119029199
C2 = [
    1.0925484305920792,
    -1.0925484305920792,
    0.31539156525252005,
    -1.0925484305920792,
    0.5462742152960396
]
C3 = [
    -0.5900435899266435,
    2.890611442640554,
    -0.4570457994644658,
    0.3731763325901154,
    -0.4570457994644658,
    1.445305721320277,
    -0.5900435899266435
]
C4 = [
    2.5033429417967046,
    -1.7701307697799304,
    0.9461746957575601,
    -0.6690465435572892,
    0.10578554691520431,
    -0.6690465435572892,
    0.47308734787878004,
    -1.7701307697799304,
    0.6258357354491761,
]

def eval_sh(deg, sh, dirs):
    """
    Evaluate spherical harmonics at unit directions
    using hardcoded SH polynomials.
    Works with torch/np/jnp.
    ... Can be 0 or more batch dimensions.

    Args:
        deg: int SH deg. Currently, 0-3 supported
        sh: jnp.ndarray SH coeffs [..., C, (deg + 1) ** 2]
        dirs: jnp.ndarray unit directions [..., 3]

    Returns:
        [..., C]
    """
    assert deg <= 4 and deg >= 0
    assert (deg + 1) ** 2 == sh.shape[-1]
    C = sh.shape[-2]

    result = C0 * sh[..., 0]
    if deg > 0:
        x, y, z = dirs[..., 0:1], dirs[..., 1:2], dirs[..., 2:3]
        result = (result -
                C1 * y * sh[..., 1] +
                C1 * z * sh[..., 2] -
                C1 * x * sh[..., 3])
        if deg > 1:
            xx, yy, zz = x * x, y * y, z * z
            xy, yz, xz = x * y, y * z, x * z
            result = (result +
                    C2[0] * xy * sh[..., 4] +
                    C2[1] * yz * sh[..., 5] +
                    C2[2] * (2.0 * zz - xx - yy) * sh[..., 6] +
                    C2[3] * xz * sh[..., 7] +
                    C2[4] * (xx - yy) * sh[..., 8])

            if deg > 2:
                result = (result +
                        C3[0] * y * (3 * xx - yy) * sh[..., 9] +
                        C3[1] * xy * z * sh[..., 10] +
                        C3[2] * y * (4 * zz - xx - yy)* sh[..., 11] +
                        C3[3] * z * (2 * zz - 3 * xx - 3 * yy) * sh[..., 12] +
                        C3[4] * x * (4 * zz - xx - yy) * sh[..., 13] +
                        C3[5] * z * (xx - yy) * sh[..., 14] +
                        C3[6] * x * (xx - 3 * yy) * sh[..., 15])
                if deg > 3:
                    result = (result + C4[0] * xy * (xx - yy) * sh[..., 16] +
                            C4[1] * yz * (3 * xx - yy) * sh[..., 17] +
                            C4[2] * xy * (7 * zz - 1) * sh[..., 18] +
                            C4[3] * yz * (7 * zz - 3) * sh[..., 19] +
                            C4[4] * (zz * (35 * zz - 30) + 3) * sh[..., 20] +
                            C4[5] * xz * (7 * zz - 3) * sh[..., 21] +
                            C4[6] * (xx - yy) * (7 * zz - 1) * sh[..., 22] +
                            C4[7] * xz * (xx - 3 * yy) * sh[..., 23] +
                            C4[8] * (xx * (xx - 3 * yy) - yy * (3 * xx - yy)) * sh[..., 24])
    return result
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86