基于斷裂面鄰域特征的文物碎片拼接

耿國華，張鵬飛，劉雨萌，周明全，姚文敏，李 康

（西北大學(xué)信息科學(xué)與技術(shù)學(xué)院，陜西西安710127）

1 引 言

文物是一個民族的瑰寶，通過文物可以進(jìn)一步解讀歷史、還原真相。隨著計算機技術(shù)、數(shù)字化技術(shù)的快速發(fā)展，借助計算機進(jìn)行文物虛擬修復(fù)已成為研究熱點。借助計算機技術(shù)實現(xiàn)文物碎片拼接，不僅可以提高文物修復(fù)效率，還可以避免對文物造成二次損壞；因此，利用計算機實現(xiàn)文物虛擬修復(fù)具有重要的意義。

文物數(shù)字化虛擬拼接技術(shù)，依據(jù)文物的厚度特征，可將文物分為薄壁類文物和非薄壁類文物［1］。壁畫、陶器和瓷器等沒有厚度特征或厚度特征可以忽略的文物通常歸為薄壁類文物；針對薄壁文物的拼接，目前已有大量的研究［2-6］。石碑、秦俑碎片等具有厚度特征的文物通常歸為非薄壁類文物，本文選用秦俑碎片模型作為研究對象，對非薄壁類文物自動化碎片拼接展開研究。文物碎片拼接是三維點云配準(zhǔn)的一個重要應(yīng)用領(lǐng)域。點云配準(zhǔn)一般分為粗略配準(zhǔn)和精細(xì)配準(zhǔn)兩個階段。粗略配準(zhǔn)是將不同坐標(biāo)系下的點云通過旋轉(zhuǎn)平移操作統(tǒng)一到同一坐標(biāo)系下，但其配準(zhǔn)誤差較大，需在此基礎(chǔ)上利用迭代算法完成精確配準(zhǔn)，使配準(zhǔn)誤差達(dá)到最小。目前，迭代最近點算法（Iterative Closest Point，ICP）［7］及其各種改進(jìn)算法［8-17］已成為精細(xì)匹配的主流算法。樊少榮［18］選取三維碎片的三角網(wǎng)格模型進(jìn)行實驗，分別提取斷裂面的外輪廓線和內(nèi)輪廓線，通過對輪廓線的匹配，從而實現(xiàn)對碎片的拼合。該方法針對完整的斷裂面有較好的拼接效果，但針對斷裂面存在缺失的碎片，其魯棒性較差。李姬俊男［19］使用斷裂面點體積積分不變量和凹凸互補性得到初始匹配約束簇對，最后通過空間幾何一致性約束完成非匹配對之間的消除。李群輝［20］提出了一種基于凹陷凸區(qū)的斷面匹配算法，當(dāng)斷層被劃分為多個凹凸特征區(qū)域時，定義了相似的區(qū)域?qū)Γ鶕?jù)距離原理消除了偽區(qū)域?qū)?。計算相似區(qū)域?qū)Φ淖儞Q參數(shù)以實現(xiàn)碎片的粗匹配。最后采用ICP 算法進(jìn)行精細(xì)匹配。Zhao［21］提出了一種基于輪廓曲線和特征區(qū)域的塊匹配方法，提取斷裂面的輪廓曲線，使用改進(jìn)的迭代最近點算法完成剛性塊與質(zhì)心集的精細(xì)匹配，提高了點云中指定厚度的剛性塊匹配的精度和速度。袁潔［22］通過提取碎片表面紋理和斷裂部位輪廓線特征點，構(gòu)建碎片斷裂面輪廓線并計算雙向距離場，最后通過四元數(shù)算法和最近點迭代算法完成碎片拼合。陸超華［23］使用斷裂面厚度特征信息實現(xiàn)文物碎片匹配。根據(jù)厚度特征將片段分組。將厚度特征使用最長公共子序列算法進(jìn)行匹配，實現(xiàn)了兵馬俑碎片的初步匹配。最后利用ICP 算法進(jìn)行精細(xì)匹配。

綜上所述，現(xiàn)有方法大多要求碎片模型斷裂面信息完整，在斷裂部位受損的破損文物模型上的研究較少，而出土的兵馬俑碎片大多存在不同程度的破損和損壞，傳統(tǒng)的碎片拼接方法拼合誤差大且耗時。針對以上問題，本文選用秦俑碎片作為研究對象，首先，定義鄰域特征參數(shù)，提取特征點集，通過定義曲率特征參數(shù)對特征點集進(jìn)行二次篩選，確定最優(yōu)特征點集；然后，利用特征點間的相對距離和相對位置矢量與特征主方向的夾角作為特征描述符，依據(jù)集合相似理論對特征點進(jìn)行相似性度量，得到匹配點對集合并采用RANSAC 剔除誤匹配點對；最終采用SVD 計算剛體變化矩陣，通過基于K-D 樹改進(jìn)的ICP 算法完成精準(zhǔn)拼合。本文算法流程如圖1 所示。

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

