基于GPU加速的粒子流體動(dòng)力學(xué)流血模擬算法

羅中粟　潘一源　唐良甫　朱珂權(quán)　李永強(qiáng)

1(浙江工業(yè)大學(xué)國(guó)際學(xué)院　浙江 杭州 310023) 2(浙江工業(yè)大學(xué)信息工程學(xué)院　浙江 杭州 310023)

0　引　言

近年來(lái)，微創(chuàng)手術(shù)[1]得到了快速的發(fā)展和廣泛的應(yīng)用。但是，因?yàn)槲?chuàng)手術(shù)本身的特點(diǎn)，醫(yī)生需要透過(guò)內(nèi)窺鏡來(lái)觀察患者體內(nèi)的局部狀況，這樣就大大增加了醫(yī)生的手術(shù)難度，而且對(duì)醫(yī)生來(lái)說(shuō)需要擁有一定的經(jīng)驗(yàn)。因此，術(shù)前手術(shù)訓(xùn)練對(duì)醫(yī)生來(lái)說(shuō)就顯得非常重要。虛擬手術(shù)訓(xùn)練系統(tǒng)利用觸覺(jué)交互設(shè)備和計(jì)算機(jī)虛擬現(xiàn)實(shí)技術(shù)實(shí)現(xiàn)了手術(shù)過(guò)程的虛擬化，為醫(yī)生提供了一個(gè)沉浸式的訓(xùn)練平臺(tái)，被廣泛應(yīng)用于醫(yī)學(xué)培訓(xùn)與手術(shù)計(jì)劃中[2]。而組織形變與流血模擬等視覺(jué)效果是手術(shù)訓(xùn)練系統(tǒng)的重要組成部分。

真實(shí)的流血模擬大多采用光滑粒子流體動(dòng)力學(xué)(SPH)方法[3]。將流體離散成帶有屬性和體積的粒子，結(jié)合有限元的思想為每個(gè)粒子求解方程得到下一時(shí)刻的狀態(tài)，粒子的屬性通過(guò)對(duì)周圍鄰居粒子的屬性使用核函數(shù)加權(quán)求和得到。2003年Muller等[4]第一次將其用于流體模擬，然而，該工作未解決流體的不可壓性。為了解決這個(gè)問(wèn)題，Becker等在2007年將Tait方程引入SPH方法，實(shí)現(xiàn)了弱可壓SPH，或簡(jiǎn)稱WCSPH[6]。該方法部分解決了流體的不可壓性，然而引入Tait方程的代價(jià)也是很大的，為了維持流體的穩(wěn)定性，模擬的時(shí)間步長(zhǎng)也必須相應(yīng)縮短到10-5級(jí)，使得模擬的進(jìn)程非常緩慢。Solenthaler等[11]在2009年提出了預(yù)估校正不可壓SPH，即PCISPH，有效地解決了均勻密度單相流的不可壓?jiǎn)栴}，模擬的步長(zhǎng)也提升了兩個(gè)數(shù)量級(jí)。但是，算法的計(jì)算量越來(lái)越大，其計(jì)算實(shí)時(shí)性問(wèn)題一直是虛擬手術(shù)訓(xùn)練應(yīng)用的瓶頸。

為了解決大量粒子條件下的實(shí)時(shí)性問(wèn)題，首先我們使用PCISPH方法對(duì)粒子建模，然后對(duì)目標(biāo)區(qū)域劃分均勻網(wǎng)格并且使用最近相鄰粒子搜索法遍歷粒子，結(jié)果顯示算法的效率明顯提高。此外，本文使用CUDA技術(shù)[8-9]來(lái)提高計(jì)算的速度。同時(shí)本文采用一種改進(jìn)的移動(dòng)立方體方法對(duì)流體表面進(jìn)行實(shí)時(shí)渲染，增強(qiáng)了手術(shù)訓(xùn)練真實(shí)性。

1　光滑粒子流體動(dòng)力學(xué)

1.1　PCISPH算法

SPH是把流體離散成粒子的集合，每個(gè)粒子都攜帶著壓強(qiáng)、速度、加速度等屬性，對(duì)于空間中任意給定一點(diǎn)r的任意物理量A而言，使用SPH方法可以通過(guò)如下離散求和計(jì)算[10]：

(1)

式中：rj代表粒子j的位置，mj是指粒子j的質(zhì)量，ρj表示密度，W(r,h)表示半徑為h的光滑核函數(shù)。

SPH的一大優(yōu)勢(shì)在于物理量的空間導(dǎo)數(shù)可以直接通過(guò)平滑核的導(dǎo)數(shù)加權(quán)和得到，即：

