摘  要： 在巖心三維模型中，中軸是描述孔隙結(jié)構(gòu)特征的一種重要表示方法。針對(duì)現(xiàn)有的拓?fù)浼?xì)化和距離變換方法存在的中軸提取不準(zhǔn)確和不連續(xù)的問(wèn)題，提出了一種使用快速行進(jìn)水平集方法進(jìn)行距離變換的三維中軸提取算法。對(duì)比實(shí)驗(yàn)結(jié)果表明，該算法中提取的中軸在準(zhǔn)確性和連續(xù)性上有著較好的保證，實(shí)際應(yīng)用中效果良好。
　　關(guān)鍵詞： 中軸提??；水平集；快速行進(jìn)方法；距離場(chǎng)
0 引言
　　三維模型分析技術(shù)在分析巖心孔隙結(jié)構(gòu)和其相應(yīng)的統(tǒng)計(jì)特性中有著重要的應(yīng)用，中軸骨架是孔隙三維模型分析中的一項(xiàng)重要參數(shù)。中軸骨架表現(xiàn)出了三維模型的拓?fù)浣Y(jié)構(gòu)，使用中軸描述三維模型不僅能很好地表示三維模型的結(jié)構(gòu)信息，而且還能提高內(nèi)存的使用率和數(shù)據(jù)壓縮率。在巖心三維模型分析中，將中軸應(yīng)用于模型的迂曲度計(jì)算，有著重要的實(shí)際意義和工程應(yīng)用價(jià)值。
　　國(guó)內(nèi)外已有許多學(xué)者對(duì)三維中軸提取進(jìn)行了研究，目前的三維圖像骨架化方法主要有拓?fù)浼?xì)化方法、基于距離變換的方法、基于V圖的方法等。
　　拓?fù)浼?xì)化算法[1-2]能較好地保證拓?fù)浣Y(jié)構(gòu)，但對(duì)邊界噪聲較為敏感，不能保證中軸位置的準(zhǔn)確性。距離變換算法[3-5]能夠提取出不規(guī)則物體的骨架并且具有平移、旋轉(zhuǎn)、縮放不變性，但是很難保證模型的連通性。V圖方法[6-7]是利用V圖思想以中軸面或中軸線為基礎(chǔ)得到模型的骨架，一般用于生成多尺度的骨架。此外還有將距離變換和迭代并行細(xì)化相結(jié)合的方法[8]，得出的骨架具有良好的連通性和拓?fù)涞葍r(jià)性，但對(duì)于三維體數(shù)據(jù)需要多次迭代，計(jì)算量較大。
　　本文針對(duì)巖心三維孔隙模型進(jìn)行迂曲度計(jì)算等需求，提出了一種使用水平集進(jìn)行距離變換的方法求解三維結(jié)構(gòu)的骨架，使用水平集的距離變換算法與一般的距離變換方法相比具有更好的穩(wěn)定性和拓?fù)錈o(wú)關(guān)性，解決了距離變換存在的連通性問(wèn)題。該方法提取出的中軸骨架具有很好的連通性和拓?fù)涞葍r(jià)性。
1 水平集方法
　　水平集方法是1988年由OSHER S和SETHIAN J首次提出的[9]，用于解決遵循熱力學(xué)方程下火苗外形的變化過(guò)程。其基本思想是將曲面（用x表示）看成高一維空間中某一函數(shù)（稱(chēng)為水平集函數(shù)）的零水平集，同時(shí)曲面的演化也擴(kuò)充到高一維的空間中。曲線的演化轉(zhuǎn)化成一個(gè)純粹的求偏微分方程數(shù)值解問(wèn)題。在任何時(shí)間，波前的位置可以由水平集函數(shù)的零水平集來(lái)確定。將水平集函數(shù)按照它所滿(mǎn)足的發(fā)展方程進(jìn)行演化或迭代，由于水平集函數(shù)不斷進(jìn)行演化，所以對(duì)應(yīng)的零水平集也在不斷變化，當(dāng)水平集演化趨于平穩(wěn)時(shí)，演化停止，得到界面形狀。
　　假設(shè)以速度F沿x（t）法線的方向進(jìn)行移動(dòng)，若波前的粒子在路徑上，則粒子的水平集值必定為0：

　　其全微分方程為：

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

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

　　這就是參考文獻(xiàn)[9]給出的水平集公式。
2 快速行進(jìn)方法
　　考慮波前以速度F沿一個(gè)方向運(yùn)動(dòng)，根據(jù)速度函數(shù)的符號(hào)，波前單調(diào)遞增或遞減傳播，計(jì)算每個(gè)點(diǎn)的到達(dá)時(shí)間。運(yùn)動(dòng)公式可以表示為：

　　其中，F(xiàn)是速度函數(shù)，T是到達(dá)時(shí)間。
　　在三維空間中，求解上述等式需要計(jì)算出每個(gè)點(diǎn)（x，y，z）的到達(dá)時(shí)間。三維空間中梯度的離散化形式為：

　　為了求解這個(gè)等式，SETHIAN J[10]提出了一種簡(jiǎn)單但是低精度的公式：

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

