3DGS 其一：3D Gaussian Splatting for Real-Time Radiance Field Rendering

Reference：

1. 深蓝学院：NeRF基础与常见算法解析
2. GitHub: gaussian-splatting
3. 原文官网
4. A Survey on 3D Gaussian Splatting
5. 开始弃用NeRF？为什么Gaussian Splatting在自动驾驶场景如此受欢迎？

相关文章：

1. NeRF 其一：NeRF: Representing Scenes as Neural Radiance Fields for View Synthesis
2. NeRF 其二：Mip-NeRF
3. NeRF 其三：Instant-NGP
4. Structure-from-Motion Revisited(COLMAP) 流程介绍

系列文章：

1. 3DGS 其一：3D Gaussian Splatting for Real-Time Radiance Field Rendering

3D Gaussian Splatting 是表达三维场景的一种方式，不同于 NeRF 用一个点来表达，它是用一堆的 3D 高斯来表达。

问题定义：

1. 辐射场：辐射场是三维空间中光分布的表示，它捕捉光如何与环境中的表面和材料相互作用。从数学上讲，辐射场可以描述为函数 $L:R^5\to R^{+}$，其中 $L(x,y,z,\theta ,\varphi )$ 将空间中的一个点 $(x,y,z)$ 和由球面坐标 $(\theta ,\varphi )$ 指定的方向映射到非负辐射值。辐射场可以通过隐式或显式表示进行封装，每种表示都具有特定的场景表示和渲染优势。

2. 隐式辐射场：隐式辐射场表示场景中的光分布，而不明确定义场景的几何体。在深度学习时代，它经常使用神经网络来学习连续的体积场景表示。最突出的例子是 NeRF。在 NeRF 中，MLP 网络用于将一组空间坐标 $(x,y,z)$ 和观看方向 $(\theta ,\varphi )$ 映射到颜色和密度值。任何点的辐射度都不是明确存储的，而是通过查询神经网络实时计算的。因此，函数可以写成：
   $$L_{\text{implicit }}(x,y,z,\theta ,\varphi )=\text{ NeuralNetwork }(x,y,z,\theta ,\varphi ).$$这种格式允许对复杂场景进行可微分和紧凑的表示，尽管由于体积光线行进，渲染过程中的计算负载往往很高。

3. 显式辐射场：相反，显式辐射场直接表示离散空间结构中的光分布，例如体素网格或点集。该结构中的每个元素存储其在空间中的相应位置的辐射信息。这种方法允许更直接且通常更快地访问辐射数据，但代价是更高的内存使用率和潜在的更低分辨率。显式辐射场表示的一般形式可以写成：
   $$L_{\text{explicit }}(x,y,z,\theta ,\varphi )=\mathrm{DataStructure}[(x,y,z)]\cdot f(\theta ,\varphi ),(2)$$其中 DataStructure 可以是网格或点云，$f(\theta ,\varphi )$ 是基于观看方向修改辐射的函数。(查找表的感觉，内存相当爆炸)

4. 两全其美的3D Gaussian Splatting：3D GS 表示从隐式辐射场到显式辐射场的转变。它通过利用3D高斯作为灵活高效的表示，利用了这两种方法的优势。这些高斯系数经过优化，可以准确地表示场景，结合了基于神经网络的优化和显式结构化数据存储的优点。这种混合方法旨在通过更快的训练和实时性能实现高质量渲染，特别是对于复杂的场景和高分辨率输出。3D 高斯表示公式化为：
   $$L_{3\mathrm{DGS}}(x,y,z,\theta ,\varphi )=\sum _{i}G\left(x,y,z,{\mathbf{\mu }}_i,{\mathbf{\Sigma }}_i\right)\cdot c_i(\theta ,\varphi )$$

Implicit Neural Representation(隐式神经表达) 对场景建模导致了许多问题：

• query过程（可以理解为射线渲染）需要大量的采样，渲染方法成本很高
• 用了大型多层MLP，运算量大，占用内存高
• 不容易编辑
• 不能显式地对空间几何建模
• 导致“遗忘”问题(在不重新访问之前训练数据的情况下不断从流数据中学习时，会遇到灾难性的遗忘问题，这种情况在大场景下容易出现)

