CN102831280A - 一种基于最小移动二乘的无网格物理形变仿真方法 - Google Patents

Abstract

本发明提出一种基于最小移动二乘的无网格物理形变仿真方法，该方法包括如下几个步骤：首先，算法按照给定的采样分辨率对仿真材料几何模型进行粒子采样，并对采样得到的粒子材料属性进行初始化。然后，算法进入物理仿真循环，通过将粒子受力计算分解为常量和以雅可比矩阵为基础的变量两部分，来进行加速物理控制方程的求解。最后，在驱动物体表面网格形变时，通过将物体表面法向量计算分解为变量与以无网格形函数梯度为基础的常量两部分，来进行高效率的仿真。本发明的仿真结果更稳定，可获得更高的精仿真度，并且对包含同等规模仿真粒子的场景，并且本发明的计算效率仍然可以与SPH方法相媲美。

Description

一种基于最小移动二乘的无网格物理形变仿真方法

技术领域

本发明涉及虚拟现实领域，具体涉及一种基于最小移动二乘的无网格物理形变仿真方法。

背景技术

虚拟现实领域对虚拟场景的物理仿真的准确性与实时性要求越来越高，传统的非物理方法，如自由形模型、半物理方法弹簧质点模型，由于仿真精度不高，难以满足一些对精度要求较高的应用的需要，有限元等有网格算法虽然仿真精度高，但计算开销大，尤其在拓扑结构改变时，物理模型的更新还难以实时准确的完成。最近几年开始出现的无网格物理形变方法，非常适用于拓扑结构可变的物体形变仿真，并可保证较高的仿真精度。该类方法的一般处理流程为：（1）按照给定分辨率对仿真材料进行粒子采样，初始化粒子材料属性；（2）在每个物理仿真循环，通过无网格形函数近似出仿真材料的位移形变场，再通过应变、应力、内力、外力步骤求解出每个物理粒子的受力；（3）采用显式或隐式积分的方法求解牛顿第二定律PDE方程，计算出仿真材料在下一个物理仿真循环的位移。

无网格物理变形算法起初应用于流体仿真领域，后来逐渐应用弹性体变形仿真领域。目前主流的无网格物理变形算法有光滑粒子流体动力学SPH（Smooth ParticlesHydrodynamics）、最小移动二乘MLS（Moving Least Squares）等算法，其中的SPH方法实现简单、效率高，但精度低，而MLS算法虽然计算开销较大，但可取得更高的仿真精度，更加适合做弹性体的形变仿真。

无网格计算方法之间的区别不一般主要体现在形函数的选择上，在选取适当形函数Ф(x)的基础上，对于任意一个场函数f，可以按公式（1）计算出近似函数<f>：

$<f>\left(x\right)=\underset{i}{\mathrm{\&Sigma;}}{\mathrm{\&Phi;}}_i\left(x\right)f_i---\left(1\right)$

如果采用SPH方法，则其形函数的形式一般为：

${\mathrm{\&Phi;}}_i^{\mathrm{SPH}}\left(x\right)=V_i{w}_h\left(\right||x-x_i|\left|\right)---\left(2\right)$

而如果使用的是MLS计算方法，则其形函数的形式一般可表示为：

${\mathrm{\&Phi;}}_i^{\mathrm{MLS}}\left(x\right)=w_h\left(\right||x-x_i|\left|\right)p{\left(x\right)}^{T}M{\left(x\right)}^{-1}p\left(x_i\right)---\left(3\right)$

本发明采用MLS方法作为物理仿真的解算模型，该模型最早由Müller等人引入到计算机图形学，并逐渐被Pauly、Nealean、Gerszewski等人进一步用于弹性体物理仿真模型的解算。Müller曾指出MLS在物理粒子处于退化状态时会出现仿真不稳定的问题，Guo等人将MLS与模态Warping结合，并通过转换到频域来完成粒子积分的方法，实现了一种更稳定的伽辽金型无网格方法，但计算量比Müller提出的算法要大许多。Martin等人则通过引入弹性元（Elastons）的概念来解决MLS的稳定性问题，但是该方法也只能满足非实时的物理仿真需求。

