3DGS源码解读 - 高斯模型

高斯模型

三维高斯分布的概率密度函数定义为：
$$f(\mathbf{x})=\frac{1}{(2\pi {)}^{3/2}|\Sigma {|}^{1/2}}\exp \left(-\frac{1}{2}(\mathbf{x}-\mu {)}^T{\Sigma }^{-1}(\mathbf{x}-\mu )\right)$$三维高斯分布的概率密度函数可视化如下：

而在论文中选取的高斯分布概率密度函数的数学表达式如下：
$$f(\mathbf{x})=\exp \left(-\frac{1}{2}(\mathbf{x}{)}^T{\Sigma }^{-1}(\mathbf{x})\right)$$其中，$\Sigma$ 为 3D 高斯的协方差矩阵。
为什么 3D 高斯是椭球呢？
参考视频：【较真系列】讲人话-3d gaussian splatting全解(原理+代码+公式)【1】 捏雪球
事实上，不妨假设 $(x{)}^T{\Sigma }^{-1}(x)$ 等于常数 $d$，将公式 $d=(x{)}^T{\Sigma }^{-1}(x)$ 展开后就是一个椭球面。由于 $G(x)\in [0,1]$，所以 $(x{)}^T{\Sigma }^{-1}(x)\in [0,+\infty ]$，从 $[0,+\infty ]$ 逐个取不同的 $d$，就会形成大椭球面套小椭球面的情况，最终的形状就类似于一个椭球。

基本属性

高斯模型主要有四个基本属性：中心位置、协方差矩阵、体密度和球谐波系数，如下图所示：

def __init__(self, sh_degree, optimizer_type="default"):
    """ 初始化 3D 高斯模型的参数
    :param sh_degree: 球谐函数的最大次数，用于控制颜色表示的复杂度
    :param optimizer_type: 优化器类型
    """
    self.active_sh_degree = 0                 # 当前激活的球谐次数，初始为 0
    self.optimizer_type = optimizer_type      # 优化器类型
    self.max_sh_degree = sh_degree            # 允许的最大球谐次数
    self._xyz = torch.empty(0)                # 3D 高斯的中心位置（均值）
    self._features_dc = torch.empty(0)        # 第一个球谐系数，用于表示基础颜色
    self._features_rest = torch.empty(0)      # 其余的球谐系数，用于表示颜色的细节和变化
    self._scaling = torch.empty(0)            # 3D 高斯的尺度参数，控制高斯的宽度
    self._rotation = torch.empty(0)           # 3D 高斯的旋转参数，用四元数表示
    self._opacity = torch.empty(0)            # 3D 高斯的不透明度，控制可见性
    self.max_radii2D = torch.empty(0)         # 在 2D 投影中，每个高斯的最大半径
    self.xyz_gradient_accum = torch.empty(0)  # 用于累积 3D 高斯中心位置的梯度
    self.denom = torch.empty(0)               # 未明确用途的参数
    self.optimizer = None                     # 优化器，用于调整上述参数以改进模型
    self.percent_dense = 0
    self.spatial_lr_scale = 0
    self.setup_functions()

setup_functions 函数

self.setup_functions() 的作用主要是初始化激活函数和协方差矩阵。

We use a sigmoid activation function for α to constrain it in the [0 − 1) range and obtain smooth gradients, and an exponential activation function for the scale of the covariance for similar reasons.

属性	尺度（scale）	旋转（rotation）	不透明度（opacity）
激活函数	$e^x$	$v=\frac{v}{\max \left(|v{|}_p,ϵ\right)}$	$\frac{1}{1+e^{-x}}$
激活函数的反函数	$log(x)$		$log(\frac{x}{1-x})$

def setup_functions(self):
    def build_covariance_from_scaling_rotation(scaling, scaling_modifier, rotation):
        """ 构建一个表示 3D 高斯模型形变的 3x3 仿射变换矩阵 L, (N, 3, 3)
        :param scaling: 尺幅缩放矩阵
        :param scaling_modifier: 尺度缩放调节因子
        :param rotation: 旋转矩阵
        :return: 仿射矩阵右上角的部分，(N, 6)
        """
        # 仿射矩阵 L = R × S
        L = build_scaling_rotation(scaling_modifier * scaling, rotation)
        # 实际协方差矩阵 covariance = R × S × S^T × R^T
        actual_covariance = L @ L.transpose(1, 2)  # 计算实际的协方差矩阵
        symm = strip_symmetric(actual_covariance)  # 由于协方差矩阵的对称性，选取仿射矩阵右上角的部分，(N, 6)
        return symm

    # 初始化各个属性的激活函数
    self.scaling_activation = torch.exp
    self.scaling_inverse_activation = torch.log
    self.covariance_activation = build_covariance_from_scaling_rotation
    self.opacity_activation = torch.sigmoid
    self.inverse_opacity_activation = inverse_sigmoid
    self.rotation_activation = torch.nn.functional.normalize

