基于断裂面邻域特征的文物碎片拼接

耿国华，张鹏飞，刘雨萌，周明全，姚文敏，李 康

（西北大学信息科学与技术学院，陕西西安710127）

1 引 言

文物是一个民族的瑰宝，通过文物可以进一步解读历史、还原真相。随着计算机技术、数字化技术的快速发展，借助计算机进行文物虚拟修复已成为研究热点。借助计算机技术实现文物碎片拼接，不仅可以提高文物修复效率，还可以避免对文物造成二次损坏；因此，利用计算机实现文物虚拟修复具有重要的意义。

文物数字化虚拟拼接技术，依据文物的厚度特征，可将文物分为薄壁类文物和非薄壁类文物［1］。壁画、陶器和瓷器等没有厚度特征或厚度特征可以忽略的文物通常归为薄壁类文物；针对薄壁文物的拼接，目前已有大量的研究［2-6］。石碑、秦俑碎片等具有厚度特征的文物通常归为非薄壁类文物，本文选用秦俑碎片模型作为研究对象，对非薄壁类文物自动化碎片拼接展开研究。文物碎片拼接是三维点云配准的一个重要应用领域。点云配准一般分为粗略配准和精细配准两个阶段。粗略配准是将不同坐标系下的点云通过旋转平移操作统一到同一坐标系下，但其配准误差较大，需在此基础上利用迭代算法完成精确配准，使配准误差达到最小。目前，迭代最近点算法（Iterative Closest Point，ICP）［7］及其各种改进算法［8-17］已成为精细匹配的主流算法。樊少荣［18］选取三维碎片的三角网格模型进行实验，分别提取断裂面的外轮廓线和内轮廓线，通过对轮廓线的匹配，从而实现对碎片的拼合。该方法针对完整的断裂面有较好的拼接效果，但针对断裂面存在缺失的碎片，其鲁棒性较差。李姬俊男［19］使用断裂面点体积积分不变量和凹凸互补性得到初始匹配约束簇对，最后通过空间几何一致性约束完成非匹配对之间的消除。李群辉［20］提出了一种基于凹陷凸区的断面匹配算法，当断层被划分为多个凹凸特征区域时，定义了相似的区域对，根据距离原理消除了伪区域对。计算相似区域对的变换参数以实现碎片的粗匹配。最后采用ICP 算法进行精细匹配。Zhao［21］提出了一种基于轮廓曲线和特征区域的块匹配方法，提取断裂面的轮廓曲线，使用改进的迭代最近点算法完成刚性块与质心集的精细匹配，提高了点云中指定厚度的刚性块匹配的精度和速度。袁洁［22］通过提取碎片表面纹理和断裂部位轮廓线特征点，构建碎片断裂面轮廓线并计算双向距离场，最后通过四元数算法和最近点迭代算法完成碎片拼合。陆超华［23］使用断裂面厚度特征信息实现文物碎片匹配。根据厚度特征将片段分组。将厚度特征使用最长公共子序列算法进行匹配，实现了兵马俑碎片的初步匹配。最后利用ICP 算法进行精细匹配。

综上所述，现有方法大多要求碎片模型断裂面信息完整，在断裂部位受损的破损文物模型上的研究较少，而出土的兵马俑碎片大多存在不同程度的破损和损坏，传统的碎片拼接方法拼合误差大且耗时。针对以上问题，本文选用秦俑碎片作为研究对象，首先，定义邻域特征参数，提取特征点集，通过定义曲率特征参数对特征点集进行二次筛选，确定最优特征点集；然后，利用特征点间的相对距离和相对位置矢量与特征主方向的夹角作为特征描述符，依据集合相似理论对特征点进行相似性度量，得到匹配点对集合并采用RANSAC 剔除误匹配点对；最终采用SVD 计算刚体变化矩阵，通过基于K-D 树改进的ICP 算法完成精准拼合。本文算法流程如图1 所示。

图1 本文算法流程图Fig.1 Flow chart of the algorithm in this paper

2 碎片特征点提取

针对秦俑碎片断裂面模型，计算法向量、曲率和法向量内积均值，定义邻域特征参数，提取特征点集；再定义曲率特征参数优化特征点集，筛选出最优特征点集。

2.1 邻域特征参数

法向量和曲率是三维空间中特别重要的几何特征，具有旋转和平移不变性。本文采用主成分分析算法（PCA）进行计算，对采样点pi和其邻域内的所有点进行协方差分析：

将采样点的法向量ni和其邻域内其他点的法向量nij内积后求均值，用于描述邻域曲面的弯曲程度，内积均值δ表示为：

为了描述邻域内各点不规则分布的程度，计算邻域内各点到邻域拟合平面L之间的欧式距离方差S：

