基于改進(jìn)Marching Cubes算法的血流實(shí)時(shí)仿真研究

王　娜，陳國棟，陳　怡

(1.福建師范大學(xué)福清分校 電子與信息工程學(xué)院，福建 福清 350300；2.福州大學(xué) 物理與信息工程學(xué)院，福建 福州 350108)

?

基于改進(jìn)Marching Cubes算法的血流實(shí)時(shí)仿真研究

王娜1，陳國棟2*，陳怡2

(1.福建師范大學(xué)福清分校 電子與信息工程學(xué)院，福建 福清 350300；2.福州大學(xué) 物理與信息工程學(xué)院，福建 福州 350108)

針對傳統(tǒng)的Marching Cubes算法空體元檢測時(shí)間過多影響執(zhí)行效率的問題，設(shè)計(jì)了一種針對流體表面繪制的Marching Cubes改進(jìn)算法。在算法中，首先檢測了規(guī)則點(diǎn)陣的密度，然后通過設(shè)定閾值將粒子密度低于閾值的區(qū)域與密度高于閾值的區(qū)域分離，僅將密度較高的區(qū)域使用簡化版Marching Cubes算法繪制。仿真實(shí)驗(yàn)證明，與球形渲染算法和Marching Cubes算法相比，本文提出的算法減少了對空體元的訪問，提高了顯示的質(zhì)量，從而使得整體繪制算法符合實(shí)時(shí)渲染的要求。

Marching Cubes(MC)算法; 體元; 血流; 等值面; 閾值

血流仿真作為虛擬手術(shù)系統(tǒng)中的一個(gè)重要組成部分，隨著虛擬手術(shù)的進(jìn)行而實(shí)時(shí)產(chǎn)生，仿真的效果將直接決定虛擬手術(shù)系統(tǒng)的整體真實(shí)感。作為虛擬手術(shù)系統(tǒng)中的關(guān)鍵技術(shù)，血液流動(dòng)仿真也得到了快速的發(fā)展?？傮w而言，血液流動(dòng)仿真作為流體仿真領(lǐng)域中的一員，經(jīng)歷了早期的基于構(gòu)造的仿真和后期的基于物理的仿真。基于構(gòu)造的血流仿真計(jì)算量小，適合早期計(jì)算機(jī)硬件水平比較落后的情況，但其突出的問題在于真實(shí)感不強(qiáng)。因此，對血流仿真領(lǐng)域的研究已經(jīng)漸漸地轉(zhuǎn)向使用基于物理的方法。該方法相比基于構(gòu)造的仿真方法而言，真實(shí)感比較強(qiáng)。

近年來，國內(nèi)外眾多學(xué)者一直在努力對于基于物理的血流仿真領(lǐng)域的研究。文[1]結(jié)合立方插值擬質(zhì)點(diǎn)的方法來模擬血液流動(dòng)，該方法可以使用較為稀疏的網(wǎng)格對流體進(jìn)行精確的模擬，然而當(dāng)網(wǎng)格出現(xiàn)大的變形時(shí)會造成網(wǎng)格扭曲而導(dǎo)致精度降低。文[2]將彈簧模型與SPH(光滑粒子流體動(dòng)力學(xué))相結(jié)合實(shí)現(xiàn)了血液與血管的相互作用。文[3]利用流體動(dòng)力學(xué)工具結(jié)合3D技術(shù)對動(dòng)脈血管中的血流進(jìn)行了仿真。文[4]、[5]采用穩(wěn)定的半拉格朗日方法對血液流動(dòng)進(jìn)行模擬。文[6]提出了一種結(jié)合力反饋機(jī)制的滲血模擬仿真模型，該模型對傳統(tǒng)的拉格朗日法進(jìn)行了簡化。

綜上可知，血液流動(dòng)現(xiàn)象仿真的主要問題集中于實(shí)時(shí)性與真實(shí)感的矛盾上。因此，權(quán)衡真實(shí)感與實(shí)時(shí)性這對矛盾將是血流模擬中必須考慮的問題。

1　粒子位置信息的預(yù)處理

為了提高計(jì)算的實(shí)時(shí)性，需要將粒子的位置信息進(jìn)行預(yù)處理，具體做法為構(gòu)建一個(gè)可以包圍所有粒子的長方體，并將長方體均勻分割成大量的立方體，根據(jù)粒子的位置對其近鄰的立方體頂點(diǎn)進(jìn)行標(biāo)記，最終構(gòu)建出一個(gè)規(guī)則的三維點(diǎn)陣。

