，

(福州大學(xué)物理與信息工程學(xué)院，福建 福州 350116)

0 引 言

在手術(shù)中，流血現(xiàn)象十分的常見，血流仿真的真實感和實時性就顯得十分必要，國內(nèi)外的大量專家學(xué)者對此進行了研究。

關(guān)于流體的模擬研究，一直是虛擬手術(shù)中研究的熱點和難點，從本世紀開始，流體模擬研究的新方法和新思路就開始大量涌現(xiàn)。2015年, 肖苗苗等人[1]使用異構(gòu)計算平臺Open CL對SPH(smoothed particle hydrodynamics)方法進行加速,充分利用了Open CL與圖形渲染API的數(shù)據(jù)互通；2015年，Ladicky等人[2]提出一種用數(shù)據(jù)流驅(qū)動的方法來模擬流體的流動。

血流模擬的主要難點是兼顧血流流動的真實感和實時性，上面提及的方法側(cè)重于改善流體模擬的實時性，而忽略了流體模擬的真實感。為此，在研究中引入了基于梯度的光流法。由于一階近似在大范圍速度場中誤差較大，提出一種新的迭代方法，優(yōu)化速度場的細節(jié)，增強血流模擬的真實感。針對FLIP方法中網(wǎng)格插值粒子過程中計算誤差較大的問題，對原有的插值方法進行改進，以減小誤差，提高仿真的精度。

1 光流法

擬引入光流法對流體表面的速度場進行求解，并在前人研究的基礎(chǔ)上對求解過程的相關(guān)方法進行改進，以提高求解精度。

光流法求解速度場可分為以下幾步：離散求解、迭代計算和表面投影。下面對求解過程作一個簡單的介紹。

1.1 離散求解

(1)

接下來引入兩個常見的正則化矩陣：一個平滑矩陣，和一個支持較小速度場矢量的矩陣(也叫做Tikhonov 正則化矩陣)。要計算速度場u，需要將公式(1)的數(shù)據(jù)項和兩個正則化矩陣最小化，然后最小化能量項的加權(quán)和，方便接下來計算速度場：

(2)

在最小二乘意義上的離散情況的能量最小化由下式給出

βS∑j|uj|2+βT|u|2dΩ

(3)

在式中，βS和βT分別表示對應(yīng)項的比例因子，∑j表示四個維度上的和。Φr是光流應(yīng)用的時間導(dǎo)數(shù), 即參數(shù)r引起的變化。通過Φ2來計算空間梯度，其中，Φ2是通過計算要得到的速度場。

公式(3)的離散最小化產(chǎn)生如下一個線性方程組

Aofu=b

(4)

通過上式計算u。矩陣Aof包含離散化的能量項，見下式。

Aof=[Φ2]T[Φ2]+βS∑jLj+βTI

(5)

公式(5)中的三項分別對應(yīng)前面提到的數(shù)據(jù)項、平滑度和Tikhonov項。[Φ2]表示一個n×4n的矩陣，這個矩陣包含Φ2的離散梯度，L表示離散的拉普拉斯矩陣。公式(4)中。

b=[-Φ2]TΦr，Φr=Φ2-Φ1

由于公式(1)中的第二項是非線性項，在實際應(yīng)用中很難線性化，所以在實際應(yīng)用過程中常常對對應(yīng)的導(dǎo)數(shù)采取一階近似的方法[3]。

1.2 迭代計算

迭代方法結(jié)合了遞歸重采樣的分層和在粗糙層次改變數(shù)據(jù)使之發(fā)生變化這兩種方法[4]，有助于減少一階近似所產(chǎn)生的誤差。

如果第一個通過光流法計算得到的速度場與預(yù)期理想的速度場較為接近，那么就可以除去有誤差的部分，用剩余的部分產(chǎn)生第二個速度場u2，這個過程稱之為剩余迭代，迭代進行到無法計算或者輸出結(jié)果的質(zhì)量不理想為止。

迭代的重新線性化使輸入更符合我們的要求，它與分層相結(jié)合的解決方案能顯著改善速度場的最終質(zhì)量。

1.3 表面投影

計算好速度場后，需要把速度場從源表面投影到目標網(wǎng)格上。給定一個通過光流法產(chǎn)生的速度場，執(zhí)行一個二分法搜索每個單元沿著目標的梯度，直到在網(wǎng)格上找到與源單元表面值相匹配的點。將投影作為估計當前速度場值uk的一個更新步驟，其遞歸公式如下。

uk+1(x)=uk(x)+bisect(x+uk(x),-

n(x+uk(x)))

(6)