其中，dj为邻域内各点到邻域拟合平面L间的欧式距离，如图2 所示。

图2 邻域内各点到邻域拟合平面的距离Fig. 2 Distance from each point in the neighborhood to the neighborhood fitting plane

结合以上几何特征，定义邻域特征参数φ［24］为：

其中：N为采样点的个数，η是为平衡邻域曲率ε和欧式距离方差S相差过大而引入的的平衡参数。由公式（7）可知，欧式距离方差S(pi)和邻域曲率ε(pi)越大，内积均值δ(pi)越小，邻域特征参数φ(pi)就越大，采样点pi为特征点的概率越高。设定邻域特征参数阈值σ，若φ(pi)＞σ，即为特征点；否则为非特征点。

2.2 曲率特征参数

通过2.1 节，可提取出能够反映文物碎片断裂面几何特征的特征点集，但由于邻域曲率是受估算值的影响，提取出的特征点集中会存在非特征点，因此定义曲率特征参数对特征点集进行二次特征点提取，得到最优特征点集。曲率特征参数定义为

其中：k1和k2分别表示采样点pi处的最大主曲率和最小主曲率。根据公式（9），若采样点pi的曲率特征参数ξ(pi)大于其邻域内任意一点的曲率特征参数，则采样点pi为邻域凸点；若采样点pi的曲率特征参数ξ(pi)小于其邻域内任意一点的曲率特征参数，则采样点pi为邻域凹点；相比于邻域曲率，曲率特征参数可以很好地描述采样点邻域表面的凹凸性和弯曲程度。

本文选取兵马俑一号坑G10-18 号佣部分碎片进行实验，设定参数σ=0.4，图3（a）为碎片实物图，由图3（a）可以看出，断裂面受到一定程度的损坏，但本文算法有较好的鲁棒性，可以很好地完成碎片拼接。图3（b）为碎片特征点提取图，红色点为算法自动标识的特征点，可以看出，比较平缓的区域，特征信息较少，特征点分布较为疏散，而凹凸不平的区域，特征信息较多，特征点分布较为集中。

图3 G10-18 号佣特征点提取Fig. 3 Feature point extraction of G10-18

3 特征描述及相似性度量

针对上述过程提取出的秦俑碎片模型断裂面上的特征点，为其定义有效的特征描述，是拼接任务的关键。本文采用相对距离与方向夹角的集合作为特征描述符，并依据集合相似理论度量特征点间的相似性，确定匹配特征点对集合。

3.1 特征描述符定义

假设特征点集合Qkey={q1，q2，…，qm}，m为特征点的个数，qi表示第i个特征点，本文选用特征点间相对距离d和相对位置矢量与特征主方向的夹角θ作为特征的描述符［25］，特征主方向定义为特征点的散度矩阵的最大特征值对应的特征向量。设特征点qi与其余m-1 个特征点的相对距离为dij，夹角为θij，其中j∈{1，2，…，i-1，i+1，…，m}，假设共有六个特征点，其中q3的特征描述如图4 所示。特征点qi的特征描述符Fi可表示为：

图4 特征描述符的可视化Fig. 4 Visualization of feature descriptors

3.2 特征点匹配

选取两特征点集合Pkey和Qkey，对于每个特征点的描述，定义相似性度量指标为：

Step1 取特征点集Pkey和Qkey中的每一个特征点pi和qj；

Step2 计算pi与Qkey中的每个特征点qj的特征描述集合的相似性，得到相似性最高的特征点qj和相似性次高的特征点qk；

Step3 若满足式（12），则将匹配的特征点对(pi，qj)加入集合Φ。

3.3 剔除误匹配点对

通过3.2 小节得到特征点对集合Φ，但集合Φ中往往会存在误匹配特征点对，如图5 所示。剔除误匹配特征点对可以有效地提高拼接效率和精度。本文采用随机抽样一致性算法对误匹配点对进行剔除，算法思路如下：

首先，随机选取n个特征点对，用于评估样本集的参数模型即点云变换矩阵；将特征点对中源点云pi经过变换矩阵计算后，得目标点云t，计算目标点云t与匹配点云qi的偏差D：

若偏差D小于给定阈值d，则认为此匹配点对为样本内点并保存；以上步骤迭代m次，取样本内点最多的点对作为最优匹配集。算法示意图如图6 所示，（a）图为原始数据集，（b）图中红点表示样本内点，绿点表示样本外点（彩图见期刊电子版）。

图5 剔除误匹配点对Fig. 5 Eliminating mismatched point pairs

图6 随机抽样一致性算法原理示意图Fig. 6 Schematic diagram of random sampling consistency algorithm

4 碎片拼合

4.1 粗略匹配

构建方差矩阵E3×3：

