袁洁 周明全 耿国华 张雨禾

三维破碎物体的匹配与拼接是计算机图形学、模式识别、可视化技术等众多领域里一个颇具挑战的问题.20世纪末,随着计算机技术的快速发展,将计算机技术引入到文物碎片的虚拟拼接过程中,提高了文物复原的效率,对文物保护和复原至今有着重要的意义.

文物数字化虚拟拼接技术,根据复原过程中是否有专家参与可分为自动复原方法和交互式复原方法.其中,现有的刚体匹配算法主要结合破损文物断裂区域的几何特征,作为判断相邻碎片之间相似性的主要依据,从而实现破损文物的完整拼合[1].

根据文物的厚度信息,自动化虚拟复原技术可大致分为两类:薄壁类和非薄壁类文物碎片.早期薄壁类破损文物的虚拟拼接,主要以提取断裂面轮廓线,并将其映射成为二维平面曲线,实现碎片拼合[1−6].例如Cooper等[1]通过结合碎片轮廓线及其法向,以最大似然函数为基础进行自底向上的搜索,完成碎片模型的自动拼合.樊少荣等[2]区分三角网格曲面模型的外表面与断裂面从而正确提取出内轮廓线和外轮廓线,实现碎片的精确拼合.Oxholm等[3]将断裂面轮廓线上顶点的曲率、挠率以及颜色信息进行结合,组成特征属性串,采用最长公共子串的方法实现轮廓线间的匹配.Willi等[4]将贝叶斯分析应用在半自动拼合方法中,该算法对断裂面光滑平整且轴对称的三维模型有很好的匹配效果,其局限性是只能作用于有限形状的三维模型.Zhang等[5]基于模板匹配的思想确定颅骨碎片的位置关系,再结合相邻碎片轮廓曲线实现碎片拼接.Huang等[6]通过计算模型表面顶点的积分不变量值提取碎块表面尖锐的边缘线,基于向前搜索技术和表面一致性的约束方法实现碎片的两两拼合.空间轮廓曲线匹配方法在几何特征的数值化计算时具有简单高效的特点;但由于其采样点数量有限,因此这类方法对噪声较敏感并且易忽略厚度信息的部分文物碎片模型效果较差.

非薄壁类刚体匹配对文物重建、数字化遗产保护有着重要的意义[7].Papaioannou等[8]采用Z缓冲方法,获取三维模型断裂面投影,计算当前位置断裂区域的“位置误差”,通过最小化误差获得最优匹配,实现虚拟复原,该方法适用于雕塑、纪念碑等大体积破碎模型的修复.李姬俊男等[9]基于相邻断裂面间的凹凸互补性,通过提取约束性的特征簇完成碎片间的两两拼合.Sahner等[10]通过计算断裂区域顶点的坐标和法向,采用层次聚类方法进行破损文物碎块之间的拼合.由于断裂面中包含了较多的特征信息,能够准确地反映模型断裂区域的邻接关系,因此,针对非薄壁类文物碎片,基于空间曲面的碎片拼合,具有一定优势,但针对断裂面部位受损较严重的碎块,匹配的正确性无法保证.

因此针对破损文物断裂部位边缘受损而引起的轮廓线不能充分表示断裂面几何信息的问题,本文基于文献[11]中的特征线提取方法,提出一种基于Morse-Smale的断裂面拓扑几何特征的破碎刚体自动拼接方法.本文算法包括断裂面拓扑图的生成、四边形描述符的定义、匹配集的确定以及碎片拼合4个主要部分,步骤如图1所示.本文算法的核心思想是采用曲面四边形描述符表示断裂面几何特征信息.在传统的基于空间曲线的碎片拼合方法中,针对断裂部位边缘受损而提取的轮廓线不精确,错误的曲线匹配导致断裂面渗透;传统的基于空间曲面的碎片拼合方法,并没有考虑断裂面特征点之间的拓扑关系,可能出现断裂面的错误配准.因此,为了提高碎片拼合的效果,基于Morse-Smale复形[12]本身的四边形性质及其可控性和鲁棒性,本文借助能有效融合断裂面上凹凸信息的四边形曲面,更好地反映空间曲面特征,并且可精确地确定邻接碎片的匹配关系,同时采用凹凸互补阈值方法筛选出最优匹配集,将计算量较大的空间点对应关系转换成少数的四边形匹配,有效地提高了刚体转换的效率,最后采用穷举搜索的方法实现破损碎片精确拼接.

