CN106227922A - 在Laplace‑Beltrami形状空间基于样例的弹性材料的实时仿真方法 - Google Patents

Abstract

一种在Laplace‑Beltrami形状空间基于样例的弹性材料的实时仿真的方法，其内容是：利用Laplace‑Beltrami算子以及线性有限元离散化方法计算仿真物体和样例在三维四面体网格上的特征函数和特征值；将样例的特征函数和仿真物体的特征函数进行配准；将仿真物体的位置投影到Laplace‑Beltrami形状空间，利用非线性形状插值求解样例流行上的目标形状；在当前时刻每个顶点上建立局部坐标系添加形状细节的方式得到欧氏空间下的最终形变结果；在模态基底构建的降维子空间求解非惯性系下形变体仿真动力学方程，得到新一时刻降维空间下的位移；降维空间下的位移利用预处理过程中计算的全局与降维空间下的投影矩阵，投影得到全局中的位移，并绘制形变结果。

Description

在Laplace-Beltrami形状空间基于样例的弹性材料的实时仿 真方法

技术领域

本发明属于计算机图形学、基于物理的仿真领域，具体涉及利用给定样例来指导形变体形变行为的一种实时仿真方法，从而能够获得更加自如的仿真和动画设计效果。

背景技术

逼真的效果和多样化的固体形变体的运动方式设计是近年来关注的热点。文献[Martin S,Thomaszewski B,Grinspun E,et al.“Example-based elastic materials,”Acm Transactions on Graphics,30(4),pp.76-79,2011.]提出的基于样例的仿真方法，通过给定的一组样例来控制仿真物体的形变方式得到了艺术化的效果，给物体的仿真模拟方式提供了一种新的手段。该方法通过定义一个形变能量，将仿真过程中的形变状态引导至由样例定义的流行表面上，达到设计形变的目的。文献[S,Botsch M.“Example-Driven Deformations Based on Discrete Shells”,Computer GraphicsForum.Blackwell Publishing Ltd,pp.2246-2257,2011.]将网格形变方法和网格插值方法相结合，提出了薄壳体的基于样例驱动的仿真方法。文献[Schumacher C,ThomaszewskiB,Coros S,et al.“Efficient simulation of example-based materials,”Proceedingsof the 11th ACM SIGGRAPH/Eurographics conference on ComputerAnimation.Eurographics Association,pp.1-8,2012.]将艺术指导与非刚性材料的形变相结合，给定的样例以及样例空间中插值得到的中间状态都使用不共用顶点的单元表示，直接线性插值单元的顶点位置得到样例空间中的中间状态，而不用像Martin等人提出的方法先得到应变线性插值的中间状态，然后再求解非线性优化得到符合真实的物体形状。此外，作者还对基于样例的弹性形变进行了扩展，提出了基于样例的塑性形变的仿真方法，使得物体的塑性形变也能符合设计的要求。文献[Koyama Y,Takayama K,Umetani N,et al."Real-Time example-based elastic deformation,"In:Proceedings of the ACMSIGGRAPH/Eurographics conference on Computer Animation.Lausanne,Switzerland:Eurographics Association,pp.19-24,2012]从几何方法角度，在无网格的形状匹配方法框架下做了同样的问题。由于利用几何形状匹配的方法，且将给定样例的线性插值直接作为样例空间，因此做到了实时的形变，但效果不如Matin等人做的好。文献[Song C,ZhangH,Wang X,et al."Fast corotational simulation for example-driven deformation,"Computers&Graphics,40(5),pp.49-57,2014.]提出利用一个共旋的有限元框架进行仿真，通过改进的线性柯西应变的二次能量函数，求解的能量优化问题可以通过求解一个线性系统得到，从而加快仿真的速度。文献[Zhu F,Li S,Wang G."Example-Based Materials inLaplace-Beltrami Shape Space,"Computer Graphics Forum,34(1),pp.36-46,2015.]首次提出利用Laplace-Beltrami形状空间的特征分析方法，将物体投影到特征函数张成的形状空间，可以使得输入样例与物体具有不同的拓扑关系，同时，在这个降维的形状空间通过最小化形变能量的方式求解目标形状也达到了加速的目的，整个仿真速度也得到了一定的提高，但由于求解过程中涉及到在整个网格上求解形变能量，所以效率上还是受到网格分辨率的影响，还没有达到实时的效率。从Laplace-Beltrami形状空间到欧氏空间形状的重建问题上，文献[Dey T K,Ranjan P,Wang Y."Eigen deformation of 3D models,"VisualComputer,28(6-8),pp.585-595,2012.]中提出利用仿真物体的Laplace-Beltrami特征函数在形状空间中的投影(被称为特征骨骼)模拟形变，特征骨架在欧氏空间下的重建结果与初始形状的向量差得到细节向量，利用仿真过程中每一时刻在局部坐标系下细节向量不变的特征，通过特征骨架快速模拟形变，利用细节向量获得重建的结果，获得不错的仿真效果。

