一种软组织与刚性地面的碰撞仿真方法及装置

1.本发明涉及计算机图形学技术领域，特别涉及一种软组织与刚性地面的碰撞仿真方法及装置。


背景技术：

2.在计算机图形学中，仿真固体通常可分为刚体、柔体和衣物三类。肝脏等人体软组织属于柔体一类，与刚体不同，柔体在受力后形状会产生变化并且内部各点的之间相对位置也会随之改变；而衣物通常采用二维建模(忽略厚度)，柔体则必须进行三维建模。因此，柔体仿真属于一类较为复杂的固体仿真。传统的newton(无摩擦)方法未考虑软组织与刚性地面碰撞时产生的摩擦力，因而会造成失真。
3.具体而言，传统的方法包括：
4.1.1弹簧-质点模型
5.弹簧-质点模型是最简单的形变模型，它是一种基于弹簧力的物理建模方法。在弹簧-质点系统中，仿真物体被离散为一系列的质点，各个质点之间经由弹簧网络连接，通过调节弹簧的系数可以控制物体的形变行为。
6.该模型简单明了，易于实现，在绝大多数仿真场景中能够取得较高的计算效率。但是其简单性也随之带来了一些缺点：
7.1、物体的形变行为取决于弹簧网络的规模和建立方式。
8.2、弹簧-质点模型不能维持一些必要的体积特性，例如体积守恒，这导致了该系统不适用于固体的形变仿真。弹簧-质点系统更多地被应用于衣物仿真，因为在计算机图形学中，通常将衣物视为没有厚度的物体，即不存在体积。
9.3、由于缺乏必要的理论支撑，弹簧系数难以设置。对于给定的仿真材料，其弹簧系数只能通过不断的实验尝试来达到较接近的匹配结果。
10.1.2有限元模型
11.fem是一种物理精确的形变建模方法。在数值计算中，fem是一种近似求解偏微分方程边值问题的方法，在结构力学分析中被广泛使用。
12.fem基于连续介质假设，即假设物质连续地分布于其所占有的整个空间。为了进行系统求解，需要将模型划分为一系列有限个数的子区域，这些子区域称为单元，单元与单元的连接点称为节点。fem计算过程先对每个单元建立方程，方程由单元应力、应变和位移间的关系确立；然后将所有单元对应的方程集合起来，形成整体系统方程组，系统方程组是一组以节点位移为未知量的线性方程；最后对系统方程组进行求解，得到各节点的位移。
13.fem直接将传统力学的理论分析方法应用到计算机图形学中，因而能取得理想的物理准确性。但是，精确的仿真结果往往是以巨大的计算消耗为代价。对于一个节点个数为n的模型，系统在每一仿真步(simulation step)都需要求解一个3n
×
3n规模的稀疏矩阵线性方程，这对仿真系统而言无疑是巨大的计算负担，因此fem不适用于大型的交互式应用。
14.因而，现有需要一种在软组织与刚性地面碰撞的交互式应用场景下能够兼顾仿真
的精确性和计算效率的方法。


技术实现要素：