2 碎片特征點提取

針對秦俑碎片斷裂面模型，計算法向量、曲率和法向量內(nèi)積均值，定義鄰域特征參數(shù)，提取特征點集；再定義曲率特征參數(shù)優(yōu)化特征點集，篩選出最優(yōu)特征點集。

2.1 鄰域特征參數(shù)

法向量和曲率是三維空間中特別重要的幾何特征，具有旋轉(zhuǎn)和平移不變性。本文采用主成分分析算法（PCA）進(jìn)行計算，對采樣點pi和其鄰域內(nèi)的所有點進(jìn)行協(xié)方差分析：

將采樣點的法向量ni和其鄰域內(nèi)其他點的法向量nij內(nèi)積后求均值，用于描述鄰域曲面的彎曲程度，內(nèi)積均值δ表示為：

為了描述鄰域內(nèi)各點不規(guī)則分布的程度，計算鄰域內(nèi)各點到鄰域擬合平面L之間的歐式距離方差S：

其中，dj為鄰域內(nèi)各點到鄰域擬合平面L間的歐式距離，如圖2 所示。

圖2 鄰域內(nèi)各點到鄰域擬合平面的距離Fig. 2 Distance from each point in the neighborhood to the neighborhood fitting plane

結(jié)合以上幾何特征，定義鄰域特征參數(shù)φ［24］為：

其中：N為采樣點的個數(shù)，η是為平衡鄰域曲率ε和歐式距離方差S相差過大而引入的的平衡參數(shù)。由公式（7）可知，歐式距離方差S(pi)和鄰域曲率ε(pi)越大，內(nèi)積均值δ(pi)越小，鄰域特征參數(shù)φ(pi)就越大，采樣點pi為特征點的概率越高。設(shè)定鄰域特征參數(shù)閾值σ，若φ(pi)＞σ，即為特征點；否則為非特征點。

2.2 曲率特征參數(shù)

通過2.1 節(jié)，可提取出能夠反映文物碎片斷裂面幾何特征的特征點集，但由于鄰域曲率是受估算值的影響，提取出的特征點集中會存在非特征點，因此定義曲率特征參數(shù)對特征點集進(jìn)行二次特征點提取，得到最優(yōu)特征點集。曲率特征參數(shù)定義為

其中：k1和k2分別表示采樣點pi處的最大主曲率和最小主曲率。根據(jù)公式（9），若采樣點pi的曲率特征參數(shù)ξ(pi)大于其鄰域內(nèi)任意一點的曲率特征參數(shù)，則采樣點pi為鄰域凸點；若采樣點pi的曲率特征參數(shù)ξ(pi)小于其鄰域內(nèi)任意一點的曲率特征參數(shù)，則采樣點pi為鄰域凹點；相比于鄰域曲率，曲率特征參數(shù)可以很好地描述采樣點鄰域表面的凹凸性和彎曲程度。

本文選取兵馬俑一號坑G10-18 號傭部分碎片進(jìn)行實驗，設(shè)定參數(shù)σ=0.4，圖3（a）為碎片實物圖，由圖3（a）可以看出，斷裂面受到一定程度的損壞，但本文算法有較好的魯棒性，可以很好地完成碎片拼接。圖3（b）為碎片特征點提取圖，紅色點為算法自動標(biāo)識的特征點，可以看出，比較平緩的區(qū)域，特征信息較少，特征點分布較為疏散，而凹凸不平的區(qū)域，特征信息較多，特征點分布較為集中。

圖3 G10-18 號傭特征點提取Fig. 3 Feature point extraction of G10-18

3 特征描述及相似性度量

針對上述過程提取出的秦俑碎片模型斷裂面上的特征點，為其定義有效的特征描述，是拼接任務(wù)的關(guān)鍵。本文采用相對距離與方向夾角的集合作為特征描述符，并依據(jù)集合相似理論度量特征點間的相似性，確定匹配特征點對集合。

3.1 特征描述符定義