实时的仿真效率一直是基于物理仿真永恒不变的话题。模态分析理论作为一种能够提高仿真速度的方法，源于振动理论，通过分析物体的形变行为，得到一组能够反应物体自然形变的模态，并根据需要去除频率较高的模态，减少整个系统的自由度，达到加速的目的。文献[Pentland A,Williams J."Good vibrations:modal dynamics for graphicsand animation,"Acm Siggraph Computer Graphics,23(3),pp.207-214,1989.]首次在图形学领域引入线性的模态分析，结合有限元方法仿真固体的形变，此后图形学领域陆续有相关研究出现。文献[Choi M G,Ko H S."Modal warping:real-time simulation oflarge rotational deformation and manipulation,".IEEE Transactions onVisualization&Computer Graphics,11(1),pp.91-101,2005.]提出了一种类似“刚度扭曲”的方法，基于每个仿真结点的旋转不变性，通过分析出每个仿真结点的运动旋转分量在模态空间中的表现形式，使得线性模态能够处理大形变。文献[Barbic J,James D."Real-Time subspace integration for St.Venant-Kirchhoff deformable models,".AcmTransactions on Graphics,24(3),pp.982-990,2005.]提出了一种模态导数，使得模态分析方法能够描述大形变，利用Stvk材质的应力可以表示成位移的三次多项式形式，多项式的系数可以预计算，因此可以直接在运动方程中使用多项式形式的非线性应力，无需做线性近似。文献[An S S,Kim T,James D L."Optimizing Cubature for EfficientIntegration of Subspace Deformations,".Acm Transactions on Graphics,27(5),pp.32-39,2009.]提出了一种计算降维空间弹性力的积分优化方法，可以比降维的模态导数方法更快，并且不局限于物质材料属性。文献[Harmon D,Zorin D."Subspaceintegration with local deformations,"Acm Transactions on Graphics,32(4),pp.96-96,2013.]提出了一种改进的基底的方法能够捕捉由于碰撞产生的局部形变。文献[Hahn F,Thomaszewski B,Coros S,et al."Subspace clothing simulation usingadaptive bases,"Acm Transactions on Graphics,33(4),pp.1-9,2014.]利用一种自适应基底选择机制，针对不同姿势的基底模拟布料的褶皱细节。文献[Teng Y,Meyer M,Derose T,et al."Subspace condensation:full space adaptivity for subspacedeformations,"Acm Transactions on Graphics,34(4),2015.]为了解决降维空间中无法处理碰撞的问题，提出一种将全局方法和降维方法相结合策略捕捉局部的形变细节。

在传统的利用样例驱动固体形变体运动的仿真方法中，大多数研究是仿真物体与样例采用同样的网格拓扑结构，用相同的仿真物体的不同姿势来驱动形变体的运动。Zhu等人提出的利用Laplace-Beltrami算子对网格做特征分析，得到的特征基底构建的形状空间，利用特征函数的等距不变性建立不同形状的部件间的对应关系，突破了仿真物体与样例的拓扑一致的局限性，但这些方法还不能达到实时仿真。目前基于样例的仿真方法中还有以下不足：

(1)现有的基于样例的仿真方法多是仿真物体与样例采用同样的拓扑结构，缺少样例指导的多样性。

(2)Zhu等人提出的利用Laplace-Beltrami形状空间可以弥补上述不足，但由于在计算目标形状的过程中，求解二次能量需要在整个网格上求解，导致求解过程与网格的分辨率相关，非常耗时，整个仿真效率较低。

发明内容

针对上述存在的问题，本发明的目的在于提供一种在Laplace-Beltrami形状空间基于样例的弹性材料的实时仿真方法。本发明基于样例的弹性材料的仿真方法，以Laplace-Beltrami形状空间作为降维的样例插值空间，采用部分积分单元代替全局计算的优化方法来近似求解形变能量，加快在仿真过程中最耗时的求解当前仿真物体形状与样例流行上对应目标形状的过程，并且采用降维仿真方法加快整个动力学方程求解速度，从而使整个基于样例的仿真方法不仅可以采用不同的样例来指导仿真物体，并且在效率上达到了实时。

为实现发明目的，本发明是通过以下技术方案实现的：一种在Laplace-Beltrami形状空间基于样例的弹性材料的实时仿真方法，该方法包括预处理步骤和仿真步骤：

1、预处理步骤包括如下内容：

1.1利用Laplace-Beltrami算子以及线性有限元离散化方法计算仿真物体和样例在三维四面体网格上的特征函数和特征值，特征函数表示为特征值对应的特征向量，将这些特征函数归一化；在预处理过程中求仿真物体和样例的特征函数过程是：将物体的空间结构采用四面体网格离散化方式描述，并利用线性有限元形函数计算体网格上的特征值和特征向量，求得的特征值由小到大排列，并利用归一化后的特征向量(即Laplace-Beltrami特征函数)张成的形状空间作为形状插值的计算空间；

