巖心三維圖像修復算法

岳桂華，滕奇志，何小海，陳冬冬

（四川大學電子信息學院，610065，成都）

在石油地質分析過程中，通常需要獲得三維巖心結構來定量地研究滲流的微觀機理，但在實際工作中，有很多種因素造成巖心的缺失或缺損，導致很難獲得完整的可供分析的巖心三維圖像。如CT（computed tomography）機在巖石三維成像的應用中，其分辨率和樣本尺寸是相矛盾的，為了得到高分辨率的孔隙結構圖像，掃描的樣本尺寸受到很大限制，使巖樣的代表性有所欠缺。針對這一問題，目前的主要解決方案是以二維薄片為基礎的三維重建［1-3］，得到具有與薄片特性相近的三維圖像。但是，該方法是基于圖像局部特性的數學建模，重建結果與真實巖心的形態結構存在一定差異，且對于多相、非均質以及原始灰度的巖心三維圖像的重建存在較大的困難。將大尺寸樣本切割為多個小樣本，用CT分別對小樣本進行掃描，再將獲取的小樣本的高分辨率圖像拼接成較大的三維圖像，能有效地解決這一問題。然而，由于物理切割導致三維圖像存在較大范圍的信息缺失，在小樣本圖像拼接成大樣本圖像的過程中需要恢復缺損的信息。因此，可考慮采用圖像修復的方法修復出缺損的信息，從而保持原始巖心三維圖像的結構和細節特征。

圖像修復［4］目前多針對二維圖像，對于三維圖像修復的研究還較少，且主要是針對醫學圖像小范圍信息缺失的三維圖像修復以及視頻的修復。目前應用較多的就是基于偏微分方程的圖像修復算法和基于樣例的圖像修復算法。基于偏微分方程的三維圖像修復算法［5］只適合小范圍缺失且細節信息較少的三維圖像修復，對于較大范圍缺失圖像的修復會存在明顯的模糊現象，且不能修復圖像的紋理，因此不適合大范圍缺失的巖心三維圖像修復。基于樣例的圖像修復算法［6-10］對于存在較大范圍缺失的二維圖像和視頻圖像取得了很好的效果，其中的Criminisi算法［6-7］及其改進算法通過計算邊界點的優先級值來決定填充的順序，能夠有效地修復較大范圍缺失的圖像。目前的基于樣例的圖像修復算法較多地應用在視頻修復中，然而視頻中每一幀都是具有完整結構和細節的二維圖像，算法主要是逐幀進行，修復過程是針對二維結構和紋理的修復［9-10］。巖心三維圖像的三維空間結構與視頻序列圖像存在本質的區別，因而不能采用傳統的針對視頻的基于樣例的圖像修復方法來修復巖心三維圖像。

本文針對存在較大范圍信息缺失的巖心三維圖像，提出了一種基于樣例的三維圖像修復算法，采用三維空間結構的模板同時填補三維圖像缺失部分的結構和細節，通過計算邊界體素相應模板的優先級值來決定填充的順序。同時，通過空間中的等照度面對優先級值的影響，優先合成與傳播圖像的結構，實現三維圖像結構的連接，且避免了對于等照度線方向的討論。針對Criminisi算法中采用全局搜索方法計算量大的問題，本文使用匹配中心因子和最大距離來限制匹配塊搜索范圍，結合了局部搜索和全局搜索的優點，大大減少了計算量。該算法不僅保持了原始巖心三維圖像的結構和細節形態特征，而且保持了巖心的微觀統計特性［11］。

1 巖心三維圖像修復算法

本文提出的巖心三維圖像修復算法的基本思想是：從待修復區域（目標部分）Ω的邊界δΩ上選取邊界體素p（p∈δΩ），以p為中心，設定固定大小的立方體模板Ψp，在p的匹配中心因子和最大距離限定的搜索范圍內搜索其最優匹配塊Ψq，并用Ψq中的已知體素填充Ψp中相應位置上的未知體素。