由于采用直接計(jì)算粒子與立方體頂點(diǎn)之間的距離尋找與粒子最近的立方體頂點(diǎn)的方法計(jì)算量大，本文通過粒子的三個(gè)維度的坐標(biāo)與立方體邊長a之間進(jìn)行比較，找到最近鄰的立方體頂點(diǎn)，具體算法如下：

由于長方體分割為立方體時(shí)是均勻分割，因此，可以計(jì)算出與xi最近的立方體頂點(diǎn)的x軸坐標(biāo)xa與xb，其中xa

(1)

(2)

(3)

(4)

圖1　尋找最近鄰的立方體頂點(diǎn)

2　Marching Cubes算法

在三維空間數(shù)據(jù)場中構(gòu)造等值面有很多種方法，其中最有代表性的就是Marching Cubes算法。Marching Cubes算法簡稱為MC算法，是W.E.Lorenson與H.E.Cline在1987年的時(shí)候提出的[7]。由于該算法具有容易實(shí)現(xiàn)、通用性強(qiáng)的優(yōu)點(diǎn)，因此得到較為廣泛的應(yīng)用。

2.1Marching Cubes算法的主要思想

Marching Cubes算法的主要思想是把離散的三維規(guī)則數(shù)據(jù)場構(gòu)建成一系列相鄰的立方體，采用線性插值計(jì)算出等值面與立方體棱邊的交點(diǎn)，將交點(diǎn)按一定方式連接成幾個(gè)相鄰的三角面片，作為等值面在局部體元內(nèi)的一個(gè)逼近表示。

Marching Cubes算法的流程如下：

(1) 將三維規(guī)則數(shù)據(jù)場讀入內(nèi)存中；

(2) 根據(jù)三維規(guī)則數(shù)據(jù)場，從前往后，從左往右構(gòu)建體元。每個(gè)體元的8個(gè)角點(diǎn)分成兩層，分別取自上下兩層的規(guī)則數(shù)據(jù)場；

(3) 根據(jù)等值面所采用的閾值，將體元的每個(gè)角點(diǎn)分成兩類，其一為當(dāng)體元角點(diǎn)的函數(shù)值大于閾值時(shí)，將該角點(diǎn)標(biāo)記為黑色，其二為當(dāng)體元角點(diǎn)的函數(shù)值小于閾值時(shí)，對角點(diǎn)不進(jìn)行標(biāo)記。對一個(gè)體元的所有角點(diǎn)進(jìn)行遍歷，即可得到該體元的狀態(tài)表；

(4) 根據(jù)狀態(tài)表，即可得到與等值面相交的體元邊界；

(5) 利用線性插值的方法計(jì)算出體元與等值面的交點(diǎn)位置；

(6) 利用中心差分算法求得體元各角點(diǎn)位置的法向，再通過線性插值計(jì)算三角形各個(gè)頂點(diǎn)位置的法向量；

(7) 遍歷所有層的數(shù)據(jù)；

(8) 根據(jù)三角面片的頂點(diǎn)坐標(biāo)和該位置上的法向量繪制出等值面圖象。

2.2Marching Cubes算法存在的缺陷

傳統(tǒng)的Marching Cubes算法需要通過在相鄰切片間構(gòu)建體元，每個(gè)體元上共有8個(gè)角點(diǎn)，因此每個(gè)體元的角點(diǎn)狀態(tài)共有28=256種，其相對應(yīng)的等值面分布也有256種。但是考慮到旋轉(zhuǎn)對稱和反轉(zhuǎn)對稱的因素，因此可以將等值面的種類縮小為15種[8]。如果對每個(gè)體元進(jìn)行處理，根據(jù)頂點(diǎn)位置與等值面閾值比對構(gòu)造出關(guān)系表，再由線性插值得到多個(gè)三角形，將得到的多個(gè)三角形網(wǎng)格形成等值面。一個(gè)醫(yī)學(xué)序列圖像所包含的三角形切片有可能數(shù)以百萬，進(jìn)行等值面繪制的操作對計(jì)算機(jī)硬件的要求比較高，由于數(shù)據(jù)量太大，耗費(fèi)時(shí)間長，處理代價(jià)高，實(shí)時(shí)性差。研究表明，在一般情況下，真正與等值面相交的體元只占總數(shù)量的很小部分(至多10%)，算法執(zhí)行過程中30%～70%的時(shí)間是花費(fèi)在空體元的檢測上。

3　本文所采用的方法

