基于sph的虛擬血管造影手術造影劑擴散過程模擬方法
【技術領域】
[0001 ] 本發(fā)明涉及一種基于SPH(Smoothed Particle Hydrodynamics)的虛擬血管造影 手術造影劑擴散過程模擬方法，屬于計算機仿真建模和虛擬手術領域。
【背景技術】
[0002] 心血管疾病目前是世界上死亡人數(shù)第一的疾病，微創(chuàng)介入手術是現(xiàn)代醫(yī)學手術的 一個巨大突破，其能夠有效的治療各類心血管疾病并且具有出血少、創(chuàng)傷低、恢復快、并發(fā) 癥少等優(yōu)點。其中血管造影術是微創(chuàng)介入手術中一個非常重要的部分。在進行微創(chuàng)介入手 術時醫(yī)生通過觀看X-ray產(chǎn)生的醫(yī)學圖像來進行手術，而X-ray會直接穿透血管使得最后在 醫(yī)學成像設備上無法清晰呈現(xiàn)血管的結構與狀態(tài)。由于造影劑原子量高、比重大在X-ray下 能夠清晰的成像，所以醫(yī)生通過外部注射器注射造影劑血管進行顯影來對血管病變處的血 管結構和病變位置進行診斷。目前有許多研究者對造影劑在血管中的擴散提出了許多方 法：
[0003] Cotin S等在"ICTS，an interventional cardiology training system" (Studies in healthtechnology and informatics ,vol. 70 ,pp. 59-65,2000)一文中米用 層流對血流進行建模，采用對流-擴散模型對造影劑混合與擴撒進行模擬。
[0004] Duriez C等在 "New approaches to catheter navigation for interventional radiology simulation，(Computer Aided Surgery,vol. 11 ,pp.300-308,2006)一文中米 用泊肅葉定律對血流進行建模并采用對流-擴散方程來模擬造影劑在血流中的擴散。
[0005] ffu X 等在 "Real-time modeling of vascular flow for angiography simulation"（Proc.Medicallmage Computing and Computer-Assisted Intervention (MICCAI 07)，pp. 557-565，2007)-文中同樣采用泊肅葉定律對血流進行建模，其次采用對 流-擴散方程來模擬造影劑在血流中的擴散。
[0006] F.Wang等在"A computer-based real-time simulation of interventional radiology"（Proceedings of the 29th Annual International Conference of the IEEE EMBS，pp. 108-115,2007)-文中使用血管中心線并采用恒定的速度對造影劑在血管 中進行傳播擴散。
[0007] 以上這些方法都采用了簡單化的Navier-stokes(NS)方程即泊肅葉定律或者層流 的方式對血流在血管中的流動進行建模，并對造影劑在血管中的擴撒采用對流擴散模型進 行擴散仿真模擬。以上方法忽略了"血液-造影劑"二相混合流體間的動力學交互對造影結 果的影響。由于造影劑具有良好水溶性，當造影注入到血管內后與血液迅速在血液中擴散 開來并與血液混合形成"血液-造影劑"二相混合流體，所以我們采用基于SPH的多相混合流 模型對"血液-造影劑"二相混合流體進行建模并模擬造影劑在血液中擴散造影的過程。

【發(fā)明內容】