對于原始三維圖像中的每一個體素q都存在一個灰度值（當體素沒有被填充時為空）和置信度C（q）。C（q）表示q的灰度值的可信度，初始化為C（q）=0，?q∈Ω和C（q）=1，?q∈Φ。每一個邊界體素p的模板（邊界模板）Ψp都有一個優先級值P（p）。本文提出的算法使用了Criminisi算法的關鍵思想，即由邊界模板的優先級值來決定邊界體素的填充順序。文中提及的體素指三維數據場中的采樣點。Ψp的最優匹配塊Ψq就是與Ψp最為相似的樣本塊，即

式中兩圖像塊的距離d（Ψp，Ψ^q）為邊界模板Ψp與樣本塊Ψ^q對應位置的已知體素的灰度值之差的平方和。

1.1 填充順序

與Criminisi算法利用等照度線的傳播來引導紋理合成的過程［6-7］類似，為了向目標區域傳播大尺度的結構，本文的算法根據三維圖像的等照度面給邊界體素設定相應的填充順序。等照度面指三維圖像中相同亮度或強度的面，等照度線指相同亮度或強度的線。邊界體素的填充順序由其相應的邊界模板的優先級值大小決定，每次都搜索優先級值最大的邊界模板進行填充。邊界模板Ψp的優先級值定義為

式中：D（p）為數據項，表示等照度面在邊界面處的強度；C（p）為邊界模板的置信度，用于衡量體素p周圍模塊范圍內可靠信息的數量，表示為

其中q∈Ψp∩表示Ψp中已被填充的體素，|Ψp|表示Ψp的體積。

在二維圖像修復中，Criminisi算法取等照度線在邊界垂直方向上的強度作為數據項D（p），等照度線的方向被定義為梯度方向的垂直方向。在三維圖像中，等照度線的軌跡是模糊不清的，且垂直于梯度的是一個平面，存在無限多個可能的方向［12］。如圖1所示，在某一小塊的三維空間中存在一連續的等照度面S，取以體素q為中心的任一等照度面塊ds，當區域足夠小時，ds為平面或存在切面T 與梯度向量▽Ip垂直，等照度面切線方向具有不確定性，如都為ds切線，且都與▽Ip正交。在連續的邊界面上，取足夠小的區域得到邊界平面δΩ，如圖2所示，其與邊界體素p所在的等照度面存在一定的夾角θ，則δΩ的法線np與梯度向量▽Ip的夾角也為θ。同時，等照度面的強度||等于梯度▽Ip的大小|▽Ip|。因此，等照度面在邊界面法線方向上的強度可表示為

式中：α為歸一化因子，取灰度等級為255，這樣就避免了等照度線方向的討論。當等照度面越強且與邊界面越接近垂直狀態，數據項D（p）越大，算法會優先選擇這樣的邊界體素進行填充，即算法對結構性較強的地方會優先修復。數據項D（p）在算法中起著根本性的作用，通過D（p）實現了等照度面對填充順序的影響，使得圖像的結構（也伴隨細節）優先得到合成與傳播，從而將破損的結構連接起來，達到視覺上的連續性。

圖1 三維空間中的連續等照度面

圖2 D（p）的說明圖

圖3 匹配中心因子的原理和設置說明圖

1.2 匹配中心因子和最大距離

