基于SURF特征描述符和杰卡德距離的文物碎片拼接

劉曉寧，狄宏璋，楊 穩，林芃樾，王世雄

(西北大學 信息科學與技術學院，陜西 西安 710127)

1 引 言

破碎文物的拼接在文物虛擬復原工作中的意義不言而喻,它是虛擬復原中至為關鍵的一步。由于信息技術的不斷更新，文物的數字化虛擬拼接也不斷有新的突破。目前文物的虛擬拼接技術大致可以分為斷裂面受損和斷裂面完整兩種。

對于斷裂面完整類文物，一般是采用基于斷裂面輪廓線以及模型表面幾何紋飾特征來進行匹配。比如劉軍等[1]提出一種將斷裂面幾何特征與輪廓線特征相融合的拼接方法。袁潔[2]通過定義表面約束點來提取斷裂面上的特征點,實現交互式文物拼接。趙夫群等[3]通過計算最長公共子序列實現斷裂面的細匹配，從而確定其鄰接關系。Schmid C等[4]提出一種高效的三維輪廓曲線來實現文物的拼接。

對于斷裂面受損類文物，一般是從斷裂面厚度出發，提取出魯棒性能較好的斷裂面顯著特征來實現拼接。比如李珊珊[5]提出基于表面鄰接約束并借助于領域專家先驗知識的交互式文物拼接。袁潔等[6]提出基于Morse-smale拓撲特征的文物拼接算法，通過Morse-smale復形理論來構建斷裂面的幾何四邊形拓撲圖，計算4個頂點到基準面的高度差來構建拓撲圖的特征描述符，來實現相似性檢索。

目前現有的文物拼接算法，對于局部碎片缺失、數據量大的碎片模型，拼合效果不是特別的準確。因此，本文基于SURF(Speeded Up Robust Features)特征描述符和Jaccard距離來實現文物的精準拼接。首先采用Canny邊緣檢測算子[7-10]提取出文物模型周圍的邊沿線和表面幾何紋飾線條，作為下一步特征提取的基礎。其次,通過構造多尺度空間，在模型邊緣線表面上提取出特征點，對其構造一系列具有強魯棒性更顯著的SURF特征點描述符[11-18]。然后,用Jaccard距離[19-22]計算特征點集相似度并篩選最優匹配點集，確定其鄰接關系。最后，采用ICP(Iterative Closest Point)算法[23]多次迭代計算出剛體變化的參數,從而實現碎片拼接。實驗表明在保證較精確拼接的同時，本文方法具有高效率、低冗余、強魯棒等特點。

2 Canny邊緣檢測

2.1 邊緣檢測

對于三維模型來說，邊緣是最直觀、最能體現其幾何紋飾信息的一種表示。Canny邊緣檢測滿足以下要求:(1)檢測出來的邊緣線是基于人的視覺感受；(2)提取出來的邊緣線不會影響模型表面的紋飾特征；(3)邊緣線較細，并且無縫連接，不會出現局部斷開。所以本文用Canny算子提取碎片模型邊緣線。

2.2 高斯濾波平滑圖像

首先，進行模型預處理操作,消除噪聲等因素對文物邊緣線檢測產生的影響。高斯濾波器計算如公式(1)：

i≤2k+1,j≤2k+1，

(1)

然后要與8鄰域像素矩陣做卷積運算，為了計算的高效性與實效性，故令k=1,ε=1.2,經過數據歸一化后的高斯濾波器H如公式(2)所示:

(2)

設模型中某一點像素值為e，其8鄰域像素值矩陣為B。將其與3×3的高斯濾波器模板的中心位置對齊從而進行卷積操作，即將該像素點8鄰域乘積權值疊加，作為該點高斯模糊后的亮度值w。則有:

(3)

2.3 邊緣梯度檢測

由于碎片原始模型大多紋理豐富、特征復雜，若在其上直接提取特征點會導致特征點集不夠泛化，無法更進一步匹配。因此，用Canny算子檢測出其水平、垂直以及斜對角線方向的邊緣線。Canny算子x方向和y方向的3×3模板分別為:

(4)

梯度計算公式如式(5)所示，H代表計算得到的梯度值，θ代表其梯度方向。

(5)

2.4 邊緣線生成

生成的邊緣線要最大程度反應碎片的邊緣信息和表面紋飾特征，即先進一步縮小范圍確定當前局部區域內最優點，對任一像素點P,若其梯度值大于我們設定的高閾值，則認為P為強邊緣像素，將其保留并納入邊緣像素點集合M={mi|i=1,2,3,...,n}。若其梯度值大于低閾值并且小于高閾值，則認為其為弱邊緣像素，再去觀察其8鄰接像素點梯度值，若存在梯度值大于高閾值的像素，則將其也納入到邊緣像素點集合M里，否則將其刪除。若任一像素點的梯度值小于低閾值，則直接將其刪除。其偽代碼描述如下:

Algorithm: Progress of picking up four-element pixel1 Input: P: One of a pixel;2Hthreshold: high threshold;3Lthreshold: low threshold;4 Output: P is edge pixel or not;5 if P>hthreshold then6 P must be on edge7 end if8 else if P>Lthreshold and P>Hthreshold then9 for each ∑p' do10 if p' Hthreshold then11 P must be an edge

12 break13 end if14 end for15 end if16 else17 P must be renamed

其中，P代表當前像素點，Hthreshold=0.8,Lthreshold=0.4分別為設定的高低閾值，它是通過多次實驗計算得出，當高閾值與低閾值的比值為2∶1時，提取出的邊緣線效果更優。

該算法提取的模型邊緣線更加符合人類視覺觀感，可以清晰的表示模型表面幾何紋飾特征以及邊緣輪廓信息，為下一步碎片特征點提取打下了良好的基礎，如圖1所示分別C2-4-065，C2-4-064，C2-4-033號碎片模型及其對應的邊緣線圖。

圖1 模型邊緣線提取Fig.1 Model edge line extraction

3 特征提取

提取碎片顯著特征點大都是通過曲率、極大值及極小值的計算來確定。本文結合以往先驗知識，通過構造多尺度空間使提取出來的特征點更加泛化，具有一般性。其主要步驟如下:

Step 1.構造多尺度空間。為了更進一步縮小特征點的范圍，使ε分別取值為1.1、1.3、1.4，得到的高斯濾波器分別為h1,h2和h3，如公式(6)所示。將高斯濾波器與碎片模型卷積操作會得到不同尺度空間下的模型。

(6)

Step 2.空間特征點檢測。為了保證關鍵點的尺度不變性以及高魯棒性,需要對模型鄰域和尺度空間鄰域點進行大小比較。如下圖所示,取中間的檢測點p和其鄰域8個點進行大小比較。在其尺度空間內,和其上下鄰域共18個檢測點進行大小比較。若檢測點p在其鄰域和空間鄰域內是最大值或最小值,就認為點p為候選特征點。

圖2 不同尺度空間的圖像Fig.2 Images in different scale spaces

Step 3.特征點精確定位。Step2中得到的特征點其實是屬于多尺度空間中候選特征點，其具有不穩定性與發散性，并不一定是真實的特征點。還需要對其進行曲線擬合來使其線性化,精確的確定特征點的位置，剔除掉低分辨率的特征點,增強匹配的穩定性與完整性。

圖3 特征點定位Fig.3 Feature point positioning

這里曲線擬合函數的展開式為:

(7)

對其求導取0，得到:

(8)

將該值代入展開式里得到:

(9)

4 構建SURF特征描述符

構建特征描述符大都以斷面特征點為基礎，這樣如果提取的特征點存在誤差，也會導致特征描述符失真。本文基于SURF特征描述符的魯棒性與一般性，構建出不僅包括該特征點的相關信息,還會包含其鄰域范圍內有較強光譜映射的特征點信息的SURF特征描述符。

圖4 構建特征描述符Fig.4 Construct feature description factor

如圖4中心紅點Q為所選擇的特征點，以Q為中心點結合直方圖映射構造一個4×4模板，在這個模板內每一小塊都認為是點Q的一個鄰域生長點(彩圖見期刊電子版)。箭頭方向表示每個鄰域生長點的8個梯度方向。再用卷積核與每一塊進行卷積加權運算,可得到每個鄰域生長點的8個梯度描述符。故對于任一特征點Q來說,其特征描述符可以用4×4×8維的梯度描述向量W來表示,W=(w1,w2,w3,....,w128)。在對其進行歸一化得到L=(l1,l2,...,l128)。歸一化公式如式(10):

(10)

5 相似性匹配

由于碎片模型拓撲特征的復雜性以及特征點的多維度特性,若用歐氏距離計算特征點集的匹配關系,其運算結果時間復雜度會特別的高并且會存在一定的誤差。基于此,為了對待拼接文物碎片斷裂面及模型表面紋飾特征點的相似型進行精確檢測,本文采用Jaccard距離來構造最優特征點匹配集合,從而確定其鄰接關系,大大提高了匹配檢測的效率。其步驟如下:

Step 1.構建Jaccard系數?。設碎片SA和碎片SB的特征點集合分別為A和B：

A={Ai|i=1,2,...,n}

B={Bj|j=1,2,...,k}

(11)

其中：∩表示特征點集合A與B的交集,∪表示特征點集合A與B的并集。

Step 2.構建Jaccard距離d。其公式為:

(12)

對碎片SA中任一特征點Ai,遍歷碎片SB的特征點集合B,若存在特征點Bj滿足Jaccard距離d(Ai,Bj)(j=1,2,...,n)納入匹配集合Q內。此時還不是最優,因為在B中會存在多個特征點。所以需要遍歷匹配集合Q，找出滿足Jaccard距離最小的特征點Bj,從而構成最佳匹配組。

6 剛體拼合

