测量3D点云数据中的几何体尺寸

本报告记录了本人使用3D数据常用的处理算法，对用3D点云数据代表的几何体进行尺寸测量。

1 引言

1.1 3D点云

3D点云数据是由大量的点组成的集合，每个点在空间中都有一个三维坐标 $(X, Y, Z)$ 。这些点通常用于描述物体的表面或环境的形状和结构。点云数据广泛应用于计算机视觉、三维建模、地图测绘、建筑设计等领域。

1.2 要进行几何尺寸测量的物体

一共有7个点云测量任务，包含9个不同的物体。每个任务都有来自不同视角拍摄的4个点云数据集。原始数据集格式为.npz，为在MATLAB上操作，已将它们全都转换为.mat文件。

在该工作中值得注意的是：

• 第5个任务中要被测量的两个物体叠放在一起。
• 第6个任务中的点云数据集中同时存在两个物体。
• 第7个任务中要被测量的物体是手持物体拍摄的。

如下图：

2 算法介绍

本节将对在该工作中所使用的算法进行介绍。

2.1 RANSAC

随机抽样一致（RANdom SAmple Consensus，RANSAC）算法。它采用迭代的方式从一组包含离群的被观测数据中估算出数学模型的参数。 RANSAC是一个非确定性算法，在某种意义上说，它会产生一个在一定概率下合理的结果，而更多次的迭代会使这一概率增加。此RANSAC算法在1981年由Fischler和Bolles首次提出。

RANSAC的基本假设是：

1. “内群”(inlier, 似乎译为内点群更加妥当，即正常数据，正确数据)数据可以通过几组模型的参数来叙述其分布，而“离群”(outlier,似乎译为外点群更加妥当，异常数​​据)数据则是不适合模型化的数据。
2. 数据会受噪声影响，噪声指的是离群，例如从极端的噪声或错误解释有关数据的测量或不正确的假设。
3. RANSAC假定，给定一组（通常很小）的内点群，存在一个程序，这个程序可以估算最佳解释或最适用于这一数据模型的参数。

2.1.1 范例

这里用一个简单的例子来说明，在一组数据点中找到一条最适合的线。假设，此有一组集合包含了内群点以及离群点，其中内群点包含可以被拟合到线段上的点，而离群点则是无法被拟合的点。如果我们用简单的最小二乘（Least Square，LS）法来寻找此线，我们将无法得到一条适合于内点群的直线，因为最小二乘法会受离群点影响而影响其结果。而用RANSAC，可以只由内群点来计算出模型，而且概率还够高。然而，RANSAC无法保证结果一定最好，所以必须小心选择参数，使其能有足够的概率得到正确的拟合直线。

{height=500}

2.1.2 算法概述

算法的具体流程如下：

1. 在数据中随机选择若干个点设定为内点群
2. 计算拟合内点群的模型
3. 把其它刚才没选到的点带入刚才建立的模型中，根据一定的评价标准，计算该点是否属于内点群
4. 记下内点群数量
5. 重复以上步骤
6. 比较哪次计算中内点群数量最多，内点群最多的那次所建的模型就是我们所要求的解

这里有两个问题：

1. 一开始的时候我们要随机选择多少点（ $n$ ）
2. 以及要重复做多少次（ $k$ ）

2.1.3 参数决定

假设每个点是真正内点群的几率是 $w$ ，则： $$w = 真正内点群的数目 / 数据总数$$ 通常我们不知道 $w$ 是多少， $w^𝑛$ 是所选择的 $𝑛$ 个点都是内点群的几率， $1−w^n$ 是所选择的 $n$ 个点至少有一个不是内点群的几率， $(1-w^{n})^{k}$ 是表示重复 $k$ 次都没有全部的 $n$ 个点都是内点群的几率，假设算法跑 $k$ 次以后成功的几率是 $p$ ，那么：

$$ \begin{array}{l} 1-p=\left(1-w^{n}\right)^{k} \\ p=1-\left(1-w^{n}\right)^{k} \\ k=\frac{\log (1-p)}{\log \left(1-w^{n}\right)} \end{array} $$

