FIB-SEM圖像多策略層間配準與插值算法

張 耀 滕奇志 卿粼波 吳曉紅

(四川大學電子信息學院圖像信息研究所 四川 成都 610065)

0 引 言

隨著世界常規油氣產量的持續下降與油氣勘探技術的發展，以頁巖氣、致密油、煤層氣等為代表的非常規油氣能源逐漸引起各國的重視[1]。與常規油氣資源不同，非常規油氣儲集層多以納米級孔喉為主，局部發育為微米-毫米級，因此微納米孔隙空間的結構特征是影響儲集層的重要因素，準確全面地表征微納米孔隙結構已成為非常規油氣研究的重要內容[1]。三維圖像重建技術可用于表征巖石內部的微觀結構，為研究非常規油氣提供了有力的技術支持。

聚焦離子束(FIB)與掃描電鏡(SEM)雙束系統，通過離子束的納米級表面切割和電子束的超高分辨率成像[2]，可達到精確重建微納米級孔隙的目的。但由于雙束系統的工作特性，其重建后的三維圖像縱向分辨率與橫向分辨率不一致，且大多表現為縱向分辨率偏低，這將對微納米孔隙后續研究產生較大影響，因此提高FIB-SEM三維圖像縱向分辨率顯得尤為重要。

提高三維圖像縱向分辨率的方法為層間插值，即通過向空間相鄰的兩幅圖像間插入新圖像的方式提高縱向分辨率。目前三維圖像層間插值方法主要有[3-4]：基于灰度插值、基于小波插值、基于配準插值。

基于灰度的插值方法計算復雜度低且易于實現，如線性插值與三次樣條插值，但插值精度較低。

基于小波的插值方法通過對小波系數的操作來實現層間插值，生成的圖像質量較高，但小波系數最多可表示一維點奇異信息，對圖像中重要的線段、輪廓邊界等二維信息的描述不能達到最優[5-6]，對于圖像高頻信息的插值結果存在較大誤差[6]。Cunha A L D等[7]針對小波變換的局限性，提出了非下采樣輪廓波變換方法NSCT，可準確捕獲圖像各方向的二維輪廓信息。

由于基于小波的插值方法的局限，王本峰等[8]將Curvelet變換與凸集投影算法相結合，較好地改善了 FIB-SEM序列圖層間距過大情況，但此算法僅可用于頁巖切片，通用性不佳。

基于配準的插值方法[9-12]在醫學圖像領域應用較為廣泛，其利用形態學相關知識，于插值前先對圖像進行形態學分析[9]，以獲得待插值圖像像素點的形變信息。此類方法依賴于相鄰層間解剖學結構的形狀變化及灰度變化為連續的假設。對于FIB-SEM圖像，由于FIB對巖心樣品縱深切割間隔為納米級，相鄰序列圖的特征結構形狀與灰度的相對變化具有較強的連續性，因此基于配準的插值方法適用于FIB-SEM圖像。

當前基于配準的插值算法均存在配準過程中出現異常點的問題[12]，且由于配準算法本身未考慮圖像全局相關性，因此插值后還會存在圖像高頻細節信息丟失現象[13-15]。為解決高頻細節丟失現象及運算耗時較長的問題，文獻[15]提出了融合多分辨率分析與曲率配準的層間插值算法，一定程度上提升了插值效率。但由于每層分辨率的配準結果均基于上一分辨率的配準結果，因此插值效果較大地依賴于最低級分辨率圖像形變場，算法魯棒性較低；且算法須對修改原圖并手動調整分辨率級數，但即使在最優分辨率級數下，圖像重要細節依然會丟失。

對于異常點問題，文獻[12]提出了基于曲面擬合的配準插值算法，在一定程度上抑制了配準過程中異常點的生成，但圖像高頻信息丟失現象依然存在，重建的圖像在一定程度上依然存在模糊現象。

針對當前基于配準的插值算法中存在的高頻細節丟失現象，本文提出一種針對FIB-SEM圖像的多策略層間配準與插值算法。該算法首先對原始FIB-SEM序列圖進行NSCT分解，在不破壞原圖信息的基礎上將原圖重要高頻特征與低頻概貌剝離，對二者采用不同的插值算法策略進行針對性處理。結合像素點空間位置信息與NSCT系數，對NSCT低頻子圖采用基于配準的層間插值算法進行插值，對NSCT高頻方向子帶采用多項式擬合算法進行插值，最終將新插值圖像的低頻子圖及高頻方向子帶利用逆NSCT生成新切片圖。該算法將NSCT變換與基于配準的插值算法有機結合，并融入多項式擬合方法，充分利用連續多張序列圖的高頻信息，較好地抑制了配準過程中出現異常點以及圖像高頻信息缺失的情況，提高了插值效率與生成圖片的質量。