本文所使用的方法是一種基于Marching Cubes和球形繪制的混合算法，使Marching Cubes算法中檢測的空體元數(shù)量減少，從而提高算法的效率。算法主要思想為，將預(yù)處理得到的立方體頂點(diǎn)序列進(jìn)行分層檢測，從x，y，z三個(gè)維度進(jìn)行檢測，根據(jù)設(shè)定的閾值，分離出標(biāo)記為“1”的頂點(diǎn)密度大與密度小的區(qū)域，并針對這兩種區(qū)域進(jìn)行獨(dú)立處理，密度小的區(qū)域中標(biāo)記為“1”的頂點(diǎn)直接使用圓球進(jìn)行繪制，對于密度大的區(qū)域使用Marching Cubes算法進(jìn)行繪制，在繪制過程中，通過簡化等值面與立方體邊界角點(diǎn)的計(jì)算公式以達(dá)到降低計(jì)算量的目的。

(1) 區(qū)域分離

先從x方向進(jìn)行區(qū)域分離。首先設(shè)置閾值Nmin，在后面的檢測中，將以該閾值作為x方向分割的依據(jù)。x方向的區(qū)域分離分成兩個(gè)步驟，即x軸坐標(biāo)從最小值增加的方向和x軸坐標(biāo)從最大值減少的方向。

首先從x增加方向進(jìn)行檢測與分離。讓x的值從最小值xmin開始增加，每增加一層，x的值增加Δx，計(jì)算每層中被標(biāo)注為1的立方體頂點(diǎn)數(shù)n。如果n小于被設(shè)定的閾值Nmin，那么將該層分割出來，作為球形繪制的區(qū)域。當(dāng)x增加到某個(gè)值xmin+k·Δx時(shí)該層標(biāo)注為1的立方體頂點(diǎn)數(shù)n的值不小于閾值Nmin，這時(shí)停止該方向的檢測。

然后，從x減少的方向進(jìn)行檢測與分離。讓x的值從最大值xmax開始減少，每減少一層，x的值減少Δx，計(jì)算每層中被標(biāo)注為1的立方體頂點(diǎn)數(shù)n。如果n小于被設(shè)定的閾值Nmin，那么將該層分割出來，作為球形繪制的區(qū)域。當(dāng)x減少到某個(gè)值xmax-l·Δx時(shí)該層標(biāo)注為1的立方體頂點(diǎn)數(shù)n的值不小于閾值Nmin，這時(shí)停止該方向的檢測。

對每次得到的長方體經(jīng)過三個(gè)維度上的分離，每次都要設(shè)置新的閾值，按上述步驟進(jìn)行處理，得到的分離結(jié)果如圖2所示，將原長方體劃分成以下七個(gè)長方體。其中前六個(gè)長方體中的點(diǎn)將使用圓球直接對其進(jìn)行繪制。長方體九中的點(diǎn)將使用Marching Cubes算法進(jìn)行繪制。

圖2　分離結(jié)果

(2)使用球形進(jìn)行繪制

用OpenGL繪制球形只要直接調(diào)用GLUT工具包中的glutSolidSphere()函數(shù)即可。該函數(shù)的原型為void glutSolidSphere(GLdouble radius , GLint slices , GLint stacks)。其中，參數(shù)radius代表球形的半徑，slices代表以豎直軸上線段為直徑分布的圓周線的條數(shù)，stacks代表圍繞在豎直軸周圍的線的條數(shù)。slices與stacks這兩個(gè)參數(shù)決定了圓球繪制的精度，這兩個(gè)值越大那么繪制出來的圓球越精細(xì)，但同時(shí)帶來的是計(jì)算量的上升。使用glutSolidSphere( )函數(shù)繪制圓球時(shí)，需要將其球心的位置平移至要繪制的位置，使用glTranslatef( )函數(shù)進(jìn)行繪制位置的平移。

(3)構(gòu)建體元

在(1)中，將長方體進(jìn)行了一系列的分離，最終確定了使用Marching Cubes算法進(jìn)行繪制的區(qū)域，即(1)中所得到的長方體九，將長方體中的立方體當(dāng)成體元即可。

(4)角點(diǎn)的標(biāo)記

將粒子預(yù)處理后轉(zhuǎn)化為立方體頂點(diǎn)點(diǎn)陣，在立方體頂點(diǎn)中添加一位標(biāo)志位。立方體頂點(diǎn)的標(biāo)志位為“1”時(shí)，該角點(diǎn)標(biāo)記為黑色；立方體頂點(diǎn)的標(biāo)志位為“0”時(shí)，該角點(diǎn)不標(biāo)記。

(5)確定與等值面相交的體元邊界