2.1.4 实现

MATLAB已有对RANSAC算法的封装好的实现，即：

[model,inlierIdx] = ransac(data,fitFcn,distFcn,sampleSize,maxDistance)

此外，针对各种点云数据拟合的问题，MATLAB也提供了相应的函数，且使用了RANSAC算法的优化变种：MSAC算法。即：

% 拟合平面
model = pcfitplane(ptCloudIn,maxDistance)
% 拟合球面
model = pcfitsphere(ptCloudIn,maxDistance)

2.2 PCA降维

2.2.1 主要概念

主成分分析（Principal Component Analysis, PCA）是一种广泛使用的降维技术，主要用于简化数据集的维数，同时保留数据中最重要的变化信息。从几何角度上看，PCA是找到一组正交的基向量，使得将数据旋转到新的基向量构成的空间后（这也意味着在不丢失维度信息的前提下，PCA变换是一种保距变换），能集中分布在少数几个基向量（对应着数据方差最大的方向）附近，此时就可以对数据进行降维处理。如下图：

{height=300}

2.2.2 PCA用于3D点云降维

在点云数据处理中，经常会有这样的问题：点云数据原本分布在同一平面上，但受噪声的影响，它们偏离了平面，如何找回点云数据原本所在的斜面，并给出点云数据在原始二维平面上的分布？这样的问题可以通过RANSAC方法去找到一个拟合平面来解决，这里介绍PCA方法在处理该问题上的能力和优势。 可以看到，上述问题实际上包含了两个要求：

1. 去噪。
2. 对数据进行降维。

由于点云数据原本分布在一个平面上，点云数据在该二维平面上位置的变化构成了数据方差的主要来源，又因为PCA能识别出数据的方差最大的方向，而噪声带来的方差比较小，且其方向无规律可言。因此，通过PCA找到数据的前两个主成分对应的特征子空间，就是点云原本分布的平面。然后，把点云数据投影到该平面上，就可以解决上文中的问题。而这一切，都仅仅是通过几个简单的矩阵运算来完成的。可以看到，这种方法相对于RANSAC方法而言，运算上更为简便和快捷，且能直接得到3D点云的二维分布，即，只有两个维度的点云数据，这极大地方便了后续的处理。

2.2.3 实现

以下给出PCA降维的MATLAB实现代码作为参考：

function [projected_points] = myPCA(points)
    % 中心化数据
    centered_points = points - mean(points);
    % 计算协方差矩阵
    cov_matrix = cov(centered_points);
    % 求解特征值和特征向量
    [V, D] = eig(cov_matrix);
    % 按特征值大小排序
    [~, idx] = sort(diag(D), 'descend');
    V = V(:, idx);
    % 将数据投影至前两个最大的特征值构成的特征空间
    projected_points = centered_points * V(:, 1:2);
end

2.3 最小圆覆盖算法

2.3.1 点云数据处理中的最小圆覆盖问题

在对点云的圆平面进行检测时，往往会遇到圆平面的边界难以检测的问题，如下图：

{height=300}

我们显然对这个点云数据集代表的是一个圆平面这个猜测十分有把握，但是由于部分点云没有被检测到或被处理掉，导致圆的边界有所缺失，此时，我们很难通过精确的数学方程来解决这个问题。

2.3.2 最小圆覆盖算法

上述问题可以抽象为：给定 $n$ 个点的平面坐标，求一个半径最小的圆，把 $n$ 个点全部包围。常用算法是 Welzl's algorithm。这里仅给出该方法的伪代码（如果从递归的角度去思考，把原本的要求得的包含 $|P|$ 个点的最小圆覆盖问题，分解为包含 $|R| <= 3$ 个点的子问题，然后从这些子问题递归回原问题，也许有助于理解该方法的证明。另外值得注意的是，下面的递归算法中，永远有 $|R| <= 3$ ，因为一开始 $|R| == 0$ ，而当 $|R|$ 在递推增加至 $|R| == 3$ 时，算法立即返回 ）：

