为了提高点云数据处理和应用的效率，需要对海量点云数据进行精简^[1-4]。近年来，国内外学者对点云精简进行了大量研究，并取得了大量的研究成果。经典的点云算法有包围盒法^[5]、曲率采样^[6]、保留边界法^[7]、聚类法等^[8]。Weir等^[9]利用体包围盒法对点云进行精简，该算法具有体包围盒高效的特性，但由于其精简过程是以一个或者一部分点来代替区域点，所以容易造成点云显著特征丢失。Wang等^[10]利用特征参数法提取特征，并利用均匀球面采样法对非特征点进行精简，该算法具有特征保留的特性，但易造成点云数据空洞现象。Lee等^[11]利用非均匀网格精简点云数据，该方法较好地保留了点云模型的特征点，但在平坦区域保留点分布不均匀。Yuan等^[12]为了改进传统算法经常容易丢失几何信息的问题，利用球体树构造和k-means法对点云进行精简，该算法不仅对局部细节有很好的效果，而且对模型完整性具有良好的保留效果，但其构型和聚类过程比较复杂，效率较低。杨秋翔等^[13]利用Hausdorff距离算法提取特征点，实现了点云精简，优点是几何特征保留较好，精简效率高，但容易造成数据空洞。陈璋雯等^[14]提出了基于模糊熵迭代的三维点云精简算法，该算法能很好地保留细节特征，但没有考虑非特征区域精简问题。目前，点云精简算法优点突出，精简效果较好，但也存在点云数据特征保留不完整和对小曲率点云精简易造成数据空洞等不足。

针对上述算法精简散乱点云时出现几何特征保留较少和过精简的问题，本文提出了一种融合k-means聚类^[15]和Hausdorff距离^[16-17]的点云精简算法。为了提高点数据的索引效率，该算法首先在八叉树算法的基础上构建点云数据的拓扑关系；然后利用最小二乘法进行局部曲面拟合，分别计算目标点与其邻域点集内各点连线和目标点切平面的夹角，以夹角的平均值估算目标点曲率，通过局部二次曲面拟合方程分别计算点云数据点的主曲率；接着根据主曲率计算Hausdorff距离，通过估计曲率和主曲率的Hausdorff距离来确定每个八叉树小长方体块的特征点和特征区域，根据精简目标要求实现点云特征提取；最后对非特征区域进行k-means聚类提取特征点。本文算法利用Hausdorff距离具有相似性度量的特性提取特征点，并利用k-means聚类算法对非特征区域进行再精简，能够保证点云数据整体的完整性，避免过精简造成数据空洞。

1 点云精简算法 1.1 局部点云曲面拟合

本文在最小二乘法曲面拟合的基础上对点云数据点的邻域点集进行拟合。针对实验点云散乱无序的特性，在八叉树算法的基础上构建点云数据的拓扑关系，以提高点云数据索引的效率。将目标点$ {p}_{i}$到相邻点的欧氏距离存入数组中，并对数组数据进行降序排序，取数组中前k个点$ {\left\{\left({x}_{i}, {y}_{i}, {z}_{i}\right)\right\}}_{i=1}^{k} $，邻域点集记为$ k\left({q}_{i}\right)$。点云数据目标点的k邻域点集只与目标点$ {p}_{i} $有关，与其他点及其邻域点不相关。许多学者对k值的选取做了大量实验^[6]，证明k取15~30效果较好，本文k取17。运用最小二乘二次多项式对经过坐标转换的目标点及其邻域点集$ k\left({q}_{i}\right)$拟合得到局部曲面，二次多项式曲面拟合方程$ f({x}_{i}, {y}_{i}) $ ^[18]为：

$ f({x}_{i}, {y}_{i})={a}_{1}+{a}_{2}{x}_{i}+{a}_{3}{x}_{i}^{2}+{a}_{4}{x}_{i}{y}_{i}+{a}_{5}{y}_{i}+{a}_{6}{{y}_{i}}^{2} $

(1)

曲面的一般方程可以表示为：

$ {z}_{i}=f({x}_{i}, {y}_{i}) $

(2)

对目标点及其邻域点集的坐标(x，y，z)先进行坐标转换，即转换成局部坐标系，再将坐标代入到方程(2)中求解二次多项式曲面拟合系数。

由于方程个数超过了拟合曲面方程系数，因此将拟合曲面系数的求解过程转换成式(3)，即根据拟合曲面在目标及其邻域点集处的值与目标及其邻域点集的真实值之差的平方和达到最小来求解拟合曲面系数。