由2.2節(jié)可知，等值面的種類可縮小為15種[8]。根據(jù)這十五種等值面分布類型可以構(gòu)建出一個(gè)查找表，查找表的長度為256，記錄著所有情況下等值面的連接方式，因此根據(jù)體元的角點(diǎn)狀態(tài)就可以從查找表中尋找到與等值面相交的體元邊界以及等值面的連接方式。

(6)計(jì)算體元和等值面交點(diǎn)位置

Marching Cubes算法中，計(jì)算體元與等值面的交點(diǎn)位置是其計(jì)算量的重要組成部分之一，因此這一步的計(jì)算量將很大程度上影響整體的繪制速度。本文選用計(jì)算量較小的中點(diǎn)插值計(jì)算公式：

(5)

其中，P代表等值點(diǎn)坐標(biāo)，P1、P2代表兩個(gè)端點(diǎn)的坐標(biāo)。

(7)計(jì)算體元和等值面交點(diǎn)處的法向量

與上一步的計(jì)算公式相對應(yīng)，本文計(jì)算等值點(diǎn)的法向量也是采用中點(diǎn)插值的計(jì)算公式：

(6)

其中，N代表等值點(diǎn)處的法向量，N1、N2代表兩個(gè)端點(diǎn)處的法向量。

(8)遍歷所有層的數(shù)據(jù)

(9)根據(jù)三角面片的頂點(diǎn)坐標(biāo)和該位置上的法向量繪制出等值面圖象

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

上述算法已經(jīng)在移動(dòng)工作站(2.8 GHz Intel處理器，4.0G內(nèi)存，NVIDIA QuadroFX 770M顯卡，顯存512M)上利用標(biāo)準(zhǔn)OpenGL圖形庫實(shí)現(xiàn)。我們將粒子分布設(shè)置為長方體的形狀，長方體分為兩個(gè)部分，其中一個(gè)部分的粒子密度比較大，一個(gè)部分的粒子密度比較小。使用三種方法進(jìn)行繪制，這三種方法分別是球形繪制、Marching Cubes算法、本文所用的方法。每種方法繪制的結(jié)果分成兩種形式，一種為俯視圖，一種為整體圖。球形繪制的結(jié)果如圖3(b)與4(a)所示，從其俯視圖可以發(fā)現(xiàn)，使用球形繪制時(shí)，密度比較小的區(qū)域繪制結(jié)果比較好，但是密度較大的區(qū)域會出現(xiàn)明顯的凹凸，表面不夠平整，效果欠佳。而從其整體圖也可以發(fā)現(xiàn)這個(gè)問題，即表面的光滑度不好。使用傳統(tǒng)的Marching Cubes算法進(jìn)行繪制的結(jié)果如圖3(a)和4(b)所示，從圖中可以發(fā)現(xiàn)，Marching Cubes算法繪制出來的結(jié)果中密度較大的區(qū)域繪制結(jié)果比較光滑，但是邊緣與真實(shí)情況差異較大。本文所使用的方法如圖3(c)和4(c)所示，該方法對粒子密度較大的區(qū)域采用Marching Cubes算法進(jìn)行繪制，而對邊緣的粒子采用球形繪制方法進(jìn)行繪制。

圖3　繪制結(jié)果的俯視圖

圖4　繪制結(jié)果的整體效果圖

最后，我們將該方法應(yīng)用于血流繪制，如圖5所示。其中，圖5(a)是利用本算法繪制出來的全局視角的血流圖。我們可以看到在血流最靠近邊緣的位置上有部分離散點(diǎn)，這些點(diǎn)就是使用球形繪制的；中間的部分比較平整，是使用Marching Cubes算法進(jìn)行繪制的。圖5(b)是同一次實(shí)驗(yàn)中將視角改變后所看到的結(jié)果。由于在改變視角的同時(shí)，粒子還處于運(yùn)動(dòng)情況，因此與圖5(a)并不是同一幀。將視角轉(zhuǎn)變?yōu)楦┮暤那闆r，從圖中依然可以看出邊緣是球形繪制，中間區(qū)域是Marching Cubes算法繪制。圖6分別表示血流在第40幀、第50幀、第60幀時(shí)的繪制結(jié)果。從圖中可以看到，本文的繪制算法結(jié)合了球形繪制和Marching Cubes算法，繪制出來的結(jié)果比較真實(shí)可靠。

(a)正常視圖　　　　　　(b)俯視圖(不同幀)圖5　不同視角的血液流動(dòng)現(xiàn)象繪制

(a)第40幀　　　　(b)第50幀　　(c)第60幀圖6　不同幀的繪制結(jié)果