/*
 * Calculates the sed of a set of Points. Call initially with R = empty set.
 * P is the set of points in the plane. R is the set of points lying on the boundary of the current circle.
 */

function sed(P,R)
{
    if (P is empty or |R| = 3) then
         D := calcDiskDirectly(R)
    else
        choose a p from P randomly;
        D := sed(P - {p}, R);
        if (p lies NOT inside D) then
            D := sed(P - {p}, R u {p});
    return D;
}

算法的效果如下图所示：

{height=300}

2.4 最小外接矩形算法

2.4.1 点云数据处理中的矩形边界确定问题

在对点云的矩形平面进行检测时，往往会遇到矩形平面的边界确定的问题，如下图的左侧所示：

{height=300}

对于这个问题，一种直观的想法是，先使用PCA，将点云数据旋转到变化最大的两个方向上，然后找出上述点云的在横坐标轴上的最大差值和在纵坐标轴上的最大差值，将其视为矩形的长度和宽度，并以此将矩形平面区分出来。这在矩形的长度显著大于其宽度时较为有效。但是，如果矩形的长度和宽度相差不大，那么PCA就会倾向于将原始点云数据旋转到矩形的对角线上。如下图：

{height=300}

显然，这时候再沿用上述思路，就会出现较大的误差。

2.4.2 最小外接矩形算法

最小外接矩形（Minimum Bounding Rectangle，MBR）算法可以解决上述问题，它用来找到包围给定点集的面积最小的矩形。这个矩形的边可以是任意方向，并不是限定于与坐标轴平行。该算法利用了最小外接矩形的一个性质：

点集的最小外接矩形的一条边与其凸包（凸包是一组点，凸包围成的多边形内的任何点的连线都会完全落在该多边形内。凸包这就像把橡皮筋绕在点上，然后取下与橡皮筋接触的点。）的一条边平行的特性。

利用这个性质，算法可以这样设计：

1. 将凸包的所有点连接在一起形成一个线段列表。
2. 遍历该线段列表，将点集的 $x$ 轴旋转到与当前遍历到的线段的方向上。
3. 找到点集的最顶部和最底部的 $y$ 值以及最左侧和最右侧的 $x$ 值。利用这四个值可以算出一个边界框面积。
4. 为了找到最小边界框，算法现在只需取面积最小的边界框。该边界框即为所求。

算法的效果如下图所示：

{height=300}

3 尺寸测量具体步骤与结果

3.1 预处理

这个项目提供的点云数据均如下左图所示：在 $z=0$ 的位置处会广泛地遍布着噪声点云，而目标点云仅仅是高于 $z=0$ 平面的一小部分，为了排除这些地面回波的干扰，使用pcfitplane函数（MATLAB的一个内置的基于MSAC算法的平面拟合函数）提取出平面（如下中图），然后扣除该平面，并中心化，得到目标点云（如下右图）。这样处理过后的点云将成为后续处理的原始点云。

3.2 Task1

在这个项目中，Task1被用于校对测量程序和以及得出尺度变换因子（即将点云数据的尺度映射到真实尺度）。完成这一任务的步骤如下所示：

1. 观察四组点云数据，发现其中三组点云数据比较可靠，将它们用作几何尺寸测量的原始数据。
2. 使用pcfitplane函数检测出平面。
3. 将该平面的最高点视作长方体的高。
4. 使用PCA算法进一步去除噪声，并降维至二维平面上。
5. 使用最小外接矩形算法求出平面的长和宽，并将其视为长方体的长和宽。
6. 将通过这三组数据计算出来的长宽高分别求平均，将其作为最终的测量结果。
7. 参考已经标定好的数据，计算出尺度变换因子。

结果如下：

{height=300}

The length is: 60.387119
The width is: 29.811153
The height is: 29.997723
The RATIO is: 1.0092505e+03

3.3 Task2

完成这一任务的步骤如下所示：