1. 预备知识

1.1 球谐函数

球谐函数这一块请阅读 Instant-NGP 一文内的相关片段：NeRF 其三：Instant-NGP，此处不再做赘述。

1.2 Splatting

Splatting(抛雪球) 是一种用于 光栅化(rasterizer) 3D 对象(如椭球)的技术。这些 3D 对象被映射到投影平面后得到的 2D 图形称为 splat，类似于一个点、圆、矩形或其他形状，就像雪球打在墙上留下的印记，能量从中心向外扩散并减弱(就像抛出一个雪球砸在墙上一样)。

换句话说，三维场景不用三维点表达，而用一个个高斯球来表达。将一个个高斯球投影到二维图像上的过程就称为 Splatting 算法。
该光栅化过程可以在 GPU 上并行处理，因为每个 Splat 之间是独立的。

1.3 $\alpha$ blending

$\alpha$ blending 算法主要解决“图层”叠加问题。

两幅图融合，其中图像 $I_1$ 的透明度为 ${\alpha }_1$(前景图像)，图像 $I_{BK}$ 为背景图像，融合公式计算如下：
$$I_{\mathrm{res}ult}=I_1\times {\alpha }_1+I_{BK}\times (1-{\alpha }_1)\text{\hspace{1em}(1)}$$那么按照深度由近到远的排序四张图像，其中，图像 $I_1$ 的透明度为 ${\alpha }_1$， 图像 $I_2$ 的透明度为 ${\alpha }_2$，图像 $I_3$ 的透明度为 ${\alpha }_3$，图像 $I_{BK}$ 为背景图像，融合公式计算如下：
$$\begin{array}{rl}{I}_{\text{result}}& =I_1\times {\alpha }_1+(1-{\alpha }_1)(I_2\times {\alpha }_2+(1-{\alpha }_2)(I_3\times {\alpha }_3+I_{BK}\times (1-{\alpha }_3)))\\ & ={\alpha }_1{I}_1+(1-{\alpha }_1){\alpha }_2{I}_2+(1-{\alpha }_1)(1-{\alpha }_2){\alpha }_3{I}_3+(1-{\alpha }_1)(1-{\alpha }_2)\left(1-{\alpha }_3\right)I_{BK}\end{array}\text{\hspace{1em}(2)}$$可以从底往上看，

• $I_{BK}$ 看成 $I$；
• $I_3\times {\alpha }_3+I_{BK}\times (1-{\alpha }_3)$ 合并成 $I^{\prime }$，将透明度 ${\alpha }_3$ 套入公式可得；
• $I_2\times {\alpha }_2+(1-{\alpha }_2)(I_3\times {\alpha }_3+I_{BK}\times (1-{\alpha }_3))$ 合并成 $I^{\prime \prime }$，将透明度 ${\alpha }_2$ 套入公式可得。

也可以从上往下看，直接得到 Eq.2 底部公式。

综上，$\alpha$ blending 公式可以写成：
$$C=\sum _{i\in \mathcal{N}}{c}_i{\alpha }_i\prod _{j=1}^{i-1}(1-{\alpha }_j)$$回顾体渲染公式：
$$C=\sum _{i=1}^N{T}_i(1-\exp (-{\sigma }_i{\delta }_i)){\mathbf{c}}_i\mathrm{with}{T}_i=\exp \left(-\sum _{j=1}^{i-1}{\sigma }_j{\delta }_j\right)$$使用 ${\alpha }_i$ 来替代 $1-\exp (-{\sigma }_i{\delta }_i)$，可得：
$$C=\sum _{i=1}^N{T}_i{\alpha }_i{\mathbf{c}}_i{\alpha }_i=(1-\exp (-{\sigma }_i{\delta }_i))\text{and}{T}_i=\prod _{j=1}^{i-1}(1-{\alpha }_i)$$可以发现，该公式和 $\alpha$ blending 公式一模一样。

