3DGS 的原论文中，每个高斯是不存法线信息的，然而后续工作中，也有研究人员将法线和 BRDF 信息进行训练，这样一来，就可以在后续的光栅化过程中输出 GBuffer，完成 relighting 的任务。生成初始法线，可以仅考虑高斯的均值，基于局部平面假设进行估计。

0x00 To begin with

这篇文章将会包含以下内容：

  • 什么是 PLY
  • 利用最小二乘法估计点云法线
  • 使用 PCL 的实现
  • 在 CloudCompare 中实现

For reference👇：

0x01 PLY 格式

PLY（Polygon File Format, also known as the Stanford Triangle Format），是一种用于存储多边形集合的图形对象的格式。有两种子格式：

  • ASCII：ASCII 编码表示，直观易读
  • Binary：二进制格式，存储空间小，能够快速保存和加载

PLY 文件的经典结构如下：

Header  // description field, containing element type, property and count
Vertex List  // list of data for each vertex
Face List  // list of data for each face
(List of other elements)

Header 描述文件格式，评论（comment）以及每个元素的描述，并在其后指定该元素属性。例如一个以 ASCII 子格式保存的 PLY 文件，其 Header 如下：

这里的数据不正常是因为中途改了原本二进制格式的编码，似乎导致了解析错误，后文所用的文件是正常的

在该 PLY 文件中，存储了 24077 个顶点的信息，每个顶点由三个浮点分量 x，y，z 描述位置。在 Header 结束（end_header）后，即列出 Vertex List 的内容。

此外，PLY 格式允许用户自定义元素内容，例如 3DGS 就专门设计了一个 PLY 格式。

可以在一些 DCC 软件原生地（如 Blender）或借助插件导入 PLY 文件，也可以使用专业分析点云的软件（如 CloudCompare）。

0x02 从点云估计法线

在传统的图形学中，每个顶点的法向由与该顶点邻接的面法向按某种形式加权计算得出，因此需要存储 Mesh 的文件至少提供每个顶点的位置数据以及组成各面的顶点索引数据。

对于点云数据，由于其只是在原始表面上离散采样的一系列点，不存在表面信息，因此无法直接通过由面到点的方式求解顶点法线。然而，可以考虑如下两种思路：

  • 利用表面网格化技术从点云中构建表面，然后从网格中计算顶点法线
  • 使用近似方法直接从点云中估计法线

第一种方法通过某种表面网格化技术（例如 Marching Cubes）从点云转换的 SDF 中重建表面，然后再使用传统图形学的做法计算顶点法线；第二种方法直接从点云中估计法线，其基于一个 local planar assumption：在足够小的一个局部区域内，点云的原始表面可以被一个平面（或者其它结构）近似。因此，确定表面上某一点的法线的问题就转化为估计与表面相切的平面的法线的问题。

本篇文章中，将实现第二种方法。

0x03 最小二乘法估计法线

在高中物理，经常使用最小二乘法用一条直线拟合实验数据，实际上，最小二乘法的思想可以扩展至更高维度，以拟合超平面或更复杂的几何结构

要拟合一个平面，本质上是一个三元线性回归问题，相应地可以转化为最小二乘法的平面估计问题（least-square plane fitting estimation）。

基本思路

对于给定点云的局部集合 $\mathcal{P}^k$（$k$ 为集合中点的个数），其潜在近似描述的平面可以表示为点 $x$ 和面法向 $\vec{n}$，对于集合中任意点 $p_{i}\in\mathcal{P}^k$，其到平面的距离为 $d_{i}=(p_{i}-x)\cdot \vec{n}$。最小二乘法的目标是求解 $x,n$，以使所有点到这个平面的距离平方和最小。令 $x$ 为集合点的质心：

$$ x=\frac{1}{k}\cdot \sum_{i=1}^kp_{i} $$

求解 $\vec{n}$ 的过程可以通过分析 $\mathcal{P}^k$ 的协方差矩阵 $C\in \mathbb{R}^{3\times 3}$：

$$ C=\frac{1}{k}\sum_{i=1}^{k}(p_{i}-\bar{p})(p_{i}-\bar{p})^T,~C\cdot \vec{v_{j}}=\lambda_{j}\cdot \vec{v_{j}},~j\in\{0,1,2\} $$

其中，$\bar{p}=\frac{1}{k}\sum_{j=1}^{k}p_{ij}$ ，$\lambda_{j}\in \mathbb{R}$ 为其特征值。定义去中心化的点坐标 $p_{i}'=p_{i}-\bar{p}=(x_{i}',y_{i}',z_{i}')$，将 $C$ 展开：

$$ C = \frac{1}{k} \sum_{i=1}^{k} \begin{bmatrix} x_i'^2 & x_i' y_i' & x_i' z_i' \\ x_i' y_i' & y_i'^2 & y_i' z_i' \\ x_i' z_i' & y_i' z_i' & z_i'^2 \end{bmatrix} =\begin{bmatrix} \sigma_{xx} & \sigma_{xy} & \sigma_{xz} \\ \sigma_{yx} & \sigma_{yy} & \sigma_{yz} \\ \sigma_{zx} & \sigma_{zy} & \sigma_{zz} \\ \end{bmatrix} $$

