摘 要： 在巖心三維模型中，中軸是描述孔隙結(jié)構(gòu)特征的一種重要表示方法。針對現(xiàn)有的拓?fù)浼?xì)化和距離變換方法存在的中軸提取不準(zhǔn)確和不連續(xù)的問題，提出了一種使用快速行進(jìn)水平集方法進(jìn)行距離變換的三維中軸提取算法。對比實(shí)驗(yàn)結(jié)果表明，該算法中提取的中軸在準(zhǔn)確性和連續(xù)性上有著較好的保證，實(shí)際應(yīng)用中效果良好。

0 引言

三維模型分析技術(shù)在分析巖心孔隙結(jié)構(gòu)和其相應(yīng)的統(tǒng)計特性中有著重要的應(yīng)用，中軸骨架是孔隙三維模型分析中的一項(xiàng)重要參數(shù)。中軸骨架表現(xiàn)出了三維模型的拓?fù)浣Y(jié)構(gòu)，使用中軸描述三維模型不僅能很好地表示三維模型的結(jié)構(gòu)信息，而且還能提高內(nèi)存的使用率和數(shù)據(jù)壓縮率。在巖心三維模型分析中，將中軸應(yīng)用于模型的迂曲度計算，有著重要的實(shí)際意義和工程應(yīng)用價值。

國內(nèi)外已有許多學(xué)者對三維中軸提取進(jìn)行了研究，目前的三維圖像骨架化方法主要有拓?fù)浼?xì)化方法、基于距離變換的方法、基于V圖的方法等。

拓?fù)浼?xì)化算法［1-2］能較好地保證拓?fù)浣Y(jié)構(gòu)，但對邊界噪聲較為敏感，不能保證中軸位置的準(zhǔn)確性。距離變換算法［3-5］能夠提取出不規(guī)則物體的骨架并且具有平移、旋轉(zhuǎn)、縮放不變性，但是很難保證模型的連通性。V圖方法［6-7］是利用V圖思想以中軸面或中軸線為基礎(chǔ)得到模型的骨架，一般用于生成多尺度的骨架。此外還有將距離變換和迭代并行細(xì)化相結(jié)合的方法［8］，得出的骨架具有良好的連通性和拓?fù)涞葍r性，但對于三維體數(shù)據(jù)需要多次迭代，計算量較大。

本文針對巖心三維孔隙模型進(jìn)行迂曲度計算等需求，提出了一種使用水平集進(jìn)行距離變換的方法求解三維結(jié)構(gòu)的骨架，使用水平集的距離變換算法與一般的距離變換方法相比具有更好的穩(wěn)定性和拓?fù)錈o關(guān)性，解決了距離變換存在的連通性問題。該方法提取出的中軸骨架具有很好的連通性和拓?fù)涞葍r性。

1 水平集方法

水平集方法是1988年由OSHER S和SETHIAN J首次提出的［9］，用于解決遵循熱力學(xué)方程下火苗外形的變化過程。其基本思想是將曲面

（用x表示）看成高一維空間中某一函數(shù)

（稱為水平集函數(shù)）的零水平集，同時曲面的演化也擴(kuò)充到高一維的空間中。曲線的演化轉(zhuǎn)化成一個純粹的求偏微分方程數(shù)值解問題。在任何時間，波前的位置可以由水平集函數(shù)

的零水平集來確定。將水平集函數(shù)按照它所滿足的發(fā)展方程進(jìn)行演化或迭代，由于水平集函數(shù)不斷進(jìn)行演化，所以對應(yīng)的零水平集也在不斷變化，當(dāng)水平集演化趨于平穩(wěn)時，演化停止，得到界面形狀。

假設(shè)

以速度F沿x（t）法線的方向進(jìn)行移動，若波前的粒子在路徑上，則粒子的水平集值必定為0：

其全微分方程為：

由于曲線沿法線方向移動，F(xiàn)可以表示為：

結(jié)合以上各式可以得出：

這就是參考文獻(xiàn)［9］給出的水平集公式。

2 快速行進(jìn)方法

考慮波前以速度F沿一個方向運(yùn)動，根據(jù)速度函數(shù)的符號，波前單調(diào)遞增或遞減傳播，計算每個點(diǎn)的到達(dá)時間。運(yùn)動公式可以表示為：

其中，F(xiàn)是速度函數(shù)，T是到達(dá)時間。

在三維空間中，求解上述等式需要計算出每個點(diǎn)（x，y，z）的到達(dá)時間。三維空間中梯度的離散化形式為：

為了求解這個等式，SETHIAN J［10］提出了一種簡單但是低精度的公式：

式中前向差分和后向差分的表現(xiàn)形式如下：

快速行進(jìn)算法的更新過程如下：

（1）從初始點(diǎn)開始推進(jìn)算法。將初始點(diǎn)加入到已知點(diǎn)集合，標(biāo)記為Known。

（2）找出Known點(diǎn)的鄰域，對鄰域點(diǎn)計算到達(dá)時間，將求解點(diǎn)標(biāo)記為trial。

（3）選取trial中到達(dá)時間最小的點(diǎn)，標(biāo)記為Known。

（4）查找所有Known點(diǎn)鄰域內(nèi)所有的未標(biāo)記點(diǎn)，對其計算到達(dá)時間，將求解點(diǎn)標(biāo)記為trial。

（5）轉(zhuǎn)至步驟（3），直到遍歷所有點(diǎn)，退出循環(huán)，算法結(jié)束。

3 使用水平集方法進(jìn)行距離變換求解骨架

