龐戰(zhàn)1,2，滕奇志1，何海波2

　?。?. 四川大學(xué) 電子信息學(xué)院圖像信息研究所，四川 成都 610064；2. 成都西圖科技有限公司，四川 成都 610064）

       摘要：為了解決傳統(tǒng)的基于尺度不變特征變換（SIFT）算法的圖像拼接方法無法使用SIFT特征匹配算法對巖石薄片圖像中的空白區(qū)域圖像進行特征提取，進而影響整個巖石薄片圖像拼接的問題，在尺度不變特征變換匹配算法的基礎(chǔ)上，提出了一種將巖石薄片圖像的行列中心位置作為拼接基準點的圖像拼接方法。該方法改進了特征點的提取方式、變換矩陣的計算順序以及矩陣優(yōu)化的順序，并利用相鄰圖像之間的位置關(guān)系估計出空白區(qū)域圖像的變換矩陣，實現(xiàn)了整個巖石薄片圖像的拼接。實驗結(jié)果表明，該方法可以較好地完成巖石薄片圖像中空白圖像的拼接，能夠較為完整地保留整個巖石薄片的紋理信息，具有一定的實際應(yīng)用價值。

　　關(guān)鍵詞：尺度不變特征變換；空白區(qū)域；巖石薄片；估計；圖像拼接

　　中圖分類號：TP391.4文獻標識碼：ADOI： 10.19358/j.issn.1674-7720.2017.06.015

　　引用格式：龐戰(zhàn)，滕奇志，何海波. 基于SIFT的巖石薄片圖像拼接［J］.微型機與應(yīng)用，2017,36（6）：46-50.

0引言

　　*基金項目：國家自然科學(xué)基金項目（61372174）圖像拼接是指將在同一環(huán)境或者條件下拍攝的一組相互之間具有一定重合量的圖像，通過提取特征點、圖像配準等處理后，得到一幅包含巖石薄片中各個圖像紋理信息的全景圖。目前，圖像拼接技術(shù)在很多方面都有廣泛的應(yīng)用價值［1］。

　　在石油地質(zhì)分析中，由于顯微鏡視域的限制，無法一次性采集得到寬視域巖石薄片的全景圖［2］。為了從整體上分析巖石薄片，只能先對巖石薄片分區(qū)域采集得到整個薄片的圖像，然后將圖像拼接成一幅大視域圖像。目前，圖像拼接主要采用SIFT算法，通過圖像預(yù)處理、圖像配準以及融合完成圖像拼接［35］ 。SIFT算法主要應(yīng)用在相鄰的具有明顯特征點和一定重疊量的圖片拼接中。然而，由于巖石薄片的不規(guī)則性，巖石薄片通常包括空白區(qū)域(如圖1所示)，由于無法進行特征提取，導(dǎo)致整個巖石薄片圖像的拼接失敗。

　　本文在SIFT算法的基礎(chǔ)上，通過改變SIFT特征提取的順序與方式、配準的基準點位置和變換矩陣的計算以及優(yōu)化方式，實現(xiàn)整個巖石薄片圖像的拼接。實驗結(jié)果表明，該方法能夠得到比較完整并且包含了各個圖像信息的巖石薄片全景圖。

1相關(guān)理論

　　1.1SIFT算法的基本原理

　　SIFT算法是由哥倫比亞大學(xué)的LOWE D G教授于1999年提出，并在2004年完善的［6］。該算法的具體流程如下：

　?。?）構(gòu)建尺度空間，檢測極值點

　　尺度空間是模擬圖像數(shù)據(jù)的多尺度特征，檢測圖像中具有尺度不變性的位置。1994年，LINDEBERG T發(fā)現(xiàn)唯一有可能的尺度空間內(nèi)核是高斯函數(shù)［7］。假設(shè)I(x,y)表示一幅圖像，尺度可變的高斯函數(shù)用G(x,y,σ)表示，則其尺度空間的定義如式（1）所示。

　　L(x，y,σ)=G(x,y,σ)*I(x,y)（1）

　　其中，σ代表尺度因子，其大小表示平滑程度，*為卷積運算符，高斯函數(shù)的定義用式（2）表示：