[0008] 鑒于以上方法所述現(xiàn)有技術存在的不足，本發(fā)明的目的在于提供一種基于SPH的 虛擬血管造影手術造影劑擴散過程模擬方法，能夠實時準確的模擬造影劑在血液中的造影 擴散過程。
[0009] 為達到上述目的，本發(fā)明采用以下技術方案：
[0010] 一種基于SPH的虛擬血管造影手術造影劑擴散過程模擬方法，包括模擬初始化和 模擬循環(huán)兩個步驟，具體如下：
[0011] (1)模擬初始化:對血管邊界條件進行構建，設置血管入口與出口邊界條件并對 "血液-造影劑"二相混合流粒子位置進行初始化；
[0012] (2)模擬循環(huán):通過"血液-造影劑"二相混合流體模型并采用邊界粒子修正邊界處 的混合流密度和應用無滑移邊界條件進行循環(huán)模擬的步驟依次為：更新空間網(wǎng)格鄰居信 息，更新邊界粒子屬性，更新混合流體粒子每一相的體積分數(shù)并對粒子每一相體積分數(shù)及 壓強值進行修正，計算混合流體粒子的加速度，進行碰撞檢測與粒子速度、位置的更新，造 影渲染。
[0013] ⑴模擬初始化
[0014] 首先，對邊界粒子的初始化步驟如下：
[0015] 1.以血管模型的每一個三角形面片的邊作為單位采樣向量，即^和為采 樣間隔。因此對每個三角形采樣點的位置為：
[0016]
[0017] c\ =/7/,(,(., < 1
[0018] < l
[0019] 其中血和肥為大于〇的正整數(shù)。
[0020] 2.對步驟1中采樣得到的點沿著采樣點對應三角形向外的法向量從1到η倍進行外 推，其中每一次外推的單位長度為采樣間隔|，η的最大值滿足和為SPH 模型中的光滑核半徑。外推得到的點的位置即為邊界粒子的位置，外推點對應三角形的法 向量即為外推向量，為距離最近三角形的距離。
[0021] 3.優(yōu)化通過步驟2得到的邊界粒子使其均勻覆蓋在均勻覆蓋在血管外壁表面。首 先，去除通過外推后位置在血管內部的邊界粒子;其次，通過兩個相鄰粒子的外推向量的內 積?)判斷相鄰粒子所處的位置，?./ζ 20說明粒子處于相鄰三角形間，反之處于血管 的分叉處，再分別對處于分叉和相鄰三角形間密集的邊界粒子分別采用在半徑閾值1和丫2 內進行優(yōu)化，其中Τ, <€，Τ2 <￡ β通過中間位置的粒子來替換需要進行優(yōu)化的兩個相鄰的 粒子，設置新粒子的外推向量為前兩個相鄰粒子外推向量和的規(guī)范化的值，并設置新粒子 的離最近三角形的距離的值為前兩個相鄰粒子的均值。
[0022] 其次，設置血管入口與出口邊界條件血管入口采用入口壓強Pir^P血流入口初速 度，出口采用出口壓強P?t;初始化"血液-造影劑"二相混合流粒子的位置使其充滿血管。 [0023] (2)模擬循環(huán)
[0024] 1.更新空間網(wǎng)格鄰居信息
[0025]采用空間稀疏哈希的方法對空間模擬粒子的鄰居信息進行更新。首先對每個模擬 空間內的粒子根據(jù)其位置哈希映射得到其對應的空間網(wǎng)格索引。將該空間網(wǎng)格索引作為 key，模擬的粒子索引作為value進行快速計數(shù)排序，排序后在同一網(wǎng)格下的粒子索引即為 連續(xù)存儲，最后再得到排序后的同一網(wǎng)格下粒子索引的起始位置和結束位置。這樣當在進 行鄰居查找時，對當前網(wǎng)格和其周圍的26個網(wǎng)格，通過起始和結束位置訪問排序后的結果 來得到同一網(wǎng)格下的所有粒子索引，即可快速完成鄰居查找。
[0026] 2.更新邊界粒子屬性
[0027] 邊界粒子的位置和質量是固定的，邊界粒子的密度和壓強值為距離最近的混合流 體粒子對應的屬性值，邊界粒子的速度為：
[0028]
[0029]
[0030] 其中設置β為約束ζ的大小，db為邊界粒子距離最近三角形的距離，pb為邊界粒子 的位置，Pf為流體粒子的位置，?為邊界粒子的外推向量。
[0031] 3.更新混合流體粒子每一相的體積分數(shù)并對粒子每一相體積分數(shù)及壓強值進行 修正
[0032] (a)混合流體中的每一相k的每一個粒子i的連續(xù)方程求解體積分數(shù)的變化，該連 續(xù)方程為：
[0033]
[0034] 其中，ak為粒子中相k的體積分數(shù)，而且需要滿足1% = ，um為混合流體 t 速度，Ui為相k漂移速度，漂移速度Umk表示為：
[0035]
[0036]其中，τ為控制慣性和壓力大小的參數(shù)，τ的值越大慣性和壓力效果越大，〇為控制 ahpk 擴散速度的參數(shù)，σ的值越大其擴散速度越快，％ = 為粒子i相k的質量分數(shù)，加速 Pmi = g-((U", ·ν)ιι," + 粒子i相k壓強梯度VA與體積分數(shù)梯度^%分別表示為： dt ，
[0039] pki = akipmi[0040] !>m = ^P： - P,m)
[0037]
[0038]
[0041]
[0042] 其中，Pkl為粒子i相k的壓強，pmi為粒子i混合流體壓強，光滑核函數(shù)，κ為氣態(tài) 常數(shù)，為粒子j的質量，Α為粒子i的對所有粒子j通過光滑核函數(shù)插值得到的插值密 度，其中粒子j包括邊界粒子。
[0043] (b)通過(a)中求得體積分數(shù)的變化，其后通過前一時刻的體積分數(shù)即可求得當前 時刻粒子每一相的體積分數(shù)。由于當前時刻粒子i每一相a k不一定滿足 k 所以需要對體積分數(shù)進行調整。首先對于粒子每一相體積分數(shù)ak<0,設置ak = 0,其后再縮 放該粒子中的所有相體積分來滿足Σ%=1°其中對于ak2 1或者Xkak2 1的情況，在體積縮 k 放后即可完成調整。由于體積分數(shù)參與壓力項的計算，所以需要對壓力項進行調整，對于修 正后變化的Δ Σ-，所以修正的壓強為氣^ = & +Δ/\。. k
[0044] 4.計算混合流體粒子的加速度
[0045] 混合流體中的每一相k的每一個粒子i的動量方程求解流體粒子的加速度，該動量 方程為：
[0046]
[0047] 其中，Α? =_1〃^_為混合流體密度，Pk為相k的靜止密度，g為重力加速度，Tm為混 k 合流體的黏滯應力張量，TDm為相間對流的動量轉移量，F(xiàn)e為血管壁對混合流粒子的外力。其 中邊界粒子不參與流體體積分數(shù)的傳遞，邊界粒子對流體加速度的貢獻通過壓強的梯度 ▽p,"和粘滯應力張量▽ ·τΜ來滿足無滑移邊界條件。其中動量方程右邊各項分別表示為：
[0055]其中，μL= Skakyk為粒子i中各相粘度yk的聚合粘度，ri為粒子i的位置。f r為流體 粒子與血管內表面碰撞時發(fā)生的摩擦力，μ為血管內表面的靜摩擦系數(shù)，為粒子i碰撞時 的切線速度。Fn為混合粒子碰撞血管壁的法向壓力，ri為接觸剛度，δ為粒子碰撞時的穿透深 度。fs為血管內表面對流體粒子的粘性力，當它們間的距離小于閾值ds時對粒子才產(chǎn)生相應 的粘性力，△ t為模擬時間步長，vk為血管內表面對各相的粘性系數(shù)，cU為粒子i距離表面最 近的距離，η為血管內表面的法向量。fb為血壓對流體粒子產(chǎn)生的推力，通過設置入口和出 口壓強以及入口血液初速度來描述。
[0056] 5.碰撞檢測與粒子速度、位置的更新
[0057] (a)碰撞檢測：以模擬的流體粒子的當前位置p?r為原點，以速度方向￥_為射線方 向，以v A t為射線長度，對血管做AABB快速射線檢測，對發(fā)送碰撞的粒子進行標記。
[0058] (b)位置與速度更新:粒子與血管壁發(fā)生碰撞的粒子設置其位置移動到血管壁表 面上，未與血管壁發(fā)生碰撞的粒子則通過速度更新粒子至新的位置pnext = pcur+VAt。其中， 為當前時刻的位置，為當前時刻粒子速度，A t為模擬時間間隔。再將所有的粒子的 速度