可以观察到 $C$ 是一个实对称半正定矩阵，主对角线元素代表集合 $\mathcal{P}^k$ 在 $x,y,z$ 方向上的方差，非对角元素表示不同坐标轴之间的协方差。根据谱定理，实对称矩阵 $C$ 可被分解为：

$$ C=Q\Lambda Q^{T}=Q\Lambda Q^{-1},~Q=[\vec{v_{0}},\vec{v_{1}},\vec{v_{2}}],~\Lambda=\mathrm{diag}(\lambda_{0}, \lambda_{1}, \lambda_{2}) $$

实对称矩阵，其逆等于其转置

可以将其视为某个二次型的矩阵，定义二次型函数 $Q(x)$：

$$ Q(x)=x^{T}C^{-1}x $$

对 $C$ 做特征分解：

$$ Q(x)=x^{T}C^{-1}x=(x^T Q)\Lambda^{-1} (Q^T x) $$

记 $y=Q^T x$，$y$ 为切线空间下的坐标，原二次型转换为：

假定特征向量均已经过归一化

$$ Q(x)=y^{T}\Lambda^{-1}y=\sum_{j=1}^{n} \frac{y_{j}^2}{\lambda_{j}}=\frac{y_{0}^{2}}{\lambda_{0}}+\frac{y_{1}^{2}}{\lambda_{1}}+\frac{y_{2}^{2}}{\lambda_{2}}=1 $$

其中，$y_{j}$ 为 $x$ 在协方差矩阵第 $j$ 个特征向量方向的投影。所以在几何角度上，可以把协方差矩阵 $C$ 视为定义了一个拟合点集的椭球体。

$C$ 的特征向量 $\vec{v_{j}}$ 构成一个三维正交基（切线空间基），代表集合 $\mathcal{P}^k$ 的三个主成分（Principle Components）。如果 $C$ 表示的椭球接近一个平面，则集合 $\mathcal{P}^k$ 在垂直于平面的方向上变化最小，即该方向的方差最小、特征值最小。考虑 $0\leq \lambda_{0}\leq \lambda_{1}\leq \lambda_{2}$，最小的特征值 $\lambda_{0}$ 对应的特征向量 $\vec{v_{0}}$ 即为拟合平面的法线。

K-Neighborhood Selection

如上文所述，计算点云中某点的法线需要选取与该点距离最近的 $k$ 个点（K-Neighborhood）进行估计。因此，$k$ 值的选取对于估计法线的质量有着极大影响：

  • 较大的 $k$ 尺度会导致邻域点集较大，从而引入额外噪声，导致估计的点法向失真
  • 较小的 $k$ 尺度可以捕捉局部细节，但容易受局部非均匀采样的影响，导致大范围内法线方向不稳定，并且增大开销

左：较小的 k，右：较大的 k，Point Cloud Library

个人感觉是采样频率的问题，高频相对低频能够更好地捕捉细节，但也更容易受到噪声的干扰

一般地，$k$ 值的选取可以让用户凭经验选择阈值，然而固定的阈值无法应对一些复杂（频率、空间分布不均匀）的点云，一种更好地方式是自适应地选取 $k$ 值。根据当前估计点的场景及其所需细节程度选取阈值，例如，在一些高频区域（潜在表面的转角、连接处）使用较小的 $k$，而在其它情况下可以使用较大的 $k$。判断给定点的所需细节程度可以通过分析其密度、曲率等因素。

翻转法线

PCA 方法估计的法线方向如下图左、中所示，在同一个平面上出现了向内和向外两种方向的法线；下图右则为扩展高斯图像（Extended Gaussian Image），描述了点云所有法线方向，由于数据集是 2.5D 的（即从单一视点下采集），因此法线只应分布在一个半球范围内，而实际上它们分布在整个球上。

左中：使用 PCA 得到的朝向不连续的法线，右：结果的 EGI，Point Cloud Library

为了使所有法线都以各自不同的程度朝向视点 $v_{p}$，对于点 $p_{i}$ 及其法线 $\vec{n_{i}}$，应满足：

$$ \vec{n_{i}}\cdot(v_{p}-p_{i})>0 $$

对于不符合上式的法线，将其取反，即可得到正确朝向（朝向视点）的法线。下图给出了经该方法重定向后的点云法线：

左、中：经过重定向的法线，右：重定向后的法线 EGI，Point Cloud Library

该方法仅在已知视点 $v_{p}$，且一个视图仅在该视点下采集时成立。当不知道视图的视点或者视图存在多个采集视点，则上述方法不成立。

综上所述，对于每个点 $p_{i}$，估计其法线的步骤如下：

  1. 获得与其距离最近的 $k$ 个邻居
  2. 计算 $p_{i}$ 的法线
  3. 检查计算出的法线是否朝向视点，否则翻转

0x03 利用 PCL 实现