1.2将样例的特征函数和仿真物体的特征函数进行配准，通过指定对应区域，使得仿真物体与样例之间具有相同的特征函数排序；

1.3将样例形状和仿真物体投影到各自的特征函数上，得到对应的在Laplace-Beltrami形状空间中的形状描述；

1.4利用模态分析方法，求得仿真物体非线性的模态基底作为全局坐标与降维空间下的投影矩阵；并将模态基底作为输入，利用准静态仿真方法，预先求得仿真过程中计算形变能量以及内力的局部的优化积分单元及对应的权重；

2、仿真步骤包括如下内容：

2.1将仿真物体的位置投影到Laplace-Beltrami形状空间，利用非线性形状插值求解样例流行上的目标形状，求解目标形状的过程包括两个部分：

a.首先计算用于求解目标形状所需要的各个样例的权重；

b.利用计算得出的权重求解目标形状；即在Laplace-Beltrami形状空间通过非线性插值求解一个目标插值形状，优化过程是求使得仿真物体在样例流行上的形状相对于所有样例形变能量的加权之和最小的形状；为加快求解形变能量的过程，利用特征函数以及模态基底，在预处理过程中建立Laplace-Beltrami形状空间与模态基底子空间的联系，直接将Laplace-Beltrami形状空间投影到模态降维子空间，省略每个时间步从Laplace-Beltrami形状空间投影到欧氏空间，再投影到模态降维子空间的过程，并利用不局限于材质模型，局部的优化积分单元代替全局的单元近似求解能量和能量梯度的方法；

2.2获得Laplace-Beltrami形状空间中的目标形状后，利用特征函数投影得到欧氏空间下的形状，该形状只包含整体信息，缺失局部信息，根据上一时刻的形变结果在特征函数上的投影前后的向量差构建细节向量，在当前时刻每个顶点上建立局部坐标系添加形状细节的方式得到欧氏空间下的最终形变结果；

2.3在模态基底构建的降维子空间求解非惯性系下形变体仿真动力学方程，得到新一时刻降维空间下的位移；

2.4利用刚体运动方程计算刚性运动部分，弥补降维仿真方法中，降维基底不能模拟刚性运动的不足；

2.5降维空间下的位移利用预处理过程中计算的全局与降维空间下的投影矩阵，投影得到全局中的位移，并绘制形变结果。

与现有技术相比，本发明的有益效果如下：

本发明提出的一种在Laplace-Beltrami形状空间基于样例弹性材料的实时仿真方法，兼顾基于样例的仿真多样性和效率问题。

(1)采用体网格计算Laplace-Beltrami特征函数，避免了现有的方法利用表面网格计算，在仿真过程中还需要从表面网格投影到体网格，近似计算体网格上的形变，提高了仿真的效率和精度；

(2)采用积分优化方法，利用有限个积分单元计算形变能量，提高在形状空间求解目标形状时，在欧氏空间下求解能量的效率；并利用降维子空间的仿真替代全局网格上求解，加快效率，使得整个仿真达到了实时；

(3)利用优化积分单元的方式，与以往方法相比可以不局限于形变体的材质属性，应用于非线性的材质，使得本发明可以拥有更好的通用性；

(4)从Laplace-Beltrami形状空间投影到欧氏空间的过程中，涉及到重建问题，通过建立局部坐标系，利用改进的添加细节的方法能够获得更好的形变后的三维网格。

附图说明

图1是预处理步骤流程框图和仿真步骤流程框图；

图2是四面体网格计算的Laplace-Beltrami特征函数绘制的剖面图；

图3是armadillo模型的特征函数颜色图；

图4是配准后的仿真物体与样例的特征函数颜色图；

图5是从形状空间到欧氏空间网格形变重建改进方法与原方法对比图；

图6是实验结果。

具体实施方式

为使本发明的上述目的、特征和优点能够更加明显易懂，下面通过具体实施例和附图，对本发明做进一步说明。

Laplace-Beltrami算子是定义域在R²空间中的Laplace算子，推广到定义在任意黎曼流行的Laplace算子。利用Laplace-Beltrami算子对形状进行分析，得到一组正交向量，通过引入谱分析的方法，把这组正交的基底看成频谱，从低维到高维的排列表示形状从整体到局部的信息。利用Laplace-Beltrami形状空间，可以用相同的样例指导不同的物体的形变仿真。采用在四面体上构建Laplace-Beltrami特征函数，利用线性有限元作为黎曼流行上的函数，求解体上的特征函数，可以使得在仿真过程中从形状空间直接投影到体网格上获得体网格的形变结果。

本发明提出一种在Laplace-Beltrami形状空间基于样例的弹性材料的实时仿真方法，该方法包括预处理步骤和仿真步骤：预处理步骤流程框图如图1左侧图所示，仿真步骤流程框图如图1右侧图所示。

