3D Gaussian Splatting代码及论文解析

前言

上一篇文章提到了最近在看3DGS的论文，所以这里也把自己的3DGS的学习记录发出来和大家讨论，可以结合上次的翻译文章一起看，欢迎批评指正；当然我更推荐各位好好读一下原文，我在使用3DGS的过程中发现很多的缺陷和细节在论文中都有提到，认真看一下论文也许能找到更好的idea，这篇文章也不过是对论文的一点个人见解罢了，如果能帮到各位，我很开心；

首先三维高斯是三维重建的一个方法，所以首要任务就是表示三维场景，虽然表示了三维场景但是人眼接受的照片还是二维平面，所以将三维场景渲染成为二维平面是第二个任务。这篇文章的解析也是从这两个方面结合代码进行整理的；

先放上源代码链接：https://github.com/graphdeco-inria/gaussian-splatting/

三维场景构建——如何表示一个三维场景

三维重建的任务是将场景重建出来，NeRF是隐式辐射场表示，3DGS就是显式辐射场表示；3DGS 属于显示辐射场，三维高斯形状是一个椭球，用大量的3D高斯建模三维场景，同时利用了神经网络的特性进行参数优化；采用了3D高斯作为表示方法，可以用更快的训练速度和实时性能实现高质量渲染；三维高斯公式表示为公式(1)

$$N(x;\mu ,\Sigma )=\frac{1}{\sqrt{(2\pi {)}^k|\Sigma |}}\exp \left(-\frac{1}{2}(x-\mu {)}^T{\Sigma }^{-1}(x-\mu )\right)(1)$$

其中$\mu$表示三维高斯的中心坐标，$x$为某个真实坐标点，$\Sigma$是协方差矩阵

三维场景中，协方差矩阵$\Sigma$ 可以表示为公式(2)（为什么要用旋转和尺度表示？）

$$\Sigma =RS{S}^T{R}^T(2)$$

$R$：Rotation Matrix，表示三维高斯点相对于世界坐标系的旋转

$S$：Scaling Matrix，表示三维高斯三个轴的缩放，这里的缩放给三维高斯赋予了物理意义，加入了单位？

实际代码中S是一个3维Vector ，R用一个个四元数quaternion表示（注意四元数的归一化）

def build_scaling_rotation(s, r):
    L = torch.zeros((s.shape[0], 3, 3), dtype=torch.float, device="cuda")
    R = build_rotation(r) #quat to rotation matrix
		# 每一个矩阵第1维都是三维高斯的数量
    L[:,0,0] = s[:,0]
    L[:,1,1] = s[:,1]
    L[:,2,2] = s[:,2]

    L = R @ L
		return L

def strip_lowerdiag(L):
    uncertainty = torch.zeros((L.shape[0], 6), dtype=torch.float, device="cuda")

    uncertainty[:, 0] = L[:, 0, 0]
    uncertainty[:, 1] = L[:, 0, 1]
    uncertainty[:, 2] = L[:, 0, 2]
    uncertainty[:, 3] = L[:, 1, 1]
    uncertainty[:, 4] = L[:, 1, 2]
    uncertainty[:, 5] = L[:, 2, 2]
    return uncertainty

def strip_symmetric(sym):
    return strip_lowerdiag(sym) # 上三角矩阵，在光栅化渲染部分调用？
    
def build_covariance_from_scaling_rotation(scaling, scaling_modifier, rotation):
    L = build_scaling_rotation(scaling_modifier * scaling, rotation)
    actual_covariance = L @ L.transpose(1, 2)
    symm = strip_symmetric(actual_covariance)
    return symm    

三维场景渲染——三维到二维的转换

• 如何表示一个三维场景的颜色，在我们选择使用三维高斯来表示三维场景后我们应该如何表示三维场景的颜色？

NERF中是通过网络训练颜色表示，不使用球谐函数表示；

假设有n个高斯球，第$i$个高斯球的颜色表示为

$$C_i={\alpha }_i{c}_i(3)$$

其中$i=1,\dots ,n$，$\alpha$是不透明度，$c$是颜色，颜色可以用球谐函数表示

$$c_i^j=\sum _{l=0}^{l_{max}}\sum _{m=-l}^l{h}_{i,lm}{Y}_l^m(v_i^j)(4)$$

其中