值得注意的是协方差矩阵具有丰富的数学性质：
1）协方差矩阵是对称矩阵和半正定矩阵
设 $X=(X_1,X_2,...,X_n{)}^T$ 为 $n$ 维随机变量，则其协方差矩阵为
$$\Sigma ={\left(c_{ij}\right)}_{n\times n}=\left(\begin{array}{cccc}{c}_{11}& c_{12}& \cdots & c_{1n}\\ c_{21}& c_{22}& \cdots & c_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ c_{n1}& c_{n2}& \cdots & c_{nn}\end{array}\right)$$其中，$c_{ij}=Cov(X_i,X_j)$。由于 $c_{ij}=c_{ji}$，所以协方差矩阵是对称矩阵。
由于 $\Sigma$ 是对称矩阵，所以我们只需要保存矩阵中的上三角部分元素即可，代码中的 strip_symmetric 提取协方差矩阵的上三角元素的。

对任意非零向量 $x$ 有
$$\begin{array}{rl}& {\mathbf{x}}^{\mathrm{T}}C\mathbf{x}={\mathbf{x}}^{\mathrm{T}}\mathrm{E}\left[(X-\mu )(X-\mu {)}^{\mathrm{T}}\right]\mathbf{x}=\mathrm{E}\left[{\mathbf{x}}^{\mathrm{T}}(X-\mu )(X-\mu {)}^{\mathrm{T}}\mathbf{x}\right]\\ & =\mathrm{E}\left[{\left((X-\mu {)}^{\mathrm{T}}\mathbf{x}\right)}^{\mathrm{T}}\left((X-\mu {)}^{\mathrm{T}}\mathbf{x}\right)\right]=\mathrm{E}\left({\Vert (X-\mu {)}^{\mathrm{T}}\mathbf{x}\Vert }^2\right)>=0\end{array}$$因此协方差矩阵为半正定矩阵。

2）协方差矩阵可以表示为 $\Sigma =RS{S}^T{R}^T$ 的形式，其中 $R$ 为旋转矩阵，$S$ 为尺度矩阵
实对称矩阵的谱定理：任意实对称矩阵均可被正交对角化​​。这表示存在正交矩阵 $V$ 满足 $V{V}^T=I$ 和对角矩阵 $\Lambda$ 使得 $\Sigma =V\Lambda V^T$，其中 $V$ 的列是 $\Sigma$ 的特征向量，$\Lambda$ 的对角元素是 $\Sigma$ 的特征值。
几何意义​：特征向量方向​​表示数据最大方差的方向（主成分方向）；​​特征值大小​​表示数据在该方向上的方差大小，特征值越大，对应主成分包含的信息越多。
令 $S=\sqrt{\Lambda }$，则有 $\Lambda =S{S}^T$，不妨设 $R=V$，则有：$\Sigma =RS{S}^T{R}^T=RS(RS{)}^T$。
由上述两个性质可知，我们只需要维护 $R$ 和 $S$ 即可构建协方差矩阵 $\Sigma$。

问：协方差矩阵 $\Sigma$ 的物理意义
答：
1）$\Sigma$ 表示数据在三维空间中的分布形状，通常表现为一个椭球。椭球的三个主轴方向由协方差矩阵的特征向量 $v_1$，$v_2$ 和 $v_3$ 确定。椭球的轴长由特征值的平方根 $\sqrt{{\lambda }_1}$，$\sqrt{{\lambda }_2}$ 和 $\sqrt{{\lambda }_3}$ 决定，其中 ${\lambda }_i$ 是特征值，表示沿各主轴方向的方差。
2）若三个特征值相等（${\lambda }_1={\lambda }_2={\lambda }_3$），数据分布接近球体，表示各方向方差相同，无显著方向性。
3）若特征值差异较大，数据分布为拉伸的椭球，沿最大特征值对应的方向延伸最广。
4）若协方差矩阵的秩不足（如秩为 2 或 1），则数据分布退化为低维结构：
秩为 2：数据分布在一个平面上，椭球塌缩为椭圆（一个特征值为零）；
秩为 1：数据分布在一条直线上，椭球塌缩为线段（两个特征值为零）。

问：为什么要通过 $R$ 和 $S$ 两个矩阵来维护协方差矩阵 $\Sigma$，而不是直接初始化一个协方差矩阵呢？
答：

