虛擬手術(shù)中基于粒子的多相耦合方法
【專利摘要】一種虛擬手術(shù)中基于粒子的多相耦合方法，對于虛擬手術(shù)中涉及的軟組織和剛體，分別根據(jù)三維組織表面三角網(wǎng)格模型生成粒子化三維體模型，而對于液體則采用SPH方法建立粒子化液體模型；建立軟組織、液體和剛體的粒子化力學(xué)模型，包括將軟組織和液體定義為流變體，并對軟組織粒子和液體粒子的任意兩個同相或兩個不同相粒子之間定義流變粘滯系數(shù)；進行基于邊界粒子和多相耦合核函數(shù)的碰撞檢測及響應(yīng)。本發(fā)明提出一種了應(yīng)用于虛擬手術(shù)中的通用的、精確的和高效的基于粒子的多相耦合方法，能在虛擬手術(shù)仿真中實現(xiàn)精確、高效的多相耦合仿真。
【專利說明】虛擬手術(shù)中基于粒子的多相耦合方法
【技術(shù)領(lǐng)域】
[0001]本發(fā)明屬于計算機仿真及虛擬現(xiàn)實【技術(shù)領(lǐng)域】，特別是涉及一種應(yīng)用于虛擬手術(shù)中的通用的、精確的和高效的基于粒子的多相耦合方法。
【背景技術(shù)】
[0002]隨著計算機仿真計算和虛擬現(xiàn)實技術(shù)的快速發(fā)展，虛擬手術(shù)仿真研究及其相關(guān)應(yīng)用已成為國際上發(fā)展迅速的一個領(lǐng)域，以Simbionix為代表的虛擬手術(shù)仿真器已廣泛應(yīng)用于許多種手術(shù)的外科醫(yī)生訓(xùn)練和評價中。真實手術(shù)中往往涉及不同性質(zhì)物體的多相耦合(軟組織、液體和剛體之間的耦合)問題，如在腹腔手術(shù)中手術(shù)器械與體液、軟組織的接觸耦合，體液與組織、腫瘤的接觸耦合，這些耦合現(xiàn)象的仿真對精確性和效率要求較高。
[0003]近年來，研究人員提出了很多基于物理的計算方法用于多相耦合仿真?；谖锢淼臄?shù)值計算方法一般分為兩類，一類是基于網(wǎng)格的方法或Eluer方法，另一類是基于粒子的方法或Lagrangian方法。基于網(wǎng)格的方法,如有限元法(Finite Element Method,簡稱FEM方法)等使用廣泛，但在處理自由表面問題、邊界運動問題、表面運動問題和大形變及碰撞傳播問題時并不適用，并且對于復(fù)雜的幾何模型，要生成高質(zhì)量的網(wǎng)格是一個困難耗時的過程，無網(wǎng)格法則能有效地解決這些問題[I]。光滑質(zhì)點流體動力學(xué)(Smoothed ParticleHydrodynamics,簡稱SPH方法)是近年來無網(wǎng)格法中使用最為廣泛的一種。SPH方法本身的特性使其非常適合對高度可變形物體進行仿真[2]。
[0004]M.MUller提出了一個基于粒子的方法，借助SPH方法對彈性體、塑料體和融化物體進行建模[I]，能夠有效地仿真物體的形變。之后關(guān)于物體形變的研究大多基于M.MUller的工作，并且相關(guān)改進多集中于計算的簡化[3]以及穩(wěn)定性[4]的提高。一些方法簡單地把流體、固體粒子視為整體來計算鄰域[4]，但這將導(dǎo)致流體粒子進入固體的不真實現(xiàn)象。一些研究學(xué)者采用虛粒子法[5]在固體、流體周圍設(shè)置虛粒子用于耦合，但這并不能仿真雙向耦合過程，且虛粒子的設(shè)置與計算過程相當復(fù)雜。一種新穎的、通用的可以同時用于單向雙向耦合的流固耦合的邊界處理方法在文獻[6]中被提出。該模型將固體建模成剛體，因為是剛體，所以該模型只需要與固體最外層粒子進行計算，且不需要粒子法向量，同時計算的邊界粒子不需要均勻采樣，可以稀疏可以濃密，流固邊界粒子的受力是對稱的。但是該模型只適用于剛體，不能用于流體與實體的可變形的固體之間的耦合計算。
[0005]雖然近年來研究人員在多相耦合的領(lǐng)域取得了一定的成績，但是要在虛擬手術(shù)仿真中實現(xiàn)精確、高效率的多相耦合仿真，還需要解決現(xiàn)有方法所面臨的如下一些問題:
[0006](I)隨著虛擬手術(shù)仿真場景的復(fù)雜度的提高，現(xiàn)有方法處理復(fù)雜模型場景的復(fù)雜度越來越高，越來越不能適應(yīng)復(fù)雜場景的仿真，因此，亟需一種較高真實性的、適用于任意復(fù)雜模型、多種復(fù)雜場景仿真的通用多相耦合模型。且真實人體軟組織形狀大小復(fù)雜多變，現(xiàn)有方法無法對任意復(fù)雜的人體軟組織模型精確地生成內(nèi)部體數(shù)據(jù)，為了適應(yīng)復(fù)雜虛擬手術(shù)仿真場景及模型的滿足物理真實性，需要對任意復(fù)雜的人體軟組織模型精確地生成內(nèi)部體數(shù)據(jù)。[0007](2)現(xiàn)有大部分基于物理的耦合方法通常采用類似于SPH方法的無網(wǎng)格法來對流體建模，采用線彈性模型對軟組織建模，這些方法都缺乏對流體和軟組織性質(zhì)的考慮。例如在真實世界中，幾乎所有的生物軟組織都是由固態(tài)和液態(tài)兩種成分共同構(gòu)成，生物軟組織是最典型的流變體，體液具有松弛行為而軟組織具有蠕變行為，因此從虛擬手術(shù)仿真的真實性的角度考量，在建立流體和固體模型的時候需要引入流體和固體的性質(zhì)。
[0008](3)現(xiàn)有方法在處理多物體、多相的耦合過程(包括碰撞檢測和響應(yīng))比較復(fù)雜，在碰撞檢測和響應(yīng)的過程中耗費較多資源，且通常選取與運動粒子最近的粒子作為碰撞點，缺乏精確性，容易產(chǎn)生穿透現(xiàn)象，因此需要一個精確高效的碰撞檢測和響應(yīng)算法來實現(xiàn)多物體、多相的耦合過程。
[0009]相關(guān)文獻:
[0010][I] Mul I er，Μ.，Keiserj R.，Nealenj A.，Pauly, Μ.，Gross, Μ.，&Alexa, Μ.(2004，August).Point based animation of elastic, plastic and melting objects.1n Proceedings of the2004ACM SIGGRAPH/Eurographics symposium on Computeranimation(pp.141-151).Eurographics Association.[0011][2] Liu，G.G.R.，&Liu，Μ.B.(2003).Smoothed particle hydrodynamics: ameshfree particle method.World Scientific Publishing Company.[0012][3] Gerszewskij D.，Bhattacharyaj H.，&Bargtei I, A.W.(2009，August).A point-based method for animating elastoplastic solids.1n Proceedings ofthe2009ACM SIGGRAPH/Eurographics Symposium on Computer Animation(pp.133-138).ACM.[0013][4] Solenthaler, B.，SoMifli5 J_，&Pajarola, R.(2007) _ A unified particlemodel for fluid - sol id interactions.Computer Animation and VirtualWorlds, 18(1),69-82.[0014][5] Schechter，H.，&Bridson，R.(2012).Ghost sph for animating water.ACMTransactions on Graphics(TOG)，31(4)，61.[0015][6] Akincij N.，Ihmsenj M.，Akincij G.，Solenthalerj B.，&Teschner, M.(2012).Versatile rigid-fluid coupling for incompressible SPH.ACM Transactions onGraphics (TOG) ,31 (4) ,62.
【發(fā)明內(nèi)容】