• $j$表示第j个视角，$v_i^j$表示第$i$个高斯球在第$j$个视角的姿态
• $l$表示球谐函数的阶数，取值$[0,l_{max}]$，论文中最高用到4阶
• $m$表示当前阶数下的次数，取值$[-l,l]$，
• $h_{i,lm}$是待优化的球谐系数， H i = { h i , l m ∣ 0 ≤ l ≤ l m a x , − l ≤ m ≤ l } H_i = \{h_{i,lm} | 0 \le l \le l_{max}, −l \le m \le l \} Hi​={hi,lm​∣0≤l≤lmax​,−l≤m≤l}是所有系数集合
• $Y_l^m(v_i^j)$为当前阶数的球谐基函数，只和姿态有关，和距离

由于我们是在球面表示，所以姿态可以用偏航角$\varphi$和俯仰角$\theta$表示，公式(2)可以表示为

$$c_i(\varphi ,\theta )=\sum _{l=0}^{l_{max}}\sum _{m=-l}^l{h}_{i,lm}{Y}_l^m(\varphi ,\theta )(5)$$

从公式3可以知道，对于N阶球谐函数，需要求解${\sum }_{i=1}^N(i+1)=(N+1{)}^2$个球谐系数，例如三阶球谐函数需要求解从0阶到3阶的系数综合$1+3+5+7=16$个系数；

代码中定义了两个变量：_features_dc和_features_rest ，其中：

• _features_dc表示y_{0,0}对应的系数
• _features_rest表示除了y_{0,0}以外的基函对应的系数

因此假设第i个高斯球在图像上位置为(X,Y),该位置的颜色表示为

$$C=\sum _{i\in N}{c}_i{\alpha }_i\prod _{j=1}^{i-1}(1-{\alpha }_j),(6)$$

其中，第$i$个高斯球的不透明度(opacity)${\alpha }_i$表示为：

$${\alpha }_i=(1-\exp (-{\sigma }_i{\delta }_i))$$

${\delta }_i$表示沿着射线的步长，${\sigma }_i$表示采样密度；

代码中如何初始化

根据以上两节可以知道，3DGS要优化的参数有：位置信息XYZ、颜色信息RGB、透明度信息Opacity、尺度信息Scale、旋转矩阵rotation（四元数表示）；

def __init__(self, sh_degree : int):
  self.active_sh_degree = 0
  self.max_sh_degree = sh_degree  
  self._xyz = torch.empty(0)  # 每个 Gaussian 的中心位置，通常是从点云位置初始化，包括colmap计算的特征点云
  self._features_dc = torch.empty(0)  # SH基函数的常量对应的系数
  self._features_rest = torch.empty(0)  # SH基函数除了常量以外的部分
  self._scaling = torch.empty(0)  # 使用该行的缩放和下一行的旋转来表示椭球
  self._rotation = torch.empty(0)
  self._opacity = torch.empty(0)  # 不透明度
  self.max_radii2D = torch.empty(0)  # 2D半径
  self.xyz_gradient_accum = torch.empty(0)  # 梯度累积的值，用于优化过程（如剪枝）
  self.denom = torch.empty(0)  # grads = self.xyz_gradient_accum / self.denom
  self.optimizer = None
  self.percent_dense = 0  # 用于得到阈值，决定 split/clone
  self.spatial_lr_scale = 0  # 位置的 lr 和其他参数的 lr 不同
  self.setup_functions()

def create_from_pcd(self, pcd : BasicPointCloud, spatial_lr_scale : float):
    self.spatial_lr_scale = spatial_lr_scale
    fused_point_cloud = torch.tensor(np.asarray(pcd.points)).float().cuda()
    fused_color = RGB2SH(torch.tensor(np.asarray(pcd.colors)).float().cuda())
    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])

    dist2 = torch.clamp_min(distCUDA2(torch.from_numpy(np.asarray(pcd.points)).float().cuda()), 0.0000001)
    scales = torch.log(torch.sqrt(dist2))[...,None].repeat(1, 3)
    rots = torch.zeros((fused_point_cloud.shape[0], 4), device="cuda")
    rots[:, 0] = 1

    opacities = inverse_sigmoid(0.1 * torch.ones((fused_point_cloud.shape[0], 1), dtype=torch.float, device="cuda"))

    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))
    self.max_radii2D = torch.zeros((self.get_xyz.shape[0]), device="cuda")    

