三维重建【1】3D Gaussian Splatting：实时的神经场渲染

3D Gaussian Splatting Paper Reading

论文相关：

• 论文链接：3D Gaussian Splatting for Real-Time Radiance Field Rendering
• 官方代码：https://github.com/graphdeco-inria/gaussian-splatting
• 中文摘要：辐射场方法最近彻底改变了使用多张照片或视频捕捉的场景的新视角合成技术。然而，要实现高视觉质量仍然需要成本高昂的神经网络来训练和渲染，而近期较快的方法不可避免地在速度和质量之间做出妥协。对于未受限且完整的场景（而不是孤立对象）和1080p分辨率渲染，目前还没有方法能实现实时显示率。我们引入了三个关键元素，使我们能够在保持竞争力的训练时间的同时，实现行业领先的视觉质量，并且重要的是，支持高质量的实时（>= 30 fps）新视角合成，分辨率为1080p。首先，从相机校准过程中产生的稀疏点出发，我们用3D高斯分布来表示场景，它保留了连续体积辐射场对场景优化的有益属性，同时避免了在空旷空间中的不必要计算；其次，我们进行3D高斯的交错优化/密度控制，特别是优化各向异性协方差以实现场景的精确表现；第三，我们开发了一种快速的可感知可见性的渲染算法，支持各向异性溅射，加速了训练并允许实时渲染。我们在几个已建立的数据集上展示了行业领先的视觉质量和实时渲染。

Motivations

1、要实现高视觉质量，仍然需要训练和渲染成本高昂的神经网络；
2、目前没有任何一种方法能够达到实时显示速率。

Overview of 3D Gaussian Splatting

本文作了以下三个工作：

• 第一点：引入了一种各向异性的3D Gaussian分布来作为高质量非结构化的辐射场表达，对应到上面的示意图就是从SfM点云出发，以每个点为中心生成3D Gaussian分布，这里的各向异性意思是从各个方向上看过去它都长得不一样，也就是说把一个点往不同相机位姿上投影的时候会投出不一样的样子；
• 第二点：实现了使用GPU进行快速可微的渲染，允许各向异性的Splatting和快速反向传播，这里的关键词是Splatting，它是一个比较经典的计算机图形学里面用三维点进行渲染的方法，它把三维点视作雪球往图像平面上抛，雪球在图像平面上会留下扩散的痕迹，这些点的扩散痕迹叠加在一起就构成了最后的图像，它是一种很方便的针对点云的渲染方法;
• 第三点：提出了针对3D Gaussian的优化方法，并同时进行自适应的密度控制，其中存储在点里的3D Gaussian参数需要在反向传播的时候进行优化更新，本文还会根据梯度去自适应地调整点云的分布，如果一个点的3D高斯太大，不能完全拟合这一处的细节，那就对它进行分割，用两个3D Gaussian去表示，如果一个位置的点太密集，用不了那么多的3D Gaussian，那就吞并为一个3D Gaussian。

接下来是整个系统的示意图，从左边开始首先对SfM的点云进行初始化得到3D Gaussian，然后沿黑色箭头向右，也就是借助camera的外参做投影projection，然后再接着用一个可微的光栅化渲染得到图像，有了渲染图像之后，将其与真实图像对比求Loss，再沿着蓝色箭头反向传播，蓝色箭头向上更新3D Gaussian里的参数，向下送入自适应密度控制来更新点云。

3D Gaussian

这是在别处看到的对比明显的一张图，左边这张是由方块组成的世界，简称“我的世界”；中间这张图是我们比较常见的网格面，它是由三角形所组成的世界，也就是一堆三角形拼拼拼，拼成一个三维的物体，也就是我们常见的网格面；而3D Gaussian用的是椭球所组成的世界，它是由一堆椭球，然后拼接拼接拼接，把他拼接成一个大的场景，也就是最右边这张图。
一维的概率密度分布函数大家应该比较熟悉，即：

