基於SIFT的岩石薄片圖像拼接

龐戰1,2， 滕奇志1， 何海波2

（1. 四川大學 電子資訊學院圖像資訊研究所， 四川 成都 610064；2. 成都西圖科技有限公司， 四川 成都 610064）

：為了解決傳統的基於尺度不變特徵變換（SIFT）演算法的圖像拼接方法無法使用SIFT特徵匹配演算法對岩石薄片圖像中的空白區域圖像進行特徵提取， 進而影響整個岩石薄片圖像拼接的問題， 在尺度不變特徵變換匹配演算法的基礎上， 提出了一種將岩石薄片圖像的行列中心位置作為拼接基準點的圖像拼接方法。 該方法改進了特徵點的提取方式、變換矩陣的計算順序以及矩陣優化的順序，

並利用相鄰圖像之間的位置關係估計出空白區域圖像的變換矩陣， 實現了整個岩石薄片圖像的拼接。 實驗結果表明， 該方法可以較好地完成岩石薄片圖像中空白圖像的拼接， 能夠較為完整地保留整個岩石薄片的紋理資訊， 具有一定的實際應用價值。

：尺度不變特徵變換；空白區域；岩石薄片；估計；圖像拼接

：TP391.4文獻標識碼：ADOI： 10.19358/j.issn.1674-7720.2017.06.015

引用格式：龐戰， 滕奇志， 何海波. 基於SIFT的岩石薄片圖像拼接［J］.微型機與應用， 2017,36（6）：46-50.

0引言

*基金專案：國家自然科學基金專案（61372174）圖像拼接是指將在同一環境或者條件下拍攝的一組相互之間具有一定重合量的圖像， 通過提取特徵點、圖像配准等處理後， 得到一幅包含岩石薄片中各個圖像紋理資訊的全景圖。

目前， 圖像拼接技術在很多方面都有廣泛的應用價值［1］。

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

本文在SIFT演算法的基礎上， 通過改變SIFT特徵提取的順序與方式、配准的基準點位置和變換矩陣的計算以及優化方式，

實現整個岩石薄片圖像的拼接。 實驗結果表明， 該方法能夠得到比較完整並且包含了各個圖像資訊的岩石薄片全景圖。

1相關理論

1.1SIFT演算法的基本原理

SIFT演算法是由哥倫比亞大學的LOWE D G教授于1999年提出，

並在2004年完善的［6］。 該演算法的具體流程如下：

（1）構建尺度空間， 檢測極值點

尺度空間是類比圖像資料的多尺度特徵， 檢測圖像中具有尺度不變性的位置。 1994年， LINDEBERG T發現唯一有可能的尺度空間內核是高斯函數［7］。 假設I(x,y)表示一幅圖像， 尺度可變的高斯函數用G(x,y,σ)表示， 則其尺度空間的定義如式（1）所示。

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

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

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

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

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在實驗中， 通過對原始圖像進行降採樣構造若干階尺度空間， 如圖2（a）所示， 然後在每一階尺度空間內利用高斯函數卷積得到若干層高斯圖像，

再差分生成高斯差分尺度空間， 如圖2（b）所示。 在同一個尺度內， 將每一個圖元點同時與該圖元點同尺度的8個相鄰點以及上下兩個相鄰尺度對應的18個圖元點的值進行比較， 只有當該圖元點的值是最大或者最小值時， 才稱該圖元點是圖像在該尺度下的一個特徵點， 也叫關鍵點。 最後， 對其餘的每一階尺度空間都重複以上操作， 直到檢測完所有尺度空間中的關鍵點為止。

（2）定位關鍵點，去除不穩定點

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

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

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

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

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

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

如果曲率大於式（10）結果的比值，則捨棄。LOWE D G在實驗中通過對比發現，當r=10時，能夠最大程度去除不穩定的邊緣響應點［6］。

（3）確定關鍵點的大小和方向

在初步確定了關鍵點的位置之後，可以利用圖像的局部特性為每一個關鍵點分配一個基準點，使該關鍵點的特徵描述符具有旋轉不變性。對於高斯金字塔檢測出來的關鍵點，使用其高斯金字塔鄰域視窗內圖元的梯度特徵和方向特徵對描述符進行描述。(x,y)處的梯度值m(x,y)和方向θ(x,y)的運算式如式（11）和式（12）所示。

