基于点云数据的古陶瓷碎片边缘轮廓提取技术

周 强，张 敏，李 巍，ANVARJON NORMURODOV，杜丹阳

(陕西科技大学 电气与控制工程学院，陕西 西安 710021)

0 引言

越来越多的古文物被发掘出土，其中古陶瓷类文物在出土时往往是以破碎的状态呈现在世人眼前，所以古陶瓷的修复成为了考古工作人员的一项重要工作；由于出土碎片数量较多，碎片拼接就成了古文物修复工作中复杂且耗时的一个环节；并且实物拼接中，人为失误造成的损失是不可挽回的，因此数字化碎片拼接在有效的提高了文物修复效率的同时，又降低了古文物被二次破坏的几率.数字化碎片拼接依赖于三维成像技术与三维图像处理技术，根据碎片特点主要分为两类，一类为厚壁文物，以断裂面的特征进行拼接[1,2]；另一类为薄壁类文物，这一类大多采用边缘轮廓作为拼接特征进行拼接.对于薄壁类文物，边缘轮廓提取技术成为了拼接过程中的重要支撑技术之一[3,4]，边缘轮廓提取的准确性直接决定了后期拼接的效果.

碎片特征提取方法主要有以下两类，一类是基于区域聚类的方法，还有一类是基于几何特征量估算的方法[5,6].基于区域聚类方法相对较少，张雨禾等[7]先利用特征点与非特征点的反K近邻具有不同尺度这一特性将点云进行划分，然后改进密度空间聚类方法并创建新的特征提取法则达到准确获取特征点的效果.该算法在反K近邻的高效计算方面有待加强.Bazazian D等[8]通过分析由K个点构成的K近邻定义的协方差矩阵的特征值提出了一种快速而精确的边缘特征检测方法，但该方法应用在古陶瓷碎片边缘轮廓提取算法复杂度较高.特征提取更多的还是基于几何特征的方法进行[9-11]，韩玉川等[12]首先对点云中的点构建K近邻，然后计算邻域的边缘系数来选择新的种子点，通过控制搜索视角来控制点云搜索方向，最终实现点云边缘轮廓的提取，该方法适用于大型点云边缘轮廓提取，对于小型点云所提取的边缘轮廓中含有非边缘点.陈义任等[13]先对点云的每个点构建K邻域，然后在每个邻域上进行平面拟合并将邻域点投影至拟合平面上，最终通过场力的大小来对点云的边缘点进行提取.该方法由于需要对点云中所有点进行邻域平面拟合，所以时间复杂度较高，并且当点在边缘上时，邻域平面拟合的准确性下降，反而降低了场力判断的准确性.李康等[14]基于三维网格的边重数提取主轮廓线，然后通过网格曲面分割识别断裂面然后对断裂面的二级邻接生长曲面进行曲线扫描提取次轮廓线，最后从主轮廓线和次轮廓线中得到文物碎片的特征轮廓线，该方法需要对点云进行网格化，在一定程度上丢失了碎片原有的细节特征，且基于网格提取轮廓线的时间复杂度较高.

本文针对古陶瓷碎片在文献[12]和文献[13]所做的工作下进行改进，提出一种基于点云数据的古陶瓷碎片边缘轮廓提取技术，首先对点云构建K近邻，然后针对古陶瓷碎片的特点进行全局平面拟合得到辅助投影面，将邻域点投影至拟合平面上，然后使用改进的边缘系数方法进行边缘点提取，最后再对提取的边缘进行二次优化，从而实现对古陶瓷点云模型的边缘轮廓的精确提取，最后使用F1-Score对边缘提取算法准确率进行评测.

本文方法与前面所列举方法相比，具备以下创新点：

(1)基于碎片几何特征的全局化平面拟合，不需要对点云进行网格化与每个点云邻域进行平面拟合，增强了边缘点的边缘特征并减少了算法时间复杂度.

(2)改进的边缘系数方法，对每个邻域点对其结点的作用加以权重系数，减弱较远邻域点在边缘点判断中的权重，有效的提高了点云边缘轮廓提取的准确度.

(3)对提取的边缘轮廓进行二次优化，有效的放宽了算法中阈值设置的约束性，提取的边缘轮廓更精确，不含非边缘点.