假設(shè)特征點集合Qkey={q1，q2，…，qm}，m為特征點的個數(shù)，qi表示第i個特征點，本文選用特征點間相對距離d和相對位置矢量與特征主方向的夾角θ作為特征的描述符［25］，特征主方向定義為特征點的散度矩陣的最大特征值對應(yīng)的特征向量。設(shè)特征點qi與其余m-1 個特征點的相對距離為dij，夾角為θij，其中j∈{1，2，…，i-1，i+1，…，m}，假設(shè)共有六個特征點，其中q3的特征描述如圖4 所示。特征點qi的特征描述符Fi可表示為：

圖4 特征描述符的可視化Fig. 4 Visualization of feature descriptors

3.2 特征點匹配

選取兩特征點集合Pkey和Qkey，對于每個特征點的描述，定義相似性度量指標(biāo)為：

Step1 取特征點集Pkey和Qkey中的每一個特征點pi和qj；

Step2 計算pi與Qkey中的每個特征點qj的特征描述集合的相似性，得到相似性最高的特征點qj和相似性次高的特征點qk；

Step3 若滿足式（12），則將匹配的特征點對(pi，qj)加入集合Φ。

3.3 剔除誤匹配點對

通過3.2 小節(jié)得到特征點對集合Φ，但集合Φ中往往會存在誤匹配特征點對，如圖5 所示。剔除誤匹配特征點對可以有效地提高拼接效率和精度。本文采用隨機抽樣一致性算法對誤匹配點對進(jìn)行剔除，算法思路如下：

首先，隨機選取n個特征點對，用于評估樣本集的參數(shù)模型即點云變換矩陣；將特征點對中源點云pi經(jīng)過變換矩陣計算后，得目標(biāo)點云t，計算目標(biāo)點云t與匹配點云qi的偏差D：

若偏差D小于給定閾值d，則認(rèn)為此匹配點對為樣本內(nèi)點并保存；以上步驟迭代m次，取樣本內(nèi)點最多的點對作為最優(yōu)匹配集。算法示意圖如圖6 所示，（a）圖為原始數(shù)據(jù)集，（b）圖中紅點表示樣本內(nèi)點，綠點表示樣本外點（彩圖見期刊電子版）。

圖5 剔除誤匹配點對Fig. 5 Eliminating mismatched point pairs

圖6 隨機抽樣一致性算法原理示意圖Fig. 6 Schematic diagram of random sampling consistency algorithm

4 碎片拼合

4.1 粗略匹配

構(gòu)建方差矩陣E3×3：

其中：Σ為E3×3特征值組成的對角陣，旋轉(zhuǎn)矩陣R和平移矩陣T為：

通過式（14）～式（17）得到旋轉(zhuǎn)矩陣R和平移矩陣T，從而得到點云初始位置關(guān)系，完成碎片初始匹配。

4.2 精確匹配

本文基于K-D 樹對ICP 算法進(jìn)行改進(jìn)，采用改進(jìn)的ICP 算法進(jìn)行精確配準(zhǔn)，利用奇異值分解算法計算迭代過程中的剛體變換矩陣，并通過KD 樹快速搜索最近點對，提高匹配效率。

4.2.1 ICP 算法

ICP 算法是點云精確配準(zhǔn)中常用的方法，其基本原理是尋找最優(yōu)矩陣參數(shù)R和T，使得誤差函數(shù)E最?。?/p>

其中：n為匹配點對個數(shù)，pi和qi分別為兩匹配點，R和T分別為旋轉(zhuǎn)矩陣和平移向量，ICP 算法的具體步驟：

Step1 獲取待配準(zhǔn)點云集合P和Q；

Step2 搜索尋找最近點對集合；

Step3 計算剛體變換矩陣R和T，并計算誤差函數(shù)E；

Step4 將配準(zhǔn)點云Pk變化到新坐標(biāo)系下Pk+1：

Step5 計算點對間的歐式平均距離：

Step6 當(dāng)dˉk-dˉk+1＜Θ 時迭代結(jié)束；否則重新迭代。

4.2.2 ICP 算法改進(jìn)

傳統(tǒng)的ICP 算法通過歐式距離判斷最近點，對于海量點云數(shù)據(jù)耗時且效率低，因此本文對ICP 算法做一個改進(jìn)，利用K-D tree 能快速搜索最近點對，提高點云配準(zhǔn)效率。其搜索步驟如下：

Step1 將待查詢結(jié)點與當(dāng)前樹結(jié)點進(jìn)行比較，若小于等于當(dāng)前樹結(jié)點，則進(jìn)入左子樹分支；否則進(jìn)入右子樹分支，沿搜索路徑找到最近鄰近似點；

Step2 對點進(jìn)行“回溯”操作，若在搜索路徑上有更近的其他子空間點，則跳到子空間結(jié)點去搜索距離最近點；

