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

【發(fā)明內(nèi)容】

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