(4)将现有的点云边缘轮廓提取物技术应用于古文物数字化恢复领域，基于点云数据处理方法提取古陶瓷碎片边缘轮廓，保留了边缘轮廓的细节特征，有利于提取拼接特征进行后期古陶瓷碎片拼接.

1 散乱点云结构化

由于点云数据的数字化采集、存贮机制限制，通常点云数据文件中，以随机顺序将点云采集设备通过传感器得到的空间位置信息，以逐条的形式保存，形成点云数据模型.点云模型中的空间点位置信息数据条目之间，不存在特定的逻辑关系，当需要依据特定关系对模型中部分点集进行运算或操作时，需根据特定关系，采用遍历的方式对所有点信息数据进行查询计算，导致计算效率低下.为了避免这个问题，针对特定的算法需求，首先依照一定的规则，如距离、角度等空间位置关系，对点云数据构建拓扑关系，即散乱点云结构化，成为三维点云处理中的一个重要环节.

将一个散乱的三维点云表示为：P(n)={p1，p2，p3,…,pi}，i∈[1,n].其中n表示该点云中包含点的数量，pi(xi，yi，zi)表示点云中的第i个点，其中的xi，yi，zi表示该点的三维坐标数据x、y、z.

目前建立点云的拓扑关系的方法主要有Octree、KD-Tree、栅格法等[15,16]，构建点云的拓扑关系就是为点云中的每个点构建邻域，通过邻域将空间中的离散点联系在一起，从而可以进行点云的滤波、分割、特征提取等处理.

Octree在分割的过程中会出现一些空的子空间块，这些空子空间块的产生会影响点云的查询效率.栅格法虽然原理相对简单易懂，但操作起来相对较难，尤其在处理数据量较大且密度分布不均匀的点云时，栅格边长的选取是最大的难题，边长选取的过大会导致每个小栅格中的点的数量过多，增加计算时间，边长选取过小会导致栅格的数量过多，影响最终的拼接精度.

KD-Tree是二叉树的推广，是一种分割多维数据空间的数据结构，主要应用于多维空间关键数据的搜索.通过KD-Tree对空间结构的划分可以方便快捷的确定每一个样点的邻近点的搜索路径，而且不受点云密度的影响，应用面较广，针对点云的特征提取，KD-Tree在不改变原有信息的完整度的前提下保证了特征点提取的准确性，基于此本文选择KD-Tree来构建拓扑结构.

KD-Tree构建邻域有两种构建方式：

(1)搜索距离采样点pi最近的k个点，这k个点构成点pi的一个K近邻Pik={pi1，pi2，pi3,…,pik}，这种方法称为K邻域搜索.

(2)构造以采样点pi为中心，以半径为R的球形邻域，落在该球形邻域内的点共同构成点pi的邻域Pik={pi1，pi2，pi3,…,pik}，该方法称为半径搜索.

本文先对古陶瓷碎片的点云模型中的每个点构建其最近的K个点的最近邻，通过对近邻特征的统计与判断实现边缘轮廓的初步提取，最后再对初次提取的边缘轮廓的每个点构建半径为R的邻域，通过判断落在该邻域内的点的数量来去除第一次错误提取的非边缘点.

2 边缘轮廓提取

古陶瓷碎片表面主要分为两类，一类是光滑且原来就存在的内外表面，另一类是古陶瓷在破碎过程中产生的断裂面.断裂面与内外表面相交的线叫做断裂面轮廓线，在边缘轮廓上的点称作边缘点，其余落在碎片内外表面内或者断裂面内的点称为非边缘点，如图1所示.

图1 古陶瓷碎片边缘点与非边缘点

2.1 改进的边缘系数法

在二维空间中，一个平面上均匀的分布着一些点，对某一点pi构建其邻域，落在其邻域中的其它点Pik={pi1,pi2,pi3,…,pik}都会对点pi有一个作用矢量，定义该作用矢量为引力Fi.若pi点为边缘点，则其邻域的点Pik对其作用力偏向一个方向，如图2(a)所示，F1～F4分别代表邻域点对点pi的作用引力，最终F1～F4的合力在数值上表现为一个比较大的数，而非边缘点的邻域点会比较均匀分布在其四周，如图2(b)所示，同样对其构建邻域，邻域内的点对该点的作用力F1～F4向四周分布，通过矢量分解在相反方向上互相抵消，最终合力在数值上将会是一个极小的值.

