專利名稱：一種基于最小移動二乘的無網(wǎng)格物理形變仿真方法
技術(shù)領(lǐng)域：
本發(fā)明涉及虛擬現(xiàn)實領(lǐng)域，具體涉及一種基于最小移動二乘的無網(wǎng)格物理形變仿真方法。
背景技術(shù)：
虛擬現(xiàn)實領(lǐng)域?qū)μ摂M場景的物理仿真的準(zhǔn)確性與實時性要求越來越高，傳統(tǒng)的非物理方法，如自由形模型、半物理方法彈簧質(zhì)點模型，由于仿真精度不高，難以滿足一些對精度要求較高的應(yīng)用的需要，有限元等有網(wǎng)格算法雖然仿真精度高，但計算開銷大，尤其在拓?fù)浣Y(jié)構(gòu)改變時，物理模型的更新還難以實時準(zhǔn)確的完成。最近幾年開始出現(xiàn)的無網(wǎng)格物理形變方法，非常適用于拓?fù)浣Y(jié)構(gòu)可變的物體形變仿真，并可保證較高的仿真精度。該類方法的一般處理流程為(I)按照給定分辨率對仿真材料進(jìn)行粒子采樣，初始化粒子材料屬 性；(2)在每個物理仿真循環(huán)，通過無網(wǎng)格形函數(shù)近似出仿真材料的位移形變場，再通過應(yīng)變、應(yīng)力、內(nèi)力、外力步驟求解出每個物理粒子的受力；(3)采用顯式或隱式積分的方法求解牛頓第二定律PDE方程，計算出仿真材料在下一個物理仿真循環(huán)的位移。無網(wǎng)格物理變形算法起初應(yīng)用于流體仿真領(lǐng)域，后來逐漸應(yīng)用彈性體變形仿真領(lǐng)域。目前主流的無網(wǎng)格物理變形算法有光滑粒子流體動力學(xué)SPH (Smooth ParticlesHydrodynamics)、最小移動二乘 MLS (Moving Least Squares)等算法,其中的 SPH 方法實現(xiàn)簡單、效率高，但精度低，而MLS算法雖然計算開銷較大，但可取得更高的仿真精度，更加適合做彈性體的形變仿真。無網(wǎng)格計算方法之間的區(qū)別不一般主要體現(xiàn)在形函數(shù)的選擇上，在選取適當(dāng)形函數(shù)Φ (X)的基礎(chǔ)上，對于任意一個場函數(shù)f，可以按公式(I)計算出近似函數(shù)<f> </> (x) = Σ Φ (1 )
i如果采用SPH方法，則其形函數(shù)的形式一般為ΦΤΗ(X) = ViWh(||x- \ ||)Κ2)而如果使用的是MLS計算方法，則其形函數(shù)的形式一般可表示為Φ，(χ) = Wh (||χ - χ; ||)ρ(χ)ΓΜ(χ)_1 ρ(χ,)(3)本發(fā)明采用MLS方法作為物理仿真的解算模型，該模型最早由MUller等人引入到計算機(jī)圖形學(xué)，并逐漸被Pauly、Nealean、Gerszewski等人進(jìn)一步用于彈性體物理仿真模型的解算。MUller曾指出MLS在物理粒子處于退化狀態(tài)時會出現(xiàn)仿真不穩(wěn)定的問題，Guo等人將MLS與模態(tài)Warping結(jié)合,并通過轉(zhuǎn)換到頻域來完成粒子積分的方法,實現(xiàn)了一種更穩(wěn)定的伽遼金型無網(wǎng)格方法，但計算量比MUller提出的算法要大許多。Martin等人則通過引入彈性元(Elastons)的概念來解決MLS的穩(wěn)定性問題，但是該方法也只能滿足非實時的物理仿真需求。由于無網(wǎng)格算法的物理仿真過程與幾何模型無關(guān)，物理計算過程不依賴于幾何模型的頂點、邊聯(lián)通性關(guān)系，這不僅適合對基于體素的體模型進(jìn)行變形仿真，而且當(dāng)仿真對象發(fā)生切割、斷裂等拓?fù)浣Y(jié)構(gòu)改變時可保證以較小的代價來更新物理仿真模型，因而基于MLS的無網(wǎng)格方法也經(jīng)常用作撕裂、斷裂等復(fù)雜物理過程的仿真解算模型。

