CN105069826A - 弹性物体变形运动的建模方法 - Google Patents

Abstract

本发明公开了一种弹性物体变形运动的建模方法，该方法包括：采集弹性物体的静态点云和动态点云序列；建立仿真四面体网格模型；驱动仿真四面体网格模型跟踪动态点云序列，得到跟踪变形运动序列；迭代估计弹性物体的材质属性系数和对应的参考形状；每个迭代周期均执行以下操作：获得当前材质属性系数对应的参考形状；根据该系数和参考形状，驱使仿真四面体网格模型从相同的初始形变下变形运动，得到仿真变形运动序列；计算仿真变形运动序列与跟踪变形运动序列的位置偏差；沿着位置偏差减小的方向更新材质属性系数；根据位置偏差最小时的材质属性系数和参考形状，建立弹性物体变形运动模型。通过上述建模方法，可以建立逼真的弹性物体变形运动模型。

Description

弹性物体变形运动的建模方法

技术领域

本发明涉及仿真建模技术领域，特别涉及一种弹性物体变形运动的建模方法。

背景技术

目前，近几十年来，计算机图形学技术得到了长足的发展，其研究成果在影视游戏、虚拟仿真、设计制造等领域都得到了广泛的应用。除了处理精细几何模型和渲染逼真效果，采用恰当物体数学物理模型，生成与控制真实感的运动是一个待继续深入研究的问题。

具有物理真实感的运动生成与控制是影视制作、动画领域不可或缺的重要组成部分。但是传统方法中数学模型的过度简化和参数不准确严重损害了仿真结果的精确性，阻碍了这一技术在实际工业领域的广泛应用。

传统的建模方法主要是：借助力传感设备在物体的不同位置施加强度已知的外力，然后通过位置跟踪设备，获得不同外力条件下稳定时物体的形状变化，最后利用机器学习、概率统计等理论方法建立应力与应变的关系曲线，从而获得被测量对象的数学物理模型进行仿真。然而，应用现有弹性物体运动建模方法所建运动模型不够逼真。

发明内容

本发明实施例提供了一种弹性物体变形运动的建模方法，用以建立逼真的弹性物体变形运动模型，该方法包括：

采集弹性物体的静态点云和变形运动过程中的动态点云序列；

根据静态点云，建立用于仿真弹性物体变形运动的仿真四面体网格模型；

驱动仿真四面体网格模型跟踪动态点云序列，得到仿真四面体网格模型的跟踪变形运动序列；

迭代估计弹性物体的材质属性系数和对应的仿真四面体网格模型的参考形状；每个迭代周期均执行以下操作：获得弹性物体的当前材质属性系数对应的仿真四面体网格模型的参考形状；根据当前弹性物体材质属性系数和对应的仿真四面体网格模型的参考形状，驱使仿真四面体网格模型从相同的初始形变下，仿真弹性物体变形运动，得到仿真四面体网格模型的仿真变形运动序列；计算仿真变形运动序列与跟踪变形运动序列的位置偏差；沿着使得位置偏差减小的方向更新弹性物体的材质属性系数；直到找到位置偏差最小时的材质属性系数和对应的仿真四面体网格模型的参考形状；

根据位置偏差最小时的材质属性系数和对应的仿真四面体网格模型的参考形状，建立弹性物体变形运动模型。

在一个实施例中，根据静态点云，建立用于仿真弹性物体变形运动的仿真四面体网格模型，包括：

根据静态点云，建立弹性物体的静态表面网格模型；

根据静态表面网格模型，建立用于仿真弹性物体变形运动的仿真四面体网格模型；

静态表面网格模型的每个顶点与仿真四面体网格模型的每个四面体空间重心坐标是线性插值关系。

在一个实施例中，驱动仿真四面体网格模型跟踪动态点云序列，得到仿真四面体网格模型的跟踪变形运动序列，包括：

找到弹性物体变形后所有动态点云序列与仿真四面体网格模型的所有节点的最大概率对应关系；

根据最大概率对应关系，向仿真四面体网格模型的每个节点施加虚拟外力，驱使仿真四面体网格模型的每个节点跟踪对应的动态点云序列，得到变形后的仿真四面体网格模型的每个节点的位置；

根据线性插值关系和变形后的仿真四面体网格模型的每个节点的位置，找到变形后的静态表面网格模型的每个顶点位置；

计算变形后的静态表面网格模型的每个顶点位置与弹性物体变形后对应的动态点云序列位置之间的差异，直到差异小于预设值时，得到当前点云数据对应的仿真四面体网格模型变形后的跟踪变形运动序列。

在一个实施例中，差异小于预设值包括：静态表面网格模型的每个顶点与弹性物体变形后对应的动态点云序列之间的距离小于预设距离，或它们之间的吸引力小于预设吸引力。

在一个实施例中，获得弹性物体的当前材质属性系数对应的仿真四面体网格模型的参考形状，包括：

验证当前材质属性系数和对应的仿真四面体网格模型的参考形状是否满足物理静力平衡方程。

在一个实施例中，弹性物体的材质属性系数的初值估计，包括：