　　G(x,y,σ)=12πσ2exp(－(x2+y2)/2σ2)（2）

　　LOWE D G使用高斯差分算子(DoG)主要有兩個原因：一是能夠更方便地檢測出穩(wěn)定的關(guān)鍵點；二是考慮到計算的復(fù)雜程度。

　　D(x,y,σ)=(G(x,y,kσ)－G(x,y,σ))*I(x,y)

　　=L(x,y,kσ)－L(x,y,σ)（3）

　　其中，k為閾值。LOWE D G在實驗中，通過對原始圖像進行降采樣構(gòu)造若干階尺度空間，如圖2（a）所示，然后在每一階尺度空間內(nèi)利用高斯函數(shù)卷積得到若干層高斯圖像，再差分生成高斯差分尺度空間，如圖2（b）所示。在同一個尺度內(nèi)，將每一個像素點同時與該像素點同尺度的8個相鄰點以及上下兩個相鄰尺度對應(yīng)的18個像素點的值進行比較，只有當(dāng)該像素點的值是最大或者最小值時，才稱該像素點是圖像在該尺度下的一個特征點，也叫關(guān)鍵點。最后，對其余的每一階尺度空間都重復(fù)以上操作，直到檢測完所有尺度空間中的關(guān)鍵點為止。

　?。?）定位關(guān)鍵點，去除不穩(wěn)定點

　　LOWE D G使用擬合三維二次函數(shù)的方法準確定位檢測出來的關(guān)鍵點的位置和尺度。該方法可以有效去除對比度低以及穩(wěn)定性較差的關(guān)鍵點。使用的尺度空間函數(shù)的泰勒公式展開式如式（4）所示：

　　其中，X=(x,y,σ)T，在采樣點處，計算D和D的導(dǎo)數(shù)，當(dāng)其導(dǎo)數(shù)為零時，得到極值點的X1，如式（5）所示：

　　利用尺度空間函數(shù)在極值點處的函數(shù)值去除對比度較低的響應(yīng)點。將式（5）代入式（4），計算得到極值點處函數(shù)值，如式（6）。

　　通常情況下，若|D(X1)|≥0.03，該特征點就保留，否則就去除。

　　因為高斯差分算子在邊緣處會產(chǎn)生比較強烈的邊緣響應(yīng)，在該算子中，差的極值點在邊緣處具有較大的曲率，在垂直方向上具有較小的曲率。因此，可以利用該性質(zhì)去除邊緣響應(yīng)。主曲率可以通過一個2×2的Hessian矩陣H求出，如式（7）所示。

　　D的主曲率和H的特征值成正比，若假設(shè)α為最大特征值，β為最小特征值，則矩陣H主對角線的值和行列式的值用式（8）和式（9）表示。

　　如果曲率大于式（10）結(jié)果的比值，則舍棄。LOWE D G在實驗中通過對比發(fā)現(xiàn)，當(dāng)r=10時，能夠最大程度去除不穩(wěn)定的邊緣響應(yīng)點［6］。

　　（3）確定關(guān)鍵點的大小和方向

　　在初步確定了關(guān)鍵點的位置之后，可以利用圖像的局部特性為每一個關(guān)鍵點分配一個基準點，使該關(guān)鍵點的特征描述符具有旋轉(zhuǎn)不變性。對于高斯金字塔檢測出來的關(guān)鍵點，使用其高斯金字塔鄰域窗口內(nèi)像素的梯度特征和方向特征對描述符進行描述。(x,y)處的梯度值m(x,y)和方向θ(x,y)的表達式如式（11）和式（12）所示。

　　其中，L所用的尺度為每個關(guān)鍵點各自所在的尺度。為了進一步準確給出關(guān)鍵點的方向，一般采用梯度直方圖統(tǒng)計的方法，以關(guān)鍵點為中心的鄰域窗口內(nèi)的像素的梯度和方向。直方圖峰值所代表的方向即為關(guān)鍵點的方向。通過上述各個步驟之后，檢測出的包含位置、所在尺度以及方向的關(guān)鍵點稱為圖像的特征點。

　?。?）生成特征描述子

