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

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

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

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

　　中圖分類號：TP391.4文獻(xiàn)標(biāo)識碼：ADOI： 10.19358/j.issn.1674-7720.2017.06.015

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

0引言

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

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

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

1相關(guān)理論

　　1.1SIFT算法的基本原理

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

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

　　尺度空間是模擬圖像數(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）

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

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

　　LOWE D G使用高斯差分算子(DoG)主要有兩個原因：一是能夠更方便地檢測出穩(wěn)定的關(guān)鍵點(diǎn)；二是考慮到計(jì)算的復(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在實(shí)驗(yàn)中，通過對原始圖像進(jìn)行降采樣構(gòu)造若干階尺度空間，如圖2（a）所示，然后在每一階尺度空間內(nèi)利用高斯函數(shù)卷積得到若干層高斯圖像，再差分生成高斯差分尺度空間，如圖2（b）所示。在同一個尺度內(nèi)，將每一個像素點(diǎn)同時與該像素點(diǎn)同尺度的8個相鄰點(diǎn)以及上下兩個相鄰尺度對應(yīng)的18個像素點(diǎn)的值進(jìn)行比較，只有當(dāng)該像素點(diǎn)的值是最大或者最小值時，才稱該像素點(diǎn)是圖像在該尺度下的一個特征點(diǎn)，也叫關(guān)鍵點(diǎn)。最后，對其余的每一階尺度空間都重復(fù)以上操作，直到檢測完所有尺度空間中的關(guān)鍵點(diǎn)為止。

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

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

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

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

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

　　因?yàn)楦咚共罘炙阕釉谶吘壧帟a(chǎn)生比較強(qiáng)烈的邊緣響應(yīng)，在該算子中，差的極值點(diǎn)在邊緣處具有較大的曲率，在垂直方向上具有較小的曲率。因此，可以利用該性質(zhì)去除邊緣響應(yīng)。主曲率可以通過一個2×2的Hessian矩陣H求出，如式（7）所示。

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

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

　?。?）確定關(guān)鍵點(diǎn)的大小和方向

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

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

　?。?）生成特征描述子

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

　　（5）特征點(diǎn)的匹配

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

　　1.2變換矩陣計(jì)算的基本原理

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

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

　　A1=MA2（14）

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

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

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

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

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

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

　　2.1特征點(diǎn)提取

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

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

　　（2）分別對巖石薄片圖像中的每一行和中間列進(jìn)行特征提取，并對不能正常匹配的圖片進(jìn)行標(biāo)記。

　　2.2計(jì)算變換矩陣

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3實(shí)驗(yàn)對比結(jié)果及分析

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

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

4結(jié)論

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

　　參考文獻(xiàn)

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

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

　?。?］ 鄒承明, 侯小碧, 馬靜. 基于幾何學(xué)圖像配準(zhǔn)的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.

　　［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.

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

　?。?］ 曲天偉, 安波, 陳桂蘭. 改進(jìn)的RANSAC算法在圖像配準(zhǔn)中的應(yīng)用［J］. 計(jì)算機(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．

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