从弹性物体变形运动时的多个变形模态中选择第一个主要的模态；

确定第一个主要的模态对应的振动频率；

根据第一个主要的模态对应的振动频率与实际采集数据振动频率的匹配程度，确定弹性物体的材质属性系数的初值。

在一个实施例中，求解弹性物体的仿真四面体网格模型的参考形状时，采用如下目标方程求解合力残差最小时对应的形状作为仿真四面体网格模型的参考形状：

$\min F\left(X\right)=\underset{X}{min}\left|\right|RK\left(R^T{x}^s-X\right)-Mg|{|}^2;$

其中，R是旋转矩阵，K是刚度矩阵，x^s是弹性物体的静态形状，X是仿真四面体网络模型的参考形状，M为弹性物体的质量，g为重力加速度；

驱使仿真四面体网格模型变形时，计算施加给每个节点的虚拟弹性力相对于参考形状的雅克比矩阵，雅克比矩阵为：

$\begin{array}{c}\frac{\∂f}{\∂X_{ij}}=V^{\′}{\mathrm{RB}}^{T}EB\left({\mathrm{Rx}}^s-X\right)+{\mathrm{VR}}^{\′}{B}^{T}EB\left({\mathrm{Rx}}^s-X\right)+{{\mathrm{VRB}}^{\′}}^{T}EB\left({\mathrm{Rx}}^s-X\right)+{\mathrm{VRB}}^T{\mathrm{EB}}^{\′}\left({\mathrm{Rx}}^s-X\right)\\ +{\mathrm{VRB}}^{T}EB\left(R^{\′}{x}^s-e_{3\left(j-1\right)+i}\right);\end{array}$

其中，f为施加给仿真四面体网格模型的每个节点的虚拟弹性力，X_ij为仿真四面体网格模型的参考形状中第i个节点在j方向上的位置，x^s是弹性物体的静态形状，X是仿真四面体网格模型的参考形状，V是仿真四面体网格模型的每个四面体元素的体积，E是与弹性物体的材质属性相关的常量矩阵，R是代表仿真四面体网格模型的每个四面体元素的刚性旋转矩阵。

在一个实施例中，估计弹性物体的材质属性系数，包括：

从仿真四面体网格模型的不同位置选出多个节点作为控制点；

根据弹性物体的材质分布属性，估计出不同位置的控制点的不同材质属性系数；

根据线性插值算法和不同位置的控制点的不同材质属性系数，求出仿真四面体网格模型的其他节点的材质属性系数。

与传统方法中，借助力传感设备在物体的不同位置施加强度已知的外力，然后通过位置跟踪设备获得不同外力条件下稳定时物体的形状变化，最后利用机器学习、概率统计等理论方法建立应力与应变的关系曲线，从而获得被测量对象的数学物理模型，进行建模的方法相比较，本发明实施例依赖数据驱动的方法对弹性物体变形运动的点云序列进行跟踪分析建模，具体地，首先，采集弹性物体的静态点云和动态点云序列；其次，根据所述静态点云，建立用于仿真弹性物体变形运动的仿真四面体网格模型；接着，驱动所述仿真四面体网格模型跟踪所述动态点云序列，得到所述仿真四面体网格模型的跟踪变形运动序列；然后，迭代估计弹性物体的材质属性系数和对应的仿真四面体网格模型的参考形状；每个迭代周期均执行以下操作：获得弹性物体的当前材质属性系数对应的仿真四面体网格模型的参考形状；根据当前弹性物体材质属性系数和对应的仿真四面体网格模型的参考形状，驱使所述仿真四面体网格模型从相同的初始形变下，仿真弹性物体变形运动，得到所述仿真四面体网格模型的仿真变形运动序列；计算所述仿真变形运动序列与跟踪变形运动序列的位置偏差；沿着使得位置偏差减小的方向更新弹性物体的材质属性系数；直到找到位置偏差最小时的材质属性系数和对应的仿真四面体网格模型的参考形状；最后，根据位置偏差最小时的材质属性系数和对应的仿真四面体网格模型的参考形状，建立弹性物体变形运动模型。通过上述建模方法，可以建立逼真的弹性物体变形运动模型。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，并不构成对本发明的限定。在附图中：

图1是本发明实施例中弹性物体变形运动的建模方法的流程示意图；

图2是本发明实施例中采集到的弹性物体的点云序列的示意图；

图3是本发明实施例中建立的弹性物体静态表面网格模型的示意图；

图4是本发明实施例中仿真四面体网格模型的示意图；

图5是本发明实施例中对植物模型控制点安放在轴方向位置的示意图；

图6是本发明一个实施例中弹性物体变形运动的建模实施时的示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，下面结合实施方式和附图，对本发明做进一步详细说明。在此，本发明的示意性实施方式及其说明用于解释本发明，但并不作为对本发明的限定。

图1是本发明实施例中弹性物体变形运动的建模方法的流程示意图；如图1所示，该方法包括如下步骤：

步骤101：采集弹性物体的静态点云和变形运动过程中的动态点云序列；