当变量的维度从一维扩展到三维，也就是说现在的变量是$v=[a,b,c]$，包含三个变量，假设这三个变量都服从均值为0、方差为1的正态分布，且这三个变量相互独立，那么就可以得到$p(v)=p(a)p(b)p(c)$，那么指数部分就会是$a^2+b^2+c^2$，如果用变量$v$来表示这三个向量的话，那么等式中的指数部分可以表示为$-\frac{1}{2}v{v}^T$。如果转换到一般情况，也就是说再给定一个$3\times 1$的向量定义为$x$，向量$x$包含变量$x$,$y$,$z$，这三个变量实际上可能并不是独立分布，那么做一个线性变换，将均值不一定为0、方差不一定为1的变换至均值为$\mu$、方差为1的$a$,$b$,$c$。
用矩阵形式表示即$v=A(x-\mu )$，它的均值向量即为各个变量的均值，它们各自减去均值的话就是将数据进行中心化，然后用一个矩阵$A$对向量$x$,$y$,$z$进行线性组合，让它们等于$[a,b,c]$，这个矩阵$A$称之为变换矩阵，变换之后再将其代入$p(v)$的公式中，从而得到公式$p(v)==1/(2\pi {)}^3/2exp(-1/2(x-\mu {)}^T{A}^{T}A(x-\mu ))$，此时概率密度函数仍然是$p(v)$。

为了得到$p(x)$，等式两边积分，然后对$dv$进行求导化简，得到关于$dx$的等式，则可以得出关于向量$x$的概率密度函数，但此时的式子并不完整，因为一维的高斯用均值和方差表示，多维的高斯也应该是用均值向量和协方差矩阵来表示：

接下来就是将变换矩阵$A$转换为协方差矩阵，协方差矩阵作为对称矩阵可以进行特征值分解，
也就是说如果有一个协方差矩阵Σ，则可以将其分解为$U\Lambda U^T$的转置，其中矩阵$U$的每一列都是相互正交的特征向量，而且是单位向量，满足$U{U}^T$等于单位矩阵，$\Lambda$矩阵对角线上的元素是协方差矩阵的特征值，其余元素都是零.
将协方差等式进一步进行拆分，可以看到括号中即为一组正交基乘以一个对角阵，这正好表示了一个线性变换，换言之就是变换矩阵$A=U{\Lambda }^{\frac{1}{2}}$，因此协方差矩阵$\Sigma =A{A}^T$，从而可以得到完整的3D Gaussian的概率密度分布函数
这个公式从形式上看与一维的高斯分布差不多，在一维的情况下它是一条中间高两边低的曲线，二维的时候是一个中间凸起的山峰，三维的时候会变成一个椭球，椭球的形状由协方差矩阵$\Sigma$来决定

Differentiable 3D Gaussian Splatting

如图这团点云即是3D Gaussian点云，它的每一个点里面存储了这样一些数据，

• 第一个是点的位置，这是一个默认信息，也就是坐标x,y,z，对于高斯来说它同时也是这个点的高斯分布的均值向量；
• 第二个是协方差矩阵，它用于决定当前点处的绿色椭球的形状和方向；
• 第三个是不透明度，用于渲染，做Splatting的时候，用于把点往图像平面上投影的时候，它们的扩散痕迹是通过这个不透明度叠加在一起的；
• 第四个是球谐函数，用来拟合视角相关的外观，原文中的说法是directional appearance component(color) of the radiance field

至于为什么选择3D Gaussian，本文作者的说法是：他们需要一种图元能够在拥有场景表达能力的时候可微，而且显式地支持快速渲染，于是就有了这样的3D Gaussian点云，它里面的参数可以在迭代优化的过程中更新，而且能够很容易的用splat的方法来投影到2D图像上，做非常快的$\alpha$混合，也就是渲染。

本文是通过协方差矩阵来表达这个点，这个矩阵控制着3D Gaussian的形状，文章中所用的3D Gaussian的公式如上图所示，相较于完整地我们推导的公式，它去掉了均值，因为这个高斯分布以点为中心，x,y,z的均值就在点上，已经中心化了，所以把均值设为0，而原式前面的系数是为了控制整个概率密度函数的积分为1，但是在这里是不需要的，从而使得整个分布的大小可以控制，想多大有多大，没有积分为1的限制。

