非刚性三维形状逐点对应方法及人体心脏运动仿真方法

摘要

本发明公开了一种非刚性三维形状逐点对应方法，包括获取待对应的三维网格；计算三维网格的拉普拉斯矩阵；计算拉普拉斯矩阵的广义特征值分解得到对应的特征向量和特征值；选择紧框架小波滤波器并生成对应的滤波器；迭代优化函数映射和逐点映射矩阵得到待对应的三维网格之间的逐点对应。本发明还公开了包括所述非刚性三维形状逐点对应方法的人体心脏运动仿真方法。本发明使用多尺度小波的保值作为约束，约束更加简单、紧凑和有效，而且迭代收敛速度更快，鲁棒性更高；同时，本发明方法的对应效果更好，而且计算效率更高。

权利要求

1.一种人体心脏运动仿真方法，其特征在于具体包括如下步骤：(1).获取两个不同时刻的运动心脏形状 和 作为待对应的三维网格，并记 为形状 的顶点坐标， 为形状 的顶点连接关系， 为形状 的顶点坐标， 为形状 的顶点连接关系；

(2).按照非刚性三维形状逐点对应方法，计算运动心脏形状 和 之间的对应关系；

所述对应关系使用逐点映射矩阵P标记；

(3).利用运动心脏形状 和 之间的对应关系P，对运动心脏形状 重新网格化标记为 的顶点为 顶点连接关系为 使得 与 具有相同的连接关系；

(4).利用测地距离插值技术，插值得到 和 之间的一系列形状 记R为形状序列的总数；通过求解下列优化目标函数得到 和 之间的一系列形状

式中，τ为时间间隔，且 R为形状序列的总数； 为时间间隔为τ的形状序列； 为形状 和 的弹性能量；

(5).利用步骤(4)生成的 和 的一系列形状 对运动心脏形状 和 之间的心脏运动进行仿真；

所述的非刚性三维形状逐点对应方法，包括如下步骤：S1.获取待对应的三维网格；

S2.计算步骤S1获取的三维网格的拉普拉斯矩阵；

S3.计算步骤S2得到的拉普拉斯矩阵的广义特征值分解，得到对应的特征向量和特征值；

S4.选择紧框架小波滤波器，并生成对应的滤波器；

S5.根据步骤S3得到的特征向量和特征值以及步骤S4生成的滤波器，迭代优化函数映射和逐点映射矩阵，最终得到待对应的三维网格之间的逐点对应；具体包括如下步骤：a.初始化映射矩阵P；

b.根据映射矩阵P，采用如下算式计算对应的函数映射矩阵C：式中Ck为第k次迭代优化并计算得到的函数映射矩阵C； 和 为步骤S4生成得到的滤波器；Pk为第k‑1次迭代优化并计算得到的映射矩阵P； 为包含拉普拉斯矩阵 的前k个特征向量的矩阵； 为包含拉普拉斯矩阵 的前k个特征向量的矩阵；+为矩阵的Moore‑Penrose伪逆；

c.根据步骤b得到的函数映射矩阵C，采用如下算式计算并更新映射矩阵P：式中Pk+1为第k次迭代优化并计算得到的映射矩阵P；Ck为第k次迭代优化并计算得到的函数映射矩阵C；|| ||F为Frobenius范数；

d.重复步骤b～c，直至达到设定的要求，从而得到最终的函数映射矩阵C和逐点映射矩阵P。

2.根据权利要求1所述的人体心脏运动仿真 方法，其特征在于步骤S1所述的获取待对应的三维网格，具体为获取三角网格 和 其中三角网格 具有m个顶点，三角网格具有n个顶点。

3.根据权利要求1或2所述的人体心脏运动仿真 方法，其特征在于步骤S2所述的计算步骤S1获取的三维网格的拉普拉斯矩阵，具体为采用如下算式计算三角网格 和 的拉普拉斯矩阵 和式中 为三角网格 的拉普拉斯矩阵； 为三角网格 的拉普拉斯矩阵； 为三角网格 的各顶点局部面积作为对角元素形成的对角矩阵； 为三角网格 的余切权重矩阵； 为三角网格 的各顶点局部面积作为对角元素形成的对角矩阵； 为三角网格 的余切权重矩阵。