步骤102：根据静态点云，建立用于仿真弹性物体变形运动的仿真四面体网格模型；

步骤103：驱动仿真四面体网格模型跟踪动态点云序列，得到仿真四面体网格模型的跟踪变形运动序列；

步骤104：迭代估计弹性物体的材质属性系数和对应的仿真四面体网格模型的参考形状；每个迭代周期均执行以下操作：获得弹性物体的当前材质属性系数对应的仿真四面体网格模型的参考形状；根据当前弹性物体材质属性系数和对应的仿真四面体网格模型的参考形状，驱使仿真四面体网格模型从相同的初始形变下，仿真弹性物体变形运动，得到仿真四面体网格模型的仿真变形运动序列；计算仿真变形运动序列与跟踪变形运动序列的位置偏差；沿着使得位置偏差减小的方向更新弹性物体的材质属性系数；直到找到位置偏差最小时的材质属性系数和对应的仿真四面体网格模型的参考形状；

步骤105：根据位置偏差最小时的材质属性系数和对应的仿真四面体网格模型的参考形状，建立弹性物体变形运动模型。

本发明实施例提供的弹性物体变形运动的建模方法，首先，采集弹性物体的静态点云和动态点云序列；其次，根据静态点云，建立用于仿真弹性物体变形运动的仿真四面体网格模型；接着，驱动仿真四面体网格模型跟踪动态点云序列，得到仿真四面体网格模型的跟踪变形运动序列；然后，迭代估计弹性物体的材质属性系数和对应的仿真四面体网格模型的参考形状；每个迭代周期均执行以下操作：获得弹性物体的当前材质属性系数对应的仿真四面体网格模型的参考形状；根据当前弹性物体材质属性系数和对应的仿真四面体网格模型的参考形状，驱使仿真四面体网格模型从相同的初始形变下，仿真弹性物体变形运动，得到仿真四面体网格模型的仿真变形运动序列；计算仿真变形运动序列与跟踪变形运动序列的位置偏差；沿着使得位置偏差减小的方向更新弹性物体的材质属性系数；直到找到位置偏差最小时的材质属性系数和对应的仿真四面体网格模型的参考形状；最后，根据位置偏差最小时的材质属性系数和对应的仿真四面体网格模型的参考形状，建立弹性物体变形运动模型。通过上述建模方法，可以建立逼真的弹性物体变形运动模型。

下面就对本发明实施例中提到的各步骤进行详细描述：

在上述步骤101中，我们使用三个Kinect组合，来采集弹性物体的静态点云序列和变形运动过程中的动态点云序列，图2即为采集到的弹性物体的点云序列的示意图，该步骤101即为采集弹性物体的静态点云和动态点云序列的过程。

在上述步骤102中，根据静态点云，建立用于仿真弹性物体变形运动的仿真四面体网格模型，可以包括：

根据静态点云，建立弹性物体的静态表面网格模型；由于该静态表面网格模型十分精确，又可称为静态精细表面网格模型，图3所示即为静态表面网格模型；

根据静态表面网格模型，建立用于仿真弹性物体变形运动的仿真四面体网格模型；图4即为仿真四面体网格模型的示意图；

静态表面网格模型的每个顶点与仿真四面体网格模型的每个四面体空间重心坐标是线性插值关系。

图3中静态表面网格模型有15368个顶点，对应到仿真四面体网格模型上，有9594个节点。

该步骤102即为建立用于仿真弹性物体变形运动的仿真四面体网格模型的过程。

具体实施时，首先，创建物体的精细表面网格模型(如图3所示)。然后将精细表面网格模型的网格传递给体数据生成工具TETGEN，导出用来物理仿真的相对粗糙的四面体网格，它将作为模板用以追踪点云序列(如图4所示)。为了从仿真的四面体元素中得到精细网格的位置信息，这里采用一种嵌入策略，精细网格上的顶点均可用所在四面体空间重心坐标表示，两者是线性插值关系。

在上述步骤103中，驱动仿真四面体网格模型跟踪动态点云序列，得到仿真四面体网格模型的跟踪变形运动序列，可以包括如下步骤：

找到弹性物体变形后所有动态点云序列与仿真四面体网格模型的所有节点的最大概率对应关系；

根据最大概率对应关系，向仿真四面体网格模型的每个节点施加虚拟外力，驱使仿真四面体网格模型的每个节点跟踪对应的动态点云序列，得到变形后的仿真四面体网格模型的每个节点的位置；

根据线性插值关系和变形后的仿真四面体网格模型的每个节点的位置，找到变形后的静态表面网格模型的每个顶点位置；

计算变形后的静态表面网格模型的每个顶点位置与弹性物体变形后对应的动态点云序列位置之间的差异，直到差异小于预设值时，得到当前点云数据对应的仿真四面体网格模型变形后的跟踪变形运动序列。

在一个实施例中，差异小于预设值可以包括：静态表面网格模型的每个顶点与弹性物体变形后对应的动态点云序列之间的距离小于预设距离，或它们之间的吸引力小于预设吸引力。

上述步骤103即为驱动仿真四面体网格模型变形跟踪弹性物体变形运动的过程。

