王　娜，陳國棟，陳　怡

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

?

基于改進Marching Cubes算法的血流實時仿真研究

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

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

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

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

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

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

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

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

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

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

由于長方體分割為立方體時是均勻分割，因此，可以計算出與xi最近的立方體頂點的x軸坐標xa與xb，其中xa

(1)

(2)

(3)

(4)

圖1　尋找最近鄰的立方體頂點

2　Marching Cubes算法

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

2.1Marching Cubes算法的主要思想

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

Marching Cubes算法的流程如下：

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

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

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

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

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

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

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

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

2.2Marching Cubes算法存在的缺陷

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

3　本文所采用的方法

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

(1) 區(qū)域分離

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

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

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

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

圖2　分離結(jié)果

(2)使用球形進行繪制

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

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

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

(4)角點的標記

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

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

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

(6)計算體元和等值面交點位置

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

(5)

其中，P代表等值點坐標，P1、P2代表兩個端點的坐標。

(7)計算體元和等值面交點處的法向量

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

(6)

其中，N代表等值點處的法向量，N1、N2代表兩個端點處的法向量。

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

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

4　實驗結(jié)果及分析

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

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

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

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

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

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

5　結(jié)論

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

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

[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].計算機應(yīng)用與軟件,2011,28(1):65-68.

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

[6] 施鵬, 熊岳山, 徐凱, 等. 虛擬肝臟手術(shù)中實時動態(tài)滲血效果模擬[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é)基金項目(2016J01293)；福建省科技計劃重點項目(2011H0027)；福建省中青年教師教育科研項目(JA15574)

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

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

TP391.41

A