4.根据权利要求3所述的人体心脏运动仿真 方法，其特征在于步骤S3所述的计算步骤S2得到的拉普拉斯矩阵的广义特征值分解，得到对应的特征向量和特征值，具体包括如下步骤：A.计算三角网格 的拉普拉斯矩阵 的广义特征值分解；

B.计算三角网格 的拉普拉斯矩阵 的广义特征值分解；

C.令矩阵 包含拉普拉斯矩阵 的前k个特征向量，矩阵 包含拉普拉斯矩阵 的前k个特征向量；

D.标记对角矩阵 包含拉普拉斯矩阵 的前k个特征值，对角矩阵 包含拉普拉斯矩阵 的前k个特征值。

5.根据权利要求4所述的人体心脏运动仿真 方法，其特征在于步骤S4所述的选择紧框架小波滤波器，并生成对应的滤波器，具体为选择Meyer紧框架小波滤波器g(λ)，同时设置尺度参数L＝6，并生成一组滤波器 和

说明书

非刚性三维形状逐点对应方法及人体心脏运动仿真方法

技术领域

[0001] 本发明属于计算机视觉技术领域，具体涉及一种非刚性三维形状逐点对应方法及人体心脏运动仿真方法。

背景技术

[0002] 建立非刚性形状之间的对应，在形状比较、纹理迁移和形状插值等技术中具有广泛应用。因此，建立非刚性形状之间的对应，是计算机图形学、计算机视觉和模式识别等领域中重要而基础的问题。现实世界中，三维形状的变形大部分为近似等距，因此寻找模型之间逐点近似等距对应，就成为了形状对应研究领域中的核心问题。

[0003] 现阶段最具影响力的建立非刚性三维形状近似等距对应的技术，为Maks等于2012年提出的函数映射框架。不同于以往直接建立形状逐点对应的方法，该框架首先构造了一种函数映射算子，用于映射定义在不同三维形状上的平方可积函数；然后，利用特殊函数的对应求解出高质量的逐点对应。该框架最大的优势在于：将逐点对应归结为一个非常简单高效的代数问题，从而极大降低了直接建立逐点对应时的问题复杂度。

[0004] 在该函数框架中，最为关键的步骤是选取一组合适的函数保值约束，从而形成求解基函数映射矩阵的线性方程组。这些函数保值约束通常编码了逐点映射过程中三维形状保持不变的几何特性，如三维形状的几何性质或外观等。因此，这些函数保值约束的选取，就直接决定了函数映射以及逐点对应的质量。

[0005] 目前，在等距形状对应中，广泛使用的函数约束为形状描述子、拉普拉斯‑贝尔特米算子的交换性以及流形优化等。然而，实际上获得大量信息丰富且线性无关的描述子函数通常较为困难，且这些描述子函数的特征提取的准确性往往不能充分体现变形模型的等距不变性。因此，现有的这些基于描述子函数建立非刚性形状对应的主流方法，在对应的质量和计算效率上仍存在较大不足。

发明内容

[0006] 本发明的目的之一在于提供一种对应效果好且计算效率高的非刚性三维形状逐点对应方法。

[0007] 本发明的目的之二在于提供一种包括了所述非刚性三维形状逐点对应方法的人体心脏运动仿真方法。

[0008] 本发明提供的这种非刚性三维形状逐点对应方法，包括如下步骤：

[0009] S1.获取待对应的三维网格；

[0010] S2.计算步骤S1获取的三维网格的拉普拉斯矩阵；

[0011] S3.计算步骤S2得到的拉普拉斯矩阵的广义特征值分解，得到对应的特征向量和特征值；

[0012] S4.选择紧框架小波滤波器，并生成对应的滤波器；

[0013] S5.根据步骤S3得到的特征向量和特征值以及步骤S4生成的滤波器，迭代优化函数映射和逐点映射矩阵，最终得到待对应的三维网格之间的逐点对应。

[0014] 步骤S1所述的获取待对应的三维网格，具体为获取三角网格 和 其中三角网格 具有m个顶点，三角网格 具有n个顶点。

[0015] 步骤S2所述的计算步骤S1获取的三维网格的拉普拉斯矩阵，具体为采用如下算式计算三角网格 和 的拉普拉斯矩阵 和

[0016]

[0017]