　　為了確保關(guān)鍵點不受光線等因素的影響，可以為每一個關(guān)鍵點建立一個描述子，使其具有獨特性，進而可以提高特征點的匹配率。為了保證描述子具有旋轉(zhuǎn)不變性，首先將方向旋轉(zhuǎn)到關(guān)鍵點的方向。LOWE D G的實驗中，在關(guān)鍵點所在的尺度空間中，以關(guān)鍵點為中心，選取一個8×8的窗口，然后在每個4×4的窗口內(nèi)計算8個方向的梯度方向直方圖，繪制每個梯度方向的累加可形成一個種子點，最終獲得一個4×4×8共128維的向量特征。

　　（5）特征點的匹配

　　SIFT特征向量生成后，一般使用最近鄰域算法對特征點進行匹配。利用特征點的最短歐氏距離和次最短歐氏距離的比值與所設(shè)定的閾值的大小關(guān)系作為判斷兩個特征點是否為一匹配對的依據(jù)。如果結(jié)果小于閾值，則認為這一對特征點是匹配的。LOWE D G在實驗中發(fā)現(xiàn)，將閾值設(shè)定為0.8比較合適，在計算特征點的歐氏距離之后，可以先利用BBF（BestBinFirst）算法對特征向量進行快速處理，再使用RANSAC算法去除誤匹配的特征點［7］。

　　1.2變換矩陣計算的基本原理

　　描述圖像之間位置關(guān)系的矩陣稱為變換矩陣。圖像的變換矩陣可以通過特征匹配對計算得到［89］。在同一個空間下，兩幅圖之間的變換關(guān)系用式（13）表示，其中，(x1,y1)為參考圖像上的點，(x2,y2)為目標圖像上與(x1,y1)對應(yīng)的點，M稱為變換矩陣:

　　在式（13）中，線性方程組中有8個未知變量，則至少需要4對特征點匹配對才能求解方程組的解。整理并將結(jié)果進行向量化用式（14）表示:

　　A1=MA2（14）

　　若矩陣A2AT2可逆，則可以求解得到變換矩陣M，如式（15）所示：

　　M=(A1AT2)(A2AT2)-1（15）

　　只要有足夠多可用的特征點匹配對就可以通過式（15）計算得到圖像之間的變換矩陣，從而確定圖像之間的相對位置關(guān)系。

2空白區(qū)域圖像拼接方法

　　巖石薄片中空白區(qū)域的圖像（見圖1）雖然可以使用SIFT算法對相鄰圖片進行特征匹配，但由于空白區(qū)域圖像的表面特征比較單一，紋理信息不豐富，檢測到的特征點匹配對大多都是錯誤匹配，在使用RANSAC算法［10］去除錯誤匹配之后，可用的特征匹配數(shù)目很少。如圖3所示，從實驗對比結(jié)果可以看出，巖石薄片中空白區(qū)域的相鄰圖像的特征匹配數(shù)目在使用RANSAC算法去除錯誤匹配之后，沒有可用的特征點匹配對，無法利用式（15）求出圖像的變換矩陣，從而無法確定兩張圖片的相對位置，影響整個巖石薄片圖像的拼接。

　　圖4拼接流程圖在實際應(yīng)用中，為了解決該問題，在SIFT算法的基礎(chǔ)上，提出了一種選取巖石薄片圖像的中心位置為基準點的拼接方法。并利用空白區(qū)域圖像與其相鄰圖像之間的位置關(guān)系估計出空白區(qū)域圖像的變換矩陣，實現(xiàn)了整個巖石薄片的拼接。流程圖如圖4所示。

　　2.1特征點提取

　　為了解決SIFT算法無法提取空白區(qū)域圖像的特征點的問題，本文在SIFT算法的基礎(chǔ)上，改變了特征點提取的方式，并在RANSAC算法篩選時設(shè)定閾值，如篩選結(jié)果小于閾值則視為不能正常匹配，并標記。本實驗設(shè)定的閾值為30，該閾值也可以根據(jù)圖片特征點的具體情況進行設(shè)定。如圖5所示，該示意圖表示3行7列圖片。其中，橫向箭頭表示每一行圖片的特征點提取順序，豎向箭頭表示中間一列圖片的特征點提取順序。具體步驟如下：