其中：Σ为E3×3特征值组成的对角阵，旋转矩阵R和平移矩阵T为：

通过式（14）～式（17）得到旋转矩阵R和平移矩阵T，从而得到点云初始位置关系，完成碎片初始匹配。

4.2 精确匹配

本文基于K-D 树对ICP 算法进行改进，采用改进的ICP 算法进行精确配准，利用奇异值分解算法计算迭代过程中的刚体变换矩阵，并通过KD 树快速搜索最近点对，提高匹配效率。

4.2.1 ICP 算法

ICP 算法是点云精确配准中常用的方法，其基本原理是寻找最优矩阵参数R和T，使得误差函数E最小：

其中：n为匹配点对个数，pi和qi分别为两匹配点，R和T分别为旋转矩阵和平移向量，ICP 算法的具体步骤：

Step1 获取待配准点云集合P和Q；

Step2 搜索寻找最近点对集合；

Step3 计算刚体变换矩阵R和T，并计算误差函数E；

Step4 将配准点云Pk变化到新坐标系下Pk+1：

Step5 计算点对间的欧式平均距离：

Step6 当dˉk-dˉk+1＜Θ 时迭代结束；否则重新迭代。

4.2.2 ICP 算法改进

传统的ICP 算法通过欧式距离判断最近点，对于海量点云数据耗时且效率低，因此本文对ICP 算法做一个改进，利用K-D tree 能快速搜索最近点对，提高点云配准效率。其搜索步骤如下：

Step1 将待查询结点与当前树结点进行比较，若小于等于当前树结点，则进入左子树分支；否则进入右子树分支，沿搜索路径找到最近邻近似点；

Step2 对点进行“回溯”操作，若在搜索路径上有更近的其他子空间点，则跳到子空间结点去搜索距离最近点；

Step3 重复Step1 和Step2 直到搜索路径为空结束。

5 实验结果与分析

本文选取兵马俑一号坑出土的秦俑碎片模型作为实验对象来验证算法的可行性，本实验在Visual Studio 2015 和Open GL 环境下进行编程，在Intel Core i7-3770CPU/3.4 GHz，内存32 GB的PC 机上完成，本文实验所涉及到的阈值参数均由经验和数据统计分析得出σ=0.4，d=0.1 mm，τ=0.5，τ1=0.5，τ2=0.5。

5.1 本文算法实验结果展示

本文选取一号坑兵马俑G10-12 和G10-41 秦俑碎片模型进行拼接效果展示，如图7～图8 所示。因受自然因素和人为因素的影响，兵马俑碎片大多存在缺损。针对轻微损坏的模型，本文算法可以很好地实现拼接，如图7 所示；针对严重缺损的模型，如图8 所示，本文算法也有较强的鲁棒性和准确性，拼接结果基本不受其影响。

图7 G10-12 号佣部分碎片拼合结果Fig. 7 Results of fragment assembly of G10-12

图8 G10-41 号佣部分碎片拼合结果Fig. 8 Results of fragment assembly of G10-41

为了验证本文算法的有效性和可应用性，选取G10-22 和G10-57 完整秦俑模型进行拼接实验，如图9～图10 所示。结果表明，本文算法可以应用于完整佣体的拼接，并有较强的鲁棒性，如图10 所示，通过将破损局部放大可以看出，佣体存在严重破损，但本文算法仍可以很好地完成拼合。

5.2 实验对比与分析

5.2.1 不同阈值参数分析比较

图9 G10-57 号佣正面与背面拼合结果Fig. 9 Assembling results of front and back of G10-57

图10 G10-22 号佣正面拼合结果及缺损部位展示Fig. 10 Display of the results and defects of G10-22

特征点匹配是文物碎片拼接过程中的关键步骤，能否准确地找出正确匹配点对，直接影响文物拼接结果，因此τ1和τ2的取值至关重要。为了更准确地提取匹配点对，引入阈值τ1和τ2对其进行限制，τ1表示与之匹配的特征点中最相似和次相似特征点之间的差距，τ2表示与之匹配的特征点中最相似特征点的相似程度。分别令τ1，τ2=0.4，0.5，0.6，本文算法收敛的评判标准是拼合误差小于等于1 mm，为验证算法的稳定性，当拼合误差首次小于1 mm 后，再迭代3 次，若每次误差均小于1 mm，则停止迭代，认为算法稳定并已达到收敛。如图11 所示。

图11 不同阈值参数分析比较Fig. 11 Analysis and comparison of different threshold parameters

从图11 可以看出，当τ1，τ2=0.5 时，拼合误差最小，当τ1，τ2=0.4 时，拼合误差较大，这是由于最相似和次相似点的差距小，相似性高的原因导致的；当τ1，τ2=0.6 时，虽然可以更快拟合，但匹配误差也偏高，这是由于最相似和次相似点的差距大，匹配对更少的原因导致的；综上，本文选取τ1，τ2=0.5 进行实验。