具体实施时，该步骤103为基于物理的概率跟踪方法，运动跟踪需要处理带噪声的点云数据，同时还要考虑遮挡、快速运动和大幅度变形的问题。因此我们将变形运动跟踪转化成一个最大后验概率(MAP)问题，利用期望最大方法(EM)进行迭代求解：(E步)根据当前点云和节点位置寻找最优的对应关系；(M步)移动节点位置使得上述对应关系为最大似然估计。

下面详细说明上述跟踪的过程。

假定c＝c_1:N,1≤n≤N表示采集到的点云中点的坐标，s＝s_1:k,1≤k≤K表示网格节点的位置。我们的任务就是基于给定的点云c求出形状匹配的s。点云中的点和网格上点的关系并不知晓，而是作为一个隐变量Z_kn，表明观测点Cn可能来自于节点s_k，假设Cn符合正态分布，具体为：c_n～N(s_k,∑_k)，协方差矩阵为∑_k＝σ²I，其中，σ为方差，I为单位矩阵。网格变形跟踪匹配点云可以表示成一个最大后验概率的问题：

$s=\arg \underset{S}{\max }p\left(s|c\right);$

此时，运用EM算法进行求解。在E部分，基于对隐变量p(Z_kn|s,c)(p指代的是概率分布)的期望，找到总对数联合概率最低边界logp(s,c)；在M部分，通过调整四面体网格的顶点位置来最大化上一步的最低边界：

$s=\arg \underset{S}{max}\[\log p\left(c|s\right)+\log p\left(s\right)\];$

上式的第二项反映了变形物体模型的势能，因此可以采取物理仿真方式进行优化求解。我们通过向每个节点添加虚拟力：

$f_k={\mathrm{\η\Σ}}_{n}p\left(z_{kn}\right){\mathrm{\Σ}}_k^{-1}(c_n-s_k);$

其中，η为虚拟力的比例系数；

作为外力驱动网格变形去匹配点云形状。这里运动方程为：

$M\stackrel{\·\·}{x}+D\stackrel{\·}{x}+RK(R^{T}x-X)=fext;$

这里采用共旋转的线性有限元模型来仿真变形物体，M是质量矩阵，D＝αM+βK是瑞利阻尼(Rayleighdamping)矩阵，α和β是瑞利阻尼对应的两个系数，R是旋转矩阵，通过对变形梯度做极分解得到，K是刚度矩阵，x是仿真四面体网格模型变形后的形状，X仿真四面体网格模型的参考形状，fext是外力合力。物理仿真部分采用第三方库VEGAFEM封装而成。为了加速计算，这里采用嵌套策略，每次迭代计算将外力映射给节点数量较少的跟踪四面体模型来仿真跟踪运动，然后把节点位移再插值回静态表面网格模型顶点上。运动跟踪过程即为施加外力仿真直至EM迭代收敛。

对于上述提到的运动跟踪部分，传统的技术方案一般采用非刚性配准算法来逐帧匹配模板网格和点云序列，这种方法首先运算速度不够快；其次对点云数据有较高要求，无法处理大变形运动和存在较多噪声干扰项的情形；再次匹配结果因为没有融合物理约束，网格拓扑会走样。而通过上述关于步骤103的记载，本申请技术方案采用基于概率的跟踪算法，找到网格顶点与点云的概率对应关系，并在物理引擎的驱动下驱使网格变形运动，同时匹配结果融合物理约束(满足上述的各个公式方程的约束)，网格拓扑不会走样，有更快的运行速度，更好的跟踪效果，更强的鲁棒性。上述步骤103可以使得建立弹性物体变形运动模型的点云数据噪声小，同时还可以解决遮挡、快速运动和大幅度变形的问题。

在上述步骤104中，获得弹性物体的当前材质属性系数对应的仿真四面体网格模型的参考形状，包括：

验证当前材质属性系数和对应的仿真四面体网格模型的参考形状是否满足物理静力平衡方程。

具体实施时，上述步骤104包括参数估计和参考形状的优化部分，在该部分中，我们通过数据驱动的方法估计未知的静态材质属性系数p＝(E,ν,α,β)和相应的参考形状X。其中E是杨氏模量，ν是泊松系数，α和β是瑞利阻尼对应的两个系数。求解问题可以表示成时空优化问题，下面的目标方程F测量了仿真和跟踪序列的位置偏差：

$\underset{p,X}{min}F(p,X)=\underset{t}{\Σ}\left|\right|x_t-{\hat{x}}_t|{|}^2;$

其中，是跟踪的输出结果，x_t是仿真结果在第t帧的位置。这个时空优化问题高维度，非线性且目标函数是非凸函数，传统方法不能有效地解决。因此我们提出了一种新颖的分治策略交替地迭代优化X和p，即通过该方法可以寻找出最优的材质属性系数和参考形状(即位置偏差最小时的材质属性系数和对应的仿真四面体网格模型的参考形状)。