一.预处理特征函数计算与配准

1-1.四面体上构建Laplace-Beltrami特征函数

利用Laplace-Beltrami算子在物体的表面网格上创建的特征函数，在仿真过程中，每次在Laplace-Beltrami形状空间上的投影获得的是表面网格上顶点的数据，需要投影到体上，在仿真过程中带来了一定的时间消耗。本发明的方法利用Laplace-Beltrami形状空间直接在体上创建特征函数，避免了表面网格和体网格之间的数据投影，提高了效率。另外，利用体网格上的特征函数，也方便利用优化积分单元方法求解形变能量，下面具体描述一下利用Laplace-Beltrami算子在体上求解特征函数的方法。利用“Shape-DNA”方法对Laplace-Beltrami特征函数的定义及求解方法，具体地，Laplace-Beltrami算子是定义在黎曼紧致流行上的二次可微实函数：

Δf:＝div(grad f)

其中，grad和div分别为相应的黎曼度量下定义的梯度算子和散度算子。Laplace-Beltrami算子的特征函数f与特征值λ统称为该算子的特征谱，定义为Helmholtz公式(又称为Laplacian特征值问题)：

Δf＝λf

采用四面体网格计算物体的特征函数。定义f∈C²表示一个黎曼流行M上的实函数，函数g和f定义方式相同，Nabla算子和Laplace-Beltrami算子定义为：

▽(f,g):＝＜grad f,grad g＞

Δf:＝div(grad f)

其中＜,＞表示内积。给定一个局部参数化映射：ψ:Rⁿ→R^n+k是黎曼流行R^n+k上的子流形：

$g_{\mathrm{ij}}:=<{\&PartialD;}_i\mathrm{\&psi;},{\&PartialD;}_j\mathrm{\&psi;}>,G:=\left(g_{\mathrm{ij}}\right)$

$W:=\sqrt{\det G},\left(g^{\mathrm{ij}}\right):=G^{-1}$

其中，G表示g的矩阵形式，det表示G的行列式。

在求解特征函数的计算过程中，首先要把Laplace-Beltrami算子的特征值求解问题转化为变分问题，利用Green公式，边界条件采用其中Nabla算子表示为：

用试函数乘以Helmholtz公式，得到3D空间下的表示：

其中，dσ＝Wdudvdω。

对于任意一个四面体T上顶点p₁,p₂,p₃,p₄，形函数选择一阶有限元线性形函数，将任意四面体所在的欧氏坐标系转换为正四面体的标准坐标系，则四面体内的向量表示为：

P(ξ,η,ζ)＝p₁+(p₂-p₁)ξ+(p₃-p₁)η+(p₄-p₁)ξ＝p₁N₁+p₂N₂+p₃N₃+p₄N₄

其中，ξ,η,ζ表示转换到标准坐标系的三个坐标方向,N_i表示形函数。

所述Laplace-Beltrami算子在体上的特征函数求解转化为求解一个广义特征值问题，AU＝λBU。其中，

$A=\left(a_{lm}\right):=\&Integral;\&Integral;\&Integral;\underset{j,k}{\&Sigma;}\left({\&part;}_j{N}_l\right)\left({\&part;}_k{N}_m\right)g^{ij}d\mathrm{\&sigma;}$

B＝(b_lm):＝∫∫∫N_lN_mdσ

由此，所求得的特征函数表示成特征向量空间的特征基底，在Laplace-Beltrami形状空间上的向量可以表示成这些基底的线性组合。求解的广义特征值按照从小到大排列，每个特征值对应的特征向量描述了物体的形状信息。特征值小的表示的是物体形状的全局特性，特征值大的描述的是物体形状的细节信息。求解的流行是封闭的，采用诺依曼边界条件，得到的第一个特征值为0，对应的特征向量是常函数，在基于样例的仿真过程中，将第一个特征值去掉，选取所求特征值的低维部分对应的特征向量作为基底张成的特征空间作为Laplace-Beltrami形状空间。计算得到的特征函数是一个n×m的矩阵，n表示物体的顶点个数，m表示所取低维特征函数的个数，首先把特征函数归一化，保证每个顶点在每个特征函数上的表示都是一个0到1之间的值，将特征向量的值利用颜色信息表征物体形状特点，根据上述方法获得的特征函数的绘制效果，即四面体网格计算的Laplace-Beltrami特征函数绘制的剖面图如图2所示，这是armadillo模型从中间剖开，展示模型第2个特征函数内部的绘制效果图。图3所示是该模型的第3个到第8个特征函数的绘制结果。

1-2.仿真物体与样例特征函数配准