發(fā)明內(nèi)容
本發(fā)明要解決的技術(shù)問題是以有限元方法為代表的常規(guī)的有網(wǎng)格物理形變方法雖然仿真精度高，但仿真效率較低，難以滿足實時引用的需要；本發(fā)明從物理模型求解與物體表面網(wǎng)格形變驅(qū)動兩個方面進(jìn)行加速，給出了一種高效的最小移動二乘物理形變仿真方法，可將軟組織材料的物理屬性離散到物理粒子進(jìn)行基于MLS的近似計算，再通過物理粒子來實時模擬仿真軟組織在外力作用下的形變。本發(fā)明的技術(shù)方案為一種基于最小移動二乘的無網(wǎng)格物理形變仿真方法，其包括以下步驟步驟(I)、初始化階段按照給定的采樣分辨率對仿真材料幾何模型進(jìn)行物理粒 子采樣，并對采樣得到的粒子進(jìn)行材料屬性初始化；步驟(2)、在每個物理仿真循環(huán)中動態(tài)更新每個粒子的物理狀態(tài)首先通過最小移動二乘形函數(shù)計算仿真材料的物理形變場，再計算仿真材料的位移梯度，然后基于位移場求解應(yīng)變、應(yīng)力，再進(jìn)一步計算出仿真材料的彈性體力以及體積保持力；步驟(3)、根據(jù)牛頓第二定律構(gòu)建偏微分控制方程，采用隱式或顯式積分計算出仿真材料下一個循環(huán)的起始位移；步驟(4)、繪制仿真結(jié)果，并判斷物理仿真循環(huán)是否需要繼續(xù)進(jìn)行，如果是，回到步驟(2)繼續(xù)進(jìn)行物理仿真，如果否，則仿真結(jié)束。進(jìn)一步的，步驟(2)中所述的每個粒子的物理狀態(tài)更新計算，將具體計算過程分解為變量與常量兩個部分來進(jìn)行計算，其中常量部分包括位從移梯度、應(yīng)變張量、應(yīng)力張量、受力的變量中所分解出來的不變項，并可以對其進(jìn)行預(yù)計算；變量部分則包括位移梯度、應(yīng)變張量、應(yīng)力張量、受力的變量中的取值可動態(tài)變化的項，這些變量與定義的雅可比矩陣密切相關(guān)，在每個物理計算循環(huán)，只要計算出雅可比矩陣，便可結(jié)合常量部分高效的計算出物體的受力。進(jìn)一步的，步驟(2)中所述的位移梯度的計算，首先將位移梯度的計算分解為常量和基于形函數(shù)梯度的變量的計算，然后通過將位移梯度的無網(wǎng)格計算過程進(jìn)行分解，使得位移梯度與當(dāng)前法向量結(jié)合后的乘積為仿真過程中的不變量，因而只需計算一次，可有效避免大量的重復(fù)計算。本發(fā)明與現(xiàn)有技術(shù)相比的優(yōu)點在于(I)、與基于SPH的無網(wǎng)格物理仿真方法相比，盡管SPH方法計算高效、實現(xiàn)簡單，但仿真精度較低，且不穩(wěn)定，本發(fā)明的仿真結(jié)果更穩(wěn)定，可獲得更高的精仿真度，并且對包含同等規(guī)模仿真粒子的場景，本發(fā)明的計算效率仍然可以與SPH方法相媲美。(2)、與以有限元為代表的有網(wǎng)格仿真方法相比，對于材料連續(xù)性遭到破壞的撕裂、切割等操作，本發(fā)明由于不依賴物體的拓?fù)浣Y(jié)構(gòu)，因而可更高效、更靈活地處理網(wǎng)格拓?fù)浣Y(jié)構(gòu)改變的情況，可有效降低因拓?fù)浣Y(jié)構(gòu)改變而導(dǎo)致的物理模型更新所帶來的代價。