[0018] 式中 为三角网格 的拉普拉斯矩阵； 为三角网格 的拉普拉斯矩阵；为三角网格 的各顶点局部面积作为对角元素形成的对角矩阵； 为三角网格的余切权重矩阵； 为三角网格 的各顶点局部面积作为对角元素形成的对角矩阵；
为三角网格 的余切权重矩阵。

[0019] 步骤S3所述的计算步骤S2得到的拉普拉斯矩阵的广义特征值分解，得到对应的特征向量和特征值，具体包括如下步骤：

[0020] A.计算三角网格 的拉普拉斯矩阵 的广义特征值分解；

[0021] B.计算三角网格 的拉普拉斯矩阵 的广义特征值分解；

[0022] C.令矩阵 包含拉普拉斯矩阵 的前k个特征向量，矩阵包含拉普拉斯矩阵 的前k个特征向量；

[0023] D.标记对角矩阵 包含拉普拉斯矩阵 的前k个特征值，对角矩阵 包含拉普拉斯矩阵 的前k个特征值。

[0024] 步骤S4所述的选择紧框架小波滤波器，并生成对应的滤波器，具体为选择Meyer紧框架小波滤波器g(λ)，同时设置尺度参数L＝6，并生成一组滤波器 和

[0025] 步骤S5所述的根据步骤S3得到的特征向量和特征值以及步骤S4生成的滤波器，迭代优化函数映射和逐点映射矩阵，最终得到待对应的三维网格之间的逐点对应，具体包括如下步骤：

[0026] a.初始化映射矩阵P；

[0027] b.根据映射矩阵P，采用如下算式计算对应的函数映射矩阵C：

[0028]

[0029] 式中Ck为第k次迭代优化并计算得到的函数映射矩阵C； 和 为步骤S4生成得到的滤波器；Pk为第k‑1次迭代优化并计算得到的映射矩阵P； 为包含拉普拉斯矩阵 的前k个特征向量的矩阵； 为包含拉普拉斯矩阵 的前k个特征向量的矩阵；+为矩阵的Moore‑Penrose伪逆；

[0030] c.根据步骤b得到的函数映射矩阵C，采用如下算式计算并更新映射矩阵P：

[0031]

[0032] 式中Pk+1为第k次迭代优化并计算得到的映射矩阵P；Ck为第k次迭代优化并计算得到的函数映射矩阵C；|| ||F为Frobenius范数；

[0033] d.重复步骤b～c，直至达到设定的要求，从而得到最终的函数映射矩阵C和逐点映射矩阵P。

[0034] 本发明还公开了一种包括所述非刚性三维形状逐点对应方法的人体心脏运动仿真方法，具体包括如下步骤：

[0035] (1).获取两个不同时刻的运动心脏形状 和 作为待对应的三维网格，并记为形状 的顶点坐标， 为形状 的顶点连接关系， 为形状 的顶点坐标， 为形状 的顶点连接关系；

[0036] (2).按照上述非刚性三维形状逐点对应方法，计算运动心脏形状 和 之间的对应关系；所述对应关系使用上述逐点映射矩阵P标记；

[0037] (3).利用运动心脏形状 和 之间的对应关系P，对运动心脏形状 重新网格化标记为 的顶点为 顶点连接关系为 使得 与 具有相同的连接关系；

[0038] (4).利用测地距离插值技术，插值得到 和 之间的一系列形状 记R为形状序列的总数；通过求解下列优化目标函数得到：

[0039]

[0040] 式中，τ为时间间隔，且 R为形状序列的总数； 为时间间隔为τ的形状序列； 为形状 和 的弹性能量；

[0041] (5).利用上述步骤生成的 和 的一系列插值形状 对运动心脏形状和 之间的心脏运动进行仿真。

[0042] 本发明提供的这种非刚性三维形状逐点对应方法及人体心脏运动仿真方法，使用多尺度小波的保值作为约束，约束更加简单、紧凑和有效，而且迭代收敛速度更快，鲁棒性更高；同时，本发明方法的对应效果更好，而且计算效率更高。

附图说明

[0043] 图1为本发明的对应方法的方法流程示意图。

[0044] 图2为本发明的对应方法的平均测地误差随着迭代次数的变化曲线示意图。

[0045] 图3为本发明的对应方法在不同数据库上建立对应的性能比较示意图。

[0046] 图4为本发明的对应方法在SHREC'16Partial数据集上与其他算法的CQC曲线对比示意图。