基于Laplace-Beltrami形状空间的仿真，能够使仿真物体与样例具有不同的拓扑关系，即仿真物体的四面体网格与样例的四面体网格表示不同，这使得仿真物体与样例计算得到不同的特征函数，而不同形状的特征函数需要配准，才能作为一个近似共同的形状空间的基底，对于多个满足或近似满足等距条件且具有简单特征谱的形状而言，特征函数是相似的，因此在预处理过程中通过对样例的特征函数调整特征向量的符号和顺序来配准，配准前后的仿真物体与样例的特征函数颜色图如图4所示，图中第一行表示仿真物体的特征函数，第二行是配准前样例的特征函数表示图，第三行是配准后样例的特征函数表示图。

二.基于样例材料仿真方法

2-1.基于样例材料的动力学方程

基于样例材料的仿真主要是在传统的基于有限元的形变体仿真中，引入一个基于形变的能量，将仿真物体的形变引导到样例流行空间的形状上。离散后的运动方程表示为：

$M\stackrel{\&CenterDot;\&CenterDot;}{x}+D\stackrel{\&CenterDot;}{x}+f_{int}+f_e=f_{ext}$

其中，M表示物体的质量矩阵，x,分别表示形变位移和加速度，D是阻尼矩阵，这里采用Reign阻尼，D＝αI+βK，α,β表示Reign阻尼系数，I表示单位矩阵，K表示刚度矩阵，f_int是内力，表示物体形变后的形状x与初始形状X之间的形变能量的梯度，f_e是样例引导力，表示物体形变后的形状x与样例流行上的目标形状X_target之间的形变能量的梯度，f_ext表示重力、摩擦力等外力的合力。

2-2.降维空间下求解目标形状

整个查找目标形状的过程分两个步骤。

首先，计算样例插值权重，计算方法如下：

在仿真过程中根据仿真物体当前形状c与每个样例e的接近程度，动态决定每个样例的控制强度,即样例插值权重：

$\underset{t}{min}\frac{1}{2}\left|\right|\underset{i=1}{\overset{k}{\&Sigma;}}{\mathrm{\&omega;}}_i{e}_i-c|{|}_F^2$

其中，||·||_F表示Frobenius范式。

其次，利用插值权重，求解一个能量优化问题。

对样例流行上的形状，定义其相对于所有样例的形变能量的加权之和，用于指导的目标样例t是最小化该加权能量和的对应形状：

$\underset{t}{min}\underset{i=1}{\overset{k}{\&Sigma;}}{\mathrm{\&omega;}}_{i}E(t,e_i)$

其中，E(·,·)表示两个形状之间的非线性形变能量，ω_i(1≤i≤k)衡量了每个样例的指导强度，且通过上述定义的目标形状在以形变能量作为距离度量的标准下，距离所有的输入样例最近。

传统的利用Laplace-Beltrami形状空间的基于样例材料的仿真中，求解上述最小化能量问题时，由于定义的形变能量是在欧氏空间下求解的，因此在每次迭代求解形变能量的过程中，需要将求得的Laplace-Beltrami形状空间的表示投影到欧氏空间，特征函数是定义在表面网格上的，因此投影的过程首先是投影到表面网格上，再利用线性插值方法计算对应的体网格上的形变位移，利用体网格上的形变位移才能获得形变能量,因此求解形变能量是与网格的拓扑分辨率相关的，尽管需要比较少的迭代次数，但在全局上求解形变能量也是相当耗时的，为了进一步提高仿真的效率，本发明提出在降维空间下求解形变能量，利用局部的优化积分单元代替全局单元加速求解过程。

2-3.降维空间下的形变能量及内力求解方法

2-3-1.降维基底求解方法

线性的模态分析方法通过对形变前物体的刚度矩阵进行特征值分析，得到一组特征向量作为形变的模态基底。物体的位移可以表示成模态基底的线性组合，将这种表示代入到物体的运动方程，物体的运动方程被解耦合成一组无关的常微分方程，且方程的个数远小于原方程的维度，因此求解速度大大加快。在有限元方法中，为了简单计算，通常把内力显性化，此时物体的运动方程表示为一个3n维的线性方程组：

$M\stackrel{\&CenterDot;\&CenterDot;}{u}+D\stackrel{\&CenterDot;}{u}+Ku=f_{ext}$

其中，M∈R^3n×3n表示质量矩阵，D∈R^3n×3n是阻尼矩阵，f_ext∈R³ⁿ表示系统的外力向量，K∈R^3n×3n是刚度矩阵，等于内力对位移的梯度。u,分别表示位移，速度，加速度向量。模态分析通过求解如下式的广义特征值问题将运动方程解耦合成3n个线性无关的常微分方程：MΦΛ＝KΦ，Λ∈R^3n×3n是对角矩阵，对角线上的每个值是矩阵M^-1K的特征值，Φ∈R^{3n ×3n}的每一列是对应的特征向量。Φ矩阵的列构成了物体形变空间的一组基底，称为模态基底。因此，物体的位移u(t)可以表示成基底的线性组合：u(t)＝Φq(t)，将其带入到运动方程中，并将等式两边左乘矩阵Φ^T，得到：