1 NSCT

NSCT繼承自Contourlet變換，解決了Contourlet變換中由于下采樣操作而導致的圖像Gibbs效應以及低頻頻譜泄漏等問題[8]。NSCT內部采用非下采樣塔式濾波器組NSPFB(Nonsubsampled Pyramid Filter Bank)及非下采樣方向濾波器組NSDPFB(Nonsubsampled Directional Filter Bank)來捕捉二維圖像幾何結構與紋理細節，最終以NSCT系數來表征各方向的圖像輪廓信息。NSCT首先利用NSPFB對圖像進行金字塔分解以便捕捉圖像奇異點信息，分解將獲得一個低頻圖和一個帶通圖；然后利用NSDFB將帶通圖沿一系列指定方向進行分解，提取二維方向性高頻信息(即圖像輪廓信息)[8]。

每一層NSFDB分解將產生一個非下采樣低通圖b和b與原始圖的差值圖(帶通圖)，對b繼續降采樣分解可得下一層的低頻圖及帶通圖，如圖1(a)所示。LP分解后可獲得一系列帶通圖與一個低頻子圖，NSDFB將對每個帶通圖進行處理以提取二維方向性特征，如圖1(b)所示。對于某一張帶通圖，NSDFB可將其分解為2的任意次冪個方向，即生成2的任意次冪個高頻方向子帶。

圖1 NSCT

2 結合NSCT的多策略層間配準與插值

2.1 基于配準的NSCT低頻系數插值方法

針對當前基于配準的插值算法普遍存在的高頻細節丟失以及出現異常點的情況，本文算法將NSCT與配準算法結合，對圖像低頻概貌與圖像高頻信息“解耦剝離”，以提升后續步驟中圖像配準質量及輪廓信息的插值精確度。令FIB-SEM切片圖的低頻子圖為αj0，NSCT第j層第k方向高頻子帶圖為βj,k，則原始FIB-SEM圖像素點的NSCT系數可表示為：

(1)

(2)

式中：k=?x/δx」-1;l=?y/δy」-1;u=x/δx-?x/δx」;v=y/δy?y/δy」;Bi表示第i階B樣條基函數；δx與δy分別為x方向與y方向的網格間距。如圖2所示，(a)與(c)為相鄰原始圖像，(b)為以(a)為浮動圖像并以(c)為參考圖像獲得的形變場；(d)與(f)分別為(a)與(c)的NSCT低頻子圖，(e)為以(d)為浮動圖像并以(c)為參考圖像獲得的形變場。由于?i,j決定非剛體形變程度，因此須對Φd進行優化。基于對FIB-SEM圖像具有非剛性形變特性的考慮，本文采用Rueckert等[16]提出的針對非剛性形變的基于歸一化互信息測度與平滑函數的控制點自適應移動優化方法對形變場進行迭代優化：

(3)

(4)

▽C=?C(Φd)/?Φd

(5)

Φd=Φd+μ▽C/‖▽C‖

(6)

式(4)中，H代表圖像熵或圖像聯合熵。式(3)中，D為圖像面積。根據修正后的Φd重新計算梯度向量直至其模小于給定閾值結束。

圖2 相鄰FIB-SEM序列圖配準

經過優化后的形變場TΦd表征相鄰圖像中像素的位移變化，位移變化即為位移向量Kn。如圖3所示，(a)為根據圖2(e)形變場繪制的圖2(d)與圖2(f)的位移向量二維投影示意圖，為方便讀者觀察圖像配準結果，圖中位移向量經過同比例放大處理。圖3(b)為位移向量示意圖，其中虛直線即表示位移向量。根據待插值圖像在左右兩圖的空間相對位置參數η(0<η<1)以及每個像素的位移向量Kn(n為NSCT低頻子圖上像素的位置)，可確定待插值圖像的NSCT低頻系數值與位置：

iη=(1-η)id+ηid+1

(7)

(8)

(9)

式中：φ3代表三次樣條插值。

圖3 位移向量示意圖

2.2 NSCT高頻系數層間插值方法

由于NSCT將圖像低頻概貌與高頻細節剝離，對于NSCT高頻子帶圖βj,k，其表征圖像在某一方向的輪廓細節信息。由于FIB-SEM圖像對巖心內部結構或孔隙輪廓在縱深方向變化具有良好的連續性，出于對巖心結構縱深連續性的考慮，本文根據Hermite多項式擬合插值算法原理，結合NSCT高頻系數特征，對待插值圖高頻部分βj,k進行處理。由于多項式擬合插值方法存在龍格現象[17]，因此本文采用三次Hermite插值方法進行高頻系數插值。