由于无网格算法的物理仿真过程与几何模型无关，物理计算过程不依赖于几何模型的顶点、边联通性关系，这不仅适合对基于体素的体模型进行变形仿真，而且当仿真对象发生切割、断裂等拓扑结构改变时可保证以较小的代价来更新物理仿真模型，因而基于MLS的无网格方法也经常用作撕裂、断裂等复杂物理过程的仿真解算模型。

发明内容

本发明要解决的技术问题是：以有限元方法为代表的常规的有网格物理形变方法虽然仿真精度高，但仿真效率较低，难以满足实时引用的需要；本发明从物理模型求解与物体表面网格形变驱动两个方面进行加速，给出了一种高效的最小移动二乘物理形变仿真方法，可将软组织材料的物理属性离散到物理粒子进行基于MLS的近似计算，再通过物理粒子来实时模拟仿真软组织在外力作用下的形变。

本发明的技术方案为：一种基于最小移动二乘的无网格物理形变仿真方法，其包括以下步骤：

步骤（1）、初始化阶段：按照给定的采样分辨率对仿真材料几何模型进行物理粒子采样，并对采样得到的粒子进行材料属性初始化；

步骤（2）、在每个物理仿真循环中动态更新每个粒子的物理状态：首先通过最小移动二乘形函数计算仿真材料的物理形变场，再计算仿真材料的位移梯度，然后基于位移场求解应变、应力，再进一步计算出仿真材料的弹性体力以及体积保持力；

步骤（3）、根据牛顿第二定律构建偏微分控制方程，采用隐式或显式积分计算出仿真材料下一个循环的起始位移；

步骤（4）、绘制仿真结果，并判断物理仿真循环是否需要继续进行，如果是，回到步骤（2）继续进行物理仿真，如果否，则仿真结束。

进一步的，步骤（2）中所述的每个粒子的物理状态更新计算，将具体计算过程分解为变量与常量两个部分来进行计算，其中常量部分包括位从移梯度、应变张量、应力张量、受力的变量中所分解出来的不变项，并可以对其进行预计算；变量部分则包括位移梯度、应变张量、应力张量、受力的变量中的取值可动态变化的项，这些变量与定义的雅可比矩阵密切相关，在每个物理计算循环，只要计算出雅可比矩阵，便可结合常量部分高效的计算出物体的受力。

进一步的，步骤（2）中所述的位移梯度的计算，首先将位移梯度的计算分解为常量和基于形函数梯度的变量的计算，然后通过将位移梯度的无网格计算过程进行分解，使得位移梯度与当前法向量结合后的乘积为仿真过程中的不变量，因而只需计算一次，可有效避免大量的重复计算。

本发明与现有技术相比的优点在于：

（1）、与基于SPH的无网格物理仿真方法相比，尽管SPH方法计算高效、实现简单，但仿真精度较低，且不稳定，本发明的仿真结果更稳定，可获得更高的精仿真度，并且对包含同等规模仿真粒子的场景，本发明的计算效率仍然可以与SPH方法相媲美。

（2）、与以有限元为代表的有网格仿真方法相比，对于材料连续性遭到破坏的撕裂、切割等操作，本发明由于不依赖物体的拓扑结构，因而可更高效、更灵活地处理网格拓扑结构改变的情况，可有效降低因拓扑结构改变而导致的物理模型更新所带来的代价。

附图说明

图1为本发明一种基于最小移动二乘的无网格物理形变仿真方法的整体流程图；

图2为物理粒子采样示意图，其中，(a)立方体，(b)立方体采样结果I，(c)立方体采样结果II，(d)肝脏表面网格，(e)肝脏采样结果I，(f)肝脏采样结果II；

图3为邻居粒子示意图；