(2)

以及：

(3)

在歐拉公式中，等溫流體由速度v、密度ρ和壓力p描述。這些量隨時(shí)間的變化由以下兩個(gè)公式表示：

(4)

(5)

式中：g表示重力加速度，μ表示液體的粘度。在SPH方法中，粒子的質(zhì)量和數(shù)量是恒定的，所以式(4)可以省略并且式(5)可以轉(zhuǎn)換為：

f=-▽p+ρg+μ▽2v=fp+fg+fv

(6)

根據(jù)以上的規(guī)則，可以得到各個(gè)物理量的計(jì)算公式，如對(duì)于粒子i，它的密度可以通過(guò)領(lǐng)域的加權(quán)和計(jì)算：

(7)

壓力：

(8)

粘性力：

(9)

由于平滑核的選擇決定了SPH方法的穩(wěn)定性、精確性和計(jì)算速度，我們使用了平坦和歸一化的內(nèi)核以便得到二階精度的插值。我們?cè)O(shè)計(jì)了以下核函數(shù)：

(10)

(11)

(12)

為了計(jì)算簡(jiǎn)便，本文使用一種不同但等價(jià)于文獻(xiàn)[11]方式進(jìn)行闡述。已知流體在不可壓的情況下有：

(13)

(14)

若要使流體不可壓，就需要使得粒子的密度變化率為0，于是就可以設(shè)上一時(shí)間點(diǎn)的粒子密度總是等于ρ0，這樣便有：

(15)

PCISPH的算法流程如下：

(1) 對(duì)于全部的粒子建立空間索引。

(2) 根據(jù)粒子的加速度與時(shí)間步長(zhǎng)得到下一時(shí)刻的速度和位置。

(3) 初始化壓強(qiáng)和壓力pi=0，F(xiàn)p=0。

(4) 計(jì)算粒子在當(dāng)前壓力下的下一刻加速度、速度和位置。

(5) 計(jì)算粒子在這樣的位置下的領(lǐng)域密度ρi，以及密度差Δρ=ρi-ρ0。

(6) 根據(jù)Δρ對(duì)壓強(qiáng)pi進(jìn)行修正，并計(jì)算新的壓力Fp。

(7) 根據(jù)當(dāng)前壓力計(jì)算下一步速度和位置。

1.2　最近相鄰粒子搜索法

為了提高 PCISPH算法的效率，我們使用最近相鄰粒子搜索法遍歷粒子，它的主要思想就是在目標(biāo)域中構(gòu)建臨時(shí)的網(wǎng)格。我們構(gòu)建的臨時(shí)網(wǎng)格大小等于支持域的大小，在當(dāng)我們計(jì)算某一個(gè)粒子的最近相鄰粒子時(shí)，不需要遍歷目標(biāo)域中的所有粒子，而只需要遍歷某個(gè)所在的網(wǎng)格和相鄰的網(wǎng)格。此方法不僅縮小了研究問(wèn)題的區(qū)域，而且減少了參與計(jì)算的粒子數(shù)目，這樣，程序運(yùn)算的效率就提高了。假設(shè)支持域?yàn)閗h,那么我們構(gòu)建臨時(shí)網(wǎng)格的邊長(zhǎng)也為kh。在使用最近相鄰粒子搜索法搜索的過(guò)程中，如果我們要搜索某一粒子i的最近相鄰粒子，可以首先找出某一粒子i最近相鄰粒子所在的區(qū)域范圍，即他們要么處于和粒子i相同的網(wǎng)格內(nèi)，要么處于和粒子i網(wǎng)格空間相鄰的其他26個(gè)網(wǎng)格內(nèi)，所以我們只需要搜索粒子i所在網(wǎng)格和相鄰網(wǎng)格共27個(gè)網(wǎng)格內(nèi)的粒子。程序流程如圖1所示。

圖1　最近相鄰搜索法程序流程圖

1.3　Cuda加速

在模擬流血的過(guò)程中，為了保證模擬效果的真實(shí)性同時(shí)使得流血模擬更加的靈活，粒子數(shù)目不得不增加，但是這樣做會(huì)增加運(yùn)算的負(fù)擔(dān)，容易造成實(shí)時(shí)性的損害。為此，本文利用CPU和 GPU的并行運(yùn)算，將所有可以并行執(zhí)行的運(yùn)算放在CUDA上計(jì)算來(lái)提高流血模擬以及虛擬手術(shù)訓(xùn)練系統(tǒng)的實(shí)時(shí)性。整個(gè)程序的系統(tǒng)如圖2所示。

