黃賢勝，成衛(wèi)青，林 偉，熊生春

(1.南京郵電大學 計算機學院，江蘇 南京 210023;2.中國科學院 滲流流體力學研究所，河北 廊坊 065007)

0 引 言

隨著CT掃描技術和各專業(yè)學科的發(fā)展，在物質(zhì)探測方面所具有的巨大優(yōu)勢使得CT技術在非醫(yī)學領域(如工業(yè)、地球物理、工程、農(nóng)業(yè)、安全監(jiān)測、油氣勘探等)得到了廣泛的應用[1]。文中利用CT掃描技術對砂巖石巖心進行切片掃描，得到二維CT序列圖像。

計算機進行建模的目的是通過一系列二維圖像構造出一個簡單直觀的模型，這個模型能反映真實世界中具體存在的物體，是真實對象的“數(shù)學數(shù)據(jù)抽象”。在石油地質(zhì)研究領域，對巖石的切片進行分析解讀是很有必要的。通過計算機技術可以對巖石薄片圖像進行分析處理，快速準確地獲取巖石樣品參數(shù)：如孔隙度、孔隙形狀、孔喉比和吼道配位數(shù)等特征信息。基于上述二維圖像中采集的特征信息，可以將其拓展到三維空間[2]。利用計算機進行三維巖心建模與可視化技術已經(jīng)存在并且發(fā)展了好幾十年，目前的三維建模和可視化主要有兩個作用：一個是利用解析出的數(shù)據(jù)構建出三維模型，另一個是對地底石油層是否可勘探提供理論基礎。

巖石(多孔介質(zhì))三維重構和可視化的重要意義主要表現(xiàn)在以下幾個方面：(1)傳統(tǒng)的以實驗為主要研究方法的方式已經(jīng)滿足不了某些問題的基本需求，如低孔滲儲層的巖心滲透問題、復雜儲層的巖石物理屬性問題[3]。(2)重構出的孔隙網(wǎng)絡結構可以用于更好地分析整個巖石的拓撲結構和幾何特征。多孔介質(zhì)內(nèi)部結構復雜多樣，主要由孔隙空間和骨架組成，其中孔隙空間又分為大孔隙和喉道，大孔隙代表較大的孔隙空間，喉道代表相對狹長的孔隙空間。

1 研究現(xiàn)狀

1.1 主要方法

目前，國內(nèi)外實現(xiàn)巖心建模和可視化的技術主要包括以下幾種：

(1)直接在已有的三維建模軟件上進行操作。比較常用的軟件如Blender、UG、Seamless 3d、3d Canvas等；也可以在軟件上對已經(jīng)建好的模型進行再加工再改造，滿足相關的需求。這種方法比較簡單易上手，不需要考慮復雜的結構、對應的節(jié)點坐標和采用的算法等，直接將數(shù)據(jù)輸入進去，就能得到想要的重構模型[4]。

(2)在現(xiàn)有的系統(tǒng)中，運用相關的編程語言加上可視化工具也可以實現(xiàn)三維建模的操作。比如OpenGL，與平臺無關，可以在任何的軟件系統(tǒng)中實現(xiàn)三維建?？梢暬痆5]，但是這類方法可操作性差，對專業(yè)知識要求高，工作量比較大，所以并沒有得到廣泛的應用。

(3)利用VTK進行三維重構。在VTK中可以自由獲取它的源碼，世界上大量的研究人員都是使用VTK庫進行3D圖像建模和可視化操作。

文中采用VTK和C++語言相結合的方法，直接調(diào)用VTK中封裝的MarchingCubes類—vtkMarchingCubes類進行三維建模，實現(xiàn)巖心的三維可視化功能。通過對孔隙和骨架進行分割，分析了孔喉結構，立體地展示巖心內(nèi)部的空間結構。

1.2 VTK簡介

VTK(visualization toolkit)是一個開源的免費軟件系統(tǒng)，它包含了一個C++類庫和眾多的接口層，可運用在建筑學、氣象學、醫(yī)學、生物學或者航空航天學上，通過對體、面、光源等的逼真渲染，幫助人們理解那些采取錯綜復雜而又往往規(guī)模龐大的數(shù)字呈現(xiàn)形式的科學概念或結果。它具有強大的三維圖形功能，在極大地改善可視化效果的同時又可以充分利用現(xiàn)有的圖形庫和圖形硬件[6]。

1.3 STL文件簡介