图1 本文方法步骤Fig.1 Procedures for the method proposed in the paper

1 拓扑图提取及相关术语

本文通过曲面四边形表示断裂面上的几何特征,从而实现碎片拼合,因此,本节提取能有效融合断裂面几何特征的拓扑图.以三角网格模型为研究对象,首先,依据顶点的显著度指标函数,提取关键点并分类;然后,采用最大角度法构建特征线;最后,对提取的断裂面Morse-Smale复形进行简化.

1.1 Morse-Smale理论

Morse理论主要用来研究n维空间内的d维流行M(M⊆Rn)的形状,通过Morse函数,将二维流行表面分割为一系列Morse单元,而Morse-Smale复形为这些莫斯单元的集合.其主要思想是将形状的拓扑研究同映射函数集合特性的定量分析结合起来,通过定义特征点并判断其类型,按照一定的寻径方法利用特征线将其连接起来,形成关键网络或转换成其他数据模型,从而表达流形表面的拓扑结构与形态特征,是对特征提取和几何形状的分析,保存其拓扑性质及形态的强大工具.

拓扑特征是对图像中区域结构形状的总体描述,其特点是不受旋转与平移等变形影响,描述能力强、计算速度快的优点,可有效、准确地表示物体模型或图像的局部形状.常见的图像拓扑特征有基于骨架结构表示的目标检测算法和基于形态图表示的三维识别算法.

基于Morse函数的特征提取算法在由二维图像扩展到三模型表面时,依据三角网格顶点的显著度度量指标(曲率值、法向量改变值等顶点属性)提取的极值点、升降线与Morse-Smale复形做为模型表面的显著特征,避免了没有实际意义的非显著特征提取,保证了Morse-Smale复形的拓扑完整性.

1.2 特征点提取

针对三维模型表面的网格数据,特征提取是指采用指标函数来提取出具有明显特征的点集合.本文采用曲度函数f(v)作为三角网格上各顶点的显著度指标函数[13],定义如下:

其中,v为三角网格上任意一个顶点;k2max(v)和k2min(v)分别为顶点v的最大和最小主曲率.

采用曲度函数作为指标函数是为了将曲率绝对值最大的特征点提取出来.断裂面上曲率绝对值最大的点即为凹点、凸点和鞍点,均是能够表示物体模型表面几何特征的点.

1.3 特征线构建

基于上节提取的断裂面特征点,本节通过最大角度法[14]构建特征线,构造能完整表示断裂面几何关系的拓扑图.

首先,对提取的特征点进行判定,并对其进行分类.定义vi(i=1,2,···,n)为兵马俑碎块断裂面的三角网格上的任意顶点,其中,n为三角网格数,其所对应的显著度指标函数为f(vi).若其一阶邻域点的显著度函数值均小于该顶点的显著度函数值,则定义该顶点为极大值点;同理,若其一阶邻域点的显著度函数值均大于该顶点的显著度函数值,该顶点则被定义为极小值点.若目标点的邻域中有某段连续区域内的顶点,其指标函数值均大于目标点的函数值,则称该段连续区域能取得的最大范围为上升区间;若其指标函数值均小于目标点的函数值,则称该段连续区域能取得的最大范围为下降区间.定义鞍点是同时连接两个上升区间和两个下降区间的特征点.如图2所示,分别对应断裂面上曲率绝对值最大的顶点,即特征点.

图2 顶点邻域关系图Fig.2 Diagrams of vertex neighborhoods

基于以上分析,本文采用最大角度法构建特征线,构建规则如下:对任意一个鞍点都存在2条沿着指标函数值下降至极小值点和2条沿着指标函数值上升到极大值点的特征线,计算鞍点与其邻域顶点的边,选择上升或下降最快的点作为下一个路径点位,构成特征线.以升弧为例,以鞍点为出发点,方向为鞍点上升区间中指标函数值最大的邻点所指示的方向,直至到达极大值点.设vi为三角网格上的任意顶点,其方向为则有:

式中,vj为vi的一阶邻域点.降弧构造方式与升弧类似,沿着最大角度反方向搜索.

1.4 拓扑简化

针对所提取的断裂面拓扑图,本节采用注销和移动操作[15]进行简化,从而获取准确的纯四边形拓扑图.

Morse-Smale复形中的特征点以及升弧、降弧构成初始特征集合.一般情况,特征线会将模型表面划分成若干四边形区域,所有四边形都是由两个鞍点,一个极大值点和一个极小值点构成.在现实情况中,存在一些特殊问题,如图3(a)深色弧线部分生成的环状,对某个鞍点在上升空间中寻找极大值点时,两条升弧同时连接到一个极大值点,使得拓扑图中出现环状回路,以及中心区域可能出现孤立极值点,如图3(c).

图3 环状和孤立临界点Fig.3 Diagrams of cyclic and solitary critical points

综上所述,针对拓扑图中环状回路的问题,本文采用注销的方法消除多余或者不规则的关键点.注销是一种双重边缘收缩的方法,即去除鞍点与极值点之间连成的两条重复特征线,其中包含所有的弧线以及这个鞍点和两个胞腔,结果如图3(b).若鞍点v所连接的两条降弧和鞍点x的极小值点重复,而特征线均为升弧,因此将这些升弧合并为一条升弧,且删除重复的两条降线,简化结果如图3(d).

针对初始特征集合,某些鞍点在特征线构建过程中,搜索下一个路径点位时可能存在微小误差,线段间误差积累,导致构建的初始莫尔斯复形中存在一些次要特征线,如图4所示,鞍点的下一个最优路径点位可能是沿着图4(c)中虚线所示的方向,而实际在构建过程中是沿着图4(a)或4(b)中三角形边方向.

图4 路径点位选择Fig.4 Paths selected by points

文献[15]中将这些次要特征线及其所连接的临界点一并删除,但是针对模型表面的几何特征,删除误差积累导致的次要特征线可能会造成模型表面特征的不完整性,因此本文基于此方法进行改进,简化规则如下:1)通过定义特征线的显著度,确定其重要性程度;2)设定阈值δ0,当特征线的显著度低于该阈值时,定义为次要特征线;3)针对次要特征线进行移动操作从而寻找最准确的特征线.

一般地,升弧上的指标函数值与其所在区域的指标函数差值越大,则该升弧越重要,显著度越高;降弧上的函数值与其所在区域的函数值差值越小,则该降弧越重要,显著度越高[16].因此,本文定义特征线的显著度为该特征线与其邻接区域的平均指标函数值之差.

基于以上分析,每个鞍点连接2条升弧和2条降弧,本文将其合并为一条升弧R和一条降弧r,假设降弧r,其所对应的上升域为D(m1)和D(m2);升弧R,对应的下降域为A(m1)和A(m2).则升弧R相对于A(mi)的显著度为

其中,i=1,2.

最大角度法构建特征线过程中,特征线方向沿着三角网格上三角形边前进.因此,每条特征线都是由若干个特征点pi(i=1,2,···,n)连接而成的,即升弧R是由若干首尾相连的线段组成,如图5所示.

图5 特征线构成Fig.5 Feature line composition

因此,R上的积分可近似表示为

式中,pj1,pj2是特征线R上的第j条线段上的2个端点,区域A(Mi)上的积分则可近似表示为

式中,pk1,pk2以及pk3分别为区域A(Mi)中第k个三角形的顶点;Ak为该三角形的面积.

通常,对任意的特征线都有两个邻接区域,因此会相应的有两个显著度.假设特征线R的2个邻域分别为A(M1)和A(M2),相对应的2个显著度为δ(R,A(M1))和δ(R,A(M2)),取较小值为该特征线的整体显著度:

同理,降弧r的两个邻接区域为D(m1)和D(m2),则r的显著度为