公式(6)中，函數(shù)bisect表示第一個參數(shù)沿著第二個參數(shù)傳遞的方向進行的實際查找。在這種情況下，n是指目標Φ2的正定化梯度，且通常是體積內(nèi)點的反轉(zhuǎn)。

2 血流與固體的交互仿真

對血流與固體的交互過程進行仿真，首先將血液粒子化，用求得的速度場驅(qū)動血粒子運動完成血液的整體流動，

2.1 FLIP方法

血液粒子化后，引入了FLIP(Fluid Implicit Particle)[5]的方法實現(xiàn)網(wǎng)格到粒子的插值計算。FLIP方法的步驟如下。

1) 初始化血粒子位置和速度。

2) 在交錯的mac網(wǎng)格節(jié)點上，計算附近粒子速度的加權(quán)平均值，并將速度保存。

3) 用保存的節(jié)點速度減去更新后的節(jié)點速度，對每個粒子進行插值。

4) 以網(wǎng)格節(jié)點速度值驅(qū)動血粒子的運動。

圖1 網(wǎng)格到粒子的插值

針對第3步的插值方法，提出一種改進的平流算法[6]的方法，以減少插值計算的量。方法簡要概述如下：一個單參數(shù)x1和兩個帶有速度場u1→2和u2→1的輸入Ψ1和Ψ2。在這種情況下插值結(jié)果Ψ0由下式給出

(7)

(8)

其中advect是平流算法，p為插值點，r為參數(shù)。公式(7)即為常規(guī)的重心轉(zhuǎn)換。

2.2 渲 染

血流與粒子交互的過程中，需要將粒子抽取成表面并渲染。這一過程將用到Marching cube[7～8]算法，該算法的核心是抽取等值面并按一定的方式連接起來繪制一個整體的表面。

整個交互過程的流程圖如下

圖2 血流仿真整體流程圖

3 仿真過程及結(jié)果

實驗環(huán)境介紹如下。硬件環(huán)境為：Intel Core i5-3210M，主頻2.5GHz，4G內(nèi)存，顯卡為NVIDIA GEFORCE GT 635，顯存2G。軟件環(huán)境為：操作系統(tǒng)為Windows 7系統(tǒng)(64位)，開發(fā)環(huán)境為Microsoft Visual Studio 2013(Win 32)，算法語言為C++，生成manta可執(zhí)行程序，通過打開相應(yīng)的Python文件，調(diào)用OpenGL庫渲染來達到仿真的效果。

下圖是血流向下流動時與固體相遇交互的場景仿真，分別是第14、44、74、104和134幀的圖像。

圖3 血流與固體交互過程

圖3所示過程可以看出，血流自上而下流動，遇到固體時，血流分成兩股，分開后的血流在固體下方流動一段距離后又逐漸合成一股。下表是相關(guān)的統(tǒng)計數(shù)據(jù)。

表1 血流與固體交互仿真效率的分析

從表1可以看出，隨著流體粒子數(shù)的增加，更新粒子所用的時間相應(yīng)地變長了，但從74幀以后，更新流體粒子的時間間隔趨于穩(wěn)定，整體來看，改進后的算法處理效率高。

由于血流表面的投影是在歐拉網(wǎng)格上進行的，因此網(wǎng)格數(shù)對表面的仿真精度有著重要的影響，下表顯示了在不同網(wǎng)格數(shù)的情況下仿真動畫的FPS(幀率)。這里的幀率取的是平均值。并與文獻[9]中的FPS結(jié)果進行對比。

表2血流與固體交互仿真中不同網(wǎng)格數(shù)對應(yīng)的FPS

網(wǎng)格數(shù)文獻[10]的方法文中的方法16*16*1637.5959.2632*32*324.9611.4364*64*640.453.00

隨著網(wǎng)格數(shù)量的增加，計算量也相對增加，流體動畫的幀率在下降。研究中，在粒子數(shù)等變量一定的情況下，從表2可以看出，相比文獻[9]的方法，文中的方法在網(wǎng)格數(shù)相同的情況下，F(xiàn)PS都相對較高，故仿真的精度較高，血流仿真的效果較好。

4 結(jié) 語

從血流仿真的真實感的角度出發(fā)，研究中引入了基于梯度的光流法。由于一階近似在大范圍速度場中誤差較大，提出一種新的迭代方法，優(yōu)化速度場的細節(jié)，增強血流模擬的真實感。針對FLIP方法中網(wǎng)格插值粒子過程中計算誤差較大的問題，對原有的插值方法進行改進，以減小誤差，提高仿真的精度。實驗結(jié)果表明，文中的方法基本滿足血流仿真中真實感和實時性的要求。