STL文件類型是3D SYSTEM公司在1988年編訂的接口協(xié)議，是專門為快速原型制造技術服務的三維圖形文件格式。STL文件有2種類型：ASCII格式和二進制格式。STL文件由多個三角形面片的定義組成，每個三角形面片的定義包括三角形各個頂點的三維坐標及三角形面片的法矢量[7]，這些三角面片的法向量均要求指向?qū)嶓w外部的方向，同時模型表面必須布滿三角面片，不允許有遺漏(裂縫和空洞)[8]。

ASCII格式的STL文件的基本結構[9]如下：

solid filename //STL文件路徑及文件名

facet normal x y z //三角面片法向量的3個分量值

outer loop

vertex x y z //三角面片第一個頂點坐標

vertex x y z //第二個頂點坐標

vertex x y z

End loop

End facet //第一個三角面片定義結束

……

End solid filename //整個STL文件定義結束

2 MarchingCubes簡介

2.1 三維可視化的理解

近年來，隨著計算機運算能力的不斷加強和設備的不斷改進，研究人員提出了基于體素的體繪制實現(xiàn)[10]。但是面繪制仍然是進行三維可視化的主流操作，面繪制是通過先建立網(wǎng)絡模型之后再進行渲染的方式完成三維模型構建。

2.2 MarchingCubes基本原理

MarchingCubes (MC)算法是面繪制算法中的經(jīng)典算法，而且簡單容易實現(xiàn)，得到了廣泛的應用。

首先在MC算法中要明確一個叫做“體素”的概念。體素是在三維圖像中由相鄰的八個體素點(三維數(shù)據(jù)場中相鄰兩層中相鄰的八個頂點)組成的立方體，八個體素點構成一個體素，對體素中的每一個體素點進行編號[11]，圖1表示其基本結構。

圖1 立方體體素基本結構

MC算法的主要思想是通過在體素中構建出各種不同的三角面片。大量的三角面片構成了一個等值面(isosurface extraction)，這個等值面就是重建模型的表面。等值面上的點定義如下：

{(x,y,z)|f(x,y,z)=c}

其中，(x,y,z)為空間點的坐標，f(x,y,z)為對應點數(shù)據(jù)的值，常數(shù)c為三維重構過程中給定的閾值。

2.3 MC算法的步驟

(1)確定體素中等值面的剖分方式。

判斷體素中的頂點狀態(tài)，某個頂點的值是大于或者等于給定閾值時，該頂點就被標記為“1”，反之該點的狀態(tài)就標記為“0”。在整個數(shù)據(jù)場中，在每個體素內(nèi)重構一等值面片，將得到的所有近似等值面片拼合起來，故稱之為“移動立方體”法[12]。

因為每個頂點每次有“0”或“1”兩種情況，所以總共就有28(256)種三角剖分狀態(tài)。根據(jù)對稱性和旋轉對稱性，可以把256種最終簡化細分成15種。

為了確定每個體素是對應于哪種構型，在程序中加入了一個字節(jié)長參數(shù)——構型索引index指數(shù)。該index指數(shù)每字位分別表示體素八個頂點的狀態(tài)值，如圖2所示。

圖2 索引index字節(jié)

對于圖1這一基本構型，用序列來表示其三角剖分形式：{3,11,2，-1，-1，-1，-1，-1，-1，-1，-1，-1，-1，-1，-1，-1}，此序列表示的是第3條邊、第11條邊和第2條邊上有三角面片的頂點，-1表示缺省。

(2)計算等值面頂點處的坐標。

當體素內(nèi)三角面片剖分方式確定后，由于MC算法假設棱邊的像素值是呈線性變化的，所以在體素上兩個端點的狀態(tài)分別為“0”和“1”的邊上，必有一點等于閾值(等值面的值)，此點也就是三角面片的頂點。計算該邊上等值面的頂點和該頂點的法向量。對于體素某一棱邊，設其兩個端點(P1、P2)的標量數(shù)據(jù)為V1、V2，那么插值得到的交點(即三角面片的頂點)為P，P的求解表達式為：

P=P1+(isovalue-V1)(P2-P1)/(V2-V1)

其中，P為頂點坐標，P1、P2為兩個端點的坐標，V1、V2為兩個端點處的像素值，isovalue為等值面的值。

(3)計算等值面頂點處的法向量。

要把等值面繪制出來，還需要計算三角面片的法向量。利用中心差分法計算出體素中各頂點處的梯度，再利用線性插值法求出三角面片各頂點的梯度(法向量)，實現(xiàn)等值面的繪制[12]。

