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

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

（西北大學信息科學與技術學院，陜西西安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最?。?/p>

其中：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 快速搜索最近點對，提高點云配準效率。

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