通常圖像的結構具有連續性［13］，圖像的某個體素值與其周圍的體素值存在一定的相關性，即圖像的某個體素（或像素）的鄰域具有充分的信息決定該體素值（或像素值）［8］。因此，若Ψq為邊界模板Ψp的最優匹配塊，p′為Ψp填充后產生的邊界體素，則相應的邊界模板Ψp′在Ψp周圍，如圖3所示，那么相應的最優匹配塊Ψq′也通常在Ψq的周圍，即體素q′通常在q的周圍。另一方面，填充Ψp后，新的邊界體素p′的邊界模板Ψp′中有1／8至1的已知體素與Ψp重疊，若Ψp與Ψq完全匹配，那么Ψp′與Ψq鄰近的相應位置的樣本塊Ψq′的匹配度為1／8至1。當修復圖像的結構時，圖3中的p、p′1、p″，新產生的邊界模板中的所有已知體素都由剛填充的模板決定，與圖中相應位置的樣本塊完全匹配。修復過程中，圖像的結構和細節將不斷地向目標區域傳播，新的邊界模板的信息將由周圍已填充的模板決定。因此，本算法提出了采用匹配中心因子來限定匹配塊的搜索范圍，將匹配塊搜索限定為局部搜索。匹配中心因子（MCF）是邊界體素相應的匹配塊局部搜索區域的中心體素。當MCF存在時，邊界體素的匹配塊局部搜索區域為以MCF為中心、大于或等于兩倍模板大小的立方體區域。在初始化時，各邊界體素的MCF為空，需全局搜索最優匹配塊。當填充邊界模板Ψp后，Ψp的最優匹配塊Ψq的中心體素q即為新產生的邊界體素p′的匹配中心因子。

匹配中心因子的原理和設置過程如圖3所示。圖3a中的δΩ為初始邊界（白色缺損部分），算法會優先修復圖像中結構性較強的地方（如p點），p的MCF初始化為空，匹配塊的搜索區域為整個圖像的已知區域，搜索得到的最優匹配塊是以點q為中心的樣本塊Ψq。填充完Ψp后，如圖3b所示，新的邊界點p′（如p′1和p′2）的 MCF設置為q，局部搜索得到的最優匹配塊為以q′（如q′1和q′2）為中心的Ψq′。對于整個圖像而言，Ψq′為Ψp′的最優匹配塊。匹配中心因子將匹配塊搜索限定為局部搜索，大大縮小了匹配搜索范圍。隨著修復的進行，如圖3c所示，填充Ψp′后新產生的邊界體素p″的 MCF為q′1，最優匹配塊為Ψq″。

最大距離表示在邊界體素p的匹配中心因子不為空的情況下，局部搜索得到的最優匹配塊與其邊界模板的距離的最大范圍，作為跳出局部搜索的條件限制。即當局部搜索得到的最優匹配塊與邊界模板的距離大于最大距離時，表明此處的圖像結構發生了突變，該匹配塊將不能作為其最優匹配塊，需重新全局搜索其最優匹配塊。邊界體素p′的最大距離表示為

式中：μ為比例參數，通常設定為1或臨近的值；p′為填充邊界模板Ψp后新產生的邊界體素，根據Ψp和其最優匹配塊Ψq的距離，設置p′的最大距離。在設置p′的匹配中心因子的同時，設置其最大距離。該方法結合局部搜索和全局搜索的優點，大大縮小了匹配搜索范圍，同時保證了匹配結果的準確性。

1.3 三維圖像修復算法的主要步驟

步驟1：設定模板大小、局部搜索區域大小，并提取待修復區域Ω的初始邊界δΩ0。

步驟2：計算邊界輪廓上所有體素的邊界模板的優先級值P（p）。

步驟3：搜索最大優先級值的邊界模板Ψp。

步驟4：根據邊界體素p的MCF和最大距離搜索最優匹配塊Ψq，并填充相應體素的灰度值。

步驟5：更新或設置置信度、邊界、新產生的邊界體素的MCF和最大距離，以及邊界模塊的優先級值的相關信息。

步驟6：重復步驟3～5直至δΩ=?。

選取的模板具有三維空間結構，通常為固定大小的立方體，其大小設定需有一定經驗。實驗中，我們將局部搜索區域設置為2倍模板大小的立方體區域。Ψp中新填充的體素的置信度設置為Ψp模板的置信度，即