(a)边缘点受力状况 (b)非边缘点受力状况图2 边缘点与非边缘点受力示意图

图3 边缘点与非边缘点邻域分布图

(1)

2.2 基于平面投影的边缘轮廓提取

根据古陶瓷碎片的几何形状特点，其断裂后的碎片模型可近似看作是光滑且曲率变化很小的一个曲面，文献[13]中需要首先对点云中每个点进行局部邻域平面拟合然后再基于拟合平面进行边缘点提取，该方法应用于数据量较大的点云时其时间复杂度较高.针对古陶瓷碎片的几何特点与文献[13]中的方法，本文在此提出基于最大投影面的平面拟合，首先对点云进行滤波，在保证曲面几何特征不丢失的前提下对点云进行简化，然后对采样简化后的点云使用最小二乘法进行平面拟合得到拟合平面π，然后基于平面π使用边缘系数法提取边缘点，该方法能有效减少边缘轮廓提取算法的时间复杂度，并对边缘点有特征增强作用，提高边缘点提取的准确率.

由于古陶瓷碎片的大小不一，对于曲率变化较小的碎片拟合的平面具有模型原有的一些几何特征，但对于较大的碎片如图4所示，这种碎片的器型从侧面延伸至底部，拟合的平面不能作为投影平面使用，所以对与该类型的碎片首先进行点云分割，将其分割为光滑连续的碎片模型，然后再对分割后的模型进行曲面拟合，并进行边缘点提取.

图4 大型碎片分割示意图

对于较大的碎片在分割中会重新生成多余边缘，如图5(a)所示，这些边缘在点云分割过程中已知，所以只需在碎片边缘轮廓提取完成后删除这些多余边缘即可，如图5(b)所示.

(a)多余边缘示意图 (b)删除多余边缘图5 边缘删除示意图

2.2.1 最小二乘平面拟合

首先对古陶瓷碎片曲面上均匀采样的N+1个点(N≥2)(xi，yi，zi)，(i=0,1,2,…，N-1,N)使用最小二乘法进行平面拟合，拟合平面时，N的取值保证在100个点以内，可通过点云滤波对点云进行抽样简化保证N在合适的取值范围.

设平面的一般的方程表达式为

Ax+By+Cz+D=0,C≠0

(2)

公式(2)可写为

(3)

令a0=-A/C，a1=-B/C，a2=-D/C则式(3)可写为：

z=a0x+a1y+a2

(4)

最小二乘拟合平面的原理就是要使得拟合得到的平面上的点与对应的原始点之间的距离差S最小，其中：

(5)

若要满足S最小就要使得∂S/∂ak=0，其中ak=(a0，a1,a2)，即:

(6)

(7)

2.2.2 点云投影

(8)

(9)

(10)

2.3 边缘轮廓二次优化

在使用边缘系数法提取边缘轮廓过程中，由于为了保证最大数量的提取边缘点而设置边缘系数法的阈值并非最优阈值，导致使用边缘系数法提取的边缘轮廓包含了大量的非边缘点，本文称其为离群点，如图6所示.本文采用空间滤波方法去除离群点，从而获得准确的古陶瓷碎片边缘轮廓.

图6 边缘点与非边缘球形邻域点分布情况

在此定义点云密度ρ，点云密度ρ基于原始点云数据进行计算，以每个点到其最近点PN的欧氏距离D为统计值，遍历点云中的所有点，最终对其求平均值即为点云的密度ρ，具体定义如下式：

(11)

式(11)中：n为初始点云的点的数量，Di为点云中每个点与其最近点之间的距离.

空间半径滤波是点云处理中比较常用且复杂度较低的一种滤波方法，对每个点构建半径为R的邻域，统计落在邻域内的点数，以邻域内的点的数量作为是否要保留该点的评判条件.如图7所示，若设定邻域内点的数量小于4的为离群点，则删除点K3和K2,若设置邻域内点的数量小6，则删除点K1、K2和K 3.本文中空间滤波的对象为初次使用边缘系数法提取的边缘轮廓点云，所以离群点与边缘点的固定半径球形邻域内所包含的邻域点数量差异很大，使用二次空间滤波方法可有效去除这些离群点.