圖I為本發(fā)明一種基于最小移動二乘的無網(wǎng)格物理形變仿真方法的整體流程圖；圖2為物理粒子采樣示意圖，其中，(a)立方體，(b)立方體采樣結(jié)果I，(C)立方體采樣結(jié)果II，(d)肝臟表面網(wǎng)格，(e)肝臟采樣結(jié)果I，(f)肝臟采樣結(jié)果II ；圖3為鄰居粒子示意圖；圖4為雅可比矩陣和應(yīng)變張量可視化示意圖，其中，Ca)雅可比矩陣可視化顯式，
(b)應(yīng)變張量可視化表示；圖5為應(yīng)力張量和粒子合力可視化示意圖，其中，(a)應(yīng)力張量可視化表示，(b)粒子合力可視化表示；圖6為物理粒子重力作用下掉落撞擊地面過程示意圖； 圖7為表面網(wǎng)格在物理粒子驅(qū)動下發(fā)生的表面變形動畫示意圖，其中，(a)物理粒子撞擊地面變形，(b)網(wǎng)格表面發(fā)生對應(yīng)變形動畫。
具體實施例方式下面結(jié)合附圖以及具體實施方式
進(jìn)一步說明本發(fā)明。一種基于最小移動二乘的無網(wǎng)格物理形變仿真方法，包括以下步驟步驟(I)、初始化階段按照給定的采樣分辨率對仿真材料幾何模型進(jìn)行物理粒子采樣，并對采樣得到的粒子進(jìn)行材料屬性初始化；步驟(2)、在每個物理仿真循環(huán)中動態(tài)更新每個粒子的物理狀態(tài)。首先通過最小移動二乘形函數(shù)計算仿真材料的物理形變場，再計算仿真材料的位移梯度，然后基于位移場求解應(yīng)變、應(yīng)力等關(guān)鍵物理量，再進(jìn)一步計算出仿真材料的彈性體力、體積保持力；步驟(3)、根據(jù)牛頓第二定律構(gòu)建偏微分控制方程，采用隱式或顯式積分計算出仿真材料下一個循環(huán)的起始位移；步驟(4)、繪制仿真結(jié)果，并判斷物理仿真循環(huán)是否需要繼續(xù)進(jìn)行，如果是，回到步驟(2)繼續(xù)進(jìn)行物理仿真，如果否，則仿真結(jié)束；首先，在進(jìn)入物理仿真循環(huán)之前，需要創(chuàng)建物理仿真粒子，粒子按照設(shè)定的分辨率進(jìn)行采樣，由于要求在物體內(nèi)部與邊緣處進(jìn)行粒子采樣，因此需要一種方法判斷粒子與物體表面關(guān)系，這里采用深度場的方法建立空間與物體表面的關(guān)系，即將空間剖分為規(guī)則網(wǎng)格點，計算網(wǎng)格點到物體表面的距離，由網(wǎng)格點維持物體的隱式表面，網(wǎng)格點的值為到物體表面的最近距離，當(dāng)在表面外部時距離為正值，內(nèi)部時距離為負(fù)值，然后計算出物體包圍盒，從離原點最近的頂點開始采樣，采樣的距離由給定的采樣分辨率計算出，針對每個采樣點，發(fā)現(xiàn)該點深度值如果小于等于零，表示在物體內(nèi)部或表面，則在該初始位置創(chuàng)建物理粒子，如果為正值，則表示在物體表面外部，不創(chuàng)建該粒子。由于在每個物理仿真循環(huán)需要根據(jù)當(dāng)前位移依次計算出下一時刻的位移梯度、應(yīng)變張量、應(yīng)力張量、受力、位移等物理量，為了加速計算，本發(fā)明將其中的應(yīng)變張量、應(yīng)力張量和受力計算三個關(guān)鍵步驟壓縮為一個，具體計算過程由變量與常量兩個部分組成，其中常量部分包括從位移梯度、應(yīng)變張量、應(yīng)力張量、受力等變量中所分解出來的不變項，并可以對其進(jìn)行預(yù)計算；變量部分則包括位移梯度、應(yīng)變張量、應(yīng)力張量、受力等變量中的取值可動態(tài)變化的項，這些變量與定義的雅可比矩陣密切相關(guān)，在每個物理計算循環(huán)，只要計算出雅可比矩陣，便可結(jié)合常量部分高效的計算出物體的受力。
特別地，基于物理的形變需要按照物理模型的狀態(tài)變化，根據(jù)表面網(wǎng)格頂點的當(dāng)前位移與法向量來計算下一時刻的位移與法向量，其中法向量的更新需要計算位移梯度。為了提高法向量計算的效率，本發(fā)明將位移梯度的計算分解為常量和基于形函數(shù)梯度的變量的計算。通過將位移梯度的無網(wǎng)格計算公式展開，位移梯度與當(dāng)前法向量結(jié)合后的乘積則為仿真過程中的不變量，因而只需計算一次，可有效避免大量的重復(fù)計算。本發(fā)明從物理模型求解與物體表面網(wǎng)格形變驅(qū)動兩個方面進(jìn)行加速，給出了一種高效的最小移動二乘物理形變算法，方法的具體處理流程如圖I所示。(I)基于多項式基核函數(shù)的材料屬性初始化方法如圖2所示，算法進(jìn)行初始化，按照給定的采樣分辨率對仿真物體進(jìn)行物理粒子采樣，并對采樣得到的粒子進(jìn)行材料屬性初始化。需要初始化的物理粒子屬性包括質(zhì)量、密度、體積等，并按照單位體積比率對其進(jìn)行計算。同時，每類物理屬性都會按一定的規(guī)律進(jìn)行分布，本發(fā)明采用多項式基函數(shù)來刻畫這些物理屬性的空間分布情況。
(2)基于最小移動二乘的物理模型加速求解方法在每個物理仿真循環(huán)中都需要更新計算形變體的位移、位移梯度、應(yīng)變、應(yīng)力、合力，其中位移梯度的計算公式為Vu(x) = ^ VO,.(x)u,.、4；最小移動二乘形函數(shù)Φ (X)則按照公式(5)進(jìn)行定義,對于給定的η階多項式基ρ(χ) = [1χ...χη]τ，形函數(shù)Φ(χ)可表示為Φ (X) = Wi (X，Xi)pT(X) [M(x) ]-1p (Xi) (5)其中WiU, Xi)是加權(quán)函數(shù)，它定義了無網(wǎng)格形函數(shù)的支持域和物理粒子間的耦合程度，[M(X)IT1是矩量矩陣(Moment Matrix)逆矩陣,可由公示(6)計算得到Μ(χ) =Σ Wi (X，Xi) P (Xi) pT (Xi)(6)得到物理粒子位移梯度Vu(x)后，便可繼續(xù)計算應(yīng)變、應(yīng)力、合力等物理量。由于直接計算這些物理量會涉及大量中間量的重復(fù)計算，本發(fā)明將這三個物理量的計算步驟壓縮為一個，并將具體計算過程分解為變量求解和常量求解兩部分。常量的求解在仿真開始時通過預(yù)計算完成并進(jìn)行存儲，這樣在每個物理仿真循環(huán)中，只需計算以雅可比矩陣為中間量的相關(guān)變量，即可快速求解物體所受合力的大小。(3)物體表面網(wǎng)格形變的加速計算方法按照已經(jīng)建立的基于最小移動二乘的無網(wǎng)格物理形變方法，基于每個仿真循環(huán)計算得到的位移形變場，可通過公式(7)來更新體網(wǎng)格頂點的位移。u(x) = Σ Φ (X)Ui(7)而其中的法向量更新需要按照公式(4)重新計算位移梯度，這實質(zhì)上需要在仿真循中多次重復(fù)計算形函數(shù)的導(dǎo)數(shù)，由于形函數(shù)是確定的，因而可以對這部分內(nèi)容進(jìn)行預(yù)計算，并最終通過公式(8)來更新表面網(wǎng)格頂點的法向量n'(x) = f/(x) + Vu(x)7 f/(.r)(8)其中，基于多項式基核函數(shù)的材料屬性初始化具體為本發(fā)明假設(shè)每個物理粒子有固定的質(zhì)量，質(zhì)量在其周圍的分布服從多項式基核函數(shù)分布，如公式(9)所示
權(quán)利要求
1.一種基于最小移動二乘的無網(wǎng)格物理形變仿真方法，其特征在于包括以下步驟 步驟(I)、初始化階段按照給定的采樣分辨率對仿真材料幾何模型進(jìn)行物理粒子采樣，并對采樣得到的粒子進(jìn)行材料屬性初始化； 步驟(2)、在每個物理仿真循環(huán)中動態(tài)更新每個粒子的物理狀態(tài)首先通過最小移動二乘形函數(shù)計算仿真材料的物理形變場，再計算仿真材料的位移梯度，然后基于位移場求解應(yīng)變、應(yīng)力，再進(jìn)一步計算出仿真材料的彈性體力以及體積保持力； 步驟(3)、根據(jù)牛頓第二定律構(gòu)建偏微分控制方程，采用隱式或顯式積分計算出仿真材料下一個循環(huán)的起始位移； 步驟(4)、繪制仿真結(jié)果，并判斷物理仿真循環(huán)是否需要繼續(xù)進(jìn)行，如果是，回到步驟(2)繼續(xù)進(jìn)行物理仿真，如果否，則仿真結(jié)束。
2.根據(jù)權(quán)利要求I所述的基于最小移動二乘的無網(wǎng)格物理形變仿真方法，其特征在于步驟(2)中所述的每個粒子的物理狀態(tài)更新計算，將具體計算過程分解為變量與常量兩個部分來進(jìn)行計算，其中常量部分包括位從移梯度、應(yīng)變張量、應(yīng)力張量、受力的變量中所分解出來的不變項，可以對其進(jìn)行預(yù)計算；變量部分則包括位移梯度、應(yīng)變張量、應(yīng)力張量、受力的變量中的取值可動態(tài)變化的項，這些變量與定義的雅可比矩陣密切相關(guān)，在每個物理計算循環(huán)，只要計算出雅可比矩陣，便可結(jié)合常量部分高效的計算出物體的受力。
3.根據(jù)權(quán)利要求I所述的基于最小移動二乘的無網(wǎng)格物理形變仿真方法，其特征在于步驟(2)中所述的位移梯度的計算，首先將位移梯度的計算分解為常量和基于形函數(shù)梯度的變量的計算，然后通過將位移梯度的無網(wǎng)格計算過程進(jìn)行分解，使得位移梯度與當(dāng)前法向量結(jié)合后的乘積為仿真過程中的不變量，因而只需計算一次，可有效避免大量的重復(fù)計算。
全文摘要
本發(fā)明提出一種基于最小移動二乘的無網(wǎng)格物理形變仿真方法，該方法包括如下幾個步驟首先，算法按照給定的采樣分辨率對仿真材料幾何模型進(jìn)行粒子采樣，并對采樣得到的粒子材料屬性進(jìn)行初始化。然后，算法進(jìn)入物理仿真循環(huán)，通過將粒子受力計算分解為常量和以雅可比矩陣為基礎(chǔ)的變量兩部分，來進(jìn)行加速物理控制方程的求解。最后，在驅(qū)動物體表面網(wǎng)格形變時，通過將物體表面法向量計算分解為變量與以無網(wǎng)格形函數(shù)梯度為基礎(chǔ)的常量兩部分，來進(jìn)行高效率的仿真。本發(fā)明的仿真結(jié)果更穩(wěn)定，可獲得更高的精仿真度，并且對包含同等規(guī)模仿真粒子的場景，并且本發(fā)明的計算效率仍然可以與SPH方法相媲美。
文檔編號G06F17/50GK102831280SQ20121033364
公開日2012年12月19日 申請日期2012年9月10日 優(yōu)先權(quán)日2012年9月10日
發(fā)明者李帥, 黃震, 郝愛民, 楊麗鵬 申請人:北京航空航天大學(xué)