　?。?）獲取待拼接巖石薄片圖像的行和列的數(shù)目。

　?。?）分別對巖石薄片圖像中的每一行和中間列進行特征提取，并對不能正常匹配的圖片進行標記。

　　2.2計算變換矩陣

　　由于無法得到空白區(qū)域圖片的特征點匹配對，從而無法確定空白區(qū)域圖片的相對位置。在前文特征提取順序的基礎(chǔ)上，對變換矩陣的計算順序作了相應(yīng)的調(diào)整。如圖6所示，其中每一格代表一張圖片，圓形代表本文選取的中心位置，五角星代表空白圖片，箭頭表示變換矩陣的計算方向，每一行圖片變換矩陣的計算方向都與示意圖中箭頭表示的方向一致。文中，“左（右）側(cè)”是表示分別對左側(cè)和右側(cè)都進行相同的操作。具體步驟如下：

　　（1）選取巖石薄片圖像的中心位置（圖6中圓形所在位置），并初始化其變換矩陣為單位矩陣，作為變換矩陣計算的基準點。

　?。?）從中心位置所在行開始，到圖像的起始行，分別從每一行的中間位置開始，對左（右）側(cè)的圖片計算變換矩陣，如果遇到標記圖像（圖6中圓形左（右）側(cè)五角星代表的圖片）則停止計算當(dāng)前圖片的變換矩陣，并對標記圖片左（右）側(cè)的所有圖片的變換矩陣賦值為空。

　?。?）從中心位置所在行開始到圖像的最后一行，分別從每一行的中間位置開始重復(fù)步驟（2）。

　?。?）從圖像中間列的中間位置開始，對中間列圖像的上下兩部分的圖片重復(fù)步驟（1）和（2）。

　　2.3優(yōu)化變換矩陣

　　計算完所有圖片的變換矩陣之后，本文在上一節(jié)提出的變換矩陣的計算順序的基礎(chǔ)上，對變換矩陣的優(yōu)化方式作了相應(yīng)的調(diào)整。具體優(yōu)化方式如下：

　　（1）從中心位置所在行開始到圖像的起始行，分別從每一行的中間位置左（右）側(cè)的圖片進行變換矩陣優(yōu)化，如果遇到標記圖像則停止優(yōu)化。直至每一行圖像都重復(fù)該操作為止。

　?。?）從中心位置所在行開始到圖像的最后一行，分別從每一行的中間位置開始重復(fù)步驟（1）。

　?。?）從薄片圖像中間列的中間位置開始，對中間列圖像的上下兩部分的圖片重復(fù)步驟（1）和（2）。

　　2.4估計空白區(qū)域圖片的變換矩陣

　　變換矩陣優(yōu)化完成之后，為了獲得完整的拼接結(jié)果，需要估計空白區(qū)域圖片的變換矩陣。具體方法如下。

　?。?）從中心位置所在行開始到圖像的起始行，分別從每一行的中間位置開始，對其左（右）側(cè)的圖片進行標志位掃描，遇到標記圖片，則利用該標記圖片右（左）側(cè)相鄰兩張圖片的變換矩陣與標記圖像之間的相對位置關(guān)系估計出標記圖片的變換矩陣，同時位于標記圖片左（右）側(cè)的所有圖片的變換矩陣均通過該方法依次得到。直至每一行圖像都重復(fù)該操作為止。

　?。?）從中心位置所在行開始到圖像的最后一行，分別從每一行的中間位置開始重復(fù)步驟（1）。

　?。?）從圖像中間列的中間位置開始，對中間列圖像的上下兩部分的圖片則重復(fù)步驟（1）和（2）。

　　完成上述處理之后，通過相應(yīng)的融合算法對結(jié)果圖進行融合，消除拼接縫［11］，完成整個巖石薄片圖像的拼接。

3實驗對比結(jié)果及分析

　　本文使用6行6列共36張圖片進行實驗，圖片尺寸都是5 184×3 456，實驗所用的圖像包含了實驗所需要的薄片信息。實驗結(jié)果如圖7~9所示，其中，圖7是實驗所用的圖像，圖8是基于SIFT算法的拼接結(jié)果圖，圖9是本文提出的拼接方法的結(jié)果圖。