以兩節點為例，設插值節點集為{x0,x1,…,xn}，相應的節點值集為{f0,f1,…,fn},設區間[xk-1,xk]位于[x0,xn]內，且兩邊端點值分別為fk與fk-1，則可定義此區間上的三次Hermite插值函數Hk(x)如下：

(10)

式中：dk-1與dk為插值函數Hk(x)在[xk-1,xk]端點處的一階導數。對于Hk(x)，其由區間端點函數值和一階導數決定。對于區間[x0,xn]內中間節點xi(i=1,2,…,n-1),其相應導數可用左右相鄰區段的一階差商加權方式近似：

(11)

對于區間端點，無法獲得其雙側的一階差商，因此無法使用上式方法計算導數值。針對此情況，本文采取區間端點雙側一階差商相等的策略。

本文針對FIB-SEM圖像高頻系數插值節點數的選取問題，對大量FIB-SEM圖像進行了實驗驗證對比分析。實驗結果表明，針對FIB-SEM圖像，其高頻系數插值節點數取4為最佳(即考慮連續4張FIB-SEM圖像的輪廓變化信息)，可充分利用巖心縱向結構輪廓的連續變化信息。因此本文對NSCT高頻部分的處理均采用4節點插值(即“四圖策略”)。

結合上文所述，FIB-SEM圖像多策略層間配準與插值算法可總結為以下幾個步驟：

1) 確定原始FIB-SEM序列圖中新層間圖的空間位置xη,與新層間圖關聯的4幅原始圖組成一個關聯組Fxη。

2) 對原始序列圖執行NSCT變換，獲得低頻子圖與高頻方向子帶圖。

3 實驗結果與分析

3.1 NSCT低頻子圖配準效果驗證實驗

為驗證本文算法低頻子圖配準部分處理效果，對6組原始FIB-SEM序列圖(圖像均由中石油勘探開發研究院提供)，分別進行原始圖像配準實驗和NSCT低頻子圖配準實驗。對于此兩種實驗，其配準方式均采用以后一張圖像為參考圖并以前一張圖為浮動圖的方式進行配準，并將配準優化后的結果圖與參考圖作比較。實驗結果采用業界廣泛認可的均方差MSD(Mean Square Difference)與互相關系數CC(Cross Correlation)[18]：

(12)

式中：I1為原始切片圖，I2為層間插值后圖像，切片圖像大小為m×n，R(i,j)代表參考圖(i,j)處的像素灰度值，F(i,j)代表浮動圖(i,j)處的像素灰度值。MSD值越小，CC值越趨近于1，則圖像配準精度越高，效果越好。由表1可知，對原始FIB-SEM圖的NSCT低頻圖像配準，其配準精度與配準效果要優于直接對原始FIB-SEM圖像配準。實驗結果如表1所示。

表1 原始FIB-SEM圖配準與NSCT低頻圖配準比較

續表1

3.2 FIB-SEM圖像配準與插值實驗

為驗證本文算法的正確性和有效性，將本文算法與線性插值法、penny插值法[11]及文獻[12]提出的USFBR算法作比較。線性插值法是被廣泛使用的基于場景插值的方法，penny插值法則為被廣泛認可的基于配準的層間插值算法，而USFBR算法則為近期在層間插值算法研究領域中一種優秀的基于配準的層間插值算法。本文共進行6組實驗，實驗中涉及的圖像均為FIB-SEM圖像。實驗結果如表2所示。

表2 切片圖MSD與PSNR比較

樣本圖像大小均為512×512。為保證實驗數據結果真實性，輸入均為FIB-SEM原始掃描圖像。實驗結果評判標準參數采用業界廣泛認可的均方差MSD與峰值信噪比PSNR：

(13)

6組實驗均選取五張連續原始圖像的最前兩張與最后兩張來重建中間第三張，并將重建圖與原始圖進行比較。為使讀者更加清晰地觀察比較4種方法的插值結果的區別，本文在圖5中給出了圖4白色矩形框標識區域的放大圖。

圖4 2組實驗的結果圖

為客觀比較4種算法，實驗計算了4種算法生成的各新層間圖的MSD與PSNR，如表2所示。根據表2結果可知，本文算法在客觀評價標準上優于其他對比算法，PSNR相較于USFBR算法平均提升0.48 dB以上。