${\mathrm{\&Phi;}}^{T}M\mathrm{\&Phi;}\stackrel{\&CenterDot;\&CenterDot;}{q}+{\mathrm{\&Phi;}}^{T}D\mathrm{\&Phi;}\stackrel{\&CenterDot;}{q}+{\mathrm{\&Phi;}}^{T}K\mathrm{\&Phi;}q={\mathrm{\&Phi;}}^T{f}_{ext}$

其中，Φ^TMΦ＝I，Φ^TDΦ＝Λ。Φ^TDΦ通常是稠密矩阵，当选择Raleign阻尼D＝αM+βK时，Φ^TDΦ＝αI+βΛ是对角矩阵。根据需要把Φ中对应特征值比较大的列去掉，可以减少常微分方程的个数，得到r个特征向量构成的基底，设g＝Φ^Tf_ext，则运动方程转化为r个无关的常微分方程：

${\stackrel{\&CenterDot;\&CenterDot;}{q}}_i+(\mathrm{\&alpha;}+{\mathrm{\&beta;\&lambda;}}_i){\stackrel{\&CenterDot;}{q}}_i+{\mathrm{\&lambda;}}_i{q}_i=g_i$

由于r＜＜3n，所以求解这个降维的方程组要比求解3n个方程组快得多，达到降维的目的。利用模态分析的降维仿真方法，针对Stvk材质，通过预先计算内力系数，大幅度提高仿真效率。

2-3-2.优化积分方法

模态分析的降维仿真方法只适用于Stvk材质，因为Stvk材质的内力可以表示成位移的三次多项式形式，多项式的系数可以在预计算的过程中求解，不占用仿真时间。然而其他材质的内力求解并不具有这一特性。优化积分单元的加速仿真方法弥补了这一不足，将内力表示为优化的积分单元在降维下的力的加权求和的形式。同时，这种方式在时间复杂度上是O(r³)，而模态分析方法的时间复杂度O(r⁴)，在效率上也有所提高，本发明采用这种优化的积分方法来计算内力和形变能量提高效率。

优化的积分方法求内力，从能量的角度考虑，对于降维基底构建的子空间下的位移q，形变能量表示为：

W(q)＝∫_ΩΨ(X,q)dΩ_X

其中，Ψ(X,q)是未形变区域Ω上物质点X的非负的应变能量密度。降维空间下的内力表示成能量的梯度的积分形式：

f(q)＝-▽_qW(q)＝-∫_Ω▽_qΨ(X,q)dΩ_X＝∫_Ωg(X,q)dΩ_X

其中表示降维空间下的力密度。降维空间下的内力近似表示成n个积分单元加权求和的形式：

$f\left(q\right)={\&Integral;}_{\mathrm{\&Omega;}}g(X,q){\mathrm{d\&Omega;}}_X\&ap;\underset{i=1}{\overset{n}{\&Sigma;}}{\mathrm{\&omega;}}_{i}g(X_i;q)$

其中ω_i表示第i个积分单元上对应的权重。利用一个贪心算法逐渐的增加积分单元采样点，通过计算训练集合(q,f(q))，在外力作用下，给定物质的内力与积分单元上的力加权求和所获得结果的误差最小的近似的积分单元及其权重表示，通过此方法计算得到的积分单元和对应权重在降维子空间下计算仿真过程中的形变能量。

本发明利用优化积分单元的方法，利用局部积分单元加权求和的形式求解形变能量和内力，并且在模态基底构建的降维子空间中求解，加快了求解目标形状是计算形变能量的过程。具体地，在求解形变能量的过程中，需要从Laplace-Beltrami形状空间的描述投影到欧氏空间获得仿真物体的位移，再投影到模态降维子空间下求解形变能量，这个投影过程非常耗时，本发明提出在预计算过程中，求解Laplace-Beltrami形状空间与降维子空间的投影矩阵，使得在每次计算形变能量过程中，输入的位移可以直接获得两个空间上的对应关系，省略仿真物体在欧氏空间上的描述，使得整个求解能量的过程与网格的拓扑结构无关，从而加快求解目标形状过程中反复迭代求解能量的效率。同时，在求解r维线性方程组过程中，形变物体的内力以及刚度矩阵也采用优化积分单元代替全局积分单元求解，利用降维仿真方法求解方程。

2-4.基于特征系数的三维网格重建方法

仿真过程中获得的目标形状，需要从形状空间到欧氏空间的投影获得欧氏空间下的位移，计算由样例引导作用产生的力，由于选取的特征函数是部分特征向量张成的形状空间，因此只能表示物体形状的整体信息，缺少细节特征。为了使仿真有效的进行，必须重新添加细节。本发明扩展了利用Laplace-Beltrami特征函数构建特征骨骼，再添加形变细节的方法，具体步骤如下：

(1)计算初始细节向量