[0016]為了彌補現(xiàn)有耦合方法在虛擬手術(shù)仿真中的不足，并在虛擬手術(shù)仿真中實現(xiàn)精確、高效的多相耦合仿真，本發(fā)明提出一種虛擬手術(shù)中基于粒子的多相耦合方法。
[0017]本發(fā)明的技術(shù)方案為一種虛擬手術(shù)中基于粒子的多相耦合方法，包括以下步驟:
[0018]步驟1，對于虛擬手術(shù)中涉及的軟組織和剛體，分別根據(jù)三維組織表面三角網(wǎng)格模型生成粒子化三維組織體模型，包括以下子步驟，
[0019]①對三維組織表面三角網(wǎng)格模型建立最小長方體包圍盒；
[0020]②將步驟①中最小長方體包圍盒劃分為若干大小相同的立方體的體素，在邊界不足處采用立方體的體素補充，并以最小長方體包圍盒左下角頂點為起點對所有立方體的體素編號，得到各體素的體素序號；[0021]③在三維組織表面三角網(wǎng)格模型的每一個三角面片內(nèi)部均勻填充頂點，得到填充頂點后的三維組織表面三角網(wǎng)格模型，并為所有頂點編號得到頂點序號；
[0022]④對步驟③中填充頂點后的三維組織表面三角網(wǎng)格模型，確定每一個頂點所位于的體素，并建立體素-頂點關(guān)系表，該關(guān)系表的內(nèi)容是步驟②所得每個體素的體素序號及該體素所包含頂點的頂點序號；
[0023]⑤根據(jù)步驟④中的體素-頂點關(guān)系表，將所有包含頂點的體素定義為邊界體素；
[0024]⑥根據(jù)步驟⑤中所確定的邊界體素，對最小長方體包圍盒沿著O-XYZ坐標系的任一坐標軸方向逐層確定邊界體素所包圍的內(nèi)部體素，得到所有的內(nèi)部體素；
[0025]⑦將步驟③中填充頂點后的三維組織表面三角網(wǎng)格模型上每一個頂點定義為一個粒子，得到邊界粒子；并將步驟⑥中得到的每一個內(nèi)部體素取中心位置定義為一個粒子，得到內(nèi)部粒子；邊界粒子和內(nèi)部粒子構(gòu)成粒子化三維組織體模型；
[0026]步驟2，建立軟組織、液體和剛體的粒子化力學(xué)模型，將軟組織和液體定義為流變體，并對軟組織粒子和液體粒子的任意兩個粒子之間定義流變粘滯系數(shù)；
[0027]所述剛體的粒子化力學(xué)模型，只采用步驟I所得剛體的粒子化三維組織體模型的邊界粒子表示剛體，剛體作為一個整體移動和轉(zhuǎn)動，并計算剛體所受的合外力以及轉(zhuǎn)矩；
[0028]所述軟組織的粒子化力學(xué)模型，采用步驟I所得軟組織的粒子化三維軟組織體模型表示軟組織，根據(jù)彈性力學(xué)基本方程和SPH方法求解得到各個粒子的物理量，所述物理量包括彈力、粘滯力和重力，計算粘滯力時采用軟組織粒子之間的流變粘滯系數(shù)；
[0029]所述液體的粒子化力學(xué)模型，采用SPH方法建立，根據(jù)描述液體運動的Navier-Stokes方程組和SPH方法求解得到各個粒子的物理量，所述物理量包括粘滯力、壓力和重力，計算粘滯力時采用液體粒子之間的流變粘滯系數(shù)；
[0030]步驟3，進行基于邊界粒子和多相耦合核函數(shù)的碰撞檢測及響應(yīng)，對于任意兩相的物質(zhì)A和B，A為軟組織、液體或剛體，B為軟組織或剛體，且A和B不相同，
[0031]碰撞檢測過程如下，
[0032]①初始化A型粒子半徑和B型粒子半徑，
[0033]②對每一個A型粒子，以A型粒子當前位置為圓心，一個時間步長內(nèi)移動的距離為半徑作球形搜索區(qū)域，檢測出此球形搜索區(qū)域內(nèi)所有的B型邊界粒子；計算得到一個時間步長內(nèi)A型粒子由當前位置O移動到的位置O'，以O(shè)為起點指向O'作一條射線；
[0034]③將步驟②中所得的每個B型邊界粒子投影至步驟②中的射線所在的直線上，并求解每個B型邊界粒子到相應(yīng)投影點的投影距離；
[0035]④在步驟③中檢測出所有滿足投影距離小于等于A型粒子半徑與B型邊界粒子半徑之和的B型邊界粒；
[0036]⑤在步驟④檢測出的所有的B型邊界粒子中滿足條件的B型邊界粒子為碰撞點，所述條件為，該B型邊界粒子的投影點位于步驟②中A型粒子當前運動方向的前方，并其投影點距離步驟②中A型粒子圓心O最近；
[0037]碰撞響應(yīng)過程如下，
[0038]①根據(jù)碰撞檢測過程，判斷液體粒子、軟組織粒子和剛體粒子中任意兩相粒子是否發(fā)生碰撞，設(shè)該兩相粒子為A型粒子和B型粒子；
[0039]②若A型粒子和B型粒子發(fā)生碰撞，則定義A型粒子和B型粒子作用的耦合核函數(shù)；
[0040]③根據(jù)步驟②中的耦合核函數(shù)計算A型粒子和B型粒子之間的耦合作用力，并根據(jù)沖量定理求解下一時刻A型粒子和B型粒子的速度大小和方向；
[0041 ] ④根據(jù)步驟③的計算結(jié)果，代入A型粒子和B型粒子相應(yīng)的粒子化力學(xué)模型，計算出A型粒子和B型粒子新的位移、速度，而后繼續(xù)進行碰撞檢測過程。
[0042]而且，定義粒子化力學(xué)模型中兩個粒子i和j之間的流變粘滯系數(shù)?其中粒子i和j為軟組織粒子或液體粒子，
【權(quán)利要求】
1.一種虛擬手術(shù)中基于粒子的多相耦合方法，其特征在于，包括以下步驟: 步驟I，對于虛擬手術(shù)中涉及的軟組織和剛體，分別根據(jù)三維組織表面三角網(wǎng)格模型生成粒子化三維組織體模型，包括以下子步驟， ①對三維組織表面三角網(wǎng)格模型建立最小長方體包圍盒； ②將步驟①中最小長方體包圍盒劃分為若干大小相同的立方體的體素，在邊界不足處采用立方體的體素補充，并以最小長方體包圍盒左下角頂點為起點對所有立方體的體素編號，得到各體素的體素序號； ③在三維組織表面三角網(wǎng)格模型的每一個三角面片內(nèi)部均勻填充頂點，得到填充頂點后的三維組織表面三角網(wǎng)格模型，并為所有頂點編號得到頂點序號； ④對步驟③中填充頂點后的三維組織表面三角網(wǎng)格模型，確定每一個頂點所位于的體素，并建立體素-頂點關(guān)系表，該關(guān)系表的內(nèi)容是步驟②所得每個體素的體素序號及該體素所包含頂點的頂點序號； ⑤根據(jù)步驟④中的體素-頂點關(guān)系表，將所有包含頂點的體素定義為邊界體素； ⑥根據(jù)步驟⑤中所確定的邊界體素，對最小長方體包圍盒沿著O-XYZ坐標系的任一坐標軸方向逐層確定邊界體素所包圍的內(nèi)部體素，得到所有的內(nèi)部體素； ⑦將步驟③中填充頂點后的三維組織表面三角網(wǎng)格模型上每一個頂點定義為一個粒子，得到邊界粒子；并將步驟⑥中得到的每一個內(nèi)部體素取中心位置定義為一個粒子，得到內(nèi)部粒子；邊界粒子和內(nèi)部粒子構(gòu)成粒子化三維組織體模型」 步驟2，建立軟組織、液體和剛體的粒子化力學(xué)模型，將軟組織和液體定義為流變體，并對軟組織粒子和液體粒子的任意兩個粒子之間定義流變粘滯系數(shù)； 所述剛體的粒子化力學(xué)模型，只采用步驟I所得剛體的粒子化三維組織體模型的邊界粒子表示剛體，剛體作為一個整體移動和轉(zhuǎn)動，并計算剛體所受的合外力以及轉(zhuǎn)矩； 所述軟組織的粒子化力學(xué)模型，采用步驟I所得軟組織的粒子化三維軟組織體模型表示軟組織，根據(jù)彈性力學(xué)基本方程和SPH方法求解得到各個粒子的物理量，所述物理量包括彈力、粘滯力和重力I計算粘滯力時采用軟組織粒子之間的流變粘滯系數(shù)； 所述液體的粒子化力學(xué)模型，采用SPH方法建立,根據(jù)描述液體運動的Navier-Stokes方程組和SPH方法求解得到各個粒子的物理量，所述物理量包括粘滯力、壓力和和重力，計算粘滯力時采用液體粒子之間的流變粘滯系數(shù)； 步驟3，進行基于邊界粒子和多相耦合核函數(shù)的碰撞檢測及響應(yīng)，對于任意兩相的物質(zhì)A和B，A為軟組織、液體或剛體，B為軟組織或剛體，且A和B不相同， 碰撞檢測過程如下， ①初始化A型粒子半徑和B型粒子半徑， ②對每一個A型粒子，以A型粒子當前位置為圓心，一個時間步長內(nèi)移動的距離為半徑作球形搜索區(qū)域，檢測出此球形搜索區(qū)域內(nèi)所有的B型邊界粒子；計算得到一個時間步長內(nèi)A型粒子由當前位置O移動到的位置O'，以O(shè)為起點指向O'作一條射線； ③將步驟②中所得的每個B型邊界粒子投影至步驟②中的射線所在的直線上，并求解每個B型邊界粒子到相應(yīng)投影點的投影距離； ④在步驟③中檢測出所有滿足投影距離小于等于A型粒子半徑與B型邊界粒子半徑之和的B型邊界粒；⑤在步驟④檢測出的所有的B型邊界粒子中滿足條件的B型邊界粒子為碰撞點，所述條件為，該B型邊界粒子的投影點位于步驟②中A型粒子當前運動方向的前方，并其投影點距離步驟②中A型粒子圓心O最近； 碰撞響應(yīng)過程如下， ①根據(jù)碰撞檢測過程，判斷液體粒子、軟組織粒子和剛體粒子中任意兩相粒子是否發(fā)生碰撞，設(shè)該兩相粒子為A型粒子和B型粒子； ②若A型粒子和B型粒子發(fā)生碰撞，則定義A型粒子和B型粒子作用的耦合核函數(shù)； ③根據(jù)步驟②中的耦合核函數(shù)計算A型粒子和B型粒子之間的耦合作用力，并根據(jù)沖量定理求解下一時刻A型粒子和B型粒子的速度大小和方向； ④根據(jù)步驟③的計算結(jié)果，代入A型粒子和B型粒子相應(yīng)的粒子化力學(xué)模型，計算出A型粒子和B型粒子新的位移、速度，而后繼續(xù)進行碰撞檢測過程。
2.根據(jù)權(quán)利要求1所述虛擬手術(shù)中基于粒子的多相耦合方法，其特征在于:定義粒子化力學(xué)模型中兩個粒子i和j之間的流變粘滯系數(shù)其中粒子i和j為軟組織粒子或液體粒子，
3.根據(jù)權(quán)利要求2所述虛擬手術(shù)中基于粒子的多相耦合方法，其特征在于:Α型粒子與B型粒子之間作用的耦合核函數(shù)如下，
4.根據(jù)權(quán)利要求3所述虛擬手術(shù)中基于粒子的多相耦合方法，其特征在于:A型粒子和B型粒子的耦合作用力計算如下，
【文檔編號】G06T17/30GK103559741SQ201310602225
【公開日】2014年2月5日 申請日期:2013年11月25日 優(yōu)先權(quán)日:2013年11月25日
【發(fā)明者】袁志勇, 廖祥云, 郭甲翔, 喻思嬌 申請人:武漢大學(xué)