在算法的每次迭代中，我们首先执行参考形状的优化估计，保证当前的材质属性系数p和其对应的参考形状X满足物理静态平衡的约束，这样该步骤融合了物理约束。接着我们采用标准的向下搜索方法沿着减少轨迹偏差的方向去更新材质属性系数p。我们循环执行该策略直到足够收敛。

在一个实施例中，求解弹性物体的仿真四面体网格模型的参考形状时，采用如下目标方程求解合力残差最小时对应的形状作为仿真四面体网格模型的参考形状：

$\min F\left(X\right)=\underset{X}{min}\left|\right|RK(R^T{x}^s-X)-Mg|{|}^2;$

其中，R是旋转矩阵，K是刚度矩阵，x^s是弹性物体的静态形状，X是的仿真四面体网络模型的参考形状，M为弹性物体的质量，g为重力加速度。

具体实施时，为了重建可信的仿真运动，物体模型的参考形状和静态形状应该区分开来，参考形状应该是不受重力因素影响的，否则最明显的失真现象就是当物体从静止形状开始仿真时，会发生明显的形状改变。拿植物的变形运动举例来说，植物叶子都会受重力影响，有一个先向下耷拉一下的轻微动作，那么仿真运动必然不真实，由于发明人考虑到了这个问题，并且采用了上述目标方程，求解出的合力残差施加给仿真四面体网络模型，得到的仿真四面体网络模型的参考形状。应用该方法，可以生成真实感强的弹性物体变形运动的仿真模型。

具体的推导过程为：由于物体本身重力是唯一导致这两种形状差别的原因，我们去除加速度和速度项，简化运动方程来优化求解模型的参考形状，用x^s表示物体静态形状，那么简化后的力平衡方程为：

RK(R^Tx^s-X)＝Mg；

于是为了得到最小的合力残差，有如下优化目标：

$\min F\left(X\right)=\underset{X}{min}\left|\right|RK(R^T{x}^s-X)-Mg|{|}^2;$

我们应用当前合力残差作为虚拟外力施加给仿真模型，让其通过仿真变形来更新模型参考形状。仿真持续运行直到合力残差足够小并且仿真达到稳定状态。这种方式相比于传统优化方法的好处在于其强健和快速，并且本身变形满足物理约束。

驱使仿真四面体网格模型变形时，计算施加给每个节点的虚拟弹性力相对于参考形状的雅克比矩阵，雅克比矩阵为：

$\begin{array}{c}\frac{\∂f}{\∂X_{ij}}=V^{\′}{\mathrm{RB}}^{T}EB\left({\mathrm{Rx}}^s-X\right)+{\mathrm{VR}}^{\′}{B}^{T}EB\left({\mathrm{Rx}}^s-X\right)+{{\mathrm{VRB}}^{\′}}^{T}EB\left({\mathrm{Rx}}^s-X\right)+{\mathrm{VRB}}^T{\mathrm{EB}}^{\′}\left({\mathrm{Rx}}^s-X\right)\\ +{\mathrm{VRB}}^{T}EB\left(R^{\′}{x}^s-e_{3\left(j-1\right)+i}\right);\end{array}$

其中，f为施加给仿真四面体网格模型的每个节点的虚拟弹性力，X_ij为仿真四面体网格模型参考形状中第i个节点在j方向上的位置，x^s是弹性物体的静态形状，X是仿真四面体网格模型的参考形状，V是仿真四面体网格模型的每个四面体元素的体积，E是与弹性物体的材质属性相关的常量矩阵，R是代表仿真四面体网格模型的每个四面体元素的刚性旋转矩阵。

即，为实现快速稳定的仿真以解决上述优化问题，我们采用隐式求解方法。这是需要计算弹性力相对于参考形状的雅克比矩阵

下面给出详细的推导过程：拟定是四面体网格在未变形时节点的位置，是参考形状的节点位置。共旋转的线性有限元模型的弹性力大小为：

f＝VRB^TE(R^Tx-X)；

其中，V是四面体元素的体积，E是一个6×6的常量矩阵，跟材质的弹性属性相关。矩阵R＝diag(R,R,R,R)是一个分块对角矩阵，R是代表元素的刚性旋转，通过对模型的变形梯度做极分解得到，即F＝RS，S是对称矩阵。对于一个四面体元素，F是矩阵 $\stackrel{\‾}{F}=V_s{V}_m^{-1}=V_s{B}_m;$

左上方3×3子块矩阵，其中， $V_s=\left[\begin{array}{cccc}{x}_1& x_2& x_3& x_4\\ 1& 1& 1& 1\end{array}\right];V_m=\left[\begin{array}{cccc}{X}_1& X_2& X_3& X_4\\ 1& 1& 1& 1\end{array}\right];B_m=V_m^{-1}$ 都为4×4的矩阵。

6×12的矩阵B仅依赖于X，内部由B_m构成，这里用B_ij表示B_m中第i行第j列的元素：

$B=\left[\begin{array}{cccc}\stackrel{\‾}{B_1}& \stackrel{\‾}{B_2}& \stackrel{\‾}{B_3}& \stackrel{\‾}{B_4}\end{array}\right];\stackrel{\‾}{B_l}=\left[\begin{array}{ccc}{B}_{i1}& 0& 0\\ 0& B_{i2}& 0\\ 0& 0& B_{i3}\\ B_{i2}& B_{i1}& 0\\ 0& B_{i3}& B_{i2}\\ B_{i3}& 0& B_{i1}\end{array}\right];$