其中,i=1,2.则该条降弧的整体显著度为

对任意特征线,若其特征线显著度低于所给定的阈值,则将该特征线进行移动操作,即从该特征线的起始鞍点开始,重新搜索下一个最优的路径点位,从而获取该区域内最准确并且显著度最高的特征线,构造断裂面上的纯四边形拓扑图.

2 特征描述符定义及搜索

为了将拓扑图中曲面四边形构造成能够有效表示断裂面几何信息的四边形描述符,首先,定义基准点和0值面,采用顶点高度差的方法对曲面四边形进行编码,获得四边形描述符;然后,通过凹凸互补误差筛选出最优四边形匹配集.

2.1 四边形定义

断裂面上的Morse-Smale复形,由多个四边形曲面构成,这些曲面四边形将断裂面进行剖分,四边形本身具有各向异性的特点,使得四边形网格更易与模型特征对应,因此本文将断裂面上的曲面四边形作为特征描述符,以相邻碎片断裂面间的凹凸特性为基准进行碎块拼合.设待匹配的两个断裂面分别为ϕA和ϕB,定义ϕA为源面,ϕB为目标面,对于源面ϕA上的Morse-Smale复形中的任意曲面四边形,假设两个鞍点分别为s1和s2,极大值为M1,极小值为m1.

如图6所示,设以极大值点M1所在的切平面l0为0值面,定义极大值点M1为基准点,鞍点s1、s2,极小值点m1为目标点.用h1表示鞍点s1到该切平面的垂直距离,极小值点m1到该切平面的垂直距离为h2,鞍点s2的对应距离为h3,极大值点的对应距离为h0,将该高度差值作为每个顶点的属性值,并按照一定规则进行编码.

图6 顶点高度差值Fig.6 Vertex height difference

编码规则:本文以源面ϕA上任意曲面四边形为例,以四边形中极大值点开始,顺时针编码,如图6所示的四边形中,极大值点M1为基准点,则定义hA0=0;依次鞍点s1、极小值点m1、鞍点s2所对应的属性值为hA1,hA2,hA3,则该四边形的一个属性编码为.

类似地,目标面ϕB在选取基准点时,根据匹配碎块断裂面的凹凸互补性,基准点选择特征四边形中的极小值点m1,以其切平面作为0值面,然后求得其余三个目标点到0直面的垂直距离,以同一编码规则获得任意特征四边形的属性编码.

2.2 最优匹配集搜索

基于本文定义的四边形属性,从破损的兵马俑断裂面上提取出特征四边形集,如图7所示,但由于特征四边形的数目较多,断裂面ϕA和ϕB之间的匹配关系会比较复杂.因此,本文设置凹凸互补阈值条件,减少误匹配四边形集合个数,求取少量的最优匹配集合.

图7 断裂面拓扑图Fig.7 Topological diagram of fracture surface

假设ϕA和ϕB为两个待匹配的断裂面,和分别为上的初始特征四边形集合.设sAi和sBi是待匹配的两个断裂面ϕA和ϕB上相对应的可匹配特征四边形,则定义凹凸互补误差为

假设四边形sA1为断裂面ϕA上的任意特征曲面四边形,属性编码如下:

相对应的,断裂面ϕB上对应的特征四边形sB1为

则凹凸互补误差公式为

其中,x1=hA0−hB0,x2=hA1−hB2,x3=hA3−hB3,x4=hA2−hB1.当ε越接近0,表示该对特征四边形的匹配度越高.因此,设定凹凸互补阈值ε0＞0,对四边形sA1,若在断裂面ϕB上满足ε(sA1,sBi)＜ε0的特征四边形存在m个,则取最优匹配对 (sA1,sB1)=, 并将该最优匹配对加入到最优匹配集s=;否则,删除该匹配对.本文通过凹凸互补阈值方法约束特征四边形,筛选出匹配度较高的匹配集,提高断裂面的匹配精度并且降低刚体转换的复杂性.图8是初始特征匹配到最优特征集匹配的结果对比图.

图8 匹配集优化结果Fig.8 Optimization results of matching set

3 刚体变换