图7 空间半径滤波方法示意图

首先将使用边缘系数提取的边缘轮廓点云作为输入数据，然后对其所有点构建半径为R的球形邻域，边缘点的邻域内所包含的点数量应远大于非边缘点的球形邻域所包含点的数量，通过统计该邻域内的点的数量是否大于阈值Kr来判断该点是否为真边缘点.若该点领域内的数量大于阈值说明该点四周分布有较多的点，该点即为正确的边缘点并保留该点；若邻域内点的数量小于阈值则说明该点四周分布的点较少，不具有边缘轮廓点的特征，并在原有的边缘轮廓上删除该点.

图8(a)为使用边缘系数方法直接提取的古陶瓷碎片的边缘轮廓，图8(b)为经过二次优化后的边缘轮廓，可以看出二次优化可有效的去除首次提取过程中由于各参数选取不当而误提取的离群点.

(a)未优化边缘轮廓 (b)优化后边缘轮廓图8 未优化与优化后边缘轮廓

3 本方法在古文物拼接的应用

古陶瓷碎片拼接是一种薄壁类文物，其断裂面的特征较少难以作为拼接特征，所以碎片的边缘轮廓是古陶瓷碎片的中重要特征之一，基于断裂面轮廓的古陶瓷碎片拼接流程如图9所示.

图9 古陶瓷数字化复原流程图

首先对获取的点云模型进行预处理，然后基于点云进行古陶瓷碎片的边缘轮廓提取，并通过断裂边缘轮廓上的特征点进行碎片拼接，最后通过纹理映射完成古陶瓷的数字化复原工作[5,17-21].

边缘轮廓提取是古陶瓷数字化复原前期的一个重要的步骤，古陶瓷边缘轮廓的精确提取能提升后续数字化复原工作的成功率及复原精度，本文方法是古陶瓷数字化复原的重要一部分，通过实验验证本方法能有效应用于后续的拼接工作.

4 实验结果与分析

4.1 实验平台与数据选择

为了验证本文算法的可行性与可靠性，本算法在表1所示的平台上进行实验验证.

表1 算法实验平台

本文以图10所示瓷碗为实验数据，通过高精度三维扫描仪，获取陶瓷器皿的点云模型，使用前述方法实现对碎片的边缘轮廓的提取，最后以人工建模生成的点云模型对本算法的提取准确率进行测试.

图10 破碎陶瓷碗碎片

4.2 实验参数设置

4.2.1 边缘轮廓提取K近邻K值确定

边缘轮廓提取首先需确定最近邻K值的选取，K的取值决定了后边缘轮廓提取的准确性与高效性.本文算法提出了二次优化，所以在使用边缘系数方法提取边缘轮廓中近邻数K与的选取具有较大的弹性范围.图11是对同一个碎片取不同K值时的边缘轮廓提取结果对比图.表2为各K值提取边缘轮廓消耗的时间，通过实验可得K值的取值范围在20～25时既能准确的提取边缘轮廓，又能保证提取轮廓的时间复杂度更低.

(a)K=30 (b)K=20 (c)K=10图11 不同K值下边缘轮廓提取结果

表2 不同K值提取边缘轮廓耗费时间

由图11与表2可以得出以下结论：对邻域点数量K取不同的值，K的取值越大提取的边缘轮廓包含的非边缘点更少，但同时也丢失了一些细节特征，并且消耗的时间成倍数增长，基于本实验平台，当K的值取20左右时提取的边缘轮廓能保证后期拼接需求且消耗的时间相对更少，所以通过实验得到本文算法中K值取20.

4.2.2 二次优化中邻域半径R与阈值Kr选取

本算法中通过统计点云的平均密度来确定R的选择，半径R的选取既要保证每个点的R邻域内包含有其周围的点，同时R邻域不易构建太大，太大容易导致空间滤波时造成错误滤波，通过实验可得R的取值保证在原始碎片点云密度的3～4倍即可.