图4为雅可比矩阵和应变张量可视化示意图，其中，（a）雅可比矩阵可视化显式，（b）应变张量可视化表示；

图5为应力张量和粒子合力可视化示意图，其中，（a）应力张量可视化表示，（b）粒子合力可视化表示；

图6为物理粒子重力作用下掉落撞击地面过程示意图；

图7为表面网格在物理粒子驱动下发生的表面变形动画示意图，其中，（a）物理粒子撞击地面变形，（b）网格表面发生对应变形动画。

具体实施方式

下面结合附图以及具体实施方式进一步说明本发明。

一种基于最小移动二乘的无网格物理形变仿真方法，包括以下步骤：

步骤（1）、初始化阶段按照给定的采样分辨率对仿真材料几何模型进行物理粒子采样，并对采样得到的粒子进行材料属性初始化；

步骤（2）、在每个物理仿真循环中动态更新每个粒子的物理状态。首先通过最小移动二乘形函数计算仿真材料的物理形变场，再计算仿真材料的位移梯度，然后基于位移场求解应变、应力等关键物理量，再进一步计算出仿真材料的弹性体力、体积保持力；

步骤（3）、根据牛顿第二定律构建偏微分控制方程，采用隐式或显式积分计算出仿真材料下一个循环的起始位移；

步骤（4）、绘制仿真结果，并判断物理仿真循环是否需要继续进行，如果是，回到步骤（2）继续进行物理仿真，如果否，则仿真结束；

首先，在进入物理仿真循环之前，需要创建物理仿真粒子，粒子按照设定的分辨率进行采样，由于要求在物体内部与边缘处进行粒子采样，因此需要一种方法判断粒子与物体表面关系，这里采用深度场的方法建立空间与物体表面的关系，即将空间剖分为规则网格点，计算网格点到物体表面的距离，由网格点维持物体的隐式表面，网格点的值为到物体表面的最近距离，当在表面外部时距离为正值，内部时距离为负值，然后计算出物体包围盒，从离原点最近的顶点开始采样，采样的距离由给定的采样分辨率计算出，针对每个采样点，发现该点深度值如果小于等于零，表示在物体内部或表面，则在该初始位置创建物理粒子，如果为正值，则表示在物体表面外部，不创建该粒子。

由于在每个物理仿真循环需要根据当前位移依次计算出下一时刻的位移梯度、应变张量、应力张量、受力、位移等物理量，为了加速计算，本发明将其中的应变张量、应力张量和受力计算三个关键步骤压缩为一个，具体计算过程由变量与常量两个部分组成，其中常量部分包括从位移梯度、应变张量、应力张量、受力等变量中所分解出来的不变项，并可以对其进行预计算；变量部分则包括位移梯度、应变张量、应力张量、受力等变量中的取值可动态变化的项，这些变量与定义的雅可比矩阵密切相关，在每个物理计算循环，只要计算出雅可比矩阵，便可结合常量部分高效的计算出物体的受力。

特别地，基于物理的形变需要按照物理模型的状态变化，根据表面网格顶点的当前位移与法向量来计算下一时刻的位移与法向量，其中法向量的更新需要计算位移梯度。为了提高法向量计算的效率，本发明将位移梯度的计算分解为常量和基于形函数梯度的变量的计算。通过将位移梯度的无网格计算公式展开，位移梯度与当前法向量结合后的乘积则为仿真过程中的不变量，因而只需计算一次，可有效避免大量的重复计算。

本发明从物理模型求解与物体表面网格形变驱动两个方面进行加速，给出了一种高效的最小移动二乘物理形变算法，方法的具体处理流程如图1所示。

（1）基于多项式基核函数的材料属性初始化方法

如图2所示，算法进行初始化，按照给定的采样分辨率对仿真物体进行物理粒子采样，并对采样得到的粒子进行材料属性初始化。需要初始化的物理粒子属性包括质量、密度、体积等，并按照单位体积比率对其进行计算。同时，每类物理属性都会按一定的规律进行分布，本发明采用多项式基函数来刻画这些物理属性的空间分布情况。