Step3 重復(fù)Step1 和Step2 直到搜索路徑為空結(jié)束。

5 實驗結(jié)果與分析

本文選取兵馬俑一號坑出土的秦俑碎片模型作為實驗對象來驗證算法的可行性，本實驗在Visual Studio 2015 和Open GL 環(huán)境下進(jìn)行編程，在Intel Core i7-3770CPU/3.4 GHz，內(nèi)存32 GB的PC 機上完成，本文實驗所涉及到的閾值參數(shù)均由經(jīng)驗和數(shù)據(jù)統(tǒng)計分析得出σ=0.4，d=0.1 mm，τ=0.5，τ1=0.5，τ2=0.5。

5.1 本文算法實驗結(jié)果展示

本文選取一號坑兵馬俑G10-12 和G10-41 秦俑碎片模型進(jìn)行拼接效果展示，如圖7～圖8 所示。因受自然因素和人為因素的影響，兵馬俑碎片大多存在缺損。針對輕微損壞的模型，本文算法可以很好地實現(xiàn)拼接，如圖7 所示；針對嚴(yán)重缺損的模型，如圖8 所示，本文算法也有較強的魯棒性和準(zhǔn)確性，拼接結(jié)果基本不受其影響。

圖7 G10-12 號傭部分碎片拼合結(jié)果Fig. 7 Results of fragment assembly of G10-12

圖8 G10-41 號傭部分碎片拼合結(jié)果Fig. 8 Results of fragment assembly of G10-41

為了驗證本文算法的有效性和可應(yīng)用性，選取G10-22 和G10-57 完整秦俑模型進(jìn)行拼接實驗，如圖9～圖10 所示。結(jié)果表明，本文算法可以應(yīng)用于完整傭體的拼接，并有較強的魯棒性，如圖10 所示，通過將破損局部放大可以看出，傭體存在嚴(yán)重破損，但本文算法仍可以很好地完成拼合。

5.2 實驗對比與分析

5.2.1 不同閾值參數(shù)分析比較

圖9 G10-57 號傭正面與背面拼合結(jié)果Fig. 9 Assembling results of front and back of G10-57

圖10 G10-22 號傭正面拼合結(jié)果及缺損部位展示Fig. 10 Display of the results and defects of G10-22

特征點匹配是文物碎片拼接過程中的關(guān)鍵步驟，能否準(zhǔn)確地找出正確匹配點對，直接影響文物拼接結(jié)果，因此τ1和τ2的取值至關(guān)重要。為了更準(zhǔn)確地提取匹配點對，引入閾值τ1和τ2對其進(jìn)行限制，τ1表示與之匹配的特征點中最相似和次相似特征點之間的差距，τ2表示與之匹配的特征點中最相似特征點的相似程度。分別令τ1，τ2=0.4，0.5，0.6，本文算法收斂的評判標(biāo)準(zhǔn)是拼合誤差小于等于1 mm，為驗證算法的穩(wěn)定性，當(dāng)拼合誤差首次小于1 mm 后，再迭代3 次，若每次誤差均小于1 mm，則停止迭代，認(rèn)為算法穩(wěn)定并已達(dá)到收斂。如圖11 所示。

圖11 不同閾值參數(shù)分析比較Fig. 11 Analysis and comparison of different threshold parameters

從圖11 可以看出，當(dāng)τ1，τ2=0.5 時，拼合誤差最小，當(dāng)τ1，τ2=0.4 時，拼合誤差較大，這是由于最相似和次相似點的差距小，相似性高的原因?qū)е碌?；?dāng)τ1，τ2=0.6 時，雖然可以更快擬合，但匹配誤差也偏高，這是由于最相似和次相似點的差距大，匹配對更少的原因?qū)е碌?；綜上，本文選取τ1，τ2=0.5 進(jìn)行實驗。

5.2.2 不同方法分析比較

為了驗證本文算法的優(yōu)越性，與文獻(xiàn)［18］中的算法進(jìn)行實驗結(jié)果對比，本文與文獻(xiàn)［18］對比的原因如下：（1）文獻(xiàn)［18］首次提出了基于斷裂面輪廓線的文物碎片拼接方法，此方法具有一定的先進(jìn)性。（2）文獻(xiàn)［18］與本文均是基于文物碎片斷裂面進(jìn)行碎片拼接。結(jié)果如圖12 所示。

圖12 G10-61 部分碎片拼合結(jié)果Fig. 12 Results of partial fragment assembly of G10-61