設立方體頂點處的法向量為N1、N2(標量數(shù)據(jù)仍然為V1、V2)，利用線性插值得到的三角面片頂點處法向量為：

N'=N1+(isovalue-V1)(N2-N1)/(V2-V1)

(4)引入鄰接表。

結合相關文獻可知，在實際的三維建模中，有一部分體素并不需要去遍歷，是可以直接忽略的。文中采用對鄰面延伸的方法，這樣建模的效率又得到了提高。確定了一個立方體的所有三角面片后，可根據(jù)三角面與立方體各面相交的情況，僅遍歷相關的鄰接立方體[13]。為此，可以先構建一個鄰接表分別存儲256種立方體的三角面片構建應延伸的方向(前、后、上、下、左、右)。

3 三維建模的改進

3.1 巖心二維圖像預處理

3.1.1 去虛邊

本次實驗的數(shù)據(jù)集是對致密砂巖巖心進行CT掃描得到的二維圖像。由于巖石掃描圖像存在亮度較亮或較暗、對比度不明顯、圖片畫面模糊等缺陷，嚴重影響圖像建模的準確性。因此，在不損壞圖像有用信息的前提下需要對其進行適當?shù)娜ピ肼?、去虛邊等操作?/p>

(1)獲取閾值。

高斯濾波器是一種線性濾波器，能夠有效地抑制噪聲和平滑圖像。首先，將原圖像進行顏色空間轉換，轉換成灰度圖，然后利用高斯濾波器函數(shù)進行濾波操作，將灰度圖與指定的高斯內(nèi)核進行卷積運算。最后創(chuàng)建滾動條(createTrackbar)函數(shù)，以及對圖像取不同閾值進行二值化處理的函數(shù)，通過手動控制滾動條的方式，調(diào)整二值化的閾值，直觀觀察出合適的閾值大小。在CT高分辨率切片原圖中選取240.tif、282.tif和330.tif，經(jīng)過不斷調(diào)整后，確定這三張圖像的閾值分別為174、170和167，生成的二值化圖像分別如圖3(a)、(b)和(c)所示。

(2)形態(tài)學處理。

運用形態(tài)學中的膨脹和腐蝕操作。膨脹(dilate)就是求局部最大值的操作，使圖像部分區(qū)域亮度增加。與膨脹相反，腐蝕(erode)就是求局部最小值的操作。腐蝕可以簡單理解為消除物體所有邊界點的過程，使整體更暗。特別是在黑白交界的地方，膨脹和腐蝕的現(xiàn)象會非常直觀明顯。

為去除原圖中的虛邊(發(fā)白模糊區(qū)域)，膨脹會使白色區(qū)域膨脹，腐蝕會減少白色區(qū)域。先腐蝕，讓黑色區(qū)域中的一些白色像素點消失，同時白色區(qū)域塊也會暗一些，再通過膨脹補回來是進行開運算的過程；先膨脹會將一些散落的白色小塊通過擴增連接到一起，再腐蝕削減一點是進行閉運算的過程。靈活運用開運算和閉運算，膨脹和腐蝕不同的次數(shù)，會產(chǎn)生不同的效果。對之前的三張CT原圖進行形態(tài)學操作，使用的閾值范圍是合適閾值X到255，并進行了摳圖和找輪廓操作。實際的操作分別是：

①切片240.tif經(jīng)反復試驗后依次使用了腐蝕一次和膨脹三次，使用的閾值范圍是174到255，生成圖3(d)；

②切片282.tif經(jīng)反復試驗后依次使用了腐蝕一次和膨脹四次，使用的閾值范圍是170到255，生成圖3(e)；

③切片330.tif經(jīng)反復試驗后依次使用了膨脹兩次和腐蝕一次，使用的閾值范圍是167到255，生成圖3(f)。

圖3 去虛邊和形態(tài)學預處理結果

3.1.2 分層灰度化和裁剪

分層灰度化是一個分割過程，分割可以將感興趣部位提取并顯示出來，便于對模型重建進行判斷[8]。將背景的灰度值從原來的0(全黑)變成128，再對巖心內(nèi)部進行二值化操作。明顯地觀察出黑色部分是孔隙，不規(guī)則地分布在白色的骨架中。這時的二維切片圖片就轉換成了灰度128的背景、灰度255的骨架和灰度0的孔隙。

由于文中的CT掃描圖像是2 048*2 048的tif格式的圖片，一張圖片所占內(nèi)存就高達4 M，所以要想完成大量tif格式圖片的建模就需要進行裁剪。利用Matlab裁剪工具裁去的是背景部分，保留巖心部分。圖4顯示已去虛邊的240.tif圖像及其灰度化與裁剪結果。