（2）基于最小移动二乘的物理模型加速求解方法

在每个物理仿真循环中都需要更新计算形变体的位移、位移梯度、应变、应力、合力，其中位移梯度的计算公式为：

$\&dtri;u\left(x\right)={\mathrm{\&Sigma;}}_i\&dtri;{\mathrm{\&Phi;}}_i\left(x\right)u_i---\left(4\right)$

最小移动二乘形函数Ф(x)则按照公式（5）进行定义，对于给定的n阶多项式基p(x)=[1x...xⁿ]^T，形函数Ф(x)可表示为：

Ф_i(x)＝w_i(x,x_i)p^T(x)[M(x)]^-1p(x_i)（5）

其中w_i(x,x_i)是加权函数，它定义了无网格形函数的支持域和物理粒子间的耦合程度，[M(x)]^-1是矩量矩阵（Moment Matrix）逆矩阵，可由公示（6）计算得到：

M(x)＝∑_iw_i(x,x_i)p(x_i)p^T(x_i)                    （6）

得到物理粒子位移梯度后，便可继续计算应变、应力、合力等物理量。由于直接计算这些物理量会涉及大量中间量的重复计算，本发明将这三个物理量的计算步骤压缩为一个，并将具体计算过程分解为变量求解和常量求解两部分。常量的求解在仿真开始时通过预计算完成并进行存储，这样在每个物理仿真循环中，只需计算以雅可比矩阵为中间量的相关变量，即可快速求解物体所受合力的大小。

（3）物体表面网格形变的加速计算方法

按照已经建立的基于最小移动二乘的无网格物理形变方法，基于每个仿真循环计算得到的位移形变场，可通过公式（7）来更新体网格顶点的位移。

u(x)＝∑_iФ_i(x)u_i                            （7）

而其中的法向量更新需要按照公式（4）重新计算位移梯度，这实质上需要在仿真循中多次重复计算形函数的导数，由于形函数是确定的，因而可以对这部分内容进行预计算，并最终通过公式（8）来更新表面网格顶点的法向量：

$n^{\&prime;}\left(x\right)=n\left(x\right)+\&dtri;u{\left(x\right)}^{T}n\left(x\right)---\left(8\right)$

其中，基于多项式基核函数的材料属性初始化具体为：

本发明假设每个物理粒子有固定的质量，质量在其周围的分布服从多项式基核函数分布，如公式（9）所示：