圖2　GPU加速流程圖

2　表面重建算法

為了得到真實(shí)的渲染結(jié)果，每進(jìn)行一步的SPH模擬都要將點(diǎn)云轉(zhuǎn)化為可著色的表面。本文的表面重建步驟如下:

(3) 對(duì)于三維密度場(chǎng)中每個(gè)體素，搜索其領(lǐng)域的粒子，根據(jù)體素中心到粒子距離進(jìn)行插值得到密度，插值時(shí)將距離乘以Si進(jìn)行調(diào)整。

(4) 對(duì)三維密度場(chǎng)進(jìn)行立方體匹配，得到三角形網(wǎng)格。

(5) 對(duì)三角形網(wǎng)格進(jìn)行表面細(xì)分，生成PN三角形，得到光滑的網(wǎng)格。

2.1　移動(dòng)立方體算法

移動(dòng)立方體[13]是一種經(jīng)典的從密度場(chǎng)生成三角形等值面的技術(shù)。它把密度場(chǎng)均勻劃分為格子，對(duì)格子的胞元的每個(gè)角點(diǎn)的密度場(chǎng)進(jìn)行評(píng)估，根據(jù)其是否大于某個(gè)常量閾值來(lái)決定生成三角形的形式。它的步驟如下：

(1) 檢查胞元立方體的12條邊是否與等值面相交，如果邊的兩個(gè)端點(diǎn)中一個(gè)的密度小于閾值且另一個(gè)的密度大于閾值則判為相交，否則判為不相交。每條相交的邊都生成一個(gè)對(duì)應(yīng)的三角形頂點(diǎn)。

(2) 通過(guò)胞元的8個(gè)角點(diǎn)是否大于等值面生成一個(gè)8位的索引(大于等值面該位則為1，否則為0)。每個(gè)索引對(duì)應(yīng)一種三角化的模式，因此通過(guò)索引便可得到哪些三角形頂點(diǎn)相連(如圖3所示)。

圖3　Marching Cubes算法的15種基本情況

(3) 對(duì)密度場(chǎng)進(jìn)行差分從而生成法線[14]：

(16)

當(dāng)使用移動(dòng)立方體算法時(shí)，我們需要計(jì)算離散點(diǎn)陣的每個(gè)單元格，因此生成的網(wǎng)格不夠光滑。當(dāng)我們以高分辨率繪制表面時(shí)，這種局限性將導(dǎo)致明顯的多邊形缺性。在本文中，我們使用PN三角形方法來(lái)平滑三角形網(wǎng)格邊緣和為頂點(diǎn)著色操作提供更多的采樣點(diǎn)。僅僅根據(jù)三角形的三個(gè)頂點(diǎn)及其法線就可以生成平滑的細(xì)分三角形。見(jiàn)圖4。

圖4　輸入數(shù)據(jù)：頂點(diǎn)Pi和法線Ni

如圖4所示：給定頂點(diǎn)P1,P2,P3∈R3和各個(gè)點(diǎn)的法線N1,N2,N3∈R3。首先計(jì)算出其對(duì)應(yīng)的Bezier曲面控制點(diǎn)bijk：

(1) 將bijk均勻分布在三角形上，三角形的三個(gè)頂點(diǎn)對(duì)應(yīng)三個(gè)控制點(diǎn)b300、b030、b003。

(2) 在每個(gè)角，將其對(duì)應(yīng)的兩條邊上的最近的點(diǎn)投影到與法線垂直的平面上，稱為切點(diǎn)，如在b300可以生成兩個(gè)切點(diǎn)b201、b210。以此類推在另兩個(gè)角生成b102、b012以及b021、b120，如圖5所示。

圖5　三角的貝塞爾曲面排列而成的控制網(wǎng)的系數(shù)

(3) 將中心點(diǎn)朝六個(gè)切點(diǎn)的平均點(diǎn)移動(dòng)，并越過(guò)移動(dòng)的方向再移動(dòng)1/2。

曲線PN三角形的系數(shù)和控制點(diǎn)用公式表示如下：

(17)

2.2　雙邊法線平滑