其中，L所用的尺度為每個關鍵點各自所在的尺度。為了進一步準確給出關鍵點的方向，一般採用梯度長條圖統計的方法，以關鍵點為中心的鄰域視窗內的圖元的梯度和方向。長條圖峰值所代表的方向即為關鍵點的方向。通過上述各個步驟之後，檢測出的包含位置、所在尺度以及方向的關鍵點稱為圖像的特徵點。

（4）生成特徵描述子

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

（5）特徵點的匹配

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

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

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

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

A1=MA2（14）

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

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

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

2空白區域圖像拼接方法

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

圖4拼接流程圖在實際應用中，為了解決該問題，在SIFT演算法的基礎上，提出了一種選取岩石薄片圖像的中心位置為基準點的拼接方法。並利用空白區域圖像與其相鄰圖像之間的位置關係估計出空白區域圖像的變換矩陣，實現了整個岩石薄片的拼接。流程圖如圖4所示。

2.1特徵點提取

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

（1）獲取待拼接岩石薄片圖像的行和列的數目。

（2）分別對岩石薄片圖像中的每一行和中間列進行特徵提取，並對不能正常匹配的圖片進行標記。

2.2計算變換矩陣

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

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

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

（3）從中心位置所在行開始到圖像的最後一行，分別從每一行的中間位置開始重複步驟（2）。

（4）從圖像中間列的中間位置開始，對中間列圖像的上下兩部分的圖片重複步驟（1）和（2）。

2.3優化變換矩陣

計算完所有圖片的變換矩陣之後，本文在上一節提出的變換矩陣的計算順序的基礎上，對變換矩陣的優化方式作了相應的調整。具體優化方式如下：

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

（2）從中心位置所在行開始到圖像的最後一行，分別從每一行的中間位置開始重複步驟（1）。

（3）從薄片圖像中間列的中間位置開始，對中間列圖像的上下兩部分的圖片重複步驟（1）和（2）。

2.4估計空白區域圖片的變換矩陣

變換矩陣優化完成之後，為了獲得完整的拼接結果，需要估計空白區域圖片的變換矩陣。具體方法如下。

（1）從中心位置所在行開始到圖像的起始行，分別從每一行的中間位置開始，對其左（右）側的圖片進行標誌位元掃描，遇到標記圖片，則利用該標記圖片右（左）側相鄰兩張圖片的變換矩陣與標記圖像之間的相對位置關係估計出標記圖片的變換矩陣，同時位於標記圖片左（右）側的所有圖片的變換矩陣均通過該方法依次得到。直至每一行圖像都重複該操作為止。

（3）從圖像中間列的中間位置開始，對中間列圖像的上下兩部分的圖片則重複步驟（1）和（2）。

完成上述處理之後，通過相應的融合演算法對結果圖進行融合，消除拼接縫［11］，完成整個岩石薄片圖像的拼接。

3實驗對比結果及分析

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

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

4結論

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

參考文獻

［1］ 徐鑫, 孫韶媛, 沙鈺傑，等. 一種基於RANSAC的紅外圖像拼接方法［J］. 鐳射與光電子學進展，2014，51(11):129-134.

［2］ 嶽永娟, 苗立剛, 彭思龍. 大規模顯微圖像拼接演算法［J］. 電腦應用，2006, 5(26):1012-1014.

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

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

［5］ 白廷柱, 侯喜報. 基於SIFT運算元的圖像匹配演算法研究［J］. 北京理工大學學報，2013，33(6):622-627.

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

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

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

［9］ 曲天偉, 安波, 陳桂蘭. 改進的RANSAC演算法在圖像配准中的應用［J］. 電腦應用，2010，30（7）：1849-1851.

［10］ 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．

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

（2）定位關鍵點，去除不穩定點

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

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

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

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

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

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

如果曲率大於式（10）結果的比值，則捨棄。LOWE D G在實驗中通過對比發現，當r=10時，能夠最大程度去除不穩定的邊緣響應點［6］。

（3）確定關鍵點的大小和方向