这里可以得到一个很有意思的结论：$\alpha$ blending 和 体渲染 都是相同的公式，但是具体做法是完全不一样的：体渲染是 NeRF 的一套做法，$\alpha$ blending 是 3D Gaussian Splatting 的一套做法。很多游戏都是使用 $\alpha$ blending 方法进行光栅化的。体渲染速度很慢，而光栅化很快，而且它是图形学内很成熟得东西，有很多软件可以辅助加速，比如 OpenGL，这样就可以充分利用硬件和软件的性能。

1.4 多维高斯的协方差矩阵

如果一个随机变量 $\mathbf{x}$ 服从高斯分布高斯分布 $N(\mu ,\sigma )$，那么它的概率密度函数为：
$$p\left(x\right)=\frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\frac{1}{2}\frac{{\left(x-\mu \right)}^2}{{\sigma }^2}\right).$$它的高维形式为：
$$p\left(x\right)=\frac{1}{\sqrt{{\left(2\pi \right)}^N\det \left(\mathbf{\Sigma }\right)}}\exp \left(-\frac{1}{2}{\left(\mathbf{x}-\mathbf{\mu }\right)}^T{\mathbf{\Sigma }}^{-1}\left(\mathbf{x}-\mathbf{\mu }\right)\right)$$这里前面系数不太重要，公式可简写成：$G(\mathbf{x})=e^{-\frac{1}{2}}(\mathbf{x}{)}^T{\Sigma }^{-1}(\mathbf{x})$。

现在来看看协方差矩阵的性质：
二维的协方差矩阵可写成 $\left[\begin{array}{cc}a& c\\ c& b\end{array}\right]$ 的形式、三维的 协方差矩阵可写成 $\left[\begin{array}{ccc}a& d& e\\ d& b& f\\ e& f& c\end{array}\right]$ 的形式，这个矩阵一定是对称的。在斜对角不为零的情况下，它一定是正定的。对于这种对称矩阵，它一定能有一个正交的矩阵将它三角化，即：
$$\Sigma =P\Lambda P^T=P{\Lambda }^{\frac{1}{2}}({\Lambda }^{\frac{1}{2}}{)}^T{P}^T$$即 $\Sigma =RS{S}^T{R}^T$，也就是说，只要我们要去表达一个协方差矩阵，只要知道了 $R$ 和 $S$ 即可。而且 $R$ 和 $S$ 构建出的矩阵，一定是单位阵。因为 $R$ 是一个正交阵，只要 $S$ 不全为 $0$，构建出的就是一个正定矩阵。

那么这里的协方差矩阵有什么具体的含义呢？
假设有一二维矩阵 $\left[\begin{array}{cc}{\sigma }_1^2& 0\\ 0& {\sigma }_2^2\end{array}\right]$，该矩阵表示分布的两个维度 $(x_1,x_2{)}^T$ 间是没有相关性的，所以右上和左下矩阵的系数为 $0$。如果有相关性，这个地方就应该有系数：如果是正数，$x_1$ 增大 $x_2$ 就会增大；如果是负数，$x_1$ 增大 $x_2$ 就会减小。

在没有相关性时，将所有点投影到 $x_1$ 方向上，它符合 ${\sigma }_1$ 的分布，其分布如下图在 $x_1$ 轴上方黑线所示。

1.4.1 高斯与椭球体的关系

由于三维高斯不好表现，先来从二维高斯看三维高斯与椭球体的关系：
函数 $G(\mathbf{x})=e^{-\frac{1}{2}}(\mathbf{x}{)}^T{\Sigma }^{-1}(\mathbf{x})$，因为现在是二维，$\mathbf{x}=(x_1,x_2{)}^T$。当 $G(\mathbf{x})$ 为常数时，组成了这个常数的 $\mathbf{x}$，将是一个什么样的分布形状？
如果令 $G(\mathbf{x})$ 为常数，那么就是说 $(\mathbf{x}{)}^T{\Sigma }^{-1}(\mathbf{x})$ 应该是一个常数，因为 $e^{-\frac{1}{2}}$ 这一部分是固定计算的。还是以协方差矩阵 $\left[\begin{array}{cc}{\sigma }_1^2& 0\\ 0& {\sigma }_2^2\end{array}\right]$ 为例，则公式变成了：
$${\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]}^T\left[\begin{array}{cc}\frac{1}{{\sigma }_1}& 0\\ 0& \frac{1}{{\sigma }_2}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]=const,\text{即 }\frac{x_1^2}{{\sigma }_1^2}+\frac{x_2^2}{{\sigma }_2^2}=const$$这就是一个椭球，其长短轴与 $\sigma$ 以及 $const$ 有关。所以所有的等势面都会构成一个个的椭球，呈现上图中一圈一圈的样子。越内圈出现的等势面的概率越高，越往外等势面出现的概率越低。