前面说这个矩阵能够控制3D Gaussian的形状，这里需要解释一下协方差矩阵的物理含义，前面我们推导的时候写到协方差矩阵$\Sigma =A{A}^T$，我们只说这个$A$是线性变换，将任意一个分布变换为均值为$\mu$、方差为1的范围内，在二维空间里来看的话相当于先在$y$方向上进行压缩，然后再绕原点旋转一定的角度，因此矩阵$A$的构造往往可以看作一个旋转矩阵$R$和一个尺度变换矩阵$S$的乘积，其中尺度变换矩阵分别在对角线上有两个缩放系数，这是二维的情况，以此类推三维也是如此。因此，在文章里会有这样一个公式，$\Sigma =RS{S}^T{R}^T$，在迭代优化的过程中通过优化矩阵$A$的数据就能改变高斯椭球的形状、大小和方向等几何外观，从而能够使其在Splatting的时候投影出正确的效果。

这里涉及到矩阵优化的问题，直接采用$3\times 3$的矩阵进行优化会导致矩阵计算多且复杂，因此本文是通过旋转四元数，将旋转矩阵内的9个参数压缩到4个参数；而尺度变换矩阵也不需要记录整个矩阵的信息，只需要记录其在三个轴方向上的缩放比即可。因此，协方差矩阵参数的更新转变为旋转四元数和一个含缩放比信息的三维向量的更新.

（注：旋转四元数是一种用于表示三维空间中旋转的数学工具，它是四元数的一种特殊表现形式，由一个实部和三个虚部组成，旋转四元数的核心是通过旋转轴上的旋转角度进行编码，以及通过旋转轴的单位向量来表示旋转的方向。旋转四元数的实部用于表示旋转角度的余弦值，而虚部用于表示旋转轴在单位向量上的三个分量。）

说完协方差矩阵还有不透明度和球谐函数：
对于不透明度只包含一个参数$\alpha$，3D Gaussian其概率密度越高的地方越不透明，概率密度越低则越透明，由于高斯分布的概率密度在靠近中心的位置最高，因此每个3D Gaussian最中心是不透明度最高的地方，可以参考上图共有3个3D Gaussian，每个3D Gaussian在中心位置色彩最亮，而到边缘位置会逐渐发黑，这正是由于3D Gaussian在边缘位置的不透明度较低，露出了黑色的背景。
对于球谐函数系数，可以将其理解为用于表示颜色的一组数值。通常我们在表示颜色时使用RGB三原色，但采用球谐函数系数来表示颜色时，每个3D Gaussian多达48个数值，以至于该系数不止能够表达单一的颜色，而且能够表达物体表面不同位置的不同色彩和纹理。类比泰勒展开式、傅里叶级数，可以将球谐函数理解为多项式的拟合，都是通过一组不同阶的基函数线性组合而成，只不过傅里叶级数的基函数是三角函数，而球谐函数的基函数是球函数，同样地球谐函数的阶数越高其表达能力越强，可以参考右图随着阶数的升高，拟合的形状会越来越接近它原本的形状。本文将RGB三个颜色分别用三个不同的球谐函数来记录其值，使得点云在不同角度呈现不同颜色，原文是采用4阶，就是包含0、1、2、3共4阶，参数量为16，一共是R\G\B三个原色，故一共采用48个参数去存储RGB颜色参数，也就是图中的这一部分，然后对其进行加权求和，从而实现对颜色的表达，将离散的颜色数据转换为连续的颜色数据，以便于对参数进行迭代。