當填充完邊界模板Ψp后，只有Ψp周圍體素的缺損情況和信息發生了變化，其余邊界體素不變，因此只需搜索Ψp內及其相鄰體素產生的新的邊界體素，并去除已填充的邊界體素，從而得到新的邊界，然后計算在以p為中心、2倍模板大小的立方體鄰域內的所有邊界體素的模板的優先級值。整個三維圖像修復算法的主要過程就是對步驟3～5的迭代，只需在算法的開始設定模板大小和局部搜索區域大小，并提取出初始邊界，修復過程將會自動進行直至修復完成。

2 實驗結果

首先將算法應用于存在較大范圍信息缺失且結構較為簡單的三維圖像的修復，從而展示算法的修復過程和結果，以此驗證算法的正確性。在此基礎上，再對存在較大范圍信息缺失的巖心三維圖像進行修復。為了更加直觀、準確地評估巖心三維圖像的修復結果，將真實的巖心CT序列圖人為地去除中間部分，模擬巖心三維圖像的信息缺失，并將修復結果與真實的原始巖心CT序列圖進行直觀地對比觀察。同時，用文獻［5］中的三維圖像修復算法對巖心三維圖像進行了相應的修復實驗，將實驗結果與本文提出的算法的實驗結果進行了簡單、直觀的對比。巖心三維圖像的修復不僅要保持視覺上的連續性，達到圖像修復的目的，同時要保持原始巖心的微觀統計特征。由于實際應用中時常需要獲取巖心三維圖像的孔隙特性，且為了定量地統計分析各組實驗的巖心微觀統計特征，本文計算了修復結果的孔隙結構的兩點能量，以及對比了修復結果與原始巖心三維圖像的絕對孔隙度和有效孔隙度。

圖4 結構簡單的三維圖像的修復過程圖

2.1 結構簡單的三維圖像

圖4為簡單規則的三維圖像的修復過程圖，第1幅圖為原始缺損圖，最后1幅圖為修復結果，其大小為64×64×64體素，缺失大小為20×64×64體素，設置的模板大小為7×7×7體素。由圖4的修復結果可以看出，三維圖像的缺失部分得到了很好地修復，保證了三維圖像空間結構的連續性，視覺上幾乎看不出原有的缺損，從而證明了算法對于三維圖像結構修復的準確性。圖4中的系列圖是按照一定的時間間隔選取的，清晰地展示了算法的修復過程，可以看出等照度面強度對于填充順序的影響，算法優先選擇圖像的結構進行填充，很好地傳播了等照度面，有效地修復了圖像的結構。同樣，由于置信度的影響，在等照度面的影響相同的情況下，算法將會優先填充可靠信息較多的邊界模塊。

圖5 巖心三維圖像修復實驗

2.2 巖心三維圖像

我們對多組存在較大范圍缺失的巖心三維圖像進行了相應的修復實驗，圖5展示了其中1組比較典型的巖心三維圖像的修復。完整的原始巖心三維圖像的三維顯示和各個方向的切片如圖5a所示，其大小為128×128×128體素，該巖心圖像噪聲量較大，結構不清晰，且有一定的紋理信息，是具有一定代表性的巖心三維圖像。目標部分（缺失）大小為192 704體素，其結構不規則，缺失后的圖像如圖5b，其中的zy切片為全白色，表示信息全部缺失的情況。圖5c為本文算法得到的修復結果，使用的模板大小為9×9×9體素。圖5d為采用基于PDE的三維圖像修復算法所得到的修復結果，其中的迭代步長dt=0.001，迭代次數為104。由圖5d可以看出，基于PDE的三維圖像修復算法對大范圍缺失的巖心三維圖像的修復結果會存在明顯的模糊現象，且不能修復圖像的紋理細節，嚴重影響到巖心的紋理分析和滲流特性分析，直觀地說明了該算法不適用于大范圍缺失的巖心三維圖像修復。本文的修復結果很好地填充了缺失部分，保證了三維圖像空間結構的連續性，保持了原始三維圖像中已知部分的結構和細節特征。