1. self._xyz ：每个 Gaussian 的中心位置，通常是从点云位置初始化，包括colmap计算的特征点云；矩阵大小为$[N,3,1]$
2. self._features_dc ：球谐函数的常量系数，也就是0阶函数的系数，从点云的颜色初始化，点云颜色归一化之后进行映射，公式为：$sh=(rgb-0.5)/c_0$；矩阵大小为$[N,3,1]$
3. self._features_rest ：球谐函数的高阶系数，全部初始化为0；矩阵大小为$[N,3,(deg+1{)}^2-1]$
4. self._scaling ：每个Gaussian三个轴的尺度系数，排序计算相邻三个点的平均距离，具体计算在Simple-knn的distCUDA2函数，初值为$scales=log(sqrt(dist2))$，**优化量为$log(scales)$；**矩阵大小为$[N,3,1]$
5. self._rotation ：每个Gaussian的旋转四元数，初始化为$[1,0,0,0]$，矩阵大小为$[N,4,1]$
6. self._opacity ：不透明度，初始值为0.1，优化量为$log(\frac{x}{1-x})$；矩阵大小为$[N,1,1]$

代码中对于待优化参数配置为torch中的nn.Parameter参数，该参数用法如下：https://blog.csdn.net/zyw2002/article/details/128170764

如何训练——3DGS优化过程

首先一个真实三维位置点$x$对应的像素位置$x^{\prime }$ 的颜色是通过如下公式获取的

$$C(x^{\prime })=\sum _{i\in N}{c}_i{\sigma }_i\prod _{j=1}^{i-1}(1-{\sigma }_j)，{\sigma }_i={\alpha }_i{G}^{\prime }(x^{\prime })$$

其中：$x^{\prime }$表示当前像素位置，$G^{\prime }(x^{\prime })$是之前提到的三维高斯函数$G(x)$对应视角的二维高斯，${\alpha }_i$为不透明度，${\sigma }_i$为$x$位置对应的不透明度，$N$为覆盖当前像素点三维高斯数量；

可以理解为不透明度${\alpha }_i$为光线透过的概率，通过二维高斯可以得到当前像素点的不透明度，和当前三维高斯颜色相乘可以得到一个颜色，最后的颜色为像素点上所有三维高斯在$x^{\prime }$位置的颜色期望；

1. 数据准备工作：

   1. 数据集初始化：相机参数，图像真值，点云数据
   2. 三维高斯模型初始化：参考前两节
   3. 渲染器初始化：
   4. 训练设置初始化：学习率，优化器，