$ D={\stackrel{n}{\sum \limits_{i=1}}\left[f\left({x}_{i}, {y}_{i}\right)-{z}_{i}\right]}^{2} $

(3)

由最小二乘算法原理可知，求解方程(3)可转换为：

$ \frac{\partial D}{\partial {a}_{i}}=0, i=\mathrm{1, 2}\dots n $

(4)

1.2 曲率的估算

首先根据局部点云曲面拟合的曲面方程计算目标点$ {p}_{i}$的切平面法向量$ {\mathit{\boldsymbol{n}}}_{\mathit{\boldsymbol{i}}}$，然后利用目标点坐标和切平面法向量计算该点的切平面方程，最后通过邻域点与目标点切平面的偏离角度平均值估算曲率^[14]。如图 1所示，估算目标点$ {p}_{i}$曲率步骤如下：

1) 将目标点$ {p}_{i}$分别和邻域内的其他k个点$ k\left({q}_{j}\right) $相连，得到k个向量$ \overrightarrow{{p}_{i}{q}_{1}}$，$ \overrightarrow{{p}_{i}{q}_{2}} $…$ \overrightarrow{{p}_{i}{q}_{{}_{k}}}$。

2) 向量$ \overrightarrow{{p}_{i}{q}_{j}}$与切平面的夹角为$ {\theta }_{ij} $，它的正弦值为向量$ \overrightarrow{{p}_{i}{q}_{j}}$与切平面法向量$ \mathit{\boldsymbol{n}} $的夹角$ {\theta }_{n}$余弦的绝对值：

$ \mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{ij}=\left|\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{n}\right|=\left|\frac{\overrightarrow{{p}_{i}{q}_{1}}\cdot \mathit{\boldsymbol{n}}}{\left|\overrightarrow{{p}_{i}{q}_{1}}\right|\times \left|\mathit{\boldsymbol{n}}\right|}\right| $

(5)

由式(5)分别计算向量$\overrightarrow{{p}_{i}{q}_{1}} $，$ \overrightarrow{{p}_{i}{q}_{2}} $$ \dots $$ \overrightarrow{{p}_{i}{q}_{{}_{k}}} $与切平面夹角的正弦值$ \mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{i1}, \mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{i2}\dots \mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{ik}$，得到目标点的估计曲率：

$ Q=\left(\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{i1}+\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{i2}+\dots +\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{ik}\right)/k $

(6)

式中，$ Q $为目标点的估计曲率；$ k $为目标点的邻域点数。

曲率代表几何物体表面高低起伏不平坦的一种衡量，曲率较大的区域一般弯曲程度较大，且包含较多的特征信息，曲率较小的区域较为平坦或者平缓。基于曲率的点云精简算法以曲率作为精简原则，保留曲率较大的点，在曲率较小的地方精简大量点易造成精简结果出现数据空洞现象；曲率较大的点与周围点共同构成物体特征信息，并且曲率存在渐变过程，而基于曲率的点云精简算法随着精简率的提高会去掉部分大曲率点周围的点，导致物体特征区域出现变形，降低了数据精度。因此，本文引入Hausdorff距离计算目标点与k邻域点的相似性来确定特征点。

1.3 基于Hausdorff距离的k-means点云聚类

假设有两组点集A、B, 则两组点集之间的Hausdorff距离定义为^[19]：

$ h\left(A, B\right)=\mathrm{m}\mathrm{a}\mathrm{x}\left\{d(A, B), d(B, A)\right\} $

(7)

式中，$ d\left(A, B\right)$和$ d\left(B, A\right)$分别表示集合A到集合B和集合B到集合A的单向Hausdorff距离。本文通过§1.1中曲面拟合方程计算点云数据目标点$ {p}_{i}$与该点邻域点集某一点$ {q}_{j}$的主曲率$ {k}_{1}\mathrm{、}{k}_{2}$和$ {{k}_{1}}^{\prime }\mathrm{、}{{k}_{2}}^{\prime }$，分别记为集合$ \left\{{k}_{1}, {k}_{2}\right\} $和$ \left\{{{k}_{1}}^{\prime }, {{k}_{2}}^{\prime }\right\}$。则目标点$ {p}_{i}$与其邻域内点曲率的Hausdorff距离$ {h}_{i1}, {h}_{i2}\dots {h}_{ik} $的计算公式如下：