1. 观察四组点云数据，发现其中两组点云数据比较可靠，将它们用作几何尺寸测量的原始数据。
2. 使用pcfitplane函数检测出平面。
3. 将该平面的最高点视作圆柱体的高。
4. 使用PCA算法进一步去除噪声，并降维至二维平面上。
5. 使用最小圆覆盖算法求出圆平面的直径，并将其视为圆柱体的直径。
6. 将通过这两组数据计算出来的直径和高分别求平均，将其作为最终的测量结果。
7. 参考已经标定好的数据，计算出真实尺寸。

结果如下：

{height=300}

The diameter is: 40.990364 mm
The height is: 40.868435 mm

3.4 Task3

这一任务比较复杂，由于观察到只有一组数据比较可靠，而且为了完成四棱锥顶部的夹角的计算，需要得到 $3$ 个侧面的方程，因此这里的步骤要稍微复杂一些，参见下面的步骤 $2,3,4$ 。完成这一任务的步骤如下所示：

1. 观察四组点云数据，发现只有一组点云数据比较可靠，将其用作几何尺寸测量的原始数据。
2. 由于四棱锥点云数据中有 $3$ 个侧面比较明显，因此在四棱锥点云数据的 $3$ 个侧面上，分别人工标注出的 $3$ 个点，通过这 $3$ 个点先粗糙地计算出这 $3$ 个侧面的法向量。
3. 将过步骤 $2$ 求得的法向量，作为pcfitplane函数的参考向量，以提高各个平面的拟合准确度，并最终得到 $3$ 个平面的方程。
4. 联合步骤 $3$ 的方程与 $z=0$， 可以求出四棱锥的顶点和底部与水平面交界的 $2$ 个端点。
5. 计算底部 $2$ 个端点的距离，作为四棱锥的底边长。
6. 使用顶点和底部的 $2$ 个端点，可得到四棱锥两条侧边，并计算出顶部的夹角（交线法）。
7. 参考已经标定好的数据，计算出真实尺寸。

结果如下：

The degree of A calculated by intersection line is: 40.854599°
The length of B is: 31.021509 mm

3.5 Task4

完成这一任务的步骤如下所示：

1. 观察四组点云数据，发现四组点云数据都比较可靠，将它们用作几何尺寸测量的原始数据。
2. 对于每组数据，都可以使用两种方法来求解：
   1. 使用pcfitsphere函数检测出半径。
   2. 使用最小二乘（LS）算法测出半径，其方程如下图所示：
3. 将这两种方法的通过这四组数据计算出来的半径分别求平均，将其作为最终的测量结果。
4. 参考已经标定好的数据，计算出真实尺寸（由于在计算过程中，发现LS算法的计算结果的方差更小，因此最终选择LS算法得出的结果作为最终结果）。

结果如下：

The radius calculated by MSAC is: 21.595808 mm
The radius calculated by LS is: 20.226751 mm

3.6 Task5

这个任务中，有两个物体被叠放在一起，增加了处理的难度。因此，这个任务被分为两部分，一部分专注于长方体的测量，另一部分专注于三棱柱的测量。

3.6.1 Task5-1

这里解决的是长方体的测量。完成这一任务的步骤如下所示：

1. 观察四组点云数据，发现只有一组点云数据比较可靠，将其用作几何尺寸测量的原始数据。
2. 原始点云数据中混杂着三棱柱带来的噪声，但由于长方体与三棱柱的分界面比较明显，因此可人工标注出位于分界面的 $3$ 个点，然后通过这 $3$ 个点构建平面方程，使用该平面方程将长方体与三棱柱分开。
3. 观察到分离后的点云数据仍然残存一些噪声，使用pcfitplane函数，并只保留近似位于同一平面上的数据。
4. 将该平面的最高点视作长方体的高。
5. 使用PCA算法进一步去除噪声，并降维至二维平面上。
6. 使用最小外接矩形算法求出平面的长和宽，并将其视为长方体的长和宽。
7. 参考已经标定好的数据，计算出真实尺寸。