5.2.2 不同方法分析比较

为了验证本文算法的优越性，与文献［18］中的算法进行实验结果对比，本文与文献［18］对比的原因如下：（1）文献［18］首次提出了基于断裂面轮廓线的文物碎片拼接方法，此方法具有一定的先进性。（2）文献［18］与本文均是基于文物碎片断裂面进行碎片拼接。结果如图12 所示。

图12 G10-61 部分碎片拼合结果Fig. 12 Results of partial fragment assembly of G10-61

从视觉角度分析，由图12 可以看出，本文算法的拼接效果较好，优于文献［18］方法的结果，从图12（b）可以看到，文献［18］方法拼接过程出现了明显的错位现象。这是由于文献［18］方法对于断裂部位存在缺失的情况，提取轮廓线出现偏差，从而导致位置匹配不准确。

为了避免视觉分析的主观性，本文还从定量角度进行分析。由于在计算拼合误差过程中，点云A 中可能存在多个不同点对应点云B 中同一点产生的误差，因此本文重新定义拼合误差。假设集合A 和B 分别为拼合后的两断裂面，定义集合A 中某点ai到集合B 的距离为点ai到集合B 中所有点的距离的最小值，计算集合A 中所有点到集合B 的距离的平均值EˉA，再计算集合B 到集合A 的距离的平均值EˉB，定义断裂面A 和B 的拼合误差为E，则：

本文的算法运行时间的评价标准如下：通过多次迭代不断优化其拼合误差，当拼合误差小于1 mm 时，继续迭代3 次，若每次误差都小于1 mm，则停止迭代，计算运行时间，此时间为算法的运行时间，通过对比算法的运行时间，从而可以评判算法的效率。

根据上述拼合误差和算法耗时的定义，分别对本文算法与文献［18］算法计算了拼合误差和算法耗时，结果如表1 所示。

表1 不同算法的性能对比Tab.1 Comparison of performances of different algorithms

从表1 中可以看出，本文算法在时间和拼合误差上都优于传统方法。在时间上，比传统算法提高了26%～40%；在拼合误差上，比传统算法减少了4%～25%。

为了展示本文算法的收敛速度，选取G10-12秦俑碎片为实验对象，选用文献［18］方法作为对比方法，对算法运行过程进行对比分析，分别计算了不同时间下的拼合误差，结果如图13 所示。本文算法的收敛速度明显优于文献［18］，在性能上有较大的提升。

在精确匹配阶段，为了验证本文改进的ICP算法效率更高，我们与传统ICP 算法进行了实验对比，选取G10-40，G10-80 和G10-87 号文物碎片模型进行实验，分别计算了两种方法的拼合误差和耗时，从表3 中可以看出，与传统ICP 相比，本文改进的ICP 算法在精度和效率上都有明显的提升，并且在点云数目越大时显著性越明显；因此利用K-D tree 改进的ICP 算法可以明显的提高配准效率。结果如表2 所示。

图13 G10-22 运行时间与匹配误差分析图Fig.13 G10-22 Running time and recombination error analysis chart

表2 ICP 算法性能对比Tab.2 Comparison of performance of ICP algorithm

6 结 论

本文提出一种基于断裂面邻域特征的碎片拼接方法，针对受损文物碎片模型，有较强的鲁棒性，并且针对稀疏点云特征难以匹配的问题，本文也给出了有效的解决方法。采用邻域特征参数和曲率特征参数两个约束对特征点二次提取；定义特征点间的相对距离和相对位置矢量与特征主方向的夹角作为特征描述符，通过集合相似理论进行度量，得到匹配点对集合，并采用RANSAC 剔除误匹配点；然后采用SVD 计算得到刚体变化矩阵，最后采用基于K-D 树改进的迭代最近点算法进行精确匹配，完成拼合。实验表明，与传统算法相比，本文算法在精度和效率上都有明显提升。

与以往的文物碎片拼接方法相比，本文算法的优势在于：（1）在特征点提取阶段，利用最小二乘法原理构造曲率特征参数，对特征点集进行优化，提高配准精度。（2）在特征点对匹配阶段，引入RANSAC 算法删除误匹配点对，筛选出最优匹配集，提高配准精度。（3）在精确匹配阶段，引入K-D tree 快速搜索最近点对，提高点云配准效率。

尽管本文算法可以在破损文物碎片模型上取得良好的效果，但现有文物数据库中，部分文物断裂面大面积破损；因此，寻找一种更稳健的特征描述符进行文物碎片拼接是我们下一步的研究方向。