15.为了解决现有技术的上述问题，本发明提供一种软组织与刚性地面的碰撞仿真方法及装置，以在软组织与刚性地面碰撞的交互式应用场景下能够实兼顾仿真的精确性和计算效率。
16.为了达到上述目的，本发明采用的技术方案为：
17.第一方面，本发明提供一种软组织与刚性地面的碰撞仿真方法，包括：
18.步骤s1、获取仿真建模后的软组织模型；
19.步骤s2、在每一个仿真步中，计算并得到所述软组织模型上每一个顶点在下一时刻的预测位置，通过以所述刚性地面作为外部约束和以柔性物体的应变能约束作为内部约束来修正所述软组织模型上每一个顶点的预测位置，并根据每一个顶点的修正后位置进行下一步的模型渲染和形变计算；
20.步骤s3、重复执行下一仿真步的计算，直至完成仿真。
21.本发明的有益效果在于：本发明以位置动力学作为基础形变模型，采用应变能约束作为内部约束，采用环境碰撞约束作为外部约束，二者构成一个系统约束方程组，进行方程组的求解得到模型顶点的新位置，从而实现形变仿真，其仿真方法稳定、高效，并且有较好的物理真实感，以在软组织与刚性地面碰撞的交互式虚拟三维场景下能够实兼顾仿真的精确性和计算效率。
22.可选地，所述步骤s2包括：
23.步骤s21、在第一个仿真步时，对所述软组织模型上的每一个顶点的基本属性进行初始化：，其中每一个顶点的基本属性包括质量mi、位置xi、速度vi；
24.步骤s22、在每一个仿真步中，以固定的时间步长δt进行时间积分计算得到所述软组织模型上每一个顶点的预测速度{v0,v1,...vn}和预测位置{p0,p1,...pn}，所述预测位置作为约束求解的初始条件；
25.步骤s23、根据当前的仿真场景，以所述刚性地面的外部约束和以柔性物体的应变能约束作为内部约束来共同构成系统约束方程组；
26.步骤s24、将每一个顶点的预测位置{p0,p1,...pn}代入系统约束方程组，使用gauss-seidel迭代法进行求解，并限制迭代次数，求解得到一组位置修正量{δp0,δp1,...δpn}；
27.步骤s25、使用所述位置修正量{δp0,δp1,...δpn}修正每一个顶点的预测位置{p0,p1,...pn}，用修正后的预测位置{x0,x1,...xn}对顶点位置属性进行更新并更新所述预测速度，之后根据每一个顶点的修正后位置进行下一步的模型渲染和形变计算。
28.根据上述描述可知，相比于传统的方法，上述的仿真过程能得到更准确的结果，同时使得约束投影更快收敛，兼顾了仿真的精确性和计算效率。
29.可选地，所述步骤s23中以所述刚性地面的外部约束包括：
30.将所述软组织模型与所述刚性地面之间的碰撞约束视为距离约束，其约束函数为：
31.c(p)＝(p-qs)
·ns
≥0；
32.其中，p为所述软组织模型上的任一顶点，qs表示在碰撞接触面上离p最近的点，ns为qs处的法向量。
33.根据上述描述可知，将软组织模型与刚性地面之间的碰撞约束视为距离约束，并限定两者的距离要大于或等于0，从而在仿真物体受重力作用而下落至地面，能防止发生仿真物体出现穿透刚性地面的现象。
34.可选地，所述步骤s23中以柔性物体的应变能约束作为内部约束包括：
35.将所述软组织模型的形变示为一个连续的位移场u，定义ω为所述软组织模型所占据的体积域，x∈ω为所述软组织模型未形变状态下的任意一点，则可得到物体变形函数φ(x)的表达式：
36.φ(x)＝x u＝x；
37.其中，变形函数描述了材料空间中的点x与变形后对应的点x之间的映射关系，所述变形函数的变形梯度张量f为：
[0038][0039]
根据变形梯度张量f可求得应变张量e为：
[0040][0041]
其中i表示单位矩阵；
[0042]
根据变形梯度张量f和应变张量e求得应变能密度ψs为：
[0043][0044]
其中，e:e＝tr(e
t
e)，函数tr(a)表示矩阵a的迹，μ、λ为拉梅系数，只与两个材料属性杨氏模量k和泊松比v有关：
[0045][0046]
由于采用四面体单元对物体体积域进行离散，单元形函数采用线性拉格朗日形函数，于是每一单元对应的变形梯度进一步表示为：
[0047][0048]
其中ds为形函数矩阵，dm为参考形函数矩阵，是一个常数矩阵，ds、dm的向量表达式分别为：
[0049]ds
＝(x
0-x
3 x
1-x
3 x
2-x3)；
[0050]dm
＝(x
0-x
3 x
1-x
3 x
2-x3)；
[0051]
对于所述软组织模型的任意一个四面体单元，定义单元体积域为则储存在一个四面体单元内的应变能约束函数为：
[0052][0053]
其中，v是四面体单元未形变时的单元体积；
[0054]
根据变形梯度f得到所述软组织模型的应力张量p(f)：
[0055][0056]
则所述软组织模型的应力张量p(f)进一步化简为：
[0057]
p(f)＝f(2μe λtr(e)i)；
[0058]
根据应力张量和变形梯度的关系式得到应变能梯度为：
[0059][0060]
对于所述软组织模型的四面体单元，只需求解4个顶点对应的应变能量梯度即可，定义为i＝{1,2,3,4}，则得到所述软组织模型的每一个四面体应变能梯度的最终表达式为：
[0061][0062]
根据上述描述可知，通过提供一种适用于柔体形变仿真的应变能约束方法，该方法结合了位置动力学的计算效率优势以及传统连续介质力学成熟的理论成果，使仿真稳定、高效并且可控。
[0063]
可选地，所述杨氏模量k为[0.3,0.5]，所述泊松比v为[0.3,0.5]。
[0064]
可选地，所述步骤s24包括：
[0065]
将每一个顶点的预测位置{p0,p1,...pn}代入系统约束方程组；
[0066]
在每一次迭代时，为每个约束函数cj(p)找到一个合适的位置修正量δp，使之满足cj(p δp)＝0，则对每一个约束方程单独进行线性化为：
[0067][0068]
其中，表示为约束函数cj(p)的梯度算子求解系统将δp限制在的方向进行求解，将质点的个体质量考虑进来，可得到所述软组织模型的每一个顶点i的位置修正量δpi的推导式为:
[0069][0070][0071]
式中，s表示拉格朗日乘子，k∈[0,1]表示刚度参数，表示矢量c(p)在pi方向上的梯度算子；
[0072]
完成限制的迭代次数后，得到一组位置修正量{δp0,δp1,...δpn}。
[0073]
可选地，所述步骤s25中更新所述预测速度为通过二阶向后微分公式更新所述预测速度。
[0074]
根据上述描述可知，采用一种修改版本的迭代法，在每一次迭代时对这n个线性方
程拆分，每个方程独立线性化、独立求解，从而显著提高计算效率。
[0075]
可选地，所述步骤s1包括：
[0076]
采用圣维南-基尔霍夫模型对软组织进行建模，得到仿真建模后的软组织模型。
[0077]
可选地，所述迭代次数为[3,8]，所述时间步长为[0.001,0.01]。
[0078]
根据上述描述可知，基于位置动力学的形变仿真效果还取决于时间步长和迭代次数，本发明限定了迭代次数和时间步长的范围，保证了软组织与刚性地面的碰撞仿真效果。
[0079]
第二方面，本发明提供一种软组织与刚性地面的碰撞仿真装置，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现第一方面的一种软组织与刚性地面的碰撞仿真方法。
[0080]
其中，第二方面所提供的一种软组织与刚性地面的碰撞仿真装置所对应的技术效果参照第一方面所提供的一种软组织与刚性地面的碰撞仿真方法的相关描述。
附图说明
[0081]
图1为本发明实施例的一种软组织与刚性地面的碰撞仿真方法的主要流程示意图；
[0082]
图2为本发明实施例涉及的形变映射的示意图；
[0083]
图3为本发明实施例涉及的环境碰撞约束中距离判断的示意图；
[0084]
图4为本发明实施例涉及的肝脏模型1的建模示意图；
[0085]
图5为本发明实施例涉及的肝脏模型2的建模示意图；
[0086]
图6为本发明实施例涉及的肝脏模型3的建模示意图；
[0087]
图7为本发明实施例涉及的肝脏模型2的仿真示意图；
[0088]
图8为本发明实施例涉及的肝脏模型3的仿真示意图；
[0089]
图9为现有技术中的肝脏模型2的仿真示意图；
[0090]
图10为本发明实施例的一种软组织与刚性地面的碰撞仿真装置的结构示意图。
[0091]
【附图标记说明】
[0092]
1：一种软组织与刚性地面的碰撞仿真装置；
[0093]
2：处理器；
[0094]
3：存储器。
具体实施方式
[0095]
为了更好的理解上述技术方案，下面将参照附图更详细地描述本发明的示例性实施例。虽然附图中显示了本发明的示例性实施例，然而应当理解，可以以各种形式实现本发明而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了能够更清楚、透彻地理解本发明，并且能够将本发明的范围完整的传达给本领域的技术人员。
[0096]
实施例一
[0097]
请参照图1至图8，一种软组织与刚性地面的碰撞仿真方法，包括步骤：
[0098]
步骤s1、获取仿真建模后的软组织模型；
[0099]
在本实施例中，步骤s1包括：
[0100]
采用圣维南-基尔霍夫模型对软组织进行建模，得到仿真建模后的软组织模型。
[0101]
其中，本实施例的软组织模型为肝脏模型，弹性形变将使得势能在形变物体内积累，对于固体而言，该势能也称为应变能(strain energy)，存在弹性势能函数的材料被为超弹性(hyperelastic)材料，常见的超弹性材料包括橡胶、海绵等。考虑到现实医学实验常用橡胶作为人体模型的制作材料，因此本发明将肝脏视为超弹性材料处理，并采用圣维南-基尔霍夫(saint venant-kirchhoff，stvk)模型进行建模。
[0102]
其中，本实施例中分别使用了三个不同精度的肝脏模型进行剖分。三个模型剖分前后的相关数据如表1所示，剖分效果依次如图4至图6所示。
[0103]
表1肝脏模型四面体剖分数据
[0104][0105]
步骤s2、在每一个仿真步中，计算并得到软组织模型上每一个顶点在下一时刻的预测位置，通过以刚性地面作为外部约束和以柔性物体的应变能约束作为内部约束来修正软组织模型上每一个顶点的预测位置，并根据每一个顶点的修正后位置进行下一步的模型渲染和形变计算；
[0106]
其中，位置动力学适用于液体、刚体、柔体以及衣物仿真，将不同类型的物体的形变计算在统一的仿真框架下进行，这样一来可以方便地在同一场景中实现多类物体间的交互仿真。
[0107]
位置动力学在计算流程上有别于弹簧-质点模型与有限元模型，弹簧-质点模型与有限元模型在每一仿真步中都需要计算出每个节点的内力和外力，通过合力计算出新的节点速度，进而根据速度更新位置。而位置动力学模型省去了速度层面的计算，直接处理位置信息，其基本思想是先预估节点的下一时刻的位置，然后通过约束方程对位置进行修正，把它投影到最终的位置上，同时用这个位置去更新速度。因此，位置动力学具有明显的计算优势，特别适用于交互式应用环境。
[0108]
在物理真实性方面，位置动力学虽然不能达到像有限元模型那样足够精确的物理仿真结果，但通过建立基于连续介质力学的约束方程也能够获得相对真实的结果，足够保证视觉上的合理性。相比于弹簧-质点模型，位置动力学有处理大形变的能力，即使在显式时间积分场景下也能保证系统稳定性，并且具有可控性。
[0109]
综上所述，位置动力学是一种快速、稳定、可控的形变模型，虽然在力学真实感上有所欠缺，但考虑到本实施例研究尚未涉及对肝脏生物力学的分析，主要涉及的是肝脏与刚性地面的碰撞仿真，并且交互式环境为应用背景，因此本发明选定位置动力学作为基础形变模型。
[0110]
其中，位置动力学通过模型的顶点集和一组约束方程来对仿真物体进行数学建模，通过顶点的位移实现物体的形变，而通过约束方程限制顶点间的位置关系，从而规范了物体的形变行为。在每一仿真步中，系统都需要对约束方程组进行求解，并根据求解结果将
顶点投影到新的位置上。
[0111]
在位置动力学计算框架中，每个顶点i包含质量mi、位置xi、速度vi三个基本属性。顶点的运动方程可由牛顿第二定律得到，如式1所示。
[0112][0113]
其中，表示加速度，fi表示作用在顶点i上的所有合力。
[0114]
由此，步骤s2包括：
[0115]
步骤s21、在第一个仿真步时，对软组织模型上的每一个顶点的基本属性进行初始化：其中每一个顶点的基本属性包括质量mi、位置xi、速度vi；
[0116]
步骤s22、在每一个仿真步中，以固定的时间步长δt进行时间积分计算得到软组织模型上每一个顶点的预测速度{v0,v1,...vn}和预测位置{p0,p1,...pn}，预测位置作为约束求解的初始条件；
[0117]
其中，计算每一个顶点的预测速度{v0,v1,...vn}和预测位置{p0,p1,...pn}如式2和式3所示：
[0118]
vi(tn δt)＝vi(tn) δtw
ifi
(tn)
ꢀꢀꢀ
(2)
[0119]
pi(tn δt)＝xi(tn) δtvi(tn δt)
ꢀꢀꢀ
(3)
[0120]
步骤s23、根据当前的仿真场景，以刚性地面的外部约束和以柔性物体的应变能约束作为内部约束来共同构成系统约束方程组；
[0121]
其中，位置动力学还有一个优势在于对碰撞的处理十分方便，同样可以通过约束方程的形式实现。碰撞约束为外部约束，用于定义模型与其他物体之间的位置约束关系。碰撞约束方程在执行系统求解之前，根据顶点的当前位置生成。碰撞约束方程与应变能约束方程直接集成，共同构成系统方程组，并一起进行求解。
[0122]
考虑简单的环境碰撞，即动态的仿真物体与其他静态固体之间的碰撞，例如仿真物体在重力作用下下落，为防止发生仿真物体出现穿透墙面的现象，需要将顶点的投影位置限定在一个有效的范围内。其中，静态固体包括地板、墙壁等。
[0123]
由此，步骤s23中以刚性地面的外部约束包括：
[0124]
将软组织模型与刚性地面之间的碰撞约束视为距离约束，其约束函数为：
[0125]
c(p)＝(p-qs)
·ns
ꢀꢀꢀ
(4)
[0126]
其中，p为软组织模型上的任一顶点，qs表示在碰撞接触面上离p最近的点，ns为qs处的法向量。
[0127]
如图3所示，物体受重力作用逐渐下落至地面，在这一仿真场景中需保证式(4)不小于0，即碰撞约束表达式为：c(p)≥0。
[0128]
其中，本实施例在位置动力学基本计算框架下引入连续介质力学理论，提出一种适用于柔体形变仿真的应变能约束方法。该方法结合位置动力学的计算效率优势以及传统连续介质力学成熟的理论成果，使仿真稳定、高效并且可控。本实施例引入应变能函数es(x)，定义能量约束方程c(x)＝es(x)＝0。下面将给出能量函数es(x)的计算过程，即步骤s23中以柔性物体的应变能约束作为内部约束包括：
[0129]
将软组织模型的形变示为一个连续的位移场u，定义ω为软组织模型所占据的体
积域，x∈ω为软组织模型未形变状态下的任意一点，则可得到物体变形函数φ(x)的表达式：
[0130]
φ(x)＝x u＝x
ꢀꢀꢀ
(5)
[0131]
其中，变形函数描述了材料空间中的点x与变形后对应的点x之间的映射关系，变形函数的变形梯度张量f为：
[0132][0133]
stvk模型使用非线性格林-拉格朗日(green
–
lagrange)应变张量，根据变形梯度张量f可求得应变张量e为：
[0134][0135]
其中i表示单位矩阵；
[0136]
根据变形梯度张量f和应变张量e求得应变能密度ψs为：
[0137][0138]
其中，e:e＝tr(e
t
e)，函数tr(a)表示矩阵a的迹，μ、λ为拉梅系数，只与两个材料属性杨氏模量k和泊松比v有关：
[0139][0140]
由于采用四面体单元对物体体积域进行离散，单元形函数采用线性拉格朗日形函数，于是每一单元对应的变形梯度进一步表示为：
[0141][0142]
其中ds为形函数矩阵，dm为参考形函数矩阵，是一个常数矩阵，ds、dm的向量表达式分别为：
[0143]ds
＝(x
0-x
3 x
1-x
3 x
2-x3)
ꢀꢀꢀ
(11)
[0144]dm
＝(x
0-x
3 x
1-x
3 x
2-x3)
ꢀꢀꢀ
(12)
[0145]
对于软组织模型的任意一个四面体单元，定义单元体积域为则储存在一个四面体单元内的应变能约束函数为：
[0146][0147]
其中，v是四面体单元未形变时的单元体积；
[0148]
根据变形梯度f得到软组织模型的应力张量p(f)：
[0149][0150]
则软组织模型的应力张量p(f)进一步化简，具体为将式(8)代入式(14)得到：
[0151][0152]
由式(7)可得应变张量的微分表达式，如式(15)所示：
[0153][0154]
将式(16)代入式(15)，可得到应力张量和变形梯度的最终关系式，如式(17)所示：
[0155]
p(f)＝f(2μe λtr(e)i)
ꢀꢀꢀ
(17)
[0156]
根据应力张量和变形梯度的关系式得到应变能梯度为：
[0157][0158]
对于软组织模型的四面体单元，只需求解4个顶点对应的应变能量梯度即可，定义为i＝{1,2,3,4}，则得到软组织模型的每一个四面体应变能梯度的最终表达式为：
[0159][0160]
步骤s24、将每一个顶点的预测位置{p0,p1,...pn}代入系统约束方程组，使用gauss-seidel迭代法进行求解，并限制迭代次数，求解得到一组位置修正量{δp0,δp1,...δpn}；
[0161]
其中，若使用牛顿迭代法求解一组包含n个约束方程的方程组时，每一次的迭代都需要对一个n
×
n的线性方程组进行求解计算，计算量大，影响仿真系统。因此，求解系统选用一种修改版本的迭代法，每一次迭代时对这n个线性方程拆分，每个方程独立线性化、独立求解。bender等验证该方法计算效率要明显地快于牛顿迭代法，并且在收敛性上只略差于牛顿迭代法。
[0162]
由此，步骤s24包括：
[0163]
将每一个顶点的预测位置{p0,p1,...pn}代入系统约束方程组；
[0164]
在每一次迭代时，为每个约束函数cj(p)找到一个合适的位置修正量δp，使之满足cj(p δp)＝0，则对每一个约束方程单独进行线性化为：
[0165][0166]
其中，表示为约束函数cj(p)的梯度算子求解系统将δp限制在的方向进行求解，将质点的个体质量考虑进来，可得到软组织模型的每一个顶点i的位置修正量δpi的推导式为:
[0167][0168][0169]
式中，s表示拉格朗日乘子，k∈[0,1]表示刚度参数，表示矢量c(p)在pi方向上的梯度算子；
[0170]
完成限制的迭代次数后，得到一组位置修正量{δp0,δp1,...δpn}。步骤s25、使用位置修正量{δp0,δp1,...δpn}修正每一个顶点的预测位置{p0,p1,...pn}，用修正后的预测位置{x0,x1,...xn}对顶点位置属性进行更新并更新预测速度，之后根据每一个顶点的修正后位置进行下一步的模型渲染和形变计算。
[0171]
其中，更新预测速度为通过二阶向后微分公式更新预测速度，具体为：
[0172][0173]
步骤s3、重复执行下一仿真步的计算，直至完成仿真。
[0174]
由此，结合步骤s2和步骤s3对于肝脏模型的仿真验证具体如下：
[0175]
对于仿真材料，有两个可控参数：刚度系数k，泊松比υ。刚度系数定义为约束的强度，刚度参数越大，越不易产生形变，而泊松比是指材料在单向受拉或受压时，横向正应变与轴向正应变的绝对值的比值，也叫横向变形系数，它是反映材料横向变形的弹性常数。对于求解系统，也有两个可控参数：迭代次数和时间步长。迭代次数越多，求解结果越精确，但求解时间也会相应增加，而时间步长越短，顶点位置的更新速度越快，仿真结果越准确，但同样系统的计算负担会相应增加。为了便于分析比较，在所有验证及所有模型中均采用相同的参数：k＝0.4，v＝0.4，迭代次数取5，时间步长取0.005。在其他实施例中，杨氏模量k为[0.3,0.5]，泊松比v为[0.3,0.5]。可选地，迭代次数为[3,8]，时间步长为[0.001,0.01]。
[0176]
实验1：碰撞约束测试
[0177]
为了更好地进行结果比较和性能分析，在本实施例中，除了交互式外力算法的测试外，本实施例将只考虑重力对模型产生的形变影响。而为了达到更好地仿真效果，必须先实现环境碰撞约束。因此，在实验1首先对环境碰撞约束的实现效果进行实验分析。
[0178]
在未加环境碰撞约束之前，系统不会将地板的约束影响纳入到形变计算中，因此当肝脏模型在下降过程中撞到地板时并不会产生应有的碰撞效果，而是直接穿过地板继续下落，如图9所示。
[0179]
在加入环境碰撞约束后，系统在每一步的形变计算中都会对模型进行位置有效性判断和修正，从而将模型的运动空间限制在有限的范围内，即地板之上。肝脏模型在未接触地板之前的形变状态与未加约束时的形变状一致，而一旦与地板发生接触立即停止下落，并在地板的碰撞作用下产生合理的碰撞反弹效果，直至最后静止在地板上，完整的下落过程如图7所示。
[0180]
实验结果表明，将环境碰撞约束应用在动态仿真物体与静态环境物体的碰撞交互场景中能达到得较好的效果，可以将肝脏模型的活动空间限制在有限的范围内，并且产生必应的碰撞反应，从而使仿真结果更加合理、更为真实。
[0181]
实验2：性能分析
[0182]
根据本实施例前文的理论论述，在基于位置动力学应变能约束的形变模型中，系统约束方程组的方程个数与模型单元个数一致。为了验证模型单元数量对形变仿真的性能影响，除了使用步骤s1中的3个肝脏模型进行实验外，还选取一个具有12829个单元的模型作为参照。几个实验模型均只受到重力的影响，下落过程中也均会与地板发生碰撞，其中肝脏模型3的仿真效果如图8所示。
[0183]
3个模型形变计算的平均求解时间如表1所示:
[0184]
表1形变仿真平均求解时间
[0185]
模型单元个数平均求解时间(ms)单元平均时间(ms)肝脏模型15960.508.39
×
10-4
肝脏模型253955.209.64
×
10-4
肝脏模型373256.809.28
×
10-4
参照模型1282912.961.01
×
10-3
[0186]
实验结果表明系统求解时间随单元数量的增加而增加。为了量化说明系统求解时间与单元个数的关系，将表中平均求解时间除于单元个数得到每个单元对应的约束方程的平均求解时间，3个模型的单元平均求解时间近似相等，可见系统求解时间与单元个数近似为线性关系。同时，实验结果也验证了本实施例所提供的仿真方法在应对单元个数超过10000较大精度的模型仿真的可行性。
[0187]
由于在位置动力学中系统约束方程组中每个方程均是独立进行求解的，利用这一特点可以采用并行计算来进一步提高系统的求解速度。本实施例采用openmp实现多线程并行计算，并与单线程的求解系统进行实验对比，实验结果如表2所示。实验结果表明，采用openmp多线程并行计算可以有效地减少系统的求解时间，增加系统仿真效率。同时实验发现随着单元个数的增加，多线程并行计算的优势越显著。这也进一步提高了位置动力学在具有大规模网格的仿真场景中的可行性。
[0188]
表2单线程、多线程求解系统性能对比
[0189]
模型单线程求解时间(ms)openmp多线程求解时间(ms)加速比肝脏模型10.500.401.25肝脏模型25.201.852.81肝脏模型36.802.103.24参照模型12.963.773.44
[0190]
由此，根据本发明仿真场景的特点选取位置动力学作为基础形变模型。采用应变能约束作为内部约束，采用环境碰撞约束作为外部约束，二者构成一个系统约束方程组，进行方程组的求解得到模型顶点的新位置，从而实现形变仿真。最后，通过实验验证形变仿真方法的可行性，实验结果证明，仿真方法稳定、高效，并且有较好的物理真实感，能基本满足虚拟三维场景中软组织与刚性地面碰撞仿真的需求。
[0191]
实施例二
[0192]
请参照图10，一种软组织与刚性地面的碰撞仿真装置1，包括存储器3、处理器2及存储在存储器3上并可在处理器2上运行的计算机程序，所述处理器2执行所述计算机程序时实现上述实施例一中的步骤。
[0193]
由于本发明上述实施例所描述的装置/装置，为实施本发明上述实施例的方法所采用的装置/装置，故而基于本发明上述实施例所描述的方法，本领域所属技术人员能够了解该装置/装置的具体结构及变形，因而在此不再赘述。凡是本发明上述实施例的方法所采用的装置/装置都属于本发明所欲保护的范围。
[0194]
本领域内的技术人员应明白，本发明的实施例可提供为方法、装置或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例，或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机
可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等)上实施的计算机程序产品的形式。
[0195]
本发明是参照根据本发明实施例的方法、设备(装置)和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。
[0196]
应当注意的是，在权利要求中，不应将位于括号之间的任何附图标记理解成对权利要求的限制。词语“包含”不排除存在未列在权利要求中的部件或步骤。位于部件之前的词语“一”或“一个”不排除存在多个这样的部件。本发明可以借助于包括有若干不同部件的硬件以及借助于适当编程的计算机来实现。在列举了若干装置的权利要求中，这些装置中的若干个可以是通过同一个硬件来具体体现。词语第一、第二、第三等的使用，仅是为了表述方便，而不表示任何顺序。可将这些词语理解为部件名称的一部分。
[0197]
此外，需要说明的是，在本说明书的描述中，术语“一个实施例”、“一些实施例”、“实施例”、“示例”、“具体示例”或“一些示例”等的描述，是指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。
[0198]
尽管已描述了本发明的优选实施例，但本领域的技术人员在得知了基本创造性概念后，则可对这些实施例作出另外的变更和修改。所以，权利要求应该解释为包括优选实施例以及落入本发明范围的所有变更和修改。
[0199]
显然，本领域的技术人员可以对本发明进行各种修改和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也应该包含这些修改和变型在内。