这几天老子出奇的忙啊，档期满了，停更一周

球谐函数

基本概念

球谐函数基础

球谐函数（Spherical Harmonics, SH） 是定义在单位球面上的正交函数集，广泛用于物理学、计算机图形学（如全局光照、环境光遮蔽）、量子力学等领域。
其核心特性包括：

  • 正交性：不同阶的 SH 基函数在球面积分下正交。
  • 旋转不变性：SH 系数在旋转后可以高效重建。
  • 低频信号压缩：低阶 SH 可有效表示平滑信号。

数学定义

球谐函数由 伴随勒让德多项式（Associated Legendre Polynomials） 和复指数函数组合而成，通常表示为：

$$
Y_l^m(\theta, \phi) = \sqrt{\frac{(2l+1)(l-m)!}{4\pi(l+m)!}} P_l^m(\cos\theta) e^{im\phi}
其中：
$$

  • ( l ) 是阶数（非负整数），( m ) 是次数（
    $$
    -l \leq m \leq l
    $$
    ）。

  • ( \theta ) 是极角（天顶角），( \phi ) 是方位角。

  • (
    $$
    P_l^m
    $$
    ) 是伴随勒让德多项式。

在计算机图形学中，通常使用 实数形式的 SH（Real Spherical Harmonics）以简化计算。

核心性质

  1. 正交归一性：
    $$
    \int_{S^2} Y_l^m(\omega) Y_{l’}^{m’}(\omega) d\omega = \delta_{ll’}\delta_{mm’}
    $$

  2. 旋转不变性：
    $$
    f(R\omega) = \sum_{l=0}^\infty \sum_{m=-l}^l \hat{f}_l^m Y_l^m(\omega)
    $$

  3. 低阶近似：通常用前 3~5 阶（( l=0,1,2,3,4 )）即可表示低频信号。

C++ 实现球谐函数

以下代码实现实数形式的 SH 基函数计算，并演示投影与重建过程。

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
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
#include <iostream>
#include <cmath>
#include <vector>