在初始时刻，将物体投影到预先选取的少数特征函数张成的形状空间得到特征系数，再将此系数利用特征函数投影到欧氏空间，获得初始的投影网格，得到的投影网格只具有整体的信息，缺失大量的细节信息，投影前后形状的向量差值求得初始细节向量；利用在仿真过程中细节向量在局部坐标系保持不变的特征，在每个时间步重新投影到物体上；

(2)建立局部坐标系

在每个时间步，物体上每个点重新建立局部坐标系，建立的这个局部坐标系需要在每个时间步保持统一，本发明采用的局部坐标系构建策略是：在初始的时刻，为每个顶点选取其所在的若干个四面体中质量最好的一个四面体；四面体质量评价如下：

$6\sqrt{2}\&CenterDot;V\&CenterDot;\frac{l_{\mathrm{harm}}}{l_{\mathrm{rms}}^4}$

其中，得到的四面体质量是一个0到1之间的值，正四面体的质量最好，值为1。选择质量最好的四面体，是为了防止在仿真过程中，顶点所在的四面体形变之后变为狭长四面体，导致被选为局部坐标系的参考点的两个顶点，在形变过程中产生三点共线或接近共线的情况，导致建立的局部坐标系不准确，细节恢复不好，影响重建效果；

在某一四面体中，顶点的局部坐标系建立方式首先选择与其具有邻接关系的任意两个顶点，该顶点与这两个顶点可以计算两个向量差，任取一个向量作为局部坐标系的一个方向，这两个向量张成的平面的法向作为局部坐标系的第二个方向，第三个方向选为已标定的两个局部坐标系方向的叉乘。

(3)添加细节向量

利用Laplace-Beltrami特征函数构建特征骨骼的形变方法中，每个时间步，以初始的细节向量作为参考值，利用局部坐标系的变换添加到形变后的网格重建结果上；仿真过程中，随着时间的推进，相对于初始形状，形变越来越大，当形变体发生大形变时利用初始的细节向量重建效果不理想；实验过程中，观察到形变效果在每步仿真过程中有一个累加的过程，相对于上一时刻，当前形状的形变是小形变，利用这个特点，在每个时间步，将上一时刻的细节向量作为参考值添加到下一帧，获得重建网格。本发明利用改进的方法进行实验，图5所示是从形状空间到欧氏空间网格形变重建改进方法与原方法对比图，第一行第一张图是初始形状经过一段时间的仿真的一个中间的形变状态，第二张图片是第一张图的形状在一个时间步形变下获得的网格，第二行第一张图是将形变后的状态从Laplace-Beltrami形状空间利用7个特征函数重建到欧氏空间下的结果，缺少了很多的细节，第二张是在局部坐标系下利用网格的初始细节信息重建的结果，可以看出经过一段时间的仿真之后，形变相对于初始状态比较大，因此利用初始细节信息的添加方式重建效果不理想，不满足仿真的要求。第三张图是采用本发明中的方法将第一行第一张图的细节信息添加到第二张图在7个特征函数重建结果的网格上，重建效果获得很大改善。

三.降维仿真中刚体和形变体的分离

利用模态基底构建降维子空间的仿真方法，由于模态基底是在预计算过程中得到的，在仿真过程中不重新计算，模态基底的线性表示不能模拟如旋转平移等刚性运动，因此，在降维仿真中，要模拟形变体的刚性行为，需要将形变体仿真与形变体的刚性行为分离。

3-1降维仿真中刚性行为的分离

在降维空间中分离刚性运动的方法就是通过建立一个非惯性系，添加几个非惯性力，来抽离形变体运动中的刚性运动部分，在形变体仿真之后再单独添加刚性运动的方式计算。根据Terzopoulous等人研究，将降维空间下的非惯性系下的四个力表示为：科里奥利力，惯性力，欧拉力和离心力，具体地：

$f_i^{cor}=-2m\&CenterDot;\mathrm{\&omega;}\&times;(U_i\&CenterDot;\stackrel{\&CenterDot;}{q}),f_i^{ine}=-m\&CenterDot;\stackrel{\&CenterDot;}{v}$

$f_i^{eul}=-m\&CenterDot;\stackrel{\&CenterDot;}{\mathrm{\&omega;}}\&times;{\tilde{x}}_i\left(q\right),f_i^{cen}=-m\&CenterDot;\mathrm{\&omega;}\&times;(\mathrm{\&omega;}\&times;{\tilde{x}}_i(q\left)\right)$

其中，v,ω分别表示非惯性系的线速度和角速度，m表示顶点质量，U_i表示顶点x_i的降维基底，是非惯性系下顶点x_i的位置。根据上式计算合力f_i ^fic＝f_i ^cor+f_i ^ine+f_i ^eul+f_i ^cen，则非惯性系中基于样例的降维空间下的运动方程为：