限于篇幅，圖5列出了6組實驗中2組實驗的結果示意圖：(a)-(e)為頁巖A樣本序列中某一張圖片的重建結果；(f)-(i)為致密碳酸鹽巖B樣本序列某一張圖的重建結果(為便于讀者觀察各算法結果的差異，此樣本重建結果的放大圖像經過亮度增強)。根據圖5所示結果可以看出，線性插值生成的圖像不僅會出現圖像信息明顯丟失的現象(尤以圖5(b)最為顯著)，而且圖像存在拖影的視覺現象和模糊現象；penny算法相較于線性插值算法，拖影現象明顯減少，孔隙準確率有所提升，但模糊現象更為嚴重；USFBR算法相較于penny算法在圖片質量上有了進一步提升，塊狀異常區域減少，但模糊現象仍存在；本文算法在USFBR算法基礎上有了進一步提升，減小了亮度差異較大區域與原圖的差異，減少了模糊現象。因此，本文算法是此4種算法中最優的。

4 結 語

本文針對FIB-SEM圖像，提出多策略層間配準與插值算法。該算法利用NSCT可有效捕獲二維圖像輪廓信息的特點，將圖像高頻輪廓細節與低頻概貌進行分離。低頻概貌部分結合基于配準的插值算法，提升低頻部分配準與插值的準確度；高頻輪廓細節部分利用融入端點優化策略的分段三次Hermite多項式，提升新層間圖像輪廓結構的精確度與質量。實驗結果顯示，本文算法的實驗結果優于各對比算法的實驗結果與其他對比算法對比，PSNR平均提升0.48 dB以上。但本文算法須對圖像進行多層結構分解，并于分解過程中引入非下采樣技術以減少圖像信息的損失，算法內部涉及迭代運算，若對原始FIB-SEM序列圖采用串行運算方式進行處理，則運算時耗大于線性插值，因此本文算法在速度上需進一步優化。

參 考 文 獻

[1] 賈承造,鄭民,張永峰.中國非常規油氣資源與勘探開發前景[J].石油勘探與開發,2012,39(2):129-136.

[2] 孫亮,王曉琦,金旭,等.微納米孔隙空間三維表征與連通性定量分析[J].石油勘探與開發,2016,43(3):490-498.

[3] Ale? N,Olivier S,Oscar A,et al.Constrained reverse diffusion for thick slice Interpolation of 3D volumetric MRI images [J].Computerized Medical Imaging & Graphics,2012,36(2):130-138.

[4] Jafari-Khouzani K.MRI Upsampling Using Feature-Based Nonlocal Means Approach[J].IEEE Transactions on Medical Imaging,2014,33(10):1969-1985.

[5] Chang S G,Cvetkovic Z,Vetterli M.Locally adaptive wavelet-based image interpolation[J].IEEE Transactions on Image Processing,2006,15(6):1471-1485.

[6] 武士想,尚鵬,王立功,等.小波-Lagrange方法進行醫學圖像層間插值[J].中國圖象圖形學報,2016,21(1):78-85.

[7] Cunha A L D,Zhou J,Do M N.The nonsubsampled contourlet Transform:Theory,Design,and Applications[J].IEEE Transactions on Image Processing,2006,15(10):3089-3101.

[8] 王本鋒,李景葉,陳小宏,等.基于Curvelet變換與POCS方法的三維數字巖心重建[J].地球物理學報,2015(6):2069-2078.

[9] Grevera G J,Udupa J K.An objective comparison of 3D image interpolation methods[J].IEEE Transactions on Medical Imaging,1998,17(4):642-52.

[10] Herman G T,Zheng J,Bucholtz C A.Shape-Based Interpolation[J].IEEE Computer Graphics & Applications,1992,12(3):69-79.

[11] Penney G P,Schnabel J A,Rueckert D,et al.Registration-based interpolation[J].IEEE Transactions on Medical Imaging,2004,23(7):922-6.

[12] 張杰,何小海,卿粼波,等.基于曲面擬合及配準的醫學圖像層間插值[J].四川大學學報(工程科學版),2016(S1):142-149.

[13] Frakes D H,Dasi L P,Pekkan K,et al.A new method for registration-based medical image interpolation[J].IEEE Transactions on Medical Imaging,2008,27(3):370-377.

[14] Laursen L F,Hansen M S.Registration-based interpolation realtime volume visualization[C]//Spring Conference on Computer Graphics [s.n.],2013:15-21.

[15] 王新征,雄洙,于靖.結合多分辨率修正曲率配準的層間插值[J].光學精密工程,2016,24(5):1224-1231.

[16] Rueckert D,Sonoda L I,Hayes C,et al.Nonrigid registration using free-form deformations: application to breast MR images[J].IEEE Transactions on Medical Imaging,1999,18(8):712-21.

[17] Epperson J F.On the Runge Example[J].American Mathematical Monthly,1987,94(4):329-341.

[18] 高廣珠,李忠武,余理富,等.歸一化互相關系數在圖像序列目標檢測中的應用[J].計算機工程與科學,2005(3):38-40.