阈值Kr的选取基于点云密度与邻域半径R的倍数关系来确定，边缘轮廓更接近线的特征，一点的邻域点仅分布于同一条线上，由于点云密度分布均匀，可认为边缘点的邻域内的点的数量是该点的一侧数量的2倍.如图7所示，点K1、K4为要保留的边缘点，根据前文中点云密度的定义来指导阈值Kr的选取，Kr的值取为半径R与点云密度之间倍数的2～3倍,本文中半径R为点云密度的3～4倍，所以在此选择Kr的值为倍数的2倍数，本文设置阈值Kr的值为6.

4.3 本文算法提取结果

分别使用文献[12]、文献[13]的方法与本文方法对上述的模型进行边缘轮廓提取，提取结果如图12所示，表3为各算法提取边缘轮廓所耗费时间对比.本文在对点云进行边缘提取之前进行了基于古陶瓷碎片几何特征的全局最大化投影平面拟合，将每个邻域内的点投影至拟合平面上，然后再使用改进的边缘系数方法提取边缘轮廓，最后在第一次提取边缘轮廓上通过二次优化来删除误判的边缘点.

(a)原始点云模型

表3 本文算法与其它算法耗费时间对比

通过图12与表3可以得出以下结论：本文不需要对每个点邻域进行平面拟合，有效的减少了时间复杂度，提取的边缘轮廓不存在误判断的离群点，而使用传统基于边缘系数方法提取的边缘轮廓含有较多的离群点，并且丢失了一些不尖锐的特征，影响后期碎片拼接精度与效率.

4.4 边缘提取准确率

本文使用F1-Score作为边缘提取的评价准则：

(12)

式(12)中：Score表示提取的评价分数，分数越高表示提取越准确；Precission指被分类为边缘点中真实边缘点的比例；Recall表示被预测为边缘点总的边缘点的比例.

(13)

(14)

式(13)、(14)中：TP表示正确检测到边缘点的个数，FP表示检测到非边缘点的个数，FN表示未被检测到的边缘点的个数.

使用式(12)对本算法提取的古陶瓷碎片的边缘轮廓的准确性进行评判，为了验证本文方法的准确性，使用一组人工点云模型作为测试数据，如图12中的第三个模型所示，根据模型生成过程确定正确边缘点，然后通本文算法提取该模型的边缘轮廓，最后采用F1-Score对本文与文献[12]、文献[13]的方法进行比较，实验结果对比如表4所示.

表4 本文方与其他文献方法提取准确率

通过图12与表4表明本文算法提取的边缘轮廓准确率更高，本文方法提取的古陶瓷碎片边缘轮廓主要减少了非边缘点判定为边缘点的情况，提取后的边缘轮廓中的离群点较少.同时从表4中的Recall和FP一栏可以看出提取的正确边缘点的数量也是有一定的提升，原因在于其采用了全局最大化平面拟合投影后使用改进的边缘系数法进行边缘提取与二次优化方法，在保证最大化正确提取边缘点的前提下减少了对非边缘点的误判断.

5 结论

本文借鉴已有点云边缘轮廓提取方法，提出了一种基于点云数据的古陶瓷碎片边缘轮廓提取技术.根据古陶瓷碎片几何特征对其进行最大化平面拟合，在拟合平面上使用改进的边缘系数方法来判断一个点是否具有边缘点的特性，最后对提取的边缘轮廓点云进行二次优化.

基于点云的古陶瓷碎片边缘轮廓提取避免了由点云网格化而带来的较高的时间复杂度与局部细节特征丢失，对传统边缘系数中的每个邻域点的作用矢量加以权重系数提高了边缘点提取的准确率，二次优化保证了在完整提取古陶瓷碎片边缘轮廓的前提下，有效地放宽了算法中阈值Fθ、Kr与参数K、R的取值范围，使算法对于类似问题具有更好的适应性.实验表明，该方法能准确的提取古陶瓷碎片边缘轮廓点，并去除误判断的非边缘点，最终精确高效的提取古陶瓷碎片边缘轮廓，为后续碎片拼接工作提供拼接特征，保证古陶瓷数字化复原工作正确、高效的进行.