我们要求的雅克比矩阵为：

k_q代表K的第q列，即为X第q列的偏微分，实际上也是就是节点j第i项组成部分，这里用X_ij表示。应用链式法则，我们有：

$\begin{array}{c}\frac{\∂f}{\∂X_{ij}}=V^{\′}{\mathrm{RB}}^{T}EB\left(Rx-X\right)+{\mathrm{VR}}^{\′}{B}^{T}EB\left(Rx-X\right)+{{\mathrm{VRB}}^{\′}}^{T}EB\left(Rx-X\right)+{\mathrm{VRB}}^T{\mathrm{EB}}^{\′}\left(Rx-X\right)\\ +{\mathrm{VRB}}^{T}EB\left(R^{\′}x-e_{3\left(j-1\right)+i}\right);\end{array}$

这里单引号(')代表偏微分是第q项的标准基。下面是具体每一小项的求导结果：

$V^{\′}=\frac{\∂V}{\∂X_{ij}}=\frac{\∂V}{\∂{\left(V_m\right)}_{ij}}={\mathrm{VB}}_{ji};B^{\′}=\[{\begin{array}{cccc}{\stackrel{\‾}{B_1}}^{\′}& \stackrel{\‾}{{B_2}^{\′}}& {\stackrel{\‾}{B_3}}^{\′}& \stackrel{\‾}{B_4}\end{array}}^{\′}\];{\stackrel{\‾}{B_k}}^{\′}=\frac{\∂\stackrel{\‾}{B_k}}{\∂X_{ij}}=-\left[\begin{array}{ccc}{B}_{ki}{B}_{j1}& 0& 0\\ 0& B_{ki}{B}_{j2}& 0\\ 0& 0& B_{ki}{B}_{j3}\\ B_{ki}{B}_{j2}& B_{ki}{B}_{j1}& 0\\ 0& B_{ki}{B}_{j3}& B_{ki}{B}_{j2}\\ B_{ki}{B}_{j3}& 0& B_{ki}{B}_{j1}\end{array}\right]$

关于R'的计算更加复杂，还是应用链式法则，可以得到：

$R^{\′}=\frac{\∂R}{\∂X_{ij}}\underset{k,l}{\Σ}\frac{\∂R}{\∂F_{kl}}\frac{\∂F_{kl}}{\∂X_{ij}};$

最后，关于上式右边第二项的求导，我们有：

${\stackrel{\‾}{F}}^{\′}=\frac{\∂\stackrel{\‾}{F}}{\∂{\left(V_m\right)}_{ij}}=V_s\frac{\∂B_m}{\∂{\left(V_m\right)}_{ij}}=-V_s\left[\begin{array}{cccc}{B}_{1i}{B}_{j1}& B_{1i}{B}_{j2}& B_{1i}{B}_{j3}& B_{1i}{B}_{j4}\\ B_{2i}{B}_{j1}& B_{2i}{B}_{j2}& B_{2i}{B}_{j3}& B_{2i}{B}_{j4}\\ B_{3i}{B}_{j1}& B_{3i}{B}_{j2}& B_{3i}{B}_{j3}& B_{3i}{B}_{j4}\\ B_{4i}{B}_{j1}& B_{4i}{B}_{j2}& B_{4i}{B}_{j3}& B_{4i}{B}_{j4}\end{array}\right];$

$\frac{\∂F_{kl}}{\∂X_{ij}}=\frac{\∂{\stackrel{\‾}{F}}_{kl}}{\∂(V_{m)ij}}={\left(\stackrel{\‾}{F^{\′}}\right)}_{kl}.$

我们采用滴水观音模型的参考形状。应用该参考形状进行重力的变形仿真可以得到非常准确的静态形状，与采集到的静态形状基本一致。对于参数估计部分，现有的技术一般认为物体模型的参考形状已知或将参考形状等同于静态形状，这样分析的结果在简化后必然不精确。本技术方案通过上述技术方案，能同时优化估计出物体模型的物理参数和对应的参考形状。

在一个实施例中，估计弹性物体的材质属性系数，包括：

从仿真四面体网格模型的不同位置选出多个节点作为控制点；

根据弹性物体的材质分布属性，估计出不同位置的控制点的不同材质属性系数；

根据线性插值算法和不同位置的控制点的不同材质属性系数，求出仿真四面体网格模型的其他四面体元素的材质属性系数。

具体实施时，我们采用上述方法，针对弹性物体的不同位置的不同材质，给不同位置的材质分别估算材质属性系数，那么就可以生成非均一材质的数学物理模型，使得所建仿真模型更加的真实。

具体地，我们需要弹性物体的真实运动轨迹去估计其弹性材质参数和阻尼系数。给模型赋予合适的材质分布，并给予同样的初始位置条件，就可能通过仿真重建一样的运动。考虑到参考形状可以通过上述方法分开计算得到，之前的目标方程可以重写为：

$F\left(p\right)=\underset{t,k}{\Σ}\left|\right|x_{tk}-{\hat{x}}_{tk}|{|}^2;$

t是帧的编号，k是节点编号。为了最小化F(p)，我们采用无梯度的单纯形优化方法求解。另外，由于单一材质的模型不能很好的真实还原物体变形运动，我们引入控制点的概念来解决此问题，通过指定控制点不同的材质属性，模型其他节点通过线性插值也得到不均匀分布的材质属性。如图5所示，对植物模型，我们将控制点安放在轴方向，插值权重用归一化的轴距离表示。对于生成的恐龙模型，控制点根据设想认为指定在不同位置，插值权重采用双调和函数控制。

在一个实施例中，弹性物体的材质属性系数的初值估计，包括：

从弹性物体变形运动时的多个变形模态中选择第一个主要的模态；

确定第一个主要的模态对应的振动频率；

根据第一个主要的模态对应的振动频率与实际采集数据振动频率的匹配程度，确定弹性物体的材质属性系数的初值。

上述“主要的模态”是指最小特征值对应的振动模态。

具体实施时，由于目标函数F(p)通常包含多个局部最小值。因此提供一个合适的输入参数对最终成功优化求解至关重要，即在首次估计材质属性系数时，给出一个最佳值，那么有利于后续快速的计算，提高建模的效率。具体地，我们提出了一种新颖的策略，利用模态分析和坐标下降法来得到合适的初始材质属性系数值。

在模态分析中，小的形变置换被线性化成u＝Φz,其中Φ＝[Φ₁,Φ₂,...,Φ_k]的每一列代表了一个变形模态。其可以通过广义特征值分解求得Kφi＝λ_iMφi。一般将较小特征值对应的特征向量作为基建立降维的模态坐标系对应每个模态的自然频率。

在特征值分解过程中，杨氏模量E会影响振动频率直观理解，材质越软的物体振动频率越小。

我们将真实位移投影到相应模态下从而得到真实的振动频率。如果估计的杨氏模量越接近基准真实值，那么频率差距应该越小。这里采用坐标下降法法去顺序更新每一个材质参数。在每次迭代过程中，我们基于变量变化对目标函数值变化的敏感程度线性搜索。

我们另外还用两个真实模型(杯垫和衣架)来验证我们的发明测量精确度，我们测得的杨氏模量分别为7.0e6和5.6e6。结果和真实实物分别做了静态受力和动态变形的对比。实验结果验证了发明的可行性，测量结果精确度较高。

图6即为本发明一个实施例中弹性物体变形运动的建模实施的示意图，如图6所示，首先采集构建弹性物体的静态形状和动态点云运动序列；接着是本系统的核心，一种交替迭代的优化策略轮流执行变形运动跟踪和参数估计部分，每次迭代运行结果能显著改进另一部分效果。将静态精细模型作为跟踪的模板网格，系统采用基于物理的概率跟踪算法去驱使网格变形配准每一帧点云，并输出四面体网格每一帧节点位置。下一步，参数估计部分同时优化估计材质属性系数，阻尼系数和物体模型参考形状。这里参考形状指的是物体模型不受任何外力包括重力作用时的原始形状，因此物体参考形状和受重力影响的静止形状应有所差别。在优化部分采取分治策略，给定初始估计物理参数后，根据静力平衡方程求解出当前模型参考形状，然后使用这组材质和参考形状数据进行正向仿真，得到在相同初始变形条件下的运动序列；并计算运动序列形状与跟踪结果之间的差异，作为评价该组参数的标准。多次迭代后找到最佳值；最终，系统生成一套能进行真实感交互的仿真模型。

本发明实施例提供的弹性物体变形运动的建模方法可以达到的有益技术效果为：

1)、实验采集平台容易搭建，设备平价，无需额外成本。