那么该怎样找一个将大多概率包络的高斯呢？
那么就将高斯往外面找，比如高斯 $99$ 能量的时候，就把这个值找到，就会形成一个包络面，包络面外面点的概率是很低的，里面的总概率是 $99$，使用这个高斯球来代表高斯分布的边界。

1.4.2 世界坐标系下的三维高斯到二维像素平面投影过程

假设世界坐标系三维空间中的一个点符合以下高斯分布 ${\mathbf{x}}_{\mathbf{w}}\sim N({\mathbf{\mu }}_{\mathbf{w}},{\Sigma }_w)$，现在想知道它对应的像素坐标系对应的向量 $\mathbf{z}$ 所对应的均值 ${\mathbf{\mu }}_{\mathbf{z}}$ 和协方差 ${\Sigma }_z$。

我们知道世界坐标系要先转换到相机坐标系：已知 ${\mathbf{x}}_{\mathbf{c}}=W_{cw}{\mathbf{x}}_{\mathbf{w}}+T$，这里用 $W$ 和 $T$ 指代旋转和平移因为原文是这样写的。通过这个变换就能将世界坐标系下的点转换到相机坐标系，这时的 ${\mathbf{x}}_{\mathbf{c}}$ 符合分布 ${\mathbf{x}}_{\mathbf{c}}\sim N(W_{cw}{\mathbf{\mu }}_{\mathbf{w}}+T,W{\Sigma }_w{W}^T)$。

像素坐标 $\mathbf{z}$ 与 ${\mathbf{x}}_{\mathbf{c}}$ 间又是什么样的关系呢？
$$\left[\begin{array}{c}{\vec{z}}_{2\times 1}\\ 1\end{array}\right]=\frac{1}{x_{3c}}\left[\begin{array}{cccc}\alpha & 0& c_x& 0\\ 0& \beta & c_y& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{c}{\overrightarrow{x_c}}_{3\times 1}\\ 1\end{array}\right]$$从公式内可以看到，$\mathbf{z}$ 与 ${\mathbf{x}}_{\mathbf{c}}$ 间并不是线性关系。

这里整理整理可以得到，$\mathbf{z}$ 与 ${\mathbf{x}}_{\mathbf{c}}$ 间的关系：$\mathbf{z}=F({\mathbf{x}}_{\mathbf{c}})$，这里的 $F$ 不再是线性函数，既然不是线性函数，那么从 ${\mathbf{x}}_{\mathbf{c}}$ 转到 $\mathbf{z}$，就不再是高斯分布了。为了解决这个问题，干脆对公式在 ${\mathbf{x}}_{\mathbf{c}}$ 点处进行一阶泰勒展开，即：
$$\mathbf{z}\approx F({\mathbf{\mu }}_{\mathbf{c}}{)}_{2\times 1}+J_{2\times 3}({\mathbf{x}}_{\mathbf{c}}-{\mathbf{\mu }}_{\mathbf{c}})$$泰勒展开以后，$F({\mathbf{\mu }}_{\mathbf{c}})$ 是一个 $2\times 1$ 向量，$J$ 是一个确定值，因为在 ${\mathbf{\mu }}_{\mathbf{c}}$ 位置进行泰勒展开以后，它的雅可比是一个确定值，${\mathbf{\mu }}_{\mathbf{c}}$ 也是一个确定值。所以 $\mathbf{z}$ 与 ${\mathbf{x}}_{\mathbf{c}}$ 在这里就是线性变化关系，其协方差矩阵为 ${\Sigma }_z=J{\Sigma }_c{J}^T=JW{\Sigma }_w{W}^T{J}^T)$；均值为 ${\mathbf{\mu }}_{\mathbf{z}}=F({\mathbf{\mu }}_{\mathbf{c}})=F(W{\mathbf{\mu }}_{\mathbf{x}}+T)$。