　　快速行進(jìn)算法的更新過(guò)程如下：
　?。?）從初始點(diǎn)開(kāi)始推進(jìn)算法。將初始點(diǎn)加入到已知點(diǎn)集合，標(biāo)記為Known。
　　（2）找出Known點(diǎn)的鄰域，對(duì)鄰域點(diǎn)計(jì)算到達(dá)時(shí)間，將求解點(diǎn)標(biāo)記為trial。
　?。?）選取trial中到達(dá)時(shí)間最小的點(diǎn)，標(biāo)記為Known。
　?。?）查找所有Known點(diǎn)鄰域內(nèi)所有的未標(biāo)記點(diǎn)，對(duì)其計(jì)算到達(dá)時(shí)間，將求解點(diǎn)標(biāo)記為trial。
　　（5）轉(zhuǎn)至步驟（3），直到遍歷所有點(diǎn)，退出循環(huán)，算法結(jié)束。
3 使用水平集方法進(jìn)行距離變換求解骨架
　　在本文算法中，首先提取出三維模型中所有的連通目標(biāo)，獨(dú)立處理每個(gè)連通目標(biāo)。對(duì)每個(gè)連通目標(biāo)進(jìn)行距離變換，求解距離場(chǎng)D（x，y，z）。距離場(chǎng)是目標(biāo)中所有體素點(diǎn)到邊界距離的最小值，可以通過(guò)計(jì)算歐氏距離得到。選取距離場(chǎng)中具有最大距離的體素點(diǎn)作為全局最大距離點(diǎn)。根據(jù)距離場(chǎng)計(jì)算速度函數(shù)：

　　選取全局最大距離點(diǎn)作為目標(biāo)的中心點(diǎn)，將目標(biāo)的邊界體素點(diǎn)作為水平集曲線，使用快速行進(jìn)方法進(jìn)行演化，演化速度由速度函數(shù)決定。在每個(gè)體素更新的過(guò)程中，快速行進(jìn)算法計(jì)算出了當(dāng)前體素點(diǎn)到起點(diǎn)的到達(dá)時(shí)間。
　　對(duì)于各向同性的快速行進(jìn)方法，沿著波前的垂直方向波傳播得最快。由于每個(gè)中軸點(diǎn)梯度的方向始終垂直于波前，因此中軸體素點(diǎn)具有最快的傳播速度，可以根據(jù)這一特性來(lái)提取出中軸體素點(diǎn)，三維模型的中軸即可通過(guò)選定的初始體素沿著到達(dá)時(shí)間的梯度方向回溯求得。
　　三維模型單個(gè)孔隙目標(biāo)中軸提取流程如下：
　?。?）求解三維模型的距離場(chǎng)，根據(jù)距離場(chǎng)求解出速度函數(shù)和全局最大距離點(diǎn)。
　　（2）提取三維模型中所有獨(dú)立的連通目標(biāo)，對(duì)每個(gè)目標(biāo)，將全局最大距離點(diǎn)加入點(diǎn)源集合，利用速度圖像，用快速行進(jìn)方法求解水平集，得到到達(dá)時(shí)間，將距離點(diǎn)源集合到達(dá)時(shí)間最大的點(diǎn)作為回溯起始點(diǎn)。
　　（3）從回溯起始點(diǎn)開(kāi)始，沿著到達(dá)時(shí)間的梯度方向開(kāi)始回溯，一直回溯到點(diǎn)源集合為止?；厮萃瓿珊蟮玫揭粭l由最遠(yuǎn)點(diǎn)到點(diǎn)源集合的一條分支。
　?。?）判斷分支長(zhǎng)度。若分支長(zhǎng)度大于距離場(chǎng)中的最大值，則認(rèn)為這條分支是中軸的一部分，把分支加入中軸中，并且將該分支中的所有點(diǎn)加入點(diǎn)源集合；若分支長(zhǎng)度小于距離場(chǎng)中的最大值，則剔除掉該分支，并且該目標(biāo)的中軸提取完成。
4 實(shí)驗(yàn)結(jié)果分析
　　本文以下的測(cè)試使用圖1所示的巖心三維模型，針對(duì)拓?fù)浼?xì)化方法和本文提出的基于快速行進(jìn)水平集方法兩種不同的中軸提取算法進(jìn)行了對(duì)比。兩種方法提取中的中軸結(jié)果如圖2所示。

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

5 結(jié)論
　　本文使用快速行進(jìn)的水平集方法實(shí)現(xiàn)了巖心三維孔隙模型的中軸提取，通過(guò)對(duì)比實(shí)驗(yàn)，可以看出本文提出的算法有著明顯的優(yōu)勢(shì)，提取出的骨架有著很好的準(zhǔn)確性和連通性，有實(shí)際意義和應(yīng)用價(jià)值。
　　參考文獻(xiàn)
　　[1] 王廣壘，張維忠，宋明玉，等.基于數(shù)學(xué)形態(tài)學(xué)的鞋楦特征曲線骨架的提取方法[J].青島大學(xué)學(xué)報(bào)（自然科學(xué)版），2012，25（2）：43-46.
　　[2] 徐瑩.基于數(shù)學(xué)形態(tài)學(xué)的圖像骨架提取和復(fù)原的改進(jìn)算法[J].成都信息工程學(xué)院學(xué)報(bào)，2009，24（3）：259-263.
　　[3] 張國(guó)棟，韓佳池.基于模糊距離變換的骨架剪枝算法[J].沈陽(yáng)航空航天大學(xué)學(xué)報(bào)，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.