对特征四边形进行初步筛选后,确定了待匹配断裂面上的最优特征集合S′A和,以及由此构成的最优匹配集S=, 但由于两个碎片之间的方向与位置均不同,为了能将碎块完整的拼合,需根据两个四边形匹配集计算出对齐碎片的旋转矩阵R和平移矩阵T[17].

基于特征四边形区域的质心坐标,采用四元组方法[18]求解出最优匹配集中的3组质心不共线特征四边形匹配对的旋转矩阵R和平移矩阵T,并采用穷举搜索法[19]筛选出匹配正确率较高的刚体变换,实现碎片的精确拼合,算法步骤如下:

步骤1.设任意特征四边形S的4个顶点为{p1,p2,p3,p4},点pi的三维坐标为(xi,yi,zi),定义S的质心为.

步骤2.根据三角形理论,从匹配集中筛选出所有满足该条件且质心不共线的3组特征四边形对,构成集合,其中,,设M中的元素个数为e个.

步骤3.从上述集合M中选取任意一组特征四边形对,采用四元组法计算出断裂面ϕA和ϕB之间的刚体变换矩阵R和T,遍历集合M中所有的特征对元素,求得刚体变换数组D.

步骤4.计算D[t],t=1,···,e的匹配正确率;c为正确匹配对个数,初值设为0;集合S中任意的特征四边形对经过D[t]变换之后,若满足,则表明该特征四边形对匹配正确,c=c+1,否则为错误匹配;遍历集合S可求出D[t]的匹配正确率q=c/m;遍历整个数组D[e],计算所有特征对刚体变换的正确匹配率.

步骤5.选取匹配正确率最高的刚体变换,完成待拼接碎片的精确对齐.

4 实验结果与分析

为检验本文算法的有效性,以秦始皇陵1号坑出土的部分兵马俑碎片网格模型作为实验数据.本文实验采用Visual C++(Visual studio 2010)和 OpenGL编程,在 Intel Core i7-3770 CPU/3.4GHz,4GB内存的PC机上实现.实验阈值均是通过经验值以及数据统计获取的,δ0=0.02,ε0=5.0mm,σ1=2.0mm,碎片在断裂处均存在边界受损,几何特征缺失现象.

4.1 本文算法结果

图9～11所示为G10-26号俑、G10-19号俑和G10-23号俑采用本文算法在部分邻接碎片上的实验结果,可以看出,本文算法在边缘受损的碎片模型上可取得较好的拼合结果.图12～14所示为本文算法在一组整俑上的拼接结果,实验结果表明,本文算法在断裂部位边缘受损的整俑模型上,如图12中矩形框标注的有小部分缺失的碎片,仍可取得较好的拼合结果.本文实验数据采用的兵马俑模型均为陶俑模型,碎片断裂面具有一定厚度,碎片数量较多并且网格数据总量较大,因此在邻接碎片的搜索时复杂度较高.实验数据如表1所示,邻接碎片在断裂面提取的特征四边形随碎片网格数增加,碎片的数量决定了邻接碎片的断裂面个数,其对特征四边形个数也有直接的影响;由表1中数据可得,本文算法可取得较高的匹配精度.

图9 G10-26号俑部分邻接碎片拼接结果Fig.9 Reassembly results of some adjacent fragments of Warriors G10-26

图10 G10-19号俑部分邻接碎片拼接结果Fig.10 Reassembly results of some adjacent fragments of Warriors G10-19

图11 G10-23号俑部分邻接碎片拼接结果Fig.11 Reassembly results of some adjacent fragments of Warriors G10-23

图12 G10-18号俑部分邻接碎片拼接结果Fig.12 Reassembly results of some adjacent fragments of Warriors G10-18

本文将图7中的碎片裁剪为薄壁模型后,断裂面上的特征四边形提取结果示意图如图15,有效特征四边形定义为四边均处于碎片模型断裂面上的四边形,如图15中矩形框标注所示.由于薄壁模型断裂面上的特征点数量较少,提取的有效特征四边形数量不足从而无法完整表示断裂面的几何结构,造成特征匹配不准确的问题,从而影响后续邻接碎片位置的判断.因此,本文算法针对薄壁类文物碎片模型,拼合效果较差.