[0047] 图5为本发明的对应方法在完整模型，缺失模型Cuts和Holes上的对应结果示意图。

[0048] 图6为本发明的基于非刚性三维形状逐点对应方法的人体心脏运动仿真方法的流程示意图。

[0049] 图7为本发明的人体心脏运动仿真方法在不同时刻的人体心脏的对应结果示意图。

具体实施方式

[0050] 如图1所示为本发明的对应方法的方法流程示意图。

[0051] 函数映射的定义：

[0052] 给定两个三维形状，即流形 和 记定义在它们表面的平方可积函数空间分别为 和 令T: 表示流形 和 之间的逐点映射，在此基础上进一步定义 和 之间的函数映射算子TF: 该算子将流形 上的函数映射
至流形 上，即TF的像定义为

[0053] 分别给定函数空间 和 的一组正交基 和 则函数映射TF将具有矩阵表示形式C＝(cji),即

[0054]

[0055] 通常选择流形 和 的拉普拉斯‑贝尔特米算子的前k个特征函数分别作为函数空间 和 截断的正交基 和 由于这组基具有类似于经典欧氏空间傅里叶基的调和性质，故内积 也被称为(流形)傅里叶系数。采用截断基后的矩阵C为低维矩阵，即 极大降低了求解的复杂度。

[0056] 矩阵C的计算：

[0057] 矩阵C完全编码了未知逐点映射T的几何特性，因而其求解成为函数映射框架中的关键。实际上，在TF(f)的计算公式中，TF将函数的展开系数从基 线性变化至基下的展开系数。为此，给定流形 和 上的q对对应函数 和现有获得这些对应函数的方法主要为 和 上匹配的描述子函数
或分片信息等，将这些函数在特征基函数 和 下的系数分别作为矩阵F和G(均 )引入正则项，通过求解下列优化函数得到函数映射矩阵C：

[0058]

[0059] 分别为流形 和 的拉普拉斯‑贝尔特米算子，|| ||F为Frobenius范数。

[0060] 逐点对应的计算：

[0061] 流形 和 的原始逐点对应T可从相应的函数映射TF中恢复。具体来说，给定上的任意点x,它在流形 的对应可视为T(x)＝maxyTF(f)(y)，max为取向量的最大值操作，f＝δx为位于点x∈M处的Delta函数。因为<δx,φi>＝φi(x),在离散情形下，将矩阵 和的列分别对应形状 和 的拉普拉斯贝尔特米算子的前k个特征向量 和则这些矩阵的每一行向量即为位于相应点处的Delta函数的傅里叶系数。根据Plancherel's定理，函数傅里叶系数之间的距离等价于函数本身之间的L2距离，因此我们可以利用最近邻搜索寻找矩阵 与 最相近的列向量，从而有效的恢复出逐点对应T。

[0062] 因此，本发明将为函数映射框架构建新的函数约束，该约束基于谱流形小波的内蕴变形不变性，并在此基础上设计出高效的计算函数映射矩阵C的方法，由此获得高质量的形状对应。

[0063] 谱流形小波及其等距变形不变性：

[0064] 为了克服小波在图和流形等不规则空间上的平移和缩放的困难，Hammond等类比了经典小波在频率空间上的特性，采用相似的方法构造了谱图小波。在本发明的框架中，将其推广到流形上并称之为谱流形小波(SMWs，Spectral Manifold Wavelets)。

[0065] 给定谱滤波器g(λ): 该函数为光滑的紧支撑实值函数，则具有尺度s且位于点y处的SMWs定义为 同时，为了捕捉信号的低频信息，定义另一个实值函数h: 作为低频滤波器，该函数满足h(0)＞0，且当λ→∞，有h(λ)→0。类似的，位于点y处的尺度函数定义为

[0066] 谱流形小波具有等距不变性，因此给定具有等距变形关联的形状 和 如果它们之间的逐点对应T是等距的，则对于 我们有 此外，因为小波可以有效地表征形状局部的多尺度邻域拓扑结构，为充分恢复映射T的等距特性，将多尺度谱流形小波的保值约束应用至函数映射中，进一步可得
(为了简单起见，以下将尺度函数 表示为
其小波的缩放参数索引l＝0)，该方程表示形状局部的每个频带结构信息均在映射T中得到了充分的保留。

[0067] 保多尺度小波的函数映射：