$ d(k, k\prime )=\underset{i=\mathrm{1, 2}}{\mathrm{m}\mathrm{a}\mathrm{x}}\left(\underset{j=\mathrm{1, 2}}{\mathrm{m}\mathrm{i}\mathrm{n}}\left(\frac{‖{k}_{i}-{{k}_{j}}^{\prime }‖}{‖{k}_{i}‖+‖{{k}_{j}}^{\prime }‖}\right)\right) $

(8)

$ d(k\prime , k)=\underset{j=\mathrm{1, 2}}{\mathrm{m}\mathrm{a}\mathrm{x}}\left(\underset{i=\mathrm{1, 2}}{\mathrm{m}\mathrm{i}\mathrm{n}}\left(\frac{‖{k}_{i}-{{k}_{j}}^{\prime }‖}{‖{k}_{i}‖+‖{{k}_{j}}^{\prime }‖}\right)\right) $

(9)

$ h\left(k, k\prime \right)=\mathrm{m}\mathrm{a}\mathrm{x}\left\{d(k, k\prime ), d(k\prime , k)\right\} $

(10)

因此，目标点$ {p}_{i}$曲率的Hausdorff距离$ {H}_{{p}_{i}} $可以表示为：

$ {H}_{{p}_{i}}=\mathrm{m}\mathrm{a}\mathrm{x}\left({h}_{i1}, {h}_{i2}\dots {h}_{ik}\right) $

(11)

如果式(8)和式(9)中的分母$ ‖{k}_{i}‖+‖{{k}_{j}}^{\prime }‖$无限接近零并且分子远大于分母时，会造成计算结果出现较大偏差。为了避免分母过小所带来的误差，本文引入限制变量$ \partial $，如果分母小于$ \partial $时，用$ \frac{‖{k}_{i}‖+‖{{k}_{j}}^{\prime }‖+\partial }{2} $来表示分母。

本文通过目标点与其邻域点主曲率的Hausdorff距离来确定目标点与其邻域点的相似性，进而确定点云模型中点的曲率变化程度，结合曲率对点云进行分类。图 2为两种方法提取的点云特征点分布图，图 2(a)为原始点云；图 2(b)为估计曲率分布模型，本文根据曲率大小将曲率分成大曲率、中曲率和小曲率3类，其中红色表示曲率较大的点云，青色表示曲率较小、区域平坦的点云，黄色表示介于红色和青色之间的点云；图 2(c)为Hausdorff距离分布模型，分层设色方法与曲率估计方法相同。由图 2(b)可知，在兔子特征较为明显的地方，其曲率较大；在兔子的腿部和耳朵等地方，红色和黄色点分布较多。对比图 2(b)和图 2(c)可以发现，图 2(c)在兔子特征明显的耳朵和腿部的红色点和黄色点区域面积要大于图 2(b)。

点云k-means聚类算法的核心是利用点云数据点与点之间的欧氏距离作为聚类指标，即选择k个聚类中心，记为$ J({j}_{1}, {j}_{2}\dots {j}_{k}) $，形成k个点簇，记为$ C\left({c}_{1}\left({p}_{1i}\right), {c}_{2}\left({p}_{2i}\right)\dots {c}_{k}\left({p}_{ki}\right)\right)$。首先分别计算点云中每个点到聚类中心点集$ J({p}_{1}, {p}_{2}\dots {p}_{k}) $的欧氏距离，比较计算所得的欧氏距离，将每个点云数据点归为距离聚类中心最近的点；然后根据每个聚类中心所包含的点集重新计算聚类中心$ J\text{'}({j}_{1}, {j}_{2}\dots {j}_{k}) $，比较此次聚类中心点集和前一次聚类中心点集，如果相邻两次聚类中心点集距离之和小于预先设定值，则停止聚类，否则重复计算聚类中心；最后以聚类中心代替每一类别中的点，从而达到精简的目的^[20]。

本文k-means聚类中k值的选取方法如下：(1)设原始点云数量为S₀，八叉树立方体数为S_c，目标精简率为$ \delta $，则精简后的点云数量为$ \left(1-\delta \right){S}_{0} $，本文设定k为精简后的点云数量的10%，即$ \left(1-\delta \right){S}_{0}/10 $，假设$ \left(1-\delta \right){S}_{0}/10$ < 2S_c，则k取值为$ 2{S}_{c}$；(2)分析点云主曲率的Hausdorff距离并确定初精简点集，将初精简集中的点分别划分到八叉树的小长方体中，根据k值大小和小长方体中点的个数来确定小长方体中聚类中心点的个数，按照每个小长方体的聚类中心点的个数均匀地选取聚类中心点。

若本文选取第i个八叉树划分算法中的小长方体作为示例，设该小长方体的非特征点集为$ {D}_{i} $，点数为$ {n}_{i}$，体积为$ {V}_{i}$，则它的密度$ {\rho }_{i}$为：

$ {\rho }_{i}={n}_{i}/{V}_{i} $

(12)

k值的计算公式为：

$ {k}_{i}=\left\lceil {k{\rho _i}/\sum\limits_{j = 1}^n {_{}{\rho _j}} } \right\rceil $

(13)

式中，$ \left\lceil \cdot \right\rceil $为向上取整符号；n为小长方体的数量。随机选取$ {k}_{i} $个数据点作为它的初始聚类中心。

点云精简算法在精简过程中得到特征点和非特征点后，往往会直接删除非特征点。本文使用Hausdorff距离计算出特征点和非特征点后，得到特征区域，并利用k-means聚类算法对非特征点进行处理。由于k-means具有聚类的特性，从而保证了点云模型非特征区域能够保留部分点且分布均匀。

融合k-means聚类和Hausdorff距离的点云精简算法步骤如下:

1) 将Matlab读取到的点云数据点坐标存入矩阵O_CloudPoints中。