2. 学习率更新：update_learning_rate

   • 公式参考get_expon_lr_func 函数
   • 延迟因子定义，使用正弦函数进行平滑控制，用于初始阶段的延迟调整，论文中$\text{lr\_delay\_mult}=0.01，\text{lr\_delay\_steps}=0$，因此延迟因子并未启用;

   $$\text{delay\_rate}(s)=\{\begin{array}{ll}\text{lr\_delay\_mult}+(1-\text{lr\_delay\_mult})\cdot \sin \left(0.5\cdot \pi \cdot \min \left(\frac{s}{\text{lr\_delay\_steps}},1\right)\right),& \text{if }\text{lr\_delay\_steps}>0\\ 1.0,& \text{otherwise}\end{array}$$

   

   假设lr_delay_mult=0.01，lr_delay_steps=500，这个函数是正弦曲线在500的位置达到极值

• 定义步数 $s$对应的归一化时间 $t$：

  $t=\min \left(\frac{s}{\text{max\_steps}},1\right)$

• 学习率 $\text{lr}(s)$可以通过以下公式计算：在 $log$空间中对$\text{lr\_init}$ 和$\text{lr\_final}$进行线性插值，并乘以延迟因子$\text{delay\_rate}(s)$，下面这个函数是一直在衰减的，如下图所示

  $\text{lr}(s)=\text{delay\_rate}(s)\cdot \exp \left(\log (\text{lr\_init})\cdot (1-t)+\log (\text{lr\_final})\cdot t\right)$

  

1. 球谐函数阶数更新：代码中每1000次增加一次球谐函数阶数，最大阶数为3

2. 渲染render

3. Loss计算：$L_1$的loss和图像相似度ssim的loss（在3DGS的衍生论文中设计了很多其他方面的loss）

4. 3DGS的稠密化和剪枝：densify_and_prune ，参数优化和密度控制两个过程交叉进行

   

   • Gaussian密集化：自适应的增加3D高斯的密度，来保证整个3D高斯适应场景细节，主要专注于几何特征缺少和高斯覆盖较大区域。此处通过grad阈值和场景范围阈值来控制高斯的拆分和克隆，grad阈值表征了视角空间下位置梯度，场景阈值表示世界空间场景范围；在大于grad阈值同时小于场景范围的情况下，复制一个高斯并移动到位置梯度的方向上；在小于grad阈值并且大于场景范围情况下，用按比例缩小的2个高斯替代；（场景范围是通过**getNerfppNorm**计算得到，使用所有相机平移$t$，计算一个中心点和半径，乘上一个预设的系数作为场景范围，场景半径的1%）
   • **Clone，densify_and_clone ：**对于xyz_gradient大于grad阈值，但是当前3维高斯的scaling小于阈值的进行克隆，新克隆的3维高斯不继承原高斯的梯度，这样在优化的时候2个高斯会向不同的方向优化

   

   • Split，densify_and_split ：对于xyz_gradient大于阈值，而且scaling大于阈值的进行分割，修改高斯的xyz和scaling，其余保持不变；

     新尺度的计算：$scalin{g}_{new}={\log }_e^{\frac{exp(scaling)}{0.8\ast N}}$，同时xyz是随机采样生成的

   • Prune，densify_and_prune：1)不透明度低于阈值也就是几乎透明的高斯；2) 二维高斯半径radius过大，即当前视图内过大；3) 三维高斯的scaling过大的，即在世界空间内体积过大；三个条件取交集，剔除这部分高斯

5. 不透明度重置**reset_opacity ：**文章的优化方法会导致相机附近的高斯密度不合理增加，因此每3000次进行一次不透明度的削弱，对其进行0.01的缩放

   

6. 梯度操作：梯度更新然后清零optimizer.**step**()和optimizer.**zero_grad**(set_to_none = True)

   函数作用参考：https://blog.csdn.net/PanYHHH/article/details/107361827