PCL（Point Cloud Library）是一个开源的点云处理库，专门用于处理三维点云数据。PCL 集成了包括点云文件的读写、特征提取、分割配准、可视化等模块，可以方便地处理、分析点云数据。此次从点云估计法线的任务可以借助 PCL 完成，我的配置信息为：

  • 平台：Windows 11
  • PCL 版本：PCL 1.15.0
  • IDE：VisualStudio 2022

PCL 可在 Github 官方仓库下载：

Releases · PointCloudLibrary/pcl

配置过程可参考：

PCL点云库安装及学习（2021.7.28）

Trouble Shooting

配置过程中遇到了一些问题，这里一并提一下：

  1. error: PCL requires C++14 or above
  2. error: Potential runtime error due to aligned malloc mismatch! You likely have to compile your code with AVX enabled or define EIGEN_MAX_ALIGN_BYTES=32
  3. \[pcl::PCDReader::read\]Could not find file (乱码)
  4. Assertion failed: index > 0 && index < count, file C:\Program Files\flann\include\flann/algorithms/kdtree_single_index.h, line 492

第一个问题默认的 C++ 版本太低，在项目–>属性–>C/C++–>语言–>C++语言标准中选择更高版本（17 或 20）即可。

第二个问题，内存对齐问题，可以通过改变指令集解决，在项目–>属性–>C/C++–>代码生成–>启用增强指令集中选择高级矢量扩展 (X86/X64) (/arch:AVX) 即可。

第三个问题，下载的 PCL 在根目录以及第三方依赖的 lib 目录下会提供 Debug 和 Release 两种版本的库，其中 PCL 的 lib 文件中 Debug 版本的库在文件名末尾会以‘d’标识，而第三方依赖的库则会以‘-gd-’标识。链接器链接动态库时把 Debug 和 Release 版本的混淆了，所以在项目–>属性–>链接器–>输入–>附加依赖项中只填入对应版本的库即可。

第四个问题，表面上是断言 k-d tree 建树错误，然而检查了一下原始数据，发现是文件错了，数据包含了 e+23 和 e-20 这种极端值。重新下载原始文件就正常了，推测是中途改了一遍编码格式导致的。

Code

 1#include <pcl/io/ply_io.h>
 2#include <pcl/point_types.h>
 3#include <pcl/filters/filter.h>
 4#include <pcl/features/normal_3d.h>
 5#include <pcl/visualization/pcl_visualizer.h>
 6
 7int main() {
 8    // read ply file
 9    pcl::PointCloud<pcl::PointXYZ>::Ptr cloud(new pcl::PointCloud<pcl::PointXYZ>);
10    if (pcl::io::loadPLYFile<pcl::PointXYZ>("Test.ply", *cloud) != 0) {
11        PCL_ERROR("Failed to load point cloud file!");
12        return -1;
13    }
14    //[DEBUG]: check data
15    //std::cout << "Loaded cloud size: " << cloud->size() << std::endl;
16    //for (size_t i = 0; i < std::min<size_t>(10, cloud->size()); ++i) {
17    //    auto& pt = cloud->points[i];
18    //    std::cout << i << ": (" << pt.x << ", " << pt.y << ", " << pt.z << ")\n";
19    //}
20
21    // filter invalid points
22    pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_filtered(new pcl::PointCloud<pcl::PointXYZ>);
23    std::vector<int> nan_indices;
24    pcl::removeNaNFromPointCloud(*cloud, *cloud_filtered, nan_indices);
25    if (cloud->empty()) {
26        PCL_ERROR("Error: Point cloud is empty after removing NaNs!");
27        return -1;
28    }
29
30    // estimate normal
31    pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> ne;
32    ne.setInputCloud(cloud_filtered);
33    ne.setViewPoint(100.0, 0.0, 0.0);
34    pcl::search::KdTree<pcl::PointXYZ>::Ptr tree(new pcl::search::KdTree<pcl::PointXYZ>());
35    ne.setSearchMethod(tree);
36    ne.setKSearch(100);
37
38    pcl::PointCloud<pcl::Normal>::Ptr normals(new pcl::PointCloud<pcl::Normal>);
39    ne.compute(*normals);
40
41    // vis
42    pcl::visualization::PCLVisualizer viewer("Normals");
43    viewer.addPointCloud(cloud_filtered, "cloud");
44    viewer.setPointCloudRenderingProperties(
45        pcl::visualization::PCL_VISUALIZER_COLOR,
46        0, 1, 0,
47        "cloud"
48    );
49    viewer.addPointCloudNormals<pcl::PointXYZ, pcl::Normal>(cloud_filtered, normals, 10, 0.05, "normals");
50    
51    while (!viewer.wasStopped()) {
52        viewer.spinOnce();
53    }
54}

运行程序，不出意外就能看到法线：

viewpoint=(0,0,0), k=20

viewpoint=(0,0,0), k=100

viewpoint=(100,0,0), k=100

可以观察到由于选取的点云为平面，较大的 $k$ 值得到的法线在整个平面上朝向更加均匀，且在转角处有更为平滑的过渡。

0x04 在 CloudCompare 中实现

然后按照指示选择合适的 Local surface model、Neighbors、Orientation，即可计算法线，最后选中点云，在 Properties 栏中对法线进行可视化：