2) 在八叉树算法的基础上构建点云数据的空间拓扑关系。八叉树构建完成后，对八叉树划分结果建立索引并存入数组O_OctTree中。

3) 在八叉树算法的基础上构建每个数据点的k-邻域索引并存入K_Nears中。

4) 在最小二乘曲面拟合算法的基础上对点云数据每个点的邻域点集进行拟合，根据拟合方程计算主曲率，通过§1.3的计算方法计算目标点主曲率的Hausdorff距离。

5) 利用步骤4)中的方法依次计算点云数据所有点的Hausdorff距离。

6) 按照由大到小的顺序对点云数据所有点的主曲率的Hausdorff距离进行排序，根据精简比例以及二次精简保留的k值确定第1次保留的点数$ \left(1-\delta \right){S}_{0}-k$，以排序后的点云数据点序号为$ \left(1-\delta \right){S}_{0}-k$的点的Hausdorff距离作为阈值ε。根据构建的O_OctTree遍历每个小长方体包含的点云数据点，并将非特征点加入$ {D}_{i} $中。

7) 根据式(12)和式(13)计算第i个小长方体的k值大小($ {k}_{i} $)，随机选取小长方体内的$ {k}_{i}$个数据点作为该小长方体的初始聚类中心，记为$ {J}_{i}({p}_{1}, {p}_{2}\dots {p}_{{k}_{i}})$。根据k-means算法确定最终的聚类中心$ J\text{'}({p}_{1}, {p}_{2}\dots {p}_{{k}_{i}})$，以该聚类结果代替原始数据完成精简。

8) 将Hausdorff距离和k-means算法提取的结果进行融合，得到最终的精简数据。

融合k-means聚类和Hausdorff距离的点云精简算法流程如图 3所示。

2 实验结果与分析

将本文的融合算法分别与基于k-means聚类的点云精简算法和基于Hausdorff的点云精简算法进行对比，以验证本文算法的准确性和有效性。采用经典的兔子数据和马数据(该数据来源于斯坦福大学)作为实验数据，点云数量分别为31 606和48 485。兔子模型的特征较为集中，主要集中在兔子四肢和耳朵等处；马模型的特征主要集中在马的四肢和头部等地方。图 4是3种算法对兔子点云精简后的结果，图 5是3种算法对马点云精简后的结果，表 1是3种精简算法对兔子和马点云精简后的精简率和点云数量对比。图 4(b)和图 5(b)表示k-means聚类精简结果，其精简结果具有k-means聚类的特性，精简后的点云分布均匀，并具有一定填充空洞的效果，但易精简大量特征点，造成特征区域变形，丢失大量特征信息。图 4(c)和图 5(c)表示基于Hausdorff距离的精简结果，由于Hausdorff距离能判断相似性，所以该算法能保留特征点周围大量的点，但其精简结果出现了大量的空洞现象。图 4(d)和图 5(d)表示本文融合算法的精简结果，通过图 4、图 5和表 1可以看出，在总的保留点数差不多的情况下，本文算法在兔子和马的特征区域保留了大量点，并且本文融合算法精简结果的轮廓保留完整，未造成数据空洞。