PN三角形可以解決網(wǎng)格過(guò)于粗糙的問(wèn)題，然而，為了生成光滑的表面，僅僅用PN三角形是不夠的。PN三角形細(xì)分出來(lái)的法線會(huì)導(dǎo)致照明呈凹凸?fàn)?，這主要是由于立方體匹配生成的網(wǎng)格中有大量的窄條，而在被PN三角形細(xì)分之后，這些窄條形的三角形會(huì)生成更多窄條形的三角形，從而使網(wǎng)格質(zhì)量顯著下降，見(jiàn)圖6。由于法線插值與三角形位置有關(guān)，因此低質(zhì)量的網(wǎng)格最終會(huì)導(dǎo)致法線不平整。一般來(lái)說(shuō)，為了提高網(wǎng)格質(zhì)量，可以對(duì)模型進(jìn)行重網(wǎng)格化，然而這類操作通常非常耗時(shí)，不適用于實(shí)時(shí)應(yīng)用。

圖6　使用PN三角形細(xì)分前后的網(wǎng)格

由于法線平滑與否主要影響屏幕繪制時(shí)的渲染效果，本文選擇在屏幕空間對(duì)法線進(jìn)行雙邊濾波的方法[15]來(lái)平滑法線。這樣做的優(yōu)勢(shì)主要有：(1) 速度快，二維雙邊濾波可以使用兩個(gè)分離的易于平行的一維濾波器近似；(2) 效果好，實(shí)驗(yàn)證明這種方法可以很好地去除光滑表面法線不平整的瑕疵。

首先將細(xì)分的網(wǎng)格光柵化至兩張紋理中，分別以浮點(diǎn)向量存儲(chǔ)位置以及法線信息。然后使用本文設(shè)計(jì)的雙邊濾波器對(duì)這兩張紋理進(jìn)行處理。本文定義它的形式為：

(18)

式中：Ω是像素的領(lǐng)域。

(19)

以及

(20)

該濾波器綜合了空間位置以及法向量自身的信息對(duì)法向量進(jìn)行平滑。它的設(shè)計(jì)概念主要源自兩點(diǎn)：

(1) 由于觀察到流體的表面曲率通常不高，幾乎不會(huì)出現(xiàn)直角和銳角，因此像素的法向相似，更可能屬于同一平面，應(yīng)允許平滑，不相似的更可能屬于不同平面，不應(yīng)該平滑。

(2) 像素的空間位置相近的更可能在流體表面上測(cè)的距離相近，因此應(yīng)允許平滑；空間位置相離較遠(yuǎn)的像素其測(cè)地距離通常更遠(yuǎn)，不應(yīng)該平滑。

因此，本文定義衡量空間位置的函數(shù)為高斯函數(shù)的形式：

(21)

同時(shí)定義法線相似函數(shù)為：

g(nx,ny)=saturate(nx·ny)

(22)

式中：saturate(x)為將x截?cái)嗟絒0,1]的截?cái)嗪瘮?shù)。

3　實(shí)驗(yàn)結(jié)果

本文所做的工作是在一個(gè)完整的病人特定的腦腫瘤切除模擬器開(kāi)發(fā)的背景下提出的。該模擬器提供了兩臺(tái)力反饋機(jī)器，有限元組織模型和在切除期間進(jìn)行拓?fù)溥m應(yīng)的模擬，真實(shí)的視覺(jué)反饋包括高分辨率的紋理和照明，平滑的投射陰影、深度和流血。

表1顯示了我們的方法在CPU和GPU條件下不同粒子數(shù)時(shí)的計(jì)算速度。實(shí)驗(yàn)結(jié)果表明雖然從CPU到GPU的數(shù)據(jù)傳輸需要很多時(shí)間，但使用CUDA技術(shù)顯著提高了計(jì)算速度。隨著粒子數(shù)量的增加，它能更有效地提高手術(shù)訓(xùn)練系統(tǒng)的實(shí)時(shí)性。

表1　采用CPU和GPU的模擬時(shí)間

此外，為了提高計(jì)算速度，我們使用最近相鄰粒子搜索法來(lái)提高PCISPH算法的效率。圖7顯示了該方法相比全局搜索提高了PCISPH算法的效率。

圖7　使用全員搜索方法和最近相鄰粒子搜索法的PCISPH仿真時(shí)間

圖8顯示了SPH粒子的分布和進(jìn)行表面渲染后的粒子分布。從圖中可以看出，網(wǎng)格的邊緣不是很平滑，所以在本文中，我們使用PN三角形法來(lái)細(xì)分三角形，結(jié)果表明，使用我們的表面重建算法后，表面變得更加光滑了。因此，我們的方法大大提高了虛擬手術(shù)訓(xùn)練系統(tǒng)的真實(shí)性。為了清楚地看到虛擬手術(shù)系統(tǒng)的實(shí)際操作，我們給出了在用于腦手術(shù)的虛擬手術(shù)模擬系統(tǒng)中直接渲染粒子的結(jié)果，如圖9所示。