结果如下：

{height=300}

The length of cuboid is: 42.524376 mm
The width of cuboid is: 43.529556 mm
The height of cuboid is: 20.391687 mm

3.6.2 Task5-2

这里解决的是三棱柱的测量。完成这一任务的步骤如下所示：

1. 观察四组点云数据，发现只有一组点云数据比较可靠，将其用作几何尺寸测量的原始数据。
2. 原始点云数据中混杂着长方体带来的噪声，但由于三棱柱与长方体的分界面比较明显，因此可人工标注出位于分界面的 $3$ 个点，然后通过这 $3$ 个点构建平面方程，使用该平面方程将三棱柱与长方体分开。
3. 由于三棱柱点云数据中的 $5$ 个侧面都比较明显，因此在三棱柱点云数据的 $5$ 个侧面上，分别人工标注出的 $3$ 个点，通过这 $3$ 个点先粗糙地计算出这 $5$ 个侧面的法向量。
4. 将过步骤 $3$ 求得的法向量，作为pcfitplane函数的参考向量，以提高各个平面的拟合准确度，并最终得到 $5$ 个平面的方程。
5. 联合步骤 $4$ 的方程， 可以求出三棱柱的 $6$ 个端点。
6. 通过这 $6$ 个端点，可以多次求得三棱柱的各个个尺寸，对其求平均，视为三棱柱的尺寸。
7. 参考已经标定好的数据，计算出真实尺寸。

结果如下（本人认为结果存在很大误差）：

The length of A is: 44.597523 mm
The length of B is: 30.320036 mm
The length of C is: 47.664532 mm
The length of D is: 14.551025 mm

3.7 Task6

这个任务中，点云数据中同时包含了两个物体被叠放在一起，增加了处理的难度。因此，这个任务被分为两部分，一部分专注于三棱柱的测量，另一部分专注于半圆柱体的测量。

3.7.1 Task6-1

这里解决的是三棱柱的测量。完成这一任务的步骤如下所示：

1. 观察四组点云数据，发现其中两组点云数据比较可靠，将它们用作几何尺寸测量的原始数据。
2. 原始点云数据中混杂着半圆柱带来的噪声，但由于三棱柱与半圆柱的分界面比较明显，因此可人工标注出位于分界面的 $3$ 个点，然后通过这 $3$ 个点构建平面方程，使用该平面方程将三棱柱与长方体分开。
3. 由于三棱柱点云数据中的 $4$ 个侧面（除底面）都比较明显，因此在三棱柱点云数据的 $4$ 个侧面上，分别人工标注出的 $3$ 个点，通过这 $3$ 个点先粗糙地计算出这 $4$ 个侧面的法向量。
4. 将过步骤 $3$ 求得的法向量，作为pcfitplane函数的参考向量，以提高各个平面的拟合准确度，并最终得到 $4$ 个平面的方程。另外，用顶面方程的法向量近似底面方程的法向量，且将底面方程的截距看作是 $0$ ，由此，得到全部 $5$ 个方程。
5. 联合步骤 $4$ 的方程， 可以求出三棱柱的 $6$ 个端点。
6. 通过这 $6$ 个端点，可以多次求得三棱柱的各个个尺寸，对其求平均（底面的三角形的边长的误差太大，不参与平均），视为三棱柱的尺寸。
7. 参考已经标定好的数据，计算出真实尺寸。

结果如下（本人认为结果存在很大误差）：

The length of A is: 16.184345 mm
The length of B is: 13.783567 mm
The length of C is: 22.052763 mm
The length of D is: 57.732460 mm

3.7.2 Task6-2

这里解决的是半圆柱体的测量。完成这一任务的步骤如下所示：