表1給出了使用本文算法進行的8組修復實驗的各項統計數據，其中的有效孔隙度［14］在一定程度上表示了孔隙的連通性。兩點概率函數S（r）常用于描述巖心孔隙結構的統計特征，表示三維圖像中任意相隔距離為r體素的兩點都為同一相的概率。兩點能量［2-3］相當于修復結果和原始三維圖像的兩點概率函數的差的平方和，在一定程度上描述了修復結果與原始三維圖像的微觀統計特性的差異。圖6展示了修復結果與原始三維圖像在x、y、z各個方向上的兩點概率函數S（r）。由表1中的兩點能量和圖6的兩點概率函數比較可以看出，修復結果與原始巖心圖像的兩點概率十分相近，修復結果保持了原始完整巖心圖像的統計特征，說明此三維圖像修復算法對于巖心三維圖像的修復是準確有效的。

表1 8組實驗的各項統計數據

圖6 表1中第3組實驗兩點概率函數

3 結 論

本文提出的三維圖像修復算法能夠修復較大范圍缺失的巖心三維圖像，修復結果保持了三維圖像的空間連續性以及原始圖像的結構和細節特征，修復結果與真實完整的三維圖像的形態結構和統計特性都極為相似。

［1］ XU Zhi，TENG Qizhi，HE Xiaohai，et al.A reconstruction method for three-dimensional pore space using multiple point geology statistic based on statistical pattern recognition and microstructure characterization ［J］.International Journal for Numerical and Analytical Methods in Geomechanics，2013，37（1）：97-110.

［2］ TANG Tang，TENG Qizhi，HE Xiaohai，et al.A pixel selection rule based on the number of different-phase neighbours for the simulated annealing reconstruction of sandstone microstructure［J］.Journal of Microscopy，2009，234（3）：262-268.

［3］ YEONG C L，TORQUATO S.Reconstructing random media：II Three-dimensional media from twodimensional cuts［J］.Physical Review：E，1998，58（1）：224.

［4］ 史金鋼，齊春.一種雙約束稀疏模型圖像修復算法［J］.西安交通大學學報，2012，46（2）：6-10.

SHI Jingang，QI Chun.An image inpainting algorithm based on sparse modeling with double constraints［J］.Journal of Xi’an Jiaotong University，2012，46（2）：6-10.

［5］ MORITA S.Three dimensional image inpainting ［M］.Berlin，Germany：Springer，2006：752-761.

［6］ CRIMINISI A，PEREZ P，TOYAMA K.Object removal by exemplar-based inpainting ［C］∥Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition.Piscataway，NJ，USA：IEEE，2003：721-728

［7］ CRIMINISI A，PEREZ P，TOYAMA K.Region filling and object removal by exemplar-based image inpainting［J］.IEEE Transactions on Image Processing，2004，13（9）：1200-1212.

［8］ TANG Feng，YING Yiting，WANG Jin，et al.A novel texture synthesis based algorithm for object removal in photographs［M］.Berlin，Germany：Springer，2005：248-258.

［9］ PATWARDHAN K A，SAPIRO G，BERTALMIO M.Video inpainting under constrained camera motion［J］.IEEE Transactions on Image Processing，2007，16（2）：545-553.

［10］谷伊，韓軍.基于樣本的圖像修補方法在視頻修復中的應用 ［J］.應用科學學報，2010，28（2）：163-169.

GU Yi，HAN Jun.Video restoration using exemplarbased inpainting ［J］.Journal of Applied Sciences，2010，28（2）：163-169.

［11］TENG Qizhi，YANG Dan，XU Zhi，et al.Training image analysis for three-dimensional reconstruction of porous media ［J］.Journal of Southeast University，2012，28（4）：415-421.

［12］VAN DEN ELSEN P A，MAINTZ J B A，POL E J D，et al.Automatic registration of CT and MR brain images using correlation of geometrical features ［J］.IEEE Transactions on Medical Imaging，1995，14（2）：384-396.

［13］JARVIS R A.A perspective on range finding techniques for computer vision［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1983（2）：122-139.

［14］楊勝來，魏俊之.油層物理學 ［M］.北京：石油工業出版社，2004.