　　從實驗的對比結(jié)果中可以明顯看出，當(dāng)采集得到的巖石薄片圖像中包含空白區(qū)域圖片時，如圖7實驗使用的圖像中“C0001.jpg”和“C0033.jpg”等圖片，因為無法使用SIFT算法對巖石薄片中空白區(qū)域圖片完成配準，因此，只能在拼接前剔除這些圖片所在的行或者列的所有圖片，對比可以看到基于SIFT算法的拼接結(jié)果圖中，第一張圖片所在的列和下方最后一行都存在嚴重缺失（結(jié)果圖中的黑色缺失部分）。而本文提出的圖像拼接方法則較為完整地完成了整個巖石薄片圖像的拼接，較好地反映出整個巖石薄片中顆粒的分布情況，保證了巖石薄片的完整性。因此，本文提出的圖像拼接方法能較好地完成巖石薄片空白區(qū)域圖像的拼接，能較為完整地展示巖石薄片的全貌，具有一定的實際應(yīng)用價值。

4結(jié)論

　　本文首先介紹了SIFT算法和變換矩陣計算的基本原理，同時對巖石薄片中空白區(qū)域圖像不能使用SIFT算法進行特征點提取的原因作了簡要分析，并在此基礎(chǔ)上分別對本文提出的圖像拼接方法中的特征點的提取方式、變換矩陣的計算方式、變換矩陣的優(yōu)化方式以及估計空白區(qū)域圖片的變換矩陣的具體步驟進行了描述，在不改變巖石薄片整體性基礎(chǔ)上，完成了包含空白區(qū)域的巖石薄片圖像的拼接，保證了巖石薄片紋理信息的完整性。本文的不足之處在于，若在高倍顯微鏡下巖石薄片的中心區(qū)域存在大量空白區(qū)域，則無法具體確定拼接的基準位置，對該方法的使用具有一定的限制。

　　參考文獻

　?。?］ 徐鑫, 孫韶媛, 沙鈺杰，等. 一種基于RANSAC的紅外圖像拼接方法［J］. 激光與光電子學(xué)進展，2014，51(11):129-134.

　?。?］ 岳永娟, 苗立剛, 彭思龍. 大規(guī)模顯微圖像拼接算法［J］. 計算機應(yīng)用，2006, 5(26):1012-1014.

　?。?］ 鄒承明, 侯小碧, 馬靜. 基于幾何學(xué)圖像配準的SIFT圖像拼接算法［J］. 華中科技大學(xué)學(xué)報（自然科學(xué)版），2016，44(4)：32-36.

　?。?］ 汪松, 王俊平, 萬國挺，等. 基于SIFT算法的圖像匹配方法［J］. 吉林大學(xué)學(xué)報（工學(xué)版），2013，43(S1):279-282.

　?。?］ 白廷柱, 侯喜報. 基于SIFT算子的圖像匹配算法研究［J］. 北京理工大學(xué)學(xué)報，2013，33(6):622-627.

　?。?］ LOWE D G. Distinctive image features from scaleinvariant keypoints［J］. International Journal of Computer Vision，2004，60(2): 91110.

　?。?］ LINDEBERG T. Scalespace theory: a basic tool for analyzing structures at different scales ［J］. Journal of Applied Statistics, 1994, 21(2):225-270.

　?。?］ 郭紅玉, 王鑒. 一種基于RANSAC基本矩陣估計的圖像匹配方法［J］. 紅外，2008，29（2）: 5-8.

　?。?］ 曲天偉, 安波, 陳桂蘭. 改進的RANSAC算法在圖像配準中的應(yīng)用［J］. 計算機應(yīng)用，2010，30（7）：1849-1851.

　?。?0］ FISCHER M A, BOLLES R C. Random sample consensus: a paradigm for model fitting with applications to analysis and automated cartography ［J］. Communication of the ACM,1981, 24 (6):381-395．

　?。?1］ 苗啟廣, 王寶樹. 基于改進的拉普拉斯金字塔變換的圖像融合方法［J］.光學(xué)學(xué)報，2007，9(27): 1605-1610.