1. 观察四组点云数据，发现其中两组点云数据比较可靠，将它们用作几何尺寸测量的原始数据。
2. 观察到点云数据有一些噪声，使用pcfitplane函数，并只保留近似位于同一平面上的数据。
3. 使用PCA算法进一步去除噪声，并降维至二维平面上。
4. 使用最小外接矩形算法求出平面的长和宽，并将其视为长方体的长和宽。
5. 将通过这两组数据计算出来的长和宽分别求平均，将其作为最终的测量结果。
6. 参考已经标定好的数据，计算出真实尺寸。

结果如下：

{height=300}

The length is: 29.124096 mm
The width is: 21.622047 mm

3.8 Task7

这个任务中，圆柱体是被握持在手上测量的。由于存在两组点云数据，分别恰好完整地捕捉到该圆柱体的顶面和侧面。因此，这个任务被分为两部分，一部分专注于顶面的测量，另一部分专注于侧面的测量。

3.8.1 Task7-1

这里解决的是圆柱体顶面的测量。完成这一任务的步骤如下所示：

1. 观察四组点云数据，发现只有一组点云数据比较可靠，将其用作几何尺寸测量的原始数据。
2. 原始点云数据中混杂着人的肢体带来的噪声，但由于人的肢体与圆柱体顶面的分界面比较明显，因此可人工标注出位于分界面的 $3$ 个点，然后通过这 $3$ 个点构建平面方程，使用该平面方程将圆柱体顶面与其它点云分开。
3. 观察到分离后的点云数据仍然残存一些噪声，使用pcfitplane函数，并只保留近似位于同一平面上的数据。
4. 使用PCA算法进一步去除噪声，并降维至二维平面上。
5. 使用最小圆覆盖算法求出圆平面的直径，并将其视为圆柱体的直径。
6. 参考已经标定好的数据，计算出真实尺寸。

结果如下：

{height=300}

The diameter is: 30.542925 mm

3.8.2 Task7-2

这里解决的是圆柱体侧面的测量。完成这一任务的步骤如下所示：

1. 观察四组点云数据，发现只有一组点云数据比较可靠，将其用作几何尺寸测量的原始数据。
2. 原始点云数据中混杂着人的肢体带来的噪声，但由于人的肢体与圆柱体侧面的分界面比较明显，因此可人工标注出位于分界面的 $3$ 个点，然后通过这 $3$ 个点构建平面方程，使用该平面方程将圆柱体侧面与其它点云分开。
3. 观察到分离后的点云数据仍然残存一些噪声，使用pcfitplane函数，并只保留近似位于同一平面上的数据。
4. 使用PCA算法进一步去除噪声，并降维至二维平面上。
5. 使用最小外接矩形算法求出圆平面的直径，并将其视为圆柱体的高。
6. 参考已经标定好的数据，计算出真实尺寸。

结果如下：

{height=300}

The height is: 61.434860 mm

4 总结

• 特色

  在这个项目中，本人觉得比较有特色的算法设计是对一些并不明显的平面人工取点，然后取得这个平面的一个近似的法向量来作为pcfitplane函数的参考向量，使其能在一众噪声中准确地提取出目标平面。这样的思路也在点云数据中同时存在多个平面时取得很好的效果。

• 不足

  1. 没有对点云数据中存在的畸变进行针对性的处理。 在对点云数据进行观察时发现，部分点云数据存在一定的畸变，比如：

     • 反映同一个平面的点云数据呈现的是有规律的波浪形分布；
     • 有些点云数据在接近 $z = 0$ 平面处时会发散开；

     这些畸变实际上对测量结果产生很大的影响，而目前的工作并没有对这些畸变产生的原因进行分析，继而也不能提出一些针对性的算法来清除它们。

  2. 当物体摆放的位置、物体间的组合不相同，不能集成数据，这也降低了数据的利用率。

  3. 在预处理步骤清除地面回波时，所选取的距离容限太大，使得目标物体的靠近地面的点云数据也被一并清除掉，有可能对测量的准确性造成比较大的影响。

• 收获

  我在这个项目中的收获是学习到了很多有意思的算法，也将课程上的知识应用于实践中。另外，我也对最小二乘法、特征值、特征向量、PCA在数据处理上的作用有了更深的认识。