$W(r,h)=\left\{\begin{array}{cc}\frac{315}{64\mathrm{\&pi;}{h}^9}{(h^2-r^2)}^3& \mathrm{ifr}<h\\ 0& \mathrm{ifr}>=h\end{array}---\left(9\right)\right.$

其中r是到物理粒子的距离，h是核函数的支持域，核函数可以归一化，并满足∫_xW(|x-x₀|,h)dx＝1，这样粒子i处密度可以按公式（10）计算：

ρ_i＝∑_jm_jw_ij                                       （10）

其中w_ij＝W(|x_j-x_i|,h_i)，物理粒子i表示的体积则为：

v_i＝m_i/ρ_i                                          （11）

物体的质量是固定的，如果为塑性变形，则密度与体积则会随着相对于粒子距离的变化而变化，质量m_i与支持域h_i需要在仿真开始前按如下方式进行初始化：对每个物理粒子，先计算其最近k个物理粒子的平均距离

k由用户设置，支持域h_i可设置为平均距离的倍数，如图3所示，质量可被初始化为

这里ρ是材料密度，s是所有物理粒子的缩放尺度。

其中，物理模型的加速求解具体如下：

本发明采用更为紧凑的粒子合力f求解公式：

$\left\{\begin{array}{c}f=f_i+f_b\\ f_i=f_e+f_v=(F_e+F_v)d_i\\ f_e=F_e{d}_i,f_v=F_v{d}_i\end{array}---\left(12\right)\right.$

其中，f_i为粒子内力，由弹性力f_e与体积保持力f_v组成，f_b为外力，如拉扯、重力、撞击力等，其中的主要计算开销是对粒子内力f_i的计算，它被分解为常量与变量两个部分来进行加速求解。

常量部分d_i：如公式（13）所示，常量部分只需要在物理仿真初始化过程中计算即可，而在物体材料连续性破坏后，只需花费极小的代价来进行重新更新：

$\left\{\begin{array}{c}{d}_i=A^{-1}(-\underset{j}{\mathrm{\&Sigma;}}{x}_{\mathrm{ij}}{w}_{\mathrm{ij}})\\ d_j=A^{-1}\left(x_{\mathrm{ij}}{w}_{\mathrm{ij}}\right)\\ A=\underset{j}{\mathrm{\&Sigma;}}{x}_{\mathrm{ij}}{x}_{\mathrm{ij}}^T{w}_{\mathrm{ij}}\end{array}\right.---\left(13\right)$

其中x_ij为粒子i与粒子j之间材料距离差，w_ij为粒子i与粒子j之间距离差的加权函数值。

变量部分(F_e+F_v)：是在每个仿真循环中都需要重新计算的部分，包括物体弹性力的可变部分F_e和体积保持力的可变部分F_v，其计算按照公式（14）所示：

$\left\{\begin{array}{c}{F}_e=-2v_i{J}_i{\mathrm{\&sigma;}}_i,{\mathrm{\&epsiv;}}_i=J_i^T{J}_i-I,{\mathrm{\&sigma;}}_i={\mathrm{C\&epsiv;}}_i\\ F_v=-v_i{k}_v\left(\right|J|-1)\left[\begin{array}{c}{(J_v\&times;J_w)}^T\\ {(J_w\&times;J_u)}^T\\ {(J_u\&times;J_v)}^T\end{array}\right]\\ J=[J_x,J_y,J_z]=\left[\begin{array}{c}{J}_u^T\\ J_v^T\\ J_w^T\end{array}\right]\\ J_i=\left[\begin{array}{c}\&dtri;u{|}_{x_i}^T\\ \&dtri;u{|}_{x_i}^T\\ \&dtri;u{|}_{x_i}^T\end{array}\right]+I\end{array}\right.---\left(14\right)$

其中，σ_i、ε_i、v_i、J_i、k_v分别代表粒子i的应力张量、应变张量、体积、雅可比矩阵和材料体积保持系数，，C是四阶张量，由杨氏模量与泊松比构成。而图4和图5通过一个实例给出了计算结果的可视化示意图。

计算物体合力后，便可通过隐式或显式积分方法求解偏微分动力学控制方程，进而通过公式（15）来近似计算每个物理粒子在t+1时刻的位移u：

$\mathrm{\&rho;}\frac{{\&PartialD;}^2{x}^{\&prime;}}{\&PartialD;t^2}=\mathrm{\&rho;}\frac{{\&PartialD;}^{2}u}{\&PartialD;t^2}=f/m---\left(15\right)$

其中，ρ、f、m分别为材料密度、合力与质量。图6以粒子单元的形式给出了物体在重力作用下掉落并撞击地面发生形变过程示意图。

进一步的，物体表面网格形变的加速计算具体为：

法向量是表面网格模型绘制所需要的关键信息，只有正确计算法向量才能保证表面网格在光照模型下的正确着色。每个仿真循环中需要更新的顶点法向量n'可通过局部位移场的梯度来近似计算，计算公式如（8）所示。

由于针对每个时间步长都需要计算位移场梯度而

的计算需要用到公式（4），因此，在每次迭代顶点邻域内的物理单元时要进行大量矩阵预算，非常耗时。本发明则将

的计算展开，如公式（16）所示，由于

是个常量，因此可以将形函数场梯度

与初始法向量n(x)这两项合并，通过一次预计算得到，从而可以减少大量的数据读取操作和矩阵计算，可极大提高物体表面网格形变的计算效率，而且改进后的计算方法也更加便于并行加速。图7给出了表面网格在物理粒子驱动下所发生的变形结果。

$\left\{\begin{array}{c}{n}^{\&prime;}\left(x\right)=n\left(x\right)+{\mathrm{\&Sigma;}}_i\left(\&dtri;{\mathrm{\&Phi;}}_i^T\left(x\right)u_i\right)n\left(x\right)\\ n^{\&prime;}\left(x\right)=n\left(x\right)+{\mathrm{\&Sigma;}}_i(\&dtri;{\mathrm{\&Phi;}}_i^T\left(x\right)n\left(x\right))u_i\end{array}---\left(16\right)\right.$

表1则统计了其它同类无网格物理算法的性能数据，并将其与本发明所实现的MLS最小移动二乘形变物理方法进行了对比，其中Zhu方法（Zhu B.,Gu L.,Peng X.,et al.APoint-Based Simulation Framework for Minimally Invasive Surgery[A].in Proc.ISMBS[C],2010:130-138）与Qin方法（Qin J.,Chui Y.,Pang W.,et al.Learning Blood Management inOrthopedic Surgery through Gameplay[J],Presented at IEEE Computer Graphics and Applications,2010:45-57）均为基于SPH的无网格物理仿真方法。从表1可以看出，SPH虽然计算效率较高，但仿真精度较低，且不稳定，而本发明所给出的方法更稳定，并可获得更高阶的精度，而且对包含1000个左右仿真粒子的场景，本发明所设计方法可取得与SPH方法处于同一个数量级的计算效率。

表1本发明方法与其他无网格物理仿真方法性能对比

方法出处	物理粒子数	计算时间
			Zhu方法	1000	5.31ms
Qin方法	1000	15.70ms
			本发明方法	1000	6.25ms

本发明未详细阐述的技术内容属于本领域技术人员的公知技术。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (3)

1.一种基于最小移动二乘的无网格物理形变仿真方法，其特征在于包括以下步骤：

步骤（1）、初始化阶段：按照给定的采样分辨率对仿真材料几何模型进行物理粒子采样，并对采样得到的粒子进行材料属性初始化；

步骤（2）、在每个物理仿真循环中动态更新每个粒子的物理状态：首先通过最小移动二乘形函数计算仿真材料的物理形变场，再计算仿真材料的位移梯度，然后基于位移场求解应变、应力，再进一步计算出仿真材料的弹性体力以及体积保持力；

步骤（3）、根据牛顿第二定律构建偏微分控制方程，采用隐式或显式积分计算出仿真材料下一个循环的起始位移；

步骤（4）、绘制仿真结果，并判断物理仿真循环是否需要继续进行，如果是，回到步骤（2）继续进行物理仿真，如果否，则仿真结束。

2.根据权利要求1所述的基于最小移动二乘的无网格物理形变仿真方法，其特征在于：步骤（2）中所述的每个粒子的物理状态更新计算，将具体计算过程分解为变量与常量两个部分来进行计算，其中常量部分包括位从移梯度、应变张量、应力张量、受力的变量中所分解出来的不变项，可以对其进行预计算；变量部分则包括位移梯度、应变张量、应力张量、受力的变量中的取值可动态变化的项，这些变量与定义的雅可比矩阵密切相关，在每个物理计算循环，只要计算出雅可比矩阵，便可结合常量部分高效的计算出物体的受力。

3.根据权利要求1所述的基于最小移动二乘的无网格物理形变仿真方法，其特征在于：步骤（2）中所述的位移梯度的计算，首先将位移梯度的计算分解为常量和基于形函数梯度的变量的计算，然后通过将位移梯度的无网格计算过程进行分解，使得位移梯度与当前法向量结合后的乘积为仿真过程中的不变量，因而只需计算一次，可有效避免大量的重复计算。

Images