$\stackrel{\&CenterDot;\&CenterDot;}{q}+U^{T}DU\stackrel{\&CenterDot;}{q}+U^T{f}^{\mathrm{int}}\left(Uq\right)+U^T{f}^e\left(Uq\right)=U^T(R^T{f}^{ext}+f^{fic}(\stackrel{\&CenterDot;}{v},\mathrm{\&omega;},\stackrel{\&CenterDot;}{\mathrm{\&omega;}},q,\stackrel{\&CenterDot;}{q}\left)\right)$

处理完形变体仿真之后，需要计算刚性运动部分，根据刚体仿真方程：

$m\stackrel{\&CenterDot;}{v}+D_{v}v=f^{ext}$

$I\stackrel{\&CenterDot;}{\mathrm{\&omega;}}+D_{\mathrm{\&omega;}}\mathrm{\&omega;}={\mathrm{\&tau;}}^{ext}$

其中，m,I为刚体的质量和惯性矩阵，f^ext,τ^ext分别为外力和外力力矩，D_v,D_ω分别为刚体的线速度阻尼矩阵和角速度阻尼矩阵。

3-2Laplace-Beltrami形状空间基于样例的刚性变换

在仿真过程中，由于引入了非惯性系，整个基于样例的仿真也要投影到非惯性系下求解样例引导力，因此在仿真的每个时间步需要将形变位移投影到非惯性系下，利用局部位移计算Laplace-Beltrami形状空间的表示，在非惯性系下求解目标形状，并将重建的网格结构重新投影到全局坐标系下。

本发明采用了多个模型进行实验，这里展示一组实验结果。如图6所示，第一行表示两个样例形状，第二行表示仿真过程中物体向着样例方向形变结果。

以上实施例仅用以说明本发明的技术方案而非对其进行限制，本领域的普通技术人员可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明的精神和范围，本发明的保护范围应以权利要求书所述为准。

Claims (1)

1.一种在Laplace-Beltrami形状空间基于样例的弹性材料的实时仿真方法，其特征在于：该方法包括预处理步骤和仿真步骤：

一、预处理步骤包括如下内容：

1.1利用Laplace-Beltrami算子以及线性有限元离散化方法计算仿真物体和样例在三维四面体网格上的特征函数和特征值，特征函数表示为特征值对应的特征向量，将这些特征函数归一化；在预处理过程中求仿真物体和样例的特征函数过程是：将物体的空间结构采用四面体网格离散化方式描述，并利用线性有限元形函数计算体网格上的特征值和特征向量，求得的特征值由小到大排列，并利用归一化后的特征向量张成的形状空间作为形状插值的计算空间；

1.2将样例的特征函数和仿真物体的特征函数进行配准，通过指定对应区域，使得仿真物体与样例之间具有相同的特征函数排序；

1.3将样例形状和仿真物体投影到各自的特征函数上，得到对应的在Laplace-Beltrami形状空间中的形状描述；

1.4利用模态分析方法，求得仿真物体非线性的模态基底作为全局坐标与降维空间下的投影矩阵；并将模态基底作为输入，利用准静态仿真方法，预先求得仿真过程中计算形变能量以及内力的局部的优化积分单元及对应的权重；

二、仿真步骤包括如下内容：

2.1将仿真物体的位置投影到Laplace-Beltrami形状空间，利用非线性形状插值求解样例流行上的目标形状，求解目标形状的过程包括两个部分：

a.首先计算用于求解目标形状所需要的各个样例的权重；

b.利用计算得出的权重求解目标形状；即在Laplace-Beltrami形状空间通过非线性插值求解一个目标插值形状，优化过程是求使得仿真物体在样例流行上的形状相对于所有样例形变能量的加权之和最小的形状；为加快求解形变能量的过程，利用特征函数以及模态基底，在预处理过程中建立Laplace-Beltrami形状空间与模态基底子空间的联系，直接将Laplace-Beltrami形状空间投影到模态降维子空间，省略每个时间步从Laplace-Beltrami形状空间投影到欧氏空间，再投影到模态降维子空间的过程，并利用不局限于材质模型，局部的优化积分单元代替全局的单元近似求解能量和能量梯度的方法；

2.2获得Laplace-Beltrami形状空间中的目标形状后，利用特征函数投影得到欧氏空间下的形状，该形状只包含整体信息，缺失局部信息，根据上一时刻的形变结果在特征函数上的投影前后的向量差构建细节向量，在当前时刻每个顶点上建立局部坐标系添加形状细节的方式得到欧氏空间下的最终形变结果；

2.3在模态基底构建的降维子空间求解非惯性系下形变体仿真动力学方程，得到新一时刻降维空间下的位移；

2.4利用刚体运动方程计算刚性运动部分，弥补降维仿真方法中，降维基底不能模拟刚性运动的不足；

2.5降维空间下的位移利用预处理过程中计算的全局与降维空间下的投影矩阵，投影得到全局中的位移，并绘制形变结果。