说完论文中的参数后，接下来便是Splatting
3D Gaussian是3D物体，要生成图像就需要将其投影到2D平面上，文章中给出的方法是替换协方差矩阵为${\Sigma }^{’}$，具体的实现过程包括：
1、确定2D Gaussian的位置：通过3D Gaussian的位置向量以及投影矩阵，从而得到2D Gaussian的位置
2、确定2D Gaussian的协方差矩阵：由于协方差矩阵描述了一个3D Gaussian的分散程度，它描述的不是一个点，而是一种形状，在非线性变换的过程中，协方差矩阵可能会被拉伸，这样的拉伸会导致3D Gaussian投影到2D平面后不再是一个高斯，因此文章中采用雅可比矩阵将这种非线性变换转换为线性变换，
3、根据3D Gaussian的球谐系数和不透明度来计算2D Gaussian的颜色：其中$c_n$是球谐函数输出的颜色，而不透明度${\alpha }_n’$是不透明度${\alpha }_n$与3D Gaussian分布的乘积结果。

Optimization with Adaptive Density Control of 3D Gaussians

文章中用到的一些优化方法包括：
1、优化策略是随机梯度下降（SGD）
2、自己写了一些CUDA的核心用来加速
3、用快速光栅化来提高效率，这部分内容后面再具体解释
4、不透明度α用sigmoid函数将其限制在[0,1)之内，协方差矩阵的尺度用的是指数激活函数
5、最后它所采用的损失函数也是比较常见的L1 Loss，就是用渲染图像与真实图像求光度误差，然后再按一定的比例$\lambda$，加上一个SSIM结构相似性的误差。

除去优化外，文章中还提出了点云密度的自适应控制，即Adaptive Density Control of 3D Gaussians，这使得系统能够从稀疏的、质量不那么高的初始点云，甚至随机初始化的点云里面去拟合出比较好的模型：
首先是每隔一定的迭代次数，移除不透明度小于阈值的点；然后是对重建不充分的区域进行处理，这里的重建不充分包括两类，即under-reconstruction and over-reconstruction，一个欠重建一个过重建，可以看上面的右图，假设我们要拟合这个有点像草履虫的图案，在重建的初期，3D Gaussian都被初始化得比较小，无法表达想要的形状；然后，我们对这个3D Gaussian进行复制，将原本的一个高斯复制成两个，两个高斯同时训练，能加快训练速度，同时也能拟合出单个高斯难以拟合出的形状；到了训练后期，高斯的尺寸可能会过大，导致无法精细地拟合出这个形状，这时我们可以将这个高斯拆分成两个较小的高斯，它们在经过训练后能够更加贴近我们想要的形状。文章中判断重建是否充分的依据是梯度，在更新点的位置的时候，如果梯度太大说明这个位置的误差比较大，需要修正的量比较大，梯度超过阈值则会进行密度控制的操作；此外，在靠近相机的位置可能会学习到一些floaters，即不属于模型的错误重建，因此本文周期性地将不透明度重置为0，在后续迭代优化的过程中真正需要用到的点的不透明度将会回到正常的值，一些不会被用到的、不透明度低的点将会被剔除。

Fast Differentiable Rasterizer for Gaussians

1. 把整个图像划分为16×16的tiles，每个tile视锥内挑选可视的3D Gaussian
2. 每个视锥内只取置信度大于99%的高斯，并按深度排序
3. 并行地在每个tile上splat
4. 有像素的不透明度达到饱和就停止对应线程
5. 反向传播误差时按照tile对高斯进行索引
   本文的渲染方法称之为Tile-based Rasterizer，它采用tile来代替像素级别的精度，首先将整个要渲染的图像分为16×16的tiles，然后为每个图块挑选视锥内的3D Gaussian，具体地说就是只保留与视锥体相交的置信度区间为99%的高斯分布；然后，文章根据重叠的tile数量实例化每个Gaussian值，并为每个实例分配一个结合视图空间深度和tile ID的键，然后根据这些键快速排序（基于快速GPU Radix）；在排序之后，按这些Gaussian到图像平面的深度值的排序顺序从近到远往对应的tile上去做splat，把splat留下的痕迹做堆叠累积，直到不透明度饱和；每个tile都单独开了一个线程块，所以可以认为所有tile上的光栅化是并行运行的，这里的光栅化指的就是从堆叠的splat痕迹中去划分像素网格来生成像素值。

至此便是3D Gaussian Splatting的全貌了，本人也在学习阶段，有问题欢迎大家提出，一同探讨。下一篇准备写3D Gaussian Splatting的官方代码解读。