5　結(jié)論

本首先將粒子的位置進(jìn)行預(yù)處理，使得分布雜亂無章的粒子轉(zhuǎn)變?yōu)榱⒎襟w頂點(diǎn)點(diǎn)陣。得到立方體頂點(diǎn)點(diǎn)陣之后，根據(jù)實(shí)時(shí)性的要求，將傳統(tǒng)的面

繪制算法Marching Cubes算法與傳統(tǒng)的球形繪制方法相結(jié)合，形成一種新的繪制方法。該方法首先將繪制區(qū)域劃分成七個(gè)區(qū)域，其中六個(gè)區(qū)域采用傳統(tǒng)的球形繪制算法進(jìn)行繪制，另外一個(gè)區(qū)域采用精簡版的Marching Cubes算法進(jìn)行繪制。實(shí)驗(yàn)表明，該方法可以實(shí)時(shí)、真實(shí)地繪制出分布雜亂無章的粒子群。

[1] Rianto S, Li L. Fluid dynamic visualisations of cuttings-bleeding for virtual reality heart beating surgery simulation[C]// Australasian Computer Science Conference, Brisbane, Australia, January. Australian Computer Society, Inc. 2010:53-60.

[2]Qin J, Pang W M, Nguyen B P, et al. Particle-based simulation of blood flow and vessel wall interactions in virtual surgery[C]// Symposium on Information and Communication Technology, Soict 2010, Hanoi, Viet Nam, August. 2010:128-133.

[3]Reorowicz P, Obidowski D, Klosinski P, et al. Numerical simulations of the blood flow in the patient-specific arterial cerebral circle region.[J]. Journal of Biomechanics, 2014, 47(7):1642-51.

[4] 黃雷,肖雙九,顧力栩,等.虛擬手術(shù)訓(xùn)練系統(tǒng)的血流模擬[J].計(jì)算機(jī)應(yīng)用與軟件,2011,28(1):65-68.

[5] 鄭廣超. 虛擬手術(shù)系統(tǒng)中模擬手術(shù)場景的渲染和平臺的構(gòu)建[D]. 上海:上海交通大學(xué), 2008.

[6] 施鵬, 熊岳山, 徐凱, 等. 虛擬肝臟手術(shù)中實(shí)時(shí)動(dòng)態(tài)滲血效果模擬[J]. 計(jì)算機(jī)應(yīng)用, 2013, 33(10): 2911-2913.

[7] Lorensen W E, Cline H E. Marching cubes: A high resolution 3D surface construction algorithm[C]//ACM Siggraph Computer Graphics. ACM, 1987, 21(4): 163-169.

[8] Mor A B. Progressive cutting with minimal new element creation of soft tissue models for interactive surgical simulation[D]. Pittsburgh：Carnegie Mellon University, 2001.

(責(zé)任編輯：曾晶)

Research on Real-Time Blood Flow Simulation Based on the Improvement Marching Cubes Algorithm

WANG Na1， CHEN Guodong2*，CHEN Yi2

(1. Department of Electronics and Information Engineering, Fuqing Campus of Fujian Normal University, Fuqing 350300, China；2. Institute of Physics and Information Engineering, Fuzhou University,Fuzhou 350108, China)

The conventional Marching Cubes algorithm takes too much time in detecting null voxels dataset, which leads to the low efficiency. To resolve this issue and to conduct visual simulation of fluid surface images, an improved algorithm was proposed. The density of the rule lattices was detected and a threshold value was set. The particle density below the threshold region was isolated from the one above the threshold region. Only the regions with high density were simplified version of Marching Cubes algorithm. The simulation experiment results show that compared with Spherical Rendering algorithm and Marching Cubes, the Marching Cubes algorithm proposed in this paper decimates the access to null voxels with good rendering performances and it enables the real-time rendering of the model.

Marching Cubes (MC) algorithm; voxels; blood flow; isosurface; threshold

1000-5269(2016)02-0084-04

10.15958/j.cnki.gdxbzrb.2016.02.19

2015-10-07

國家自然科學(xué)基金(61471124)；福建省自然科學(xué)基金項(xiàng)目(2016J01293)；福建省科技計(jì)劃重點(diǎn)項(xiàng)目(2011H0027)；福建省中青年教師教育科研項(xiàng)目(JA15574)

王娜(1978-)，女，副教授，碩士，研究方向：計(jì)算機(jī)圖形學(xué)和虛擬現(xiàn)實(shí)，Email：studyres@126.com .

陳國棟，Email：fzucgd@126.com.

TP391.41

A