2)、运动重建满足物理约束，能处理大幅度变形，对噪声和数据确实有较强的容错性和鲁棒性。

3)、参数估计快速高效，利用分治策略同时计算物理参数和参考形状使得参数估计和后续仿真更加准确逼真。

以上仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明实施例可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种弹性物体变形运动的建模方法，其特征在于，包括：

采集弹性物体的静态点云和变形运动过程中的动态点云序列；

根据所述静态点云，建立用于仿真弹性物体变形运动的仿真四面体网格模型；

驱动所述仿真四面体网格模型跟踪所述动态点云序列，得到所述仿真四面体网格模型的跟踪变形运动序列；

迭代估计弹性物体的材质属性系数和对应的仿真四面体网格模型的参考形状；每个迭代周期均执行以下操作：获得弹性物体的当前材质属性系数对应的仿真四面体网格模型的参考形状；根据当前弹性物体材质属性系数和对应的仿真四面体网格模型的参考形状，驱使所述仿真四面体网格模型从相同的初始形变下，仿真弹性物体变形运动，得到所述仿真四面体网格模型的仿真变形运动序列；计算所述仿真变形运动序列与跟踪变形运动序列的位置偏差；沿着使得位置偏差减小的方向更新弹性物体的材质属性系数；直到找到位置偏差最小时的材质属性系数和对应的仿真四面体网格模型的参考形状；