在本文算法中，首先提取出三維模型中所有的連通目標(biāo)，獨(dú)立處理每個連通目標(biāo)。對每個連通目標(biāo)進(jìn)行距離變換，求解距離場D（x，y，z）。距離場是目標(biāo)中所有體素點(diǎn)到邊界距離的最小值，可以通過計算歐氏距離得到。選取距離場中具有最大距離的體素點(diǎn)作為全局最大距離點(diǎn)。根據(jù)距離場計算速度函數(shù)：

選取全局最大距離點(diǎn)作為目標(biāo)的中心點(diǎn)，將目標(biāo)的邊界體素點(diǎn)作為水平集曲線，使用快速行進(jìn)方法進(jìn)行演化，演化速度由速度函數(shù)決定。在每個體素更新的過程中，快速行進(jìn)算法計算出了當(dāng)前體素點(diǎn)到起點(diǎn)的到達(dá)時間。

對于各向同性的快速行進(jìn)方法，沿著波前的垂直方向波傳播得最快。由于每個中軸點(diǎn)梯度的方向始終垂直于波前，因此中軸體素點(diǎn)具有最快的傳播速度，可以根據(jù)這一特性來提取出中軸體素點(diǎn)，三維模型的中軸即可通過選定的初始體素沿著到達(dá)時間的梯度方向回溯求得。

三維模型單個孔隙目標(biāo)中軸提取流程如下：

（1）求解三維模型的距離場，根據(jù)距離場求解出速度函數(shù)和全局最大距離點(diǎn)。

（2）提取三維模型中所有獨(dú)立的連通目標(biāo)，對每個目標(biāo)，將全局最大距離點(diǎn)加入點(diǎn)源集合，利用速度圖像，用快速行進(jìn)方法求解水平集，得到到達(dá)時間，將距離點(diǎn)源集合到達(dá)時間最大的點(diǎn)作為回溯起始點(diǎn)。

（3）從回溯起始點(diǎn)開始，沿著到達(dá)時間的梯度方向開始回溯，一直回溯到點(diǎn)源集合為止?；厮萃瓿珊蟮玫揭粭l由最遠(yuǎn)點(diǎn)到點(diǎn)源集合的一條分支。

（4）判斷分支長度。若分支長度大于距離場中的最大值，則認(rèn)為這條分支是中軸的一部分，把分支加入中軸中，并且將該分支中的所有點(diǎn)加入點(diǎn)源集合；若分支長度小于距離場中的最大值，則剔除掉該分支，并且該目標(biāo)的中軸提取完成。

4 實(shí)驗(yàn)結(jié)果分析

本文以下的測試使用圖1所示的巖心三維模型，針對拓?fù)浼?xì)化方法和本文提出的基于快速行進(jìn)水平集方法兩種不同的中軸提取算法進(jìn)行了對比。兩種方法提取中的中軸結(jié)果如圖2所示。

從圖2可以看出，拓?fù)浼?xì)化方法保持了中軸的連通性，但是逐步剔除邊界點(diǎn)時，為了保證單像素性，產(chǎn)生了很多折線。使用這種中軸計算三維孔隙圖像的迂曲度時會產(chǎn)生較大的誤差。從圖3可以看出，本文使用的算法提取出的中軸較為平滑，在保持連通性的同時，很好地解決了折線的問題。圖4為本文提取出的中軸和三維模型的疊加顯示，可以看出本文算法提取出的中軸在連通性和正確性上得到了保證。

5 結(jié)論

本文使用快速行進(jìn)的水平集方法實(shí)現(xiàn)了巖心三維孔隙模型的中軸提取，通過對比實(shí)驗(yàn)，可以看出本文提出的算法有著明顯的優(yōu)勢，提取出的骨架有著很好的準(zhǔn)確性和連通性，有實(shí)際意義和應(yīng)用價值。

參考文獻(xiàn)

［1］ 王廣壘，張維忠，宋明玉，等?；跀?shù)學(xué)形態(tài)學(xué)的鞋楦特征曲線骨架的提取方法［J］。青島大學(xué)學(xué)報（自然科學(xué)版），2012，25（2）：43-46.

［2］ 徐瑩?；跀?shù)學(xué)形態(tài)學(xué)的圖像骨架提取和復(fù)原的改進(jìn)算法［J］。成都信息工程學(xué)院學(xué)報，2009，24（3）：259-263.

［3］ 張國棟，韓佳池?；谀：嚯x變換的骨架剪枝算法［J］。沈陽航空航天大學(xué)學(xué)報，2012，29（1）：64-69.

［4］ GAGVANI N， KENCHAMMANA H D， SILVER D. Volume animation using the skeleton tree［C］。 Proceedings of IEEE Volume Visualization， 1998：47-53.

［5］ DEY T K， SUN J. Defining and computing curve-skeletons with medial geodesic function［C］。 Proceedings of the fourth Eurographics Symposium on Geometry processing， AirelaVille， Switzerland， Eurographics Association， 2006： 143-152.

［6］ 劉輝，秦茂玲，徐海峰。基于Reeb圖的三維網(wǎng)格模型骨架提取算法［J］。信息技術(shù)與信息化，2012，24（5）：1672-9528.

［7］ 吳艷花。三維模型骨架提取算法及其在檢索中的應(yīng)用［D］。廣州：中山大學(xué)，2013.

［8］ 滕奇志，康瑕，唐棠，等。基于升序復(fù)核的并行三維圖像骨架化算法［J］。光學(xué)精密工程，2009，17（10）：2528-2534.

［9］ OSHER S， SETHIAN J A. Fronts propagating with curvature-dependent speed： algorithms based on Hamilton-Jacobi formulations［J］。 Journal of Computational Physics， 1988，79：12-49.

［10］ SETHIAN J A. Level sets methods and fast marching methods（2nd edition）［M］.Cambridge University Press，1999.

編輯：jq