[0068] 因多尺度的SMWs可以有效的编码映射的等距变形不变性，受此启发，利用多尺度谱流形小波的保值作为约束，确定相应的函数映射矩阵C，然后利用该函数映射实现高质量的对应。为了更清楚地说明算法细节，接下来将在离散情形下构建我们的算法框架。

[0069] 在离散设定中，形状 和 分别表示为具有m和n个顶点的三角网格。首先计算和 的离散拉普拉斯算子，即拉普拉斯矩阵 和 其中和 是形状 和 各顶点局部面积作为对角元素形成的对角矩阵， 和 为
相应形状的余切权重矩阵。

[0070] 给定矩阵 和 分别包含 和 的前k个特征向量，将映射T: 表示成矩阵 其中如果T(i)＝j，则P(i,j)＝1，否则P(i,j)＝0，i和j是分别为形状 和 顶点的索引。进一步，标记对角矩阵
和 分别包含 和 的前k个特征值，diag为将向量转化为
对角矩阵的操作。

[0071] 现在将函数矩阵投射到相应的傅里叶基函数，可得 其中+表示Moore‑Penrose伪逆，T表示矩阵的转置。由于拉普拉斯矩阵特征向量关于面积权重正交，可得 其中I为单位矩阵,即伪逆

[0072] 给定滤波器g(λ)，可以得到两个谱流形小波矩阵和 它们的每一列向量为具有尺度s、位于相应索引点的SMW。
在实际应用中，尺度参数s同样被采样成一组离散点 因小波矩阵 和 相应的傅里叶系数矩阵为 和 根据Plancherel's定理，我们希望求得的逐
点映射矩阵P以及相应的函数映射矩阵C满足每一个尺度sl的保小波约束条件。该策略可用公式表示为下列优化问题：

[0073]

[0074] 在上述目标函数中，因矩阵P和C均为未知矩阵且它们之间相互影响，为此，我们采用交替更新P和C的方法，近似求解该优化问题。具体来说，首先给定初始矩阵P，然后利用下述优化问题求解出C：

[0075]

[0076] 上述优化问题实际上具有解析解，如果让矩阵C满足：

[0077]

[0078] 然而直接求解满足上述所有条件的矩阵C并不容易实现，因为这种方法需要求解大型线性方程组。但是，如果使用合适的滤波器 使得基于该滤波器所形成的小波函数能够成为流形函数空间里的紧框架，则矩阵C的求解将避免解线性方程组，且仅需要计算简单的矩阵和向量乘法，从而极大降低计算复杂度。

[0079] 下面简要介绍紧框架的相关内容。对于具有n个顶点 的网格，给定一组离散尺度 和合适的滤波器g(λ)和h(λ)，则小波组 (为了简单，以下记h(λ)为g(s0λ)，尺度函数函数为小波基 )，则该小波组能构成形状函数空间的一个标架，且具有边界 和 其中 特别是，如果a＝b，则该小波称为紧框架，进一步的，如果a＝b＝1，则称为Parseval框架。

[0080] 如果由滤波器g(λ)所生成的小波函数组 和 分别能够生成 和 函数空间里的Parseval紧框架，则矩阵C可由下式计算获得：

[0081]

[0082] 上述关于C的算法仅使用了矩阵与向量的乘积，避免了以前方法需要求解大型方程组的问题，计算效率得到大幅度提升。

[0083] 得到更新的矩阵C后，我们将采用最近邻搜索方法恢复出逐点对应矩阵P。

[0084] 也就是计算形状 上第i个顶点在形状 上的映射T(i)，可通过求解下列优化问题

[0085]

[0086] 该 目 标 函 数 可 通 过 最 近 邻 搜 索 方 法 求 解 得 到 P ，并 利 用更新矩阵C。使用这种交替更新的迭代方式逐步提炼矩阵C和P，直至获得理想的对应精度。

[0087] 因此，综上所述，本发明提供的这种非刚性三维形状逐点对应方法，包括如下步骤：

[0088] S1.获取待对应的三维网格；具体为获取三角网格 和 其中三角网格 具有m个顶点，三角网格 具有n个顶点；

[0089] S2.计算步骤S1获取的三维网格的拉普拉斯矩阵；具体为采用如下算式计算三角网格 和 的拉普拉斯矩阵 和

[0090]

[0091]