圖4 切片240.tif去虛邊、分層灰度化與裁剪結果

3.2 對移動立方體算法(MC算法)的改進

(1)三角面片頂點計算優(yōu)化。

根據(jù)前面介紹的傳統(tǒng)MC算法，如果使用線性插值法直接計算三角面片的頂點(等值點)和對應的法向量，那么過程會很復雜，運算量巨大，同時也會消耗大量的運算時間，而且如果體素量過大，數(shù)據(jù)場數(shù)據(jù)量更高，則生成的三角面片和真實的三角面片就有較大的出入。為了降低算法運行時間，有人提出使用黃金分割點作為該邊的等值點[14]，由于每個體素都很小，確定每個體素每條邊的黃金分割點仍需要經(jīng)過一些計算，效率上值得再次商榷。文中采取的方式是將體素邊分為四等份，提取四等分點，該邊的四分之一處的點即為三角面片的頂點，計算上簡單方便，極大減少了計算量，同時產(chǎn)生的三角面片比較準確。

(2)頂點法向量計算優(yōu)化。

法向量也能進一步的簡化，法向量計算改為對體素棱邊上的頂點的法向量直接求平均，具有很好的效果。

經(jīng)過一系列的改進后，偽代碼如下：

輸入：預處理后的二維圖像

輸出：三維模型

Begin：

1.初始化二元標志數(shù)組List={Flag，p}，F(xiàn)lag表示立方體p是否被訪問過，初始化CubeQueue隊列(利用存放待處理的立方體)，初始化存儲256種立方體的鄰接關系的鄰接表neiborCases；

2.然后檢測數(shù)據(jù)集，選取一個構型索引index不為0或255的立方體作為p，并將其壓入隊列CubeQueue中，并置List[p].Flag=0；

while(!CubeQueue.empty()){

3.從CubeQueue中取出隊首立方體C； //編號C的立方體

4.if(List[C].Flag==1) continue； //表示該編號的立方體已經(jīng)處理過

5.令List[C].Flag=1；

6.根據(jù)當前立方體體素的頂點狀態(tài)得到其構型索引index；

7.利用構型索引index值查找三角剖分表得到當前體素的剖分形式，結合文中所述的計算方法求出面片的頂點坐標和頂點處的法向量，并將其作為輸出的數(shù)據(jù)；

8.通過index查找neiborCases鄰接表，將鄰接的立方體壓入到CubeQueue隊列中去；

9.用VTK工具繪制三維模型}

End

4 實驗結果與分析

文中實驗使用的主機有8G RAM，操作系統(tǒng)為64位 Windows 7。開發(fā)使用Cmake，Visual Studio 2013嵌入VTK類庫、Python-opencv圖像處理函數(shù)庫和Pycharm 2018編輯器。使用VC++、Python和Matlab語言編程實現(xiàn)所有功能。

(1)模型對比。

對100張預處理后的二維圖像進行三維建模，結果如圖5所示，前兩張圖是經(jīng)典MC算法生成的模型，后兩張是改進MC算法生成的模型。

圖5 改進前后算法重構模型

(2)3D實物打印。

在3D打印機上讀取文中生成的STL文件，采用9 400樹脂材料進行打印，精確的尺寸為89.253 mm*81.139 mm*49.750 mm。圖6是打印出的實物模型。

圖6 致密砂巖巖心放大3D打印實物

(3)實驗結果總結。

由于建模的圖片數(shù)量有限，而且新舊算法都是對巖心進行建模，都有較高的精度，所以改進前后的整體模型相異性較小。改進后的巖心內(nèi)部的孔隙都能較好地重構出來，能夠保證建模的準確性。表1是多次實驗結果的匯總。

表1 改進的MC算法和原始的MC算法的對比

由表1可見，改進后的MC算法無論從建模的運算時間上還是生成STL文件的大小上，都有不同程度的優(yōu)化，提高了效率。在保證生成模型準確性的基礎上，縮短了運算時間，減小了最終文件大小，有一定的優(yōu)勢，增強了MC算法繪制的實用性。

5 結束語

針對巖心的CT圖像設計了有效的圖像預處理方法，對傳統(tǒng)建模方法MC進行了改進，并利用VTK編程實現(xiàn)了一個對二維CT圖像序列進行三維建模的系統(tǒng)。實驗結果表明改進MC算法在保證生成模型準確性的基礎上提高了三維建模的效率。