在初步確定了關鍵點的位置之後，可以利用圖像的局部特性為每一個關鍵點分配一個基準點，使該關鍵點的特徵描述符具有旋轉不變性。對於高斯金字塔檢測出來的關鍵點，使用其高斯金字塔鄰域視窗內圖元的梯度特徵和方向特徵對描述符進行描述。(x,y)處的梯度值m(x,y)和方向θ(x,y)的運算式如式（11）和式（12）所示。

其中，L所用的尺度為每個關鍵點各自所在的尺度。為了進一步準確給出關鍵點的方向，一般採用梯度長條圖統計的方法，以關鍵點為中心的鄰域視窗內的圖元的梯度和方向。長條圖峰值所代表的方向即為關鍵點的方向。通過上述各個步驟之後，檢測出的包含位置、所在尺度以及方向的關鍵點稱為圖像的特徵點。

（4）生成特徵描述子

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

（5）特徵點的匹配

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

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

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

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

A1=MA2（14）

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

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

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

2空白區域圖像拼接方法

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

圖4拼接流程圖在實際應用中，為了解決該問題，在SIFT演算法的基礎上，提出了一種選取岩石薄片圖像的中心位置為基準點的拼接方法。並利用空白區域圖像與其相鄰圖像之間的位置關係估計出空白區域圖像的變換矩陣，實現了整個岩石薄片的拼接。流程圖如圖4所示。

2.1特徵點提取

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

（1）獲取待拼接岩石薄片圖像的行和列的數目。

（2）分別對岩石薄片圖像中的每一行和中間列進行特徵提取，並對不能正常匹配的圖片進行標記。

2.2計算變換矩陣

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

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

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

（3）從中心位置所在行開始到圖像的最後一行，分別從每一行的中間位置開始重複步驟（2）。

（4）從圖像中間列的中間位置開始，對中間列圖像的上下兩部分的圖片重複步驟（1）和（2）。

2.3優化變換矩陣

計算完所有圖片的變換矩陣之後，本文在上一節提出的變換矩陣的計算順序的基礎上，對變換矩陣的優化方式作了相應的調整。具體優化方式如下：

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

（2）從中心位置所在行開始到圖像的最後一行，分別從每一行的中間位置開始重複步驟（1）。

（3）從薄片圖像中間列的中間位置開始，對中間列圖像的上下兩部分的圖片重複步驟（1）和（2）。

2.4估計空白區域圖片的變換矩陣

變換矩陣優化完成之後，為了獲得完整的拼接結果，需要估計空白區域圖片的變換矩陣。具體方法如下。

（1）從中心位置所在行開始到圖像的起始行，分別從每一行的中間位置開始，對其左（右）側的圖片進行標誌位元掃描，遇到標記圖片，則利用該標記圖片右（左）側相鄰兩張圖片的變換矩陣與標記圖像之間的相對位置關係估計出標記圖片的變換矩陣，同時位於標記圖片左（右）側的所有圖片的變換矩陣均通過該方法依次得到。直至每一行圖像都重複該操作為止。

（3）從圖像中間列的中間位置開始，對中間列圖像的上下兩部分的圖片則重複步驟（1）和（2）。

完成上述處理之後，通過相應的融合演算法對結果圖進行融合，消除拼接縫［11］，完成整個岩石薄片圖像的拼接。

3實驗對比結果及分析

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

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

4結論

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

參考文獻

［1］ 徐鑫, 孫韶媛, 沙鈺傑，等. 一種基於RANSAC的紅外圖像拼接方法［J］. 鐳射與光電子學進展，2014，51(11):129-134.

［2］ 嶽永娟, 苗立剛, 彭思龍. 大規模顯微圖像拼接演算法［J］. 電腦應用，2006, 5(26):1012-1014.

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

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

［5］ 白廷柱, 侯喜報. 基於SIFT運算元的圖像匹配演算法研究［J］. 北京理工大學學報，2013，33(6):622-627.

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

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

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

［9］ 曲天偉, 安波, 陳桂蘭. 改進的RANSAC演算法在圖像配准中的應用［J］. 電腦應用，2010，30（7）：1849-1851.

［10］ 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．

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