[0092] 式中 为三角网格 的拉普拉斯矩阵； 为三角网格 的拉普拉斯矩阵；为三角网格 的各顶点局部面积作为对角元素形成的对角矩阵； 为三角网格的余切权重矩阵； 为三角网格 的各顶点局部面积作为对角元素形成的对角矩阵；
为三角网格 的余切权重矩阵；

[0093] S3.计算步骤S2得到的拉普拉斯矩阵的广义特征值分解，得到对应的特征向量和特征值；具体包括如下步骤：

[0094] A.计算三角网格 的拉普拉斯矩阵 的广义特征值分解；

[0095] B.计算三角网格 的拉普拉斯矩阵 的广义特征值分解；

[0096] C.令矩阵 包含拉普拉斯矩阵 的前k个特征向量，矩阵包含拉普拉斯矩阵 的前k个特征向量；

[0097] D.标记对角矩阵 包含拉普拉斯矩阵 的前k个特征值，对角矩阵 包含拉普拉斯矩阵 的前k个特征值；k取值
可以选定为100；

[0098] S4.选择紧框架小波滤波器，并生成对应的滤波器；具体为选择Meyer紧框架小波滤波器g(λ)，同时设置尺度参数L＝6，并生成一组滤波器 和

[0099] S5.根据步骤S3得到的特征向量和特征值以及步骤S4生成的滤波器，迭代优化函数映射和逐点映射矩阵，最终得到待对应的三维网格之间的逐点对应；具体包括如下步骤：

[0100] a.初始化映射矩阵P；

[0101] b.根据映射矩阵P，采用如下算式计算对应的函数映射矩阵C：

[0102]

[0103] 式中Ck为第k次迭代优化并计算得到的函数映射矩阵C； 和 为步骤S4生成得到的滤波器；Pk为第k‑1次迭代优化并计算得到的映射矩阵P； 为包含拉普拉斯矩阵 的前k个特征向量的矩阵； 为包含拉普拉斯矩阵 的前k个特征向量的矩阵；+为矩阵的Moore‑Penrose伪逆；

[0104] c.根据步骤b得到的函数映射矩阵C，采用如下算式计算并更新映射矩阵P：

[0105]

[0106] 式中Pk+1为第k次迭代优化并计算得到的映射矩阵P；Ck为第k次迭代优化并计算得到的函数映射矩阵C；|| ||F为Frobenius范数；

[0107] d.重复步骤b～c，直至达到设定的要求，从而得到最终的函数映射矩阵C和逐点映射矩阵P。

[0108] 图2为本发明的对应方法的平均测地误差随着迭代次数的变化曲线示意图；从图2可以看到，本发明方法得迭代收敛速度更快。

[0109] 图3为本发明的对应方法在不同数据库上建立对应的性能比较示意图；采用常用的标准数据集FAUST、SCAPE、TOSCA测试，并使用累积逐点对应测地误差曲线评价算法的准确度，曲线越高，算法越准确。同时，下表1统计了不同算法的平均测度误差和计算时间。

[0110] 表1不同算法的平均测度误差和计算时间示意表

[0111]

[0112] 从图3和表1可以看到，本发明的对应方法具有更好的准确性和高效性。

[0113] 图4为本发明的对应方法在SHREC'16Partial数据集上与其他算法的累积逐点对应测地误差曲线对比示意图；图5为本发明的对应方法在完整模型，缺失模型Cuts和Holes上的对应结果示意图。从图4和图5可以看到，本发明的对应方法在针对局部缺失模型的处理时，也具有较好的准确性和鲁棒性。

[0114] 如图6所示为本发明的人体心脏运动仿真方法的方法流程示意图：本发明还提供的这种包括所述非刚性三维形状逐点对应方法的人体心脏运动仿真方法，具体包括如下步骤：

[0115] (1).获取两个不同时刻的运动心脏形状 和 作为待对应的三维网格，并记为形状 的顶点坐标， 为形状 的顶点连接关系， 为形状 的顶点坐标， 为形状 的顶点连接关系；

[0116] (2).按照上述非刚性三维形状逐点对应方法，计算运动心脏形状 和 之间的对应关系；所述对应关系使用上述逐点映射矩阵P标记；

[0117] (3).利用运动心脏形状 和 之间的对应关系P，对运动心脏形状 重新网格化标记为 的顶点为 顶点连接关系为 使得 与 具有相同的连接关系；