有了上述这些知识，就可以学习 3D Gaussian Splatting 技术了。

2. 3D Gaussian Splatting

2.1 特点

3D Gaussian Splatting 和 NeRF 一样，所做的任务也是新视图合成。它有以下特点：

1. 使用光栅化渲染方式，而非基于射线的体渲染方式
   与 NeRF 的区别在于，光栅化的渲染方式是将三维空间的一个物体投到二维图像上形成对应的颜色；而体渲染方式是从图像上的点(视线)触发，将光线上的点进行汇聚形成一个点的颜色 $C$。所以一个是正向而另一个是逆向的过程，两个的渲染方式是完全不一样的。
   

2. 使用多个 $3D$ 高斯椭球显式的表达场景
   在 NeRF 里去表达一个三维场景的信息时，使用的是三维场景的点和它的体密度值σ、颜色C。3DGS内不再用点表达，而是使用 3D 高斯组件替代了点，所以在空间中的表达是一堆堆的3D高斯。这些3D高斯的信息没有保存在神经网络里，而是存在了硬盘上，所以它是一个显式的表达场景。

3. 推理速度快、 质量高

4. 未使用神经网络

需要注意的是，3D Gaussian Splatting 与 NeRF 是完全不同的思想。不要认为 3D Gaussian Splatting 是在改进 NeRF 的某个环节，它俩的思维方式完全不一样。

2.2 流程与关键步骤

论文大体流程如下：

2.2.1 场景表达

文中是以 3D Gaussian 的方式存储信息，每一个基本单元就是一个高斯球，用一堆高斯球来表达一个场景。
每个高斯球都有以下变量：

• 中心位置 $\mathbf{p}$；
  因为是三维高斯球，所以中心位置 $\mathbf{p}$ 也是三维的。
• 以 $R$ 和 $S$ 形式表达的协方差矩阵；
  $R$ 可以用四元数，$S$ 用三个实数表达，加起来就是一个七维变量，对应这里的协方差矩阵。
• 体密度(透明度) $\alpha$；
  一维变量。
• 球谐波系数；
  这里使用的 $J=3$，即有 $16$ 个基(系数)，那么 RGB 一共有 $48$ 个系数。当然使用的阶数越高，模型就越精确，但是要求的系数也越多。

综上，一个高斯球总共有 $59$ 个系数，只要给到这 $59$ 个系数，那么这个高斯球的性质就完全确定了。

2.2.2 整体流程

1. 基于SFM得到点云初始化 3D 高斯， 每个三维点初始化为一个高斯椭球
   初始化时的输入量使用的是 COLMAP 等 SFM 方式输出的点云，这里与 NeRF 就完全不同了，NeRF 使用的 仅是 COLMAP 等输出的相机位姿，而 3DGS 中，这些点云是有用的。3DGS 将根据这些点云进行初始化：基于这些点云的位置，会在每一个位置上放置一个高斯球，系数随机。

2. 给定摄像机内、 外参数及标答图像，将椭球 splatting 到图像上
   给定相机位姿，就可以将这一个个 3D 高斯球投影到图像上了，投影方式依照先前所描述的分布公式：
   $$\begin{array}{rl}& G(x)=e^{-\frac{1}{2}(x{)}^T{\Sigma }^{-1}(x)}\\ & \Sigma =RS{S}^T{R}^T\\ & {\Sigma }^{\prime }=JW\Sigma W^T{J}^T\end{array}$$

3. 通过 $\alpha$ blending 进行光栅化渲染
   3DGS 从近到远每个球投下来以后都形成了一个图像区，那么在重叠区域就可以进行光栅化的融合了。每个点都进行融合以后就可以得到图像。