An obvious approach would be to directly optimize the covariance matrix Σ to obtain 3D Gaussians that represent the radiance field. However, covariance matrices have physical meaning only when they are positive semi-definite. For our optimization of all our parameters, we use gradient descent that cannot be easily constrained to produce such valid matrices, and update steps and gradients can very easily create invalid covariance matrices.

协方差矩阵是半正定的情况下才具有物理意义，如果直接优化协方差矩阵 $\Sigma$ 不能保证其一定是半正定矩阵。论文作者的解决方案是维护一个尺度矩阵 $S$ 和四元数 $q$（四元数转换为旋转矩阵 $R$）。

create_from_pcd 函数

Our goal is to optimize a scene representation that allows highquality novel view synthesis, starting from a sparse set of (SfM) points without normals.

3DGS 的输入是由 SFM 得到三维点云，这些点云不需要法向量。为什么不需要法向量呢？论文中给出的解释为：

Given the extreme sparsity of SfM points it is very hard to estimate normals. Similarly, optimizing very noisy normals from such an estimation would be very challenging. Instead, we model the geometry as a set of 3D Gaussians that do not require normals.

论文作者认为从稀疏的 SFM 点云中计算法向量太困难了。
create_from_pcd 的目的是将点云数据转换为可优化的高斯参数（位置、颜色、尺度、旋转、不透明度），代码如下：

 def create_from_pcd(self, pcd: BasicPointCloud, cam_infos: int, spatial_lr_scale: float):
     # 设置空间学习率缩放因子
     self.spatial_lr_scale = spatial_lr_scale
     fused_point_cloud = torch.tensor(np.asarray(pcd.points)).float().cuda()

     # 将点云颜色（RGB）转换为球谐函数（Spherical Harmonics, SH）的 0 阶系数
     fused_color = RGB2SH(torch.tensor(np.asarray(pcd.colors)).float().cuda())
     # 创建一个全零张量，用于存储所有球谐系数（包括 0 阶和高阶）
     features = torch.zeros((fused_color.shape[0], 3, (self.max_sh_degree + 1) ** 2)).float().cuda()
     features[:, :3, 0] = fused_color
     features[:, 3:, 1:] = 0.0

     print("Number of points at initialisation : ", fused_point_cloud.shape[0])

     # 计算每个点与其最近邻点的平方距离，用于初始化高斯尺度
     # distCUDA2(...)：调用 CUDA 加速的函数，计算点云中每个点到最近邻的平方距离。
     dist2 = torch.clamp_min(distCUDA2(torch.from_numpy(np.asarray(pcd.points)).float().cuda()), 0.0000001)
     # 将距离限制在最小 1e-7，避免数值不稳定
     scales = torch.log(torch.sqrt(dist2))[..., None].repeat(1, 3)

     # 初始化旋转参数为单位四元数（无旋转），四元数格式为 (w, x, y, z)，其中 w=1 表示无旋转
     rots = torch.zeros((fused_point_cloud.shape[0], 4), device="cuda")
     rots[:, 0] = 1

     # 通过逆激活函数初始化不透明度参数
     opacities = self.inverse_opacity_activation(
         0.1 * torch.ones((fused_point_cloud.shape[0], 1), dtype=torch.float, device="cuda")
     )

     # 定义可优化参数，包括：3D高斯的位置、尺度、旋转、不透明度和球谐系数
     self._xyz = nn.Parameter(fused_point_cloud.requires_grad_(True))
     self._features_dc = nn.Parameter(features[:, :, 0:1].transpose(1, 2).contiguous().requires_grad_(True))
     self._features_rest = nn.Parameter(features[:, :, 1:].transpose(1, 2).contiguous().requires_grad_(True))
     self._scaling = nn.Parameter(scales.requires_grad_(True))
     self._rotation = nn.Parameter(rots.requires_grad_(True))
     self._opacity = nn.Parameter(opacities.requires_grad_(True))

     # 记录每个高斯点在 2D 投影中的最大半径
     self.max_radii2D = torch.zeros((self.get_xyz.shape[0]), device="cuda")

     # 创建相机名称到索引的字典，用于管理不同相机的曝光参数
     self.exposure_mapping = {cam_info.image_name: idx for idx, cam_info in enumerate(cam_infos)}
     self.pretrained_exposures = None
     # 为每个相机初始化一个 3x4 的可学习曝光矩阵
     exposure = torch.eye(3, 4, device="cuda")[None].repeat(len(cam_infos), 1, 1)
     self._exposure = nn.Parameter(exposure.requires_grad_(True))

代码中的一些疑问：
1）RGB2SH 的作用

def RGB2SH(rgb):
    return (rgb - 0.5) / C0

RGB2SH 是为了获取球谐函数的 0 阶系数。