根据位置偏差最小时的材质属性系数和对应的仿真四面体网格模型的参考形状，建立弹性物体变形运动模型。

2.如权利要求1所述的弹性物体变形运动的建模方法，其特征在于，根据所述静态点云，建立用于仿真弹性物体变形运动的仿真四面体网格模型，包括：

根据所述静态点云，建立弹性物体的静态表面网格模型；

根据所述静态表面网格模型，建立用于仿真弹性物体变形运动的仿真四面体网格模型；

所述静态表面网格模型的每个顶点与所述仿真四面体网格模型的每个四面体空间重心坐标是线性插值关系。

3.如权利要求2所述的弹性物体变形运动的建模方法，其特征在于，驱动所述仿真四面体网格模型跟踪所述动态点云序列，得到所述仿真四面体网格模型的跟踪变形运动序列，包括：

找到弹性物体变形后所有动态点云序列与所述仿真四面体网格模型的所有节点的最大概率对应关系；

根据所述最大概率对应关系，向所述仿真四面体网格模型的每个节点施加虚拟外力，驱使所述仿真四面体网格模型的每个节点跟踪对应的动态点云序列，得到变形后的所述仿真四面体网格模型的每个节点的位置；

根据所述线性插值关系和变形后的所述仿真四面体网格模型的每个节点的位置，找到变形后的静态表面网格模型的每个顶点位置；

计算变形后的静态表面网格模型的每个顶点位置与弹性物体变形后对应的动态点云序列位置之间的差异，直到差异小于预设值时，得到当前点云数据对应的仿真四面体网格模型变形后的跟踪变形运动序列。

4.如权利要求3所述的弹性物体变形运动的建模方法，其特征在于，所述差异小于预设值包括：静态表面网格模型的每个顶点与弹性物体变形后对应的动态点云序列之间的距离小于预设距离，或它们之间的吸引力小于预设吸引力。

5.如权利要求1所述的弹性物体变形运动的建模方法，其特征在于，获得弹性物体的当前材质属性系数对应的仿真四面体网格模型的参考形状，包括：

验证当前材质属性系数和对应的仿真四面体网格模型的参考形状是否满足物理静力平衡方程。

6.如权利要求1所述的弹性物体变形运动的建模方法，其特征在于，弹性物体的材质属性系数的初值估计，包括：

从弹性物体变形运动时的多个变形模态中选择第一个主要的模态；

确定所述第一个主要的模态对应的振动频率；

根据第一个主要的模态对应的振动频率与实际采集数据振动频率的匹配程度，确定弹性物体的材质属性系数的初值。

7.如权利要求1所述的弹性物体变形运动的建模方法，其特征在于，求解弹性物体的仿真四面体网格模型的参考形状时，采用如下目标方程求解合力残差最小时对应的形状作为仿真四面体网格模型的参考形状：

$\min F\left(X\right)=\underset{X}{min}\left|\right|RK\left(R^T{x}^s-X\right)-Mg|{|}^2;$

其中，R是旋转矩阵，K是刚度矩阵，x^s是弹性物体的静态形状，X是仿真四面体网络模型的参考形状，M为弹性物体的质量，g为重力加速度；

驱使仿真四面体网格模型变形时，计算施加给每个节点的虚拟弹性力相对于参考形状的雅克比矩阵，雅克比矩阵为：

$\begin{array}{c}\frac{\∂f}{\∂X_{ij}}=V^{\′}{\mathrm{RB}}^{T}EB\left({\mathrm{Rx}}^s-X\right)+{\mathrm{VR}}^{\′}{B}^{T}EB\left({\mathrm{Rx}}^s-X\right)+{{\mathrm{VRB}}^{\′}}^{T}EB\left({\mathrm{Rx}}^s-X\right)+{\mathrm{VRB}}^T{\mathrm{EB}}^{\′}\left({\mathrm{Rx}}^s-X\right)\\ +{\mathrm{VRB}}^{T}EB\left(R^{\′}{x}^s-e_{3\left(j-1\right)+i}\right);\end{array}$

其中，f为施加给仿真四面体网格模型的每个节点的虚拟弹性力，X_ij为仿真四面体网格模型的参考形状中第i个节点在j方向上的位置，x^s是弹性物体的静态形状，X是仿真四面体网格模型的参考形状，V是仿真四面体网格模型的每个四面体元素的体积，E是与弹性物体的材质属性相关的常量矩阵，R是代表仿真四面体网格模型的每个四面体元素的刚性旋转矩阵。

8.如权利要求1所述的弹性物体变形运动的建模方法，其特征在于，估计弹性物体的材质属性系数，包括：

从仿真四面体网格模型的不同位置选出多个节点作为控制点；

根据弹性物体的材质分布属性，估计出不同位置的控制点的不同材质属性系数；

根据线性插值算法和不同位置的控制点的不同材质属性系数，求出仿真四面体网格模型的其他四面体元素的材质属性系数。