從視覺角度分析，由圖12 可以看出，本文算法的拼接效果較好，優(yōu)于文獻(xiàn)［18］方法的結(jié)果，從圖12（b）可以看到，文獻(xiàn)［18］方法拼接過程出現(xiàn)了明顯的錯位現(xiàn)象。這是由于文獻(xiàn)［18］方法對于斷裂部位存在缺失的情況，提取輪廓線出現(xiàn)偏差，從而導(dǎo)致位置匹配不準(zhǔn)確。

為了避免視覺分析的主觀性，本文還從定量角度進(jìn)行分析。由于在計算拼合誤差過程中，點云A 中可能存在多個不同點對應(yīng)點云B 中同一點產(chǎn)生的誤差，因此本文重新定義拼合誤差。假設(shè)集合A 和B 分別為拼合后的兩斷裂面，定義集合A 中某點ai到集合B 的距離為點ai到集合B 中所有點的距離的最小值，計算集合A 中所有點到集合B 的距離的平均值EˉA，再計算集合B 到集合A 的距離的平均值EˉB，定義斷裂面A 和B 的拼合誤差為E，則：

本文的算法運行時間的評價標(biāo)準(zhǔn)如下：通過多次迭代不斷優(yōu)化其拼合誤差，當(dāng)拼合誤差小于1 mm 時，繼續(xù)迭代3 次，若每次誤差都小于1 mm，則停止迭代，計算運行時間，此時間為算法的運行時間，通過對比算法的運行時間，從而可以評判算法的效率。

根據(jù)上述拼合誤差和算法耗時的定義，分別對本文算法與文獻(xiàn)［18］算法計算了拼合誤差和算法耗時，結(jié)果如表1 所示。

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

從表1 中可以看出，本文算法在時間和拼合誤差上都優(yōu)于傳統(tǒng)方法。在時間上，比傳統(tǒng)算法提高了26%～40%；在拼合誤差上，比傳統(tǒng)算法減少了4%～25%。

為了展示本文算法的收斂速度，選取G10-12秦俑碎片為實驗對象，選用文獻(xiàn)［18］方法作為對比方法，對算法運行過程進(jìn)行對比分析，分別計算了不同時間下的拼合誤差，結(jié)果如圖13 所示。本文算法的收斂速度明顯優(yōu)于文獻(xiàn)［18］，在性能上有較大的提升。

在精確匹配階段，為了驗證本文改進(jìn)的ICP算法效率更高，我們與傳統(tǒng)ICP 算法進(jìn)行了實驗對比，選取G10-40，G10-80 和G10-87 號文物碎片模型進(jìn)行實驗，分別計算了兩種方法的拼合誤差和耗時，從表3 中可以看出，與傳統(tǒng)ICP 相比，本文改進(jìn)的ICP 算法在精度和效率上都有明顯的提升，并且在點云數(shù)目越大時顯著性越明顯；因此利用K-D tree 改進(jìn)的ICP 算法可以明顯的提高配準(zhǔn)效率。結(jié)果如表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 結(jié) 論

本文提出一種基于斷裂面鄰域特征的碎片拼接方法，針對受損文物碎片模型，有較強的魯棒性，并且針對稀疏點云特征難以匹配的問題，本文也給出了有效的解決方法。采用鄰域特征參數(shù)和曲率特征參數(shù)兩個約束對特征點二次提取；定義特征點間的相對距離和相對位置矢量與特征主方向的夾角作為特征描述符，通過集合相似理論進(jìn)行度量，得到匹配點對集合，并采用RANSAC 剔除誤匹配點；然后采用SVD 計算得到剛體變化矩陣，最后采用基于K-D 樹改進(jìn)的迭代最近點算法進(jìn)行精確匹配，完成拼合。實驗表明，與傳統(tǒng)算法相比，本文算法在精度和效率上都有明顯提升。

與以往的文物碎片拼接方法相比，本文算法的優(yōu)勢在于：（1）在特征點提取階段，利用最小二乘法原理構(gòu)造曲率特征參數(shù)，對特征點集進(jìn)行優(yōu)化，提高配準(zhǔn)精度。（2）在特征點對匹配階段，引入RANSAC 算法刪除誤匹配點對，篩選出最優(yōu)匹配集，提高配準(zhǔn)精度。（3）在精確匹配階段，引入K-D tree 快速搜索最近點對，提高點云配準(zhǔn)效率。

盡管本文算法可以在破損文物碎片模型上取得良好的效果，但現(xiàn)有文物數(shù)據(jù)庫中，部分文物斷裂面大面積破損；因此，尋找一種更穩(wěn)健的特征描述符進(jìn)行文物碎片拼接是我們下一步的研究方向。