def training(dataset, opt, pipe, testing_iterations, saving_iterations, checkpoint_iterations, checkpoint, debug_from):
		# 初始化
    first_iter = 0
    tb_writer = prepare_output_and_logger(dataset)
    gaussians = GaussianModel(dataset.sh_degree) #模型
    scene = Scene(dataset, gaussians) # 数据集
    gaussians.training_setup(opt) # 学习率的设置
    if checkpoint:
        (model_params, first_iter) = torch.load(checkpoint)
        gaussians.restore(model_params, opt)

    bg_color = [1, 1, 1] if dataset.white_background else [0, 0, 0]
    background = torch.tensor(bg_color, dtype=torch.float32, device="cuda")

    iter_start = torch.cuda.Event(enable_timing = True)
    iter_end = torch.cuda.Event(enable_timing = True)

    viewpoint_stack = None
    ema_loss_for_log = 0.0
    progress_bar = tqdm(range(first_iter, opt.iterations), desc="Training progress")
    first_iter += 1
    for iteration in range(first_iter, opt.iterations + 1):        
        if network_gui.conn == None:
            network_gui.try_connect()
        while network_gui.conn != None:
            try:
                net_image_bytes = None
                custom_cam, do_training, pipe.convert_SHs_python, pipe.compute_cov3D_python, keep_alive, scaling_modifer = network_gui.receive()
                if custom_cam != None:
                    net_image = render(custom_cam, gaussians, pipe, background, scaling_modifer)["render"]
                    net_image_bytes = memoryview((torch.clamp(net_image, min=0, max=1.0) * 255).byte().permute(1, 2, 0).contiguous().cpu().numpy())
                network_gui.send(net_image_bytes, dataset.source_path)
                if do_training and ((iteration < int(opt.iterations)) or not keep_alive):
                    break
            except Exception as e:
                network_gui.conn = None

        iter_start.record()

        gaussians.update_learning_rate(iteration) # 学习率更新

        # Every 1000 its we increase the levels of SH up to a maximum degree
        if iteration % 1000 == 0:
            gaussians.oneupSHdegree() # 球谐函数阶数更新

        # Pick a random Camera
        if not viewpoint_stack:
            viewpoint_stack = scene.getTrainCameras().copy()
        viewpoint_cam = viewpoint_stack.pop(randint(0, len(viewpoint_stack)-1))

        # Render
        if (iteration - 1) == debug_from:
            pipe.debug = True

        bg = torch.rand((3), device="cuda") if opt.random_background else background

        render_pkg = render(viewpoint_cam, gaussians, pipe, bg)
        image, viewspace_point_tensor, visibility_filter, radii = render_pkg["render"], render_pkg["viewspace_points"], render_pkg["visibility_filter"], render_pkg["radii"]

        # Loss
        gt_image = viewpoint_cam.original_image.cuda()
        Ll1 = l1_loss(image, gt_image)
        loss = (1.0 - opt.lambda_dssim) * Ll1 + opt.lambda_dssim * (1.0 - ssim(image, gt_image))
        loss.backward()

        iter_end.record()

        with torch.no_grad():
            # Progress bar
            ema_loss_for_log = 0.4 * loss.item() + 0.6 * ema_loss_for_log
            if iteration % 10 == 0:
                progress_bar.set_postfix({"Loss": f"{ema_loss_for_log:.{7}f}"})
                progress_bar.update(10)
            if iteration == opt.iterations:
                progress_bar.close()

            # Log and save
            training_report(tb_writer, iteration, Ll1, loss, l1_loss, iter_start.elapsed_time(iter_end), testing_iterations, scene, render, (pipe, background))
            if (iteration in saving_iterations):
                print("\n[ITER {}] Saving Gaussians".format(iteration))
                scene.save(iteration)

            # Densification
            if iteration < opt.densify_until_iter:
                # Keep track of max radii in image-space for pruning
                gaussians.max_radii2D[visibility_filter] = torch.max(gaussians.max_radii2D[visibility_filter], radii[visibility_filter])
                gaussians.add_densification_stats(viewspace_point_tensor, visibility_filter)

                if iteration > opt.densify_from_iter and iteration % opt.densification_interval == 0:
                    size_threshold = 20 if iteration > opt.opacity_reset_interval else None
                    gaussians.densify_and_prune(opt.densify_grad_threshold, 0.005, scene.cameras_extent, size_threshold)
                
                if iteration % opt.opacity_reset_interval == 0 or (dataset.white_background and iteration == opt.densify_from_iter):
                    gaussians.reset_opacity()  # 为了防止输入摄像头附近的高斯密度出现不合理的增加，在迭代一定次数后，高斯的不透明度会被设置为接近于零

            # Optimizer step
            if iteration < opt.iterations:
                gaussians.optimizer.step()
                gaussians.optimizer.zero_grad(set_to_none = True)

            if (iteration in checkpoint_iterations):
                print("\n[ITER {}] Saving Checkpoint".format(iteration))
                torch.save((gaussians.capture(), iteration), scene.model_path + "/chkpnt" + str(iteration) + ".pth")

参考链接

• 数学原理的推导：3D Gaussian Splatting
• 一文入门3DGS：https://zhuanlan.zhihu.com/p/680669616，https://zhuanlan.zhihu.com/p/678877999
• https://www.bilibili.com/video/BV1aw411E7bA/?spm_id_from=pageDriver
• https://blog.csdn.net/gwplovekimi/category_12633639.html
• 源码解读：https://blog.csdn.net/gwplovekimi/article/details/135500438
• 3DGS论文仓库：https://github.com/Awesome3DGS/3D-Gaussian-Splatting-Papers，https://github.com/MrNeRF/awesome-3D-gaussian-splatting
• EWA Splatting：https://zhuanlan.zhihu.com/p/674099667，解释了为什么用JW相乘的原理
• 几个重要论文的翻译：https://www.zhihu.com/people/he-xiao-ming-54-84/posts
• https://zhuanlan.zhihu.com/p/680669616，https://zhuanlan.zhihu.com/p/673703478
• 球谐函数的原理：https://blog.csdn.net/weixin_44518102/article/details/132676961，https://wuli.wiki/online/SphHar.html，球谐函数及Legendre多项式
• 对球谐函数的简单理解：https://zhuanlan.zhihu.com/p/676761611

最后推广一下自己的公众号