本文利用Geomagic软件对原始点云数据以及3种算法精简后的点云数据进行三维曲面重构，分别对实验数据中特征明显的区域(兔子耳朵和马四肢)进行了对比，并对兔子耳朵和马四肢的三角面片数量进行了统计分析。为了保证截取的细节部位形状、大小和位置相同，本文在裁剪原始点云重构模型和3种算法精简后的点云数据重构模型时，裁剪面的坐标和旋转角度完全相同。图 6和图 7分别是兔子耳朵和马四肢原始点云重构模型和3种算法精简后的点云数据重构模型，表 2是兔子耳朵和马四肢的三角面片数量。由表 2可知，基于Hausdorff距离的精简结果重构网格切割后，兔子耳朵处三角面片数量为5 698，马四肢处为8 490；基于k-means的精简算法在兔子耳朵处的三角面片数量为3 726，在马四肢处为4 650。由此可见，基于Hausdorff距离的精简结果在特征区域保留的信息明显多于基于k-means的精简算法，但基于Hausdorff距离的精简结果出现了大量空洞。本文融合算法融合了k-means聚类和Hausdorff距离，其重构网格在兔子耳朵处的三角面片数量为4 942，在马四肢的三角面片数量为6 908，与基于Hausdorff距离算法的结果相差不大，并且未造成数据空洞。因此，本文算法在保留更多特征信息的同时，也弥补了基于Hausdorff距离易造成数据空洞的不足，提升了算法的适用范围。

为了分析不同精简率的精简结果，利用本文算法分别对精简率目标约19.5%、34.9%、46.6%、66.2%和86.3%的兔子点云进行精简，其结果如图 8所示。其中，图 8(a)为原始兔子点云，图 8(b)~8(f)分别对应上述精简率精简后的点数。由图 8可以看出，本文融合算法在兔子特征信息多的区域保留了大量点，比如兔子四肢和耳朵等起伏较大的不平坦地方；在兔子较为平坦的地方则点比较稀疏，但对兔子的轮廓保留较好，并且在小曲率点区域的点分布均匀。

为了对精简结果质量进行评价^[21]，本文选取精简后点云重构网格的表面积和精简后兔子点云数据与原始兔子点云数据误差$ \mathrm{\Delta }(S, {S}^{\ast })=d(p, {S}^{\ast })$的标准差作为量化指标，其中$ d(p, {S}^{\ast }) $表示原始点云数据重构网格曲面S上点p到精简后点云重构网格曲面$ {S}^{\ast } $上投影点$ {p}^{\ast }$的欧氏距离^{[15, 22]}。本文分别计算了3种精简算法精简结果的质量量化指标，如图 9所示。

由图 9(a)可知，随着精简率的增大，精简结果的点云数量占比降低，3种算法精简结果的误差标准差都呈现增大的趋势，本文融合算法和基于Hausdorff距离的精简算法的误差标准差明显小于基于k-means聚类算法。由图 9(b)可知，基于Hausdorff距离算法的精简结果重构网格表面积随着精简结果点云数量占比的减少而减少，所以本文算法具备精简结果误差标准差较小和重构网格表面积与原模型表面积相差不大的优点。但由于本文算法先运用主曲率的Hausdorff距离确定初精简集，再运用k-means聚类，所以运算时间要高于另外两种算法。

综合分析图 4至图 8可以得出，基于k-means的精简算法具有高效精简的特性，且精简结果分布均匀，但其在兔子和马的特征区域丢失了大量特征信息；基于Hausdorff距离的精简算法能较完整地保留兔子和马的特征信息，但其精简结果具有大量的空洞现象；本文算法能较完整地保留兔子和马的特征信息，并且能避免形成空洞现象。

3 结语

基于Hausdorff距离的点云精简算法能较好地保留点云形状特征，提高精简效率，然而直接去除非特征区域点容易造成数据空洞；基于k-means聚类的点云精简算法有较高的精简效率, 但对细节特征保留效果较差。本文融合Hausdorff距离和k-means聚类，利用Hausdorff距离具有相似性度量的特性提取特征区域，并对非特征区域进行k-means聚类提取特征点，凸现了k-means聚类算法保留点云数据整体轮廓的特点，避免了数据空洞问题, 从而使精简结果更加完整，精简率高，且保留了较多的形状特征。实验表明，该算法不但能保持原始数据的整体轮廓且避免出现数据空洞，还能够较好地保留点云数据的局部细节形状特征。如何确定最优聚类中心以及提升融合算法的效率还有待进一步研究。