圖8　不同方法對(duì)粒子表面渲染后的效果圖

圖9　虛擬手術(shù)模擬系統(tǒng)渲染粒子結(jié)果

4　結(jié)　語(yǔ)

本文提出一種方法，有效地將流血模擬集成到神經(jīng)外科手術(shù)模擬過(guò)程當(dāng)中。本文的方法實(shí)現(xiàn)了表面出血的真實(shí)模擬，包括抽吸和燒灼。本文采用PCISPH方法對(duì)血液進(jìn)行建模并采用CUDA加速實(shí)現(xiàn)，展示了與切割仿真系統(tǒng)結(jié)合的表面流血效果。實(shí)驗(yàn)結(jié)果表明：本文模擬的流血效果具有良好的靈活性，能夠滿足不同的模擬需求，而基于GPU的實(shí)現(xiàn)大大提高了計(jì)算速度，滿足了虛擬手術(shù)中實(shí)時(shí)性的要求。此外，本文所用的表面重建算法獲得了逼真的出血效果，大大提高了手術(shù)訓(xùn)練系統(tǒng)的真實(shí)性。然而，如果我們想將此技術(shù)用于醫(yī)院的實(shí)際外科手術(shù)訓(xùn)練中，該技術(shù)仍然有很多需要改進(jìn)的地方。例如，在實(shí)際外科訓(xùn)練系統(tǒng)中，流血模擬必須與變形，吸取和切割相結(jié)合。此外，流血模擬的實(shí)時(shí)性和真實(shí)性需要進(jìn)一步提高。我們將在未來(lái)的工作中進(jìn)一步去研究。

[1] Puangmali P,Althoefer K,Seneviratne L D,et al.State-of-the-Art in Force and Tactile Sensing for Minimally Invasive Surgery[J].IEEE Sensors Journal,2008,8(4):371-381.

[2] Wang G,Zhang S,Xie H,et al.A homotopy-based sparse representation for fast and accurate shape prior modeling in liver surgical planning[J].Medical Image Analysis,2015,19(1):176.

[3] Monaghan J J.Simulating Free Surface Flows with SPH[J].Journal of Computational Physics,1994,110(2):399-406.

[4] Chen Z,Qin J,Si W,et al.Particle-based fluid simulation with small scale details[C]//IEEE International Conference on Information Science and Technology.IEEE,2014:561-564.

[5] Li Qiang,Yu Guichang,Liu Shulian,et al.Application of Computational Fluid Dynamics and Fluid Structure Interaction Techniques for Calculating the 3D Transient Flow of Journal Bearings Coupled with Rotor Systems[J].Chinese Journal of Mechanical Engineering,2012,25(5):926-932.

[6] Becker M,Teschner M.Weakly compressible SPH for free surface flows[C]//ACM Siggraph/eurographics Symposium on Computer Animation,SCA 2007,San Diego,California,Usa,August.DBLP,2007:209-217.

[7] Ihmsen M,Akinci N,Gissler M,et al.Boundary Handling and Adaptive Time-stepping for PCISPH[C]//The Workshop on Virtual Reality Interactions & Physical Simulations.DBLP,2010:79-88.

[8] Sanders J,Kandrot E.CUDA by Example:An Introduction to General-Purpose GPU Programming[M].Addison-Wesley Professional,2010.

[9] Zhou Y,Tan Y.GPU-based parallel particle swarm optimization[C]//Evolutionary Computation,2009.CEC’09.IEEE Congress on.IEEE,2009:1493-1500.

[10] Xiao C,Feng Y,Li Y,et al.Real-time and authentic blood simulation for surgical training[C]//Chinese Control and Decision Conference.2017:6832-6837.

[11] Solenthaler B,Pajarola R.Predictive-corrective incompressible SPH[J].Acm Transactions on Graphics,2009,28(3):341-352.

[12] Jin hanjun,Liu xiaoya,Liu xiaoyu.Collision Detection for Fountain Model Based on Partical System[C]//2010 3rd IEEE International Conference on Computer Science and Information Technology.IEEE,2010:502-505.

[13] Wang C,Chang T,Lin M.Reconstruction and Representation for 3D Implicit Surfaces[J].Communications in Computer & Information Science,2011,202(2):364-372.

[14] 王旭初,王贊.基于最近鄰Marching Cubes的醫(yī)學(xué)圖像三維重建[J].計(jì)算機(jī)工程與應(yīng)用,2012,48(18):154-158.

[15] 王鵬程.基于粒子方法的流體實(shí)時(shí)仿真研究[D].北京工業(yè)大學(xué),2013.