[0118] (4).利用测地距离插值技术，插值得到 和 之间的一系列形状 记R为形状序列的总数；通过求解下列优化目标函数得到：

[0119]

[0120] 式中，τ为时间间隔，且 R为形状序列的总数； 为时间间隔为τ的形状序列； 为形状 和 的弹性能量；

[0121] (5).利用上述步骤生成的 和 的一系列插值形状 对运动心脏形状和 之间的心脏运动进行仿真。

[0122] 步骤(2)所述的按照上述非刚性三维形状逐点对应方法，计算运动心脏形状 和之间的对应关系，具体包括如下步骤：

[0123] 1).获取两个不同时刻的运动心脏形状 和 作为待对应的三角网格；

[0124] 2).计算步骤1)获取的三维网格的拉普拉斯矩阵；

[0125] 3).计算步骤2)得到的拉普拉斯矩阵的广义特征值分解，得到对应的特征向量和特征值；

[0126] 4).选择紧框架小波滤波器，并生成对应的滤波器；

[0127] 5).根据步骤3)得到的特征向量和特征值以及步骤4)生成的滤波器，迭代优化函数映射和逐点映射矩阵，最终得到待对应的三维网格之间的逐点对应。

[0128] 步骤2)所述的计算步骤1)获取的三维网格的拉普拉斯矩阵，具体为采用如下算式计算三角网格 和 的拉普拉斯矩阵 和

[0129]

[0130]

[0131] 式中 为三角网格 的拉普拉斯矩阵； 为三角网格 的拉普拉斯矩阵；为三角网格 的各顶点局部面积作为对角元素形成的对角矩阵； 为三角网格的余切权重矩阵； 为三角网格 的各顶点局部面积作为对角元素形成的对角矩阵；
为三角网格 的余切权重矩阵。

[0132] 步骤3)所述的计算步骤2)得到的拉普拉斯矩阵的广义特征值分解，得到对应的特征向量和特征值，具体包括如下步骤：

[0133] A).计算三角网格 的拉普拉斯矩阵 的广义特征值分解；

[0134] B).计算三角网格 的拉普拉斯矩阵 的广义特征值分解；

[0135] C).令矩阵 包含拉普拉斯矩阵 的前k个特征向量，矩阵包含拉普拉斯矩阵 的前k个特征向量；

[0136] D).标记对角矩阵 包含拉普拉斯矩阵 的前k个特征值，对角矩阵 包含拉普拉斯矩阵 的前k个特征值。

[0137] 步骤4)所述的选择紧框架小波滤波器，并生成对应的滤波器，具体为选择Meyer紧框架小波滤波器g(λ)，同时设置尺度参数L＝6，并生成一组滤波器 和

[0138] 步骤5)所述的根据步骤3)得到的特征向量和特征值以及步骤4)生成的滤波器，迭代优化函数映射和逐点映射矩阵，最终得到待对应的三维网格之间的逐点对应，具体包括如下步骤：

[0139] a).初始化映射矩阵P；

[0140] b).根据映射矩阵P，采用如下算式计算对应的函数映射矩阵C：

[0141]

[0142] 式中Ck为第k次迭代优化并计算得到的函数映射矩阵C； 和 为步骤4)生成得到的滤波器；Pk为第k‑1次迭代优化并计算得到的映射矩阵P； 为包含拉普拉斯矩阵 的前k个特征向量的矩阵； 为包含拉普拉斯矩阵 的前k个特征向量的矩阵；+为矩阵的Moore‑Penrose伪逆；

[0143] c).根据步骤b)得到的函数映射矩阵C，采用如下算式计算并更新映射矩阵P：

[0144]

[0145] 式中Pk+1为第k次迭代优化并计算得到的映射矩阵P；Ck为第k次迭代优化并计算得到的函数映射矩阵C；|| ||F为Frobenius范数；

[0146] d).重复步骤b)～c)，直至达到设定的要求，从而得到最终的函数映射矩阵C和逐点映射矩阵P。

[0147] 如图7所示为本发明的人体心脏运动仿真方法在不同时刻的人体心脏的对应结果示意图：心脏之间的对应点使用线段连接。通过图7可以看到，本发明的人体心脏运动仿真方法，能够较好的实现不同时刻的人体心脏的运动仿真。