經過上述操作,已經得到了待拼接碎片間最優匹配特征點集合H={〈SAi,SBi〉|i=1,2,...,n}。這一步主要是將位于不同坐標系下的特征點集,轉化為同一坐標系下,來完成文物碎片的準確拼接。本文采用ICP算法來求解旋轉矩陣R和平移矩陣T，從而實現碎片模型的剛體變化。主要步驟如下:

Step 1:從最優匹配特征點集合中任意選取兩組不共線的特征點對〈SAj,SBj〉〈SAk,SBk〉(j,k∈[1,n],j≠k)。利用ICP算法計算出旋轉矩陣R和平移矩陣T,將其保存在剛體變化數組V[m]中。

Step 2:經過上述操作得到m個剛體變化參數,則剛體變化數組記為V[m]={v0,v1,v2，...,vm},其中v[i]={(Ri,Ti)|i∈[0,m]}。

Step 3: 遍歷剛體變化數組V[m],將其代入目標函數求解。則求解最優變化參數R和T就可以轉化為求滿足目標函數W最小值的最優解問題。

(13)

Step 4:求解得到最優參數R和T,完成碎片的精確拼接。

7 實驗結果及分析

本實驗采用VicualC++和Open CV編程,在Intel i7-6700k CPU/4.0 GHz,32 G內存的PC機上實現。以1號坑的C2-4號121塊右驂馬碎片模型作為實驗數據。

圖5 C2-4號右驂馬發掘現場圖Fig.5 Excavation site of No. C2-4

7.1 重組實驗結果

圖5所示為采用本文算法對c2-4號右驂馬碎片重組后的效果圖，分別為右視圖、左視圖、正視圖及其對應的邊緣線圖與提取出的三維特征點云圖。拼合時間如表1所示,其中n表示碎片塊數,t1表示提取模型邊緣線所需的時間,t2表示確定碎片及周邊碎片的幾何鄰接關系所需的時間,t3表示碎片重組拼接需要的時間。

圖6 右驂馬三視圖Fig.6 Three views of horse

表1 右驂馬拼接所需時間

Tab.1 Splicing time required for horse

n t1/s t2/st3/s12170.34657.44223.212

7.2 局部拼接實驗結果

圖6所示分別為C2-4-064,C2-4-063和C2-4-033號碎片及其拼接關系圖，其位置如圖藍色標注處(彩圖見期刊電子版)，直線標明的是其周邊鄰接碎片的位置。可以看到即使存在部分斷裂面受損的碎片，采用本文方法也能較準確地定位出其鄰接位置，取得較好的拼合效果。

圖7 拼接關系圖Fig.7 Splicing diagram

7.3 對比實驗結果

圖7(a)為采用本文方法對C2-4-18(藍色標注處)號碎片進行拼接的效果。圖7(b)為用文獻[22]方法的拼接效果圖(紅色標注處，彩圖見期刊電子版)，可以看到其幾何特征無法正確的反映出其鄰接關系，呈現出浸透現象，如黑色方框所標注。

圖8 C2-4-18號碎片拼接效果圖Fig.8 No.C2-4-18 effect drawing of fragment splicing

圖9 C2-4-13號碎片拼接效果圖Fig.9 No.C2-4-18 effect drawing of fragment splicing

圖8(a)為采用本文方法對C2-4-13號碎片(藍色標注處)進行拼接的效果。圖8(b)為用文獻[18]方法的拼接效果圖(紅色標注處，彩圖見期刊電子版)，可以看到存在明顯的縫隙過大現象，如黑色圓圈標注處所示。

表2所示為本文算法與文獻[18]、文獻[22]算法的運行時間和拼合誤差對比，拼合誤差為鄰接碎片斷裂面輪廓線上特征點集間的歐式距離之和的平均值。本文算法的運行時間與文獻[18]相比提高了16%，與文獻[22]相比提高了12%。拼合誤差相對于文獻[18]和文獻[22]來說更小，不超過0.750 mm，拼合效果更準確，適用于斷裂面受損的碎片拼接。

表2 運行時間與拼合誤差對比

8 結 論

本文主要針對文物多碎片拼接過程中出現的局部碎片缺失、紋飾幾何特征點受損而不能準確提取斷裂面紋飾特征等嚴重問題，提出基于SURF特征描述符和Jaccard距離的碎片拼接方法。先用Canny算子提取出碎片模型中能夠表示其幾何紋飾特征的邊緣線，從而進一步簡化其幾何模型。再構造多Scale空間來提取碎片模型斷裂面的特征點，構建SURF特征描述子。利用特征點集間的鄰接幾何特性，不同于歐式距離計算的高復雜度與低時性，本文用Jaccard距離計算特征點集相似度。實驗結果表明，本文方法使得拼合誤差減少到0.750 mm內，算法運行時間與文獻[18]和文獻[22]相比分別提高了12%和16%，可以有效地減少因斷裂面破損而造成的拼接縫隙過大、出現滲透等現象。

由于構建特征描述符的維數比較大，因此尋找更加魯棒、高效的特征描述符將是下一步著重的研究方向。