4. 与标答图像计算损失
   使用的损失函数为：$\mathcal{L}=(1-\lambda ){\mathcal{L}}_1+\lambda {\mathcal{L}}_{\text{D-SSIM}}$，其中

   • ${\mathcal{L}}_1$ 度量两像素间差异；
   • ${\mathcal{L}}_{\text{D-SSIM}}$ 度量两图像间结果差异。

   NeRF 内逐像素进行计算，即采集一个 batch 的像素，送进去训练，然后输出损失；而在 3DGS 内，每次采集一小批图，以图像为单位进行损失计算。

5. 梯度回传
   有了损失，就可以以梯度回传的方式更新 3DGS 球的属性，并能控制高斯球的克隆与分裂。

   • 上支更新 3D 高斯椭球体的属性；
     ① 可微分光栅化渲染器里面都是一些加法和乘法，肯定是可以回传的；
     ② 函数 $F$ 没有学习的必要，也没有系数需要去学习。所以需要更新的只有 $3D$ 高斯球的那 $59$ 个属性。

   • 下支实现 3D 高斯椭球体的克隆和分裂等。
     

     • 学习过程中，较大梯度的高斯椭球存在 欠重构(under-reconstruction) 和 过重构(over-reconstruction) 问题。
       梯度在传过来时没有更新任何参数，只是通过对这 $59$ 维导数的模值来确定当前高斯球，是否存在欠重构或过重构的问题，如果是就进行复制或分裂。这个步骤是不可导的。
       
       ① 欠重构区域的高斯椭球方差小， 进行复制操作；
       可以看到上图中的几何体，又是很难用一个高斯球去描述这个几何体的形状，所以就对高斯球进行克隆，克隆的操作是不可导的。克隆完再优化就成了右边的样子。
       ② 过重构区域的高斯椭球方差大， 进行分裂操作；
       图中可以看到方差大的高斯球太大了，拟合覆盖了全部形状，但有太多不属于这个几何形体的形状，这样描述是不准确的。
       ③ 每经过固定次数的迭代进行一次剔除操作， 剔除几乎透明的高斯椭球以及方差过大的高斯椭球。

2.3 算法伪代码

2.3.1 整体流程伪代码

2.3.2 光栅化伪代码

整体流程为：

1. 将图像分为 $16\times 16$ 个块， 为每个块筛选 视锥体(view frustum) 内的 3D高斯椭球；
2. 实例化高斯椭球，即为其分配索引值与键值；
3. 根据 键值(key) 中的深度信息对高斯椭球进行排序；
4. 将排好序的高斯椭球从近到远向对应块上做 Splatting；
5. 在每个块上做 α-blending。

One more thing，注意一个高斯球往下投影时，这一个高斯球投影到图像中每个点的 $\alpha$ 值是不一样的，每个点的透明度还考虑了点距离中心的距离，即前面所讲的 $G(z)$。所以，这个点的透明度由 $P^{\prime }\alpha$ 决定，也就是说在投影过程中，透明度是从中间向外衰减的，这也就是 Splatting 的过程。

3. 运行流程

Running

python train.py -s /home/yj/项目资料/51.3DGS/tandt_db/bicycle

使用自己数据时，最好与官方提供数据路径保持一致，或者使用命令行参数：

• --source_path / -s：
• --model_path / -m：
• --images / -i：

4. 运行可能错误

1. libcublas.so.11: undefined symbol: cublasLtHSHMatmulAlgoInit, version libcublasLt.so.11
   看样子是 conda 安装 PyTorch 造成的问题。
   解决方法：
   既然连不上动态库，那先找一下所有的 libcublasLt.so.11
   得到下面结果：
   
   找到需要的路径后，检查一下环境变量里有没有对应的 LD_LIBRARY_PATH：
   
   显然是没有的，所以我们手动添加进去，问题解决。

   export LD_LIBRARY_PATH=/home/yj/anaconda3/envs/gaussian_splatting/lib/:$LD_LIBRARY_PATH

2. 可视化时遇到：File not found: uv_mesh.vert
   解决方法：
   在路径 ./SIBR_viewers/install/shaders/core 下运行 SIBR_viewers 即可。