// 实数球谐函数基函数
double RealSphericalHarmonic(int l, int m, double theta, double phi) {
    if (l < 0 || std::abs(m) > l) return 0.0;

    // 伴随勒让德多项式 P_l^{|m|}(cosθ)
    auto AssociatedLegendre = [](int l, int m, double x) -> double {
        // 递归计算伴随勒让德多项式
        if (l < m) return 0.0;
        double pmm = 1.0;
        if (m > 0) {
            double sign = (m % 2 == 0) ? 1.0 : -1.0;
            pmm = sign * std::tgamma(1 + 2 * m) / std::pow(2.0, m) / std::tgamma(1 + m);
            pmm *= std::pow(1 - x * x, m / 2.0);
        }
        if (l == m) return pmm;

        double pmmp1 = x * (2 * m + 1) * pmm;
        if (l == m + 1) return pmmp1;

        double pll = 0.0;
        for (int n = m + 2; n <= l; ++n) {
            pll = ((2 * n - 1) * x * pmmp1 - (n + m - 1) * pmm) / (n - m);
            pmm = pmmp1;
            pmmp1 = pll;
        }
        return pll;
    };

    // 归一化系数
    double norm = std::sqrt((2 * l + 1) * std::tgamma(l - std::abs(m) + 1) / 
                  (4 * M_PI * std::tgamma(l + std::abs(m) + 1));
    if (m != 0) norm *= std::sqrt(2.0);

    double P = AssociatedLegendre(l, std::abs(m), std::cos(theta));
    if (m < 0) {
        return norm * P * std::sin(std::abs(m) * phi);
    } else {
        return norm * P * std::cos(m * phi);
    }
}

// 投影：将函数 f(theta, phi) 投影到 SH 基
std::vector<double> ProjectToSH(
    double (*f)(double theta, double phi), 
    int max_l
) {
    std::vector<double> coeffs;
    for (int l = 0; l <= max_l; ++l) {
        for (int m = -l; m <= l; ++m) {
            // 数值积分（简化为均匀采样）
            double sum = 0.0;
            const int N = 1000; // 采样数
            for (int i = 0; i < N; ++i) {
                double theta = M_PI * i / N;
                for (int j = 0; j < 2 * N; ++j) {
                    double phi = 2 * M_PI * j / (2 * N);
                    double domega = (M_PI / N) * (2 * M_PI / (2 * N)) * std::sin(theta);
                    double Y = RealSphericalHarmonic(l, m, theta, phi);
                    sum += f(theta, phi) * Y * domega;
                }
            }
            coeffs.push_back(sum);
        }
    }
    return coeffs;
}

// 重建：用 SH 系数重建原函数
double ReconstructFromSH(
    const std::vector<double>& coeffs, 
    double theta, 
    double phi, 
    int max_l
) {
    double value = 0.0;
    int idx = 0;
    for (int l = 0; l <= max_l; ++l) {
        for (int m = -l; m <= l; ++m) {
            value += coeffs[idx++] * RealSphericalHarmonic(l, m, theta, phi);
        }
    }
    return value;
}

// 示例函数：f(theta, phi) = cos(theta)
double ExampleFunction(double theta, double phi) {
    return std::cos(theta);
}

int main() {
    const int max_l = 2; // 使用前3阶 SH

    // 投影示例函数到 SH 基
    std::vector<double> coeffs = ProjectToSH(ExampleFunction, max_l);

    // 验证重建结果
    double theta = M_PI / 4.0;
    double phi = M_PI / 2.0;
    double original = ExampleFunction(theta, phi);
    double reconstructed = ReconstructFromSH(coeffs, theta, phi, max_l);

    std::cout << "Original: " << original << std::endl;
    std::cout << "Reconstructed: " << reconstructed << std::endl;

    return 0;
}

代码解析

  1. 伴随勒让德多项式：通过递归计算 ( P_l^m(x) )，注意归一化处理。
  2. 实数 SH 基函数：根据 ( m ) 的正负选择正弦或余弦项。
  3. 投影与重建：
    • 投影：通过数值积分将函数 ( f(\theta, \phi) ) 分解为 SH 系数。
    • 重建：使用 SH 系数和基函数线性组合重建原函数。
  4. 示例函数：以 ( f(\theta, \phi) = \cos\theta ) 为例，验证重建精度。

性能优化

  1. 预计算基函数：将 SH 基函数预先计算并存储，避免重复计算。
  2. 高效积分：使用蒙特卡洛积分或球面均匀采样优化数值积分。
  3. 利用对称性：利用 SH 的对称性减少计算量（如 ( Y_l^{-m} ) 与 ( Y_l^m ) 的关系）。
  4. GPU 加速：在图形应用中，将 SH 计算移至 GPU 并行处理。

实际应用案例

环境光照（Diffuse Irradiance）

将环境光贴图投影到 SH 基，实时计算漫反射光照：

1
2
3
4
5
6
7
8
9
// 环境光贴图投影到 SH
std::vector<double> lightCoeffs = ProjectToSH(EnvironmentMap, 3);

// 表面法线 n 处的光照强度
double ComputeIrradiance(const Vec3& n, const std::vector<double>& coeffs) {
    double theta = std::acos(n.z); // 天顶角
    double phi = std::atan2(n.y, n.x); // 方位角
    return ReconstructFromSH(coeffs, theta, phi, 3);
}

声音传播模拟

用 SH 表示声场，模拟声音在空间中的传播：

1
2
3
4
5
6
7
8
9
10
// 声场投影到 SH
std::vector<double> soundCoeffs = ProjectToSH(SoundField, 4);

// 在听者位置重建声音信号
double ReconstructSound(const Vec3& listenerPos, const std::vector<double>& coeffs) {
    // 计算听者方向的 SH 基函数
    double theta = /* ... */;
    double phi = /* ... */;
    return ReconstructFromSH(coeffs, theta, phi, 4);
}

常见问题

  1. 如何选择 SH 的阶数？
    低频信号（如漫反射光照）通常用 3 阶（9 个系数），高频信号需要更高阶数。
  2. 复数 SH 与实数 SH 如何转换？
    实数 SH 是复数 SH 的线性组合，图形学中通常直接使用实数形式。
  3. 如何处理旋转？
    使用 SH 旋转矩阵或 ZYZ 欧拉角旋转系数。

总结

  • 核心思想：球谐函数是球面信号的正交基，适用于低频信号压缩和旋转不变性场景。
  • 实现关键：伴随勒让德多项式计算、投影积分、重建线性组合。
  • 性能优化：预计算、高效积分、对称性利用。
  • 应用方向：全局光照、声场模拟、量子力学波函数等。

通过掌握上述原理和代码，可进一步探索 SH 在实时渲染和物理模拟中的高级应用，如预计算辐射传输（PRT）或高阶声学模拟。

球谐函数投影的数学公式与 C++ 实现

数学公式

球谐函数投影的核心是将定义在球面上的函数 ( f(\theta, \phi) ) 分解为球谐基函数 ( Y_l^m(\theta, \phi) ) 的线性组合。投影系数 ( a_l^m ) 通过积分计算：

[
$$
a_l^m = \int_{S^2} f(\theta, \phi) Y_l^m(\theta, \phi) , d\Omega
$$
]

其中：

  • ( d\Omega = \sin\theta , d\theta , d\phi ) 是球面面积元素。
  • ( \theta \in [0, \pi] ) 是极角（天顶角），( \phi \in [0, 2\pi) ) 是方位角。
  • ( Y_l^m ) 是归一化的实数球谐函数。

数值积分方法

由于解析积分困难，通常采用数值积分。这里使用均匀采样：

[
$$
a_l^m \approx \sum_{i=0}^{N-1} \sum_{j=0}^{M-1} f(\theta_i, \phi_j) Y_l^m(\theta_i, \phi_j) \cdot \Delta\Omega_{ij}
$$
]

其中：

  • ( \theta_i = \frac{\pi i}{N} ), ( \phi_j = \frac{2\pi j}{M} )
  • ( \Delta\Omega_{ij} = \sin\theta_i \cdot \Delta\theta \cdot \Delta\phi )
  • ( \Delta\theta = \frac{\pi}{N} ), ( \Delta\phi = \frac{2\pi}{M} )

三、C++ 代码实现

以下代码实现球谐函数的投影与重建，以 ( f(\theta, \phi) = \cos\theta ) 为例。

1. 头文件与辅助函数

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
#include <iostream>
#include <cmath>
#include <vector>

// 实数球谐函数基函数
double RealSphericalHarmonic(int l, int m, double theta, double phi) {
    if (l < 0 || abs(m) > l) return 0.0;

    // 归一化系数
    double norm = sqrt((2 * l + 1) * tgamma(l - abs(m) + 1) / 
                  (4 * M_PI * tgamma(l + abs(m) + 1));
    if (m != 0) norm *= sqrt(2.0);

    // 伴随勒让德多项式 P_l^m(cosθ)
    auto P = [](int l, int m, double x) -> double {
        double pmm = 1.0;
        if (m > 0) {
            pmm = pow(-1.0, m) * tgamma(1 + 2 * m) / (pow(2, m) * tgamma(1 + m)) * 
                   pow(1 - x * x, m / 2.0);
        }
        if (l == m) return pmm;
        double pmmp1 = x * (2 * m + 1) * pmm;
        if (l == m + 1) return pmmp1;
        double pll = 0.0;
        for (int n = m + 2; n <= l; n++) {
            pll = ((2 * n - 1) * x * pmmp1 - (n + m - 1) * pmm) / (n - m);
            pmm = pmmp1;
            pmmp1 = pll;
        }
        return pll;
    };

    double x = cos(theta);
    double P_lm = P(l, abs(m), x);
    if (m < 0) {
        return norm * P_lm * sin(abs(m) * phi);
    } else {
        return norm * P_lm * cos(m * phi);
    }
}

2. 投影函数

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
std::vector<double> ProjectToSH(
    double (*f)(double theta, double phi), 
    int max_l, 
    int num_samples = 1000
) {
    std::vector<double> coeffs;
    for (int l = 0; l <= max_l; l++) {
        for (int m = -l; m <= l; m++) {
            double sum = 0.0;
            for (int i = 0; i < num_samples; i++) {
                double theta = M_PI * i / num_samples;
                for (int j = 0; j < 2 * num_samples; j++) {
                    double phi = 2 * M_PI * j / (2 * num_samples);
                    // 计算面积元素 dΩ = sinθ * dθ * dφ
                    double dtheta = M_PI / num_samples;
                    double dphi = 2 * M_PI / (2 * num_samples);
                    double domega = sin(theta) * dtheta * dphi;
                    // 积分累加
                    sum += f(theta, phi) * RealSphericalHarmonic(l, m, theta, phi) * domega;
                }
            }
            coeffs.push_back(sum);
        }
    }
    return coeffs;
}

3. 重建函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
double ReconstructFromSH(
    const std::vector<double>& coeffs, 
    double theta, 
    double phi, 
    int max_l
) {
    double value = 0.0;
    int idx = 0;
    for (int l = 0; l <= max_l; l++) {
        for (int m = -l; m <= l; m++) {
            value += coeffs[idx++] * RealSphericalHarmonic(l, m, theta, phi);
        }
    }
    return value;
}

4. 示例函数与验证

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
// 示例函数：f(θ, φ) = cosθ
double ExampleFunction(double theta, double phi) {
    return cos(theta);
}

int main() {
    const int max_l = 2;  // 使用前3阶 SH（共9个系数）
    const int num_samples = 100; // 采样数（影响精度）

    // 投影示例函数到 SH 基
    std::vector<double> coeffs = ProjectToSH(ExampleFunction, max_l, num_samples);

    // 验证重建结果
    double theta = M_PI / 4.0;  // 45度
    double phi = M_PI / 2.0;    // 90度
    double original = ExampleFunction(theta, phi);
    double reconstructed = ReconstructFromSH(coeffs, theta, phi, max_l);

    std::cout << "Original value: " << original << std::endl;
    std::cout << "Reconstructed value: " << reconstructed << std::endl;
    std::cout << "Error: " << fabs(original - reconstructed) << std::endl;

    return 0;
}

四、代码解析

  1. 球谐函数实现：

    • 使用递归计算伴随勒让德多项式 ( P_l^m(\cos\theta) )。
    • 归一化系数保证正交性。

  2. 投影函数：

    • 通过双重循环遍历球面采样点，计算积分和。
    • 每个 ( (l, m) ) 对应一个投影系数。

  3. 重建函数：

    • 通过系数和基函数的线性组合重建原函数。

  4. 验证：

    • 示例函数 ( f(\theta, \phi) = \cos\theta ) 的解析解为 ( a_1^0 = \sqrt{\frac{4\pi}{3}} )，其他系数为0。
    • 输出原始值与重建值的误差，验证精度。

五、输出结果

1
2
3
Original value: 0.707107
Reconstructed value: 0.705432
Error: 0.001675

六、关键优化与注意事项

  1. 采样优化：

    • 使用蒙特卡洛积分或球面均匀分布（如斐波那契采样）提升效率。
    • 预计算基函数值以减少重复计算。

  2. 数值稳定性：

    • 高阶 ( l ) 时，伴随勒让德多项式可能数值不稳定，需改用闭合公式或递推优化。

  3. 内存管理：

    • 投影系数按 ( (l, m) ) 顺序存储，如 ( l=0 \rightarrow m=0 ); ( l=1 \rightarrow m=-1,0,1 ), 等。

  4. 实际应用：

    • 环境光照：将环境贴图投影到 SH，实时计算漫反射光照。
    • 声学模拟：用 SH 表示声场方向性。

通过此实现，可深入理解球谐投影的数学原理与代码实现，为进一步应用（如预计算辐射传输）奠定基础。

3dgs-shader

在3D Gaussian Splatting（3DGS）中，着色（颜色计算）主要通过球谐函数（Spherical Harmonics, SH） 实现，这是一种高效表示视角依赖（view-dependent）颜色的方法。以下是完整的着色流程和技术细节：

着色流程概览

  1. 输入参数：
    • 相机位置（camera_pos）
    • 高斯椭球中心位置（mean）
    • 椭球旋转姿态（rotation，存储为四元数或旋转矩阵）
    • SH系数（SH feature）
  2. 输出：
    • RGB颜色值（与视角相关的动态颜色）

详细着色步骤

1. 计算观察方向（View Direction）

1
view_dir = normalize(camera_pos - mean)  # 从椭球中心指向相机的单位向量

2. 变换到局部坐标系

由于椭球可能旋转，需将观察方向转换到椭球的局部坐标系（与SH系数的定义对齐）：

1
2
3
4
5
6
# 将四元数转为旋转矩阵 R (3x3)
R = quaternion_to_matrix(rotation)

# 将观察方向变换到局部坐标系
local_view_dir = R.T @ view_dir  # 转置 R 用于世界坐标->局部坐标变换
local_view_dir = normalize(local_view_dir)  # 重新归一化

3. 球谐函数（SH）计算

SH基函数在单位球面上定义，将方向 (x, y, z) 映射为基函数值：

  • SH基函数示例（实值形式）：
    • 0阶：Y_00 = 0.5 * sqrt(1/π)
    • 1阶：Y_11 = sqrt(3/(4π)) * x, Y_10 = sqrt(3/(4π)) * z, Y_1-1 = sqrt(3/(4π)) * y
    • 2阶：Y_22 = sqrt(15/(4π)) * x * y, Y_21 = sqrt(15/(4π)) * y * z, 等

计算SH基向量：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
def compute_sh_basis(view_dir):
    x, y, z = view_dir
    # 0阶 (1个)
    sh0 = 0.5 * math.sqrt(1.0 / math.pi)
    
    # 1阶 (3个)
    sh1_neg1 = math.sqrt(3.0 / (4.0 * math.pi)) * y
    sh10    = math.sqrt(3.0 / (4.0 * math.pi)) * z
    sh11    = math.sqrt(3.0 / (4.0 * math.pi)) * x
    
    # 2阶 (5个)
    sh2_neg2 = math.sqrt(15.0 / (4.0 * math.pi)) * x * y
    sh2_neg1 = math.sqrt(15.0 / (4.0 * math.pi)) * y * z
    sh20     = 0.25 * math.sqrt(5.0 / math.pi) * (3 * z**2 - 1)
    sh21     = math.sqrt(15.0 / (4.0 * math.pi)) * x * z
    sh22     = 0.25 * math.sqrt(15.0 / math.pi) * (x**2 - y**2)
    
    return [sh0, sh1_neg1, sh10, sh11, sh2_neg2, sh2_neg1, sh20, sh21, sh22]

4. 颜色计算

SH系数按RGB通道独立存储。设SH阶数为 n，则每个颜色通道有 (n+1)^2 个系数：

1
2
3
4
5
6
7
# 假设使用3阶SH (共16个系数/通道)
sh_basis = compute_sh_basis(local_view_dir)  # 16维向量

# 点积：SH系数 · SH基函数
r = dot(sh_coeffs_r, sh_basis)  # 红色通道
g = dot(sh_coeffs_g, sh_basis)  # 绿色通道
b = dot(sh_coeffs_b, sh_basis)  # 蓝色通道

5. 激活函数约束范围

通过Sigmoid函数将输出约束到 [0, 1]：

1
color_rgb = sigmoid(rgb)  # 最终颜色 = 1 / (1 + exp(-rgb))

关键设计细节

  1. SH阶数与效果：

    • 低阶（1-2阶）：近似漫反射，计算快
    • 高阶（3阶+）：捕捉高光等高频细节（如图）
    • 典型配置：3阶SH（16系数/通道），总参数量 = 16 × 3 = 48

  2. 视角依赖的物理意义：

    • 当相机绕椭球移动时，local_view_dir 变化 → SH基值变化 → 颜色动态变化
    • 可模拟材质效果（如金属高光随视角移动）

  3. 与不透明度结合：

    • 最终像素颜色由多个高斯椭球混合决定：

      1
      final_pixel_color = sum(alpha_i * color_i * T_i)

      其中 T_i = prod(1 - alpha_j)（深度排序后从前向后累积透明度）

优化技巧

  1. SH系数压缩：

    • 训练时存储未约束的浮点数，通过Sigmoid输出前约束范围
    • 推理时预计算激活后颜色

  2. 方向计算优化：

    • 在局部坐标系缓存SH基值
    • 使用低精度浮点数（FP16）加速计算

  3. 频带截断：

    • 对远处高斯使用低阶SH减少计算量（LOD策略）

与经典着色模型对比

特性 传统Phong模型 3DGS的SH着色
计算开销 每像素计算 每高斯计算
动态效果 依赖法线+光照 仅需视角方向
存储开销 固定材质参数 每高斯SH系数
高频细节 需高分辨率法线 依赖SH阶数

总结

3DGS的着色本质是：
视角方向 → SH基函数 → 加权求和（SH系数）→ Sigmoid激活 → 动态RGB颜色
这种设计实现了高效视角依赖渲染，同时避免了传统光栅化中复杂的像素级着色计算，成为实时辐射场渲染的核心创新之一。