表1 兵马俑碎片实验数据Table 1 Experimental datas of the Terracotta Army fragments

图13 G10-36号俑部分邻接碎片拼接结果Fig.13 Reassembly results of some adjacent fragments of Warriors G10-36

4.2 算法运行时间

本文算法针对成组的兵马俑碎片拼合时间如表2所示,其中,断裂面曲面四边形特征提取为T1,本文每次提取2个碎片的断裂面四边形特征;特征匹配集获取为T2,同理,每次获取2个碎片的特征匹配集;T3为碎片拼合.实验数据表明,随着碎片数据的增加,运行时间也会随之增加,这是因为断裂面上四边形特征描述符提取为本文算法最耗时的部分,需要遍历断裂面网格上的所有顶点.

图14 G10-22号俑部分邻接碎片拼接结果Fig.14 Reassembly results of some adjacent fragments of Warriors G10-22

图15 薄壁碎片断裂面特征提取图Fig.15 Feature extraction in thin fragments

4.3 对比实验结果

图16所示分别为采用基于形态图方法[20]及基于断裂面几何特征方法[7]的碎片拼合结果,选取3组碎片模型为实验数据,分别为:G10-18号俑,G10-19号俑及G10-36号俑部分邻接碎片.该3组邻接碎片具有不同的特点:G10-18号俑碎片模型断裂部位边缘受损较少;G10-19号俑碎片模型断裂部位边缘受损较严重且断裂面存在缺损;G10-36号俑邻接碎片断裂部位面积较大,碎片存在缺失.

表2 本文算法运行时间Table 2 Execute times of the proposed algorithm

若两个断裂面可实现正确匹配,那么当两个碎片完成精确配准之后,相邻的两个断裂面之间不应该发生明显的碰撞,存在缝隙过大以及相互重叠交错[21].图16(a)所示为G10-18号俑碎片采用文献[20]中提出的基于形态图方法的拼合结果,将碎片模型转换成点云数据,可看出,当碎片的断裂处存在缺损时,基于断裂面出轮廓线的拼合方法会出现相互渗透的现象,如图16(a)中矩形框中标注所示,难以达到较好的拼合结果;图16(b)中所示为采用文献[7]的基于断裂面特征簇的拼合方法对G10-19号俑部分碎片的拼合结果,由于碎片断裂面存在缺损,则导致几何信息无法有效表示断裂面结构,造成断裂面相互渗透现象,如图16(b)中矩形框标注所示;当碎片存在小部分的缺损时,则会造成不准确的位置匹配,出现碎片表面纹饰信息不连续的现象,如图16(c)所示.

实验结果表明,相较于其他拓扑特征的拼合方法和基于断裂面几何特征的拼合方法,本文算法通过融合断裂面几何特征的四边形曲面,增强了特征描述符的约束性,能够准确有效地判断碎片的邻接位置,针对断裂面边缘受损的碎片具有更好的拼合效果.

5 结语

针对破损文物断裂部位边缘受损而引起的轮廓线不能充分表示断裂面几何特征的问题,本文提出了一种基于Morse-Smale的断裂面拓扑特征的破损文物自动拼接算法.该算法将断裂面上特征点的曲度函数作为特征检测算法,根据断裂面几何特征提取特征点;采用最大角度法对特征点进行特征线构建,并计算特征线重要度简化拓扑图,得到断裂面的纯四边形拓扑图;将拓扑图中四边形曲面构造成为能有效表示断裂面集合特征的描述符,弥补了传统方法反映局部特征的不足,提高了匹配的精度和效率,同时扩展了匹配算法的使用范围.本文算法需要较为精确地识别出碎片的断裂面,解决薄壁类破损碎片断裂面特征平滑,无法取得精确地拼合结果的难题.

由于本文算法采用了拓扑图四边形为特征描述符,因此在提取特征描述符时较耗时,寻找一种更加鲁棒、高效地特征提取算法以减少计算量,将是我们下一步工作的重点.

图16 传统算法实验结果Fig.16 Experimental results of traditional methods