一种基于逆有限元优化分析的连续形变恢复方法与流程

本发明涉及形变仿真的技术领域，具体涉及一种基于逆有限元优化分析的连续形变恢复方法，其为一种基于逆有限元优化分析的为可形变模型重建出符合给定输入目标姿态序列的连续平滑的形变仿真序列的方法。

背景技术：

为可形变模型构造出符合物理规律并且与指定关键帧形态相吻合的形变空间已经成为了研究者们普遍关注的问题。当在可形变模型的运动序列中指定一组姿态为关键帧后，人们通常关注如何得到这个模型的形变空间，使得该形变空间能够在包含关键帧的同时对可形变模型的形变序列进行连续平滑的描述，这种连续形变恢复问题通常被称为关键帧插值或是形状插值。

现有的关键帧插值方法通常可以被分为两类：几何形状插值和基于物理的关键帧插值方法。几何形状插值方法通常基于在相邻的两个几何网格间根据特定标准(如：等距标准)寻找一致性，之后在两者互相对应的图元间创造插值。虽然几何形状插值方法可以在两个给定几何形状之间提供快速且合理的插值，但是通过这类方法得到的插值结果中，动力学效果和真实感是无法恢复的。基于物理的关键帧插值方法通常可以在很少人为干预的情况下提供真实感，这种方法试图创造一个需要最少非物理控制外力来匹配目标姿态或约束的合理运动，并且这类方法通常被转化成一个时空优化问题。然而基于物理的关键帧插值方法通常面临着以下主要问题：较高的计算复杂度，线弹性模型在大形变下的失真，对于关键帧目标姿态的表示通常是具有相同连接属性的网格或是假设关键帧位移的存在以及固定材质信息下对输入姿态幅度的限制。

为了解决上述问题，本发明提出了一种基于逆有限元优化分析的连续形变恢复方法，结合模型化简技术和旋转应变坐标方法对有限元仿真进行化简和扩展。之后通过求解两个时空优化问题分别求解出有限元模型的关键帧位移和符合输入目标姿态序列的连续平滑的形变仿真序列。

技术实现要素：

本发明要解决的技术问题是：针对基于物理的关键帧插值技术面临的问题，提供了一种符合输入目标姿态的连续形变恢复方法，该方法基于有限元方法为可形变模型建模，使用模型化简技术降低计算复杂度，并使用旋转应变坐标方法将其扩展至大形变；该方法引入符号距离场隐式地表示目标姿态以绕过对网格一致性的需求，并在求解时空优化问题中将材质信息作为其中一个参数进行优化以更好的匹配目标姿态。

本发明采用的技术方案为：一种基于逆有限元优化分析的连续形变恢复方法，包括以下四个步骤：

步骤(1)、符号距离场序列的生成：根据输入的同一物体的4d姿态数据，根据不同的输入类型，选择不同的方法生成符号距离场序列以隐式地表达姿态：若输入的是多边形网格文件序列，则使用快速行进法生成符号距离场；若输入的是4dct/mri文件，首先对其进行分割并提取出表面体数据，再利用基于径向基函数的重建方法生成符号距离场序列；

步骤(2)、有限元仿真模型的建立：对初始姿态网格进行六面体有限元离散化，为每个六面体单元建立运动方程并集成全局运动方程，利用模态分析技术进行模型化简，并采用旋转应变坐标技术将线性模型扩展至无失真的大形变；

步骤(3)、有限元模型关键帧位移的生成：定义一个目标函数序列，它可以依照时空顺序测量由有限元位移场表示的模型姿态到各个输入目标姿态的差异，并结合伴随方法和共轭梯度优化方法求解时空优化问题以得到尽可能符合输入姿态序列的关键帧位移；

步骤(4)、连续形变空间的恢复：将关键帧位移投影到旋转应变坐标空间中，并对旋转应变坐标向量进行平滑插值，之后将插值结果通过求解泊松问题重建为连续平滑的全局位移场。定义另一目标函数为有限元模型位移与重建的连续平滑位移场的差异与控制施力量的总和，并通过求解时空优化问题得到最优的形变控制信息，并应用于有限元仿真以得到符合输入姿态序列的连续形变仿真。

本发明的原理在于：

(1)使用有限元方法可以较为真实的对可形变模型进行仿真，通过模型化简可以在较好保持精确度的同时极大降低计算复杂度，旋转应变坐标方法可以在将线弹性模型扩展至大形变的同时保持全局位移与子空间位移向量间的简单映射，便于后续对目标函数梯度的计算。

(2)使用符号距离场表示目标姿态为比较当前模型姿态与目标姿态间的差异提供简单的判断依据，由此定义的目标函数在经过最小化的过程后可以得到符合目标姿态的有限元模型关键帧位移。

(3)在旋转应变坐标空间中对关键帧位移进行插值并重建可以较好地保留形变过程中的旋转信息，使得插值结果更为自然。

(4)由于目标函数具有二次形式，适合使用共轭梯度优化方法进行最小化求解，在求解过程中，使用伴随方法可以快速简单地求解目标函数相对于控制参数的梯度。

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

1、本发明提出的有限元关键帧位移生成方法具有对于输入数据格式的普遍适用性，可以作为其他基于物理的关键帧插值技术的预处理。

2、本发明采用的模型化简技术、伴随方法等对于时空优化问题的求解具有较好的速度优势。

3、本发明将材质信息作为参数进行优化提高了对输入姿态重建的鲁棒性。

附图说明

图1为基于逆有限元优化分析的连续形变恢复方法的处理流程图；

图2为由体数据序列生成对应的符号距离场序列的示意图；

图3为符号距离场中空间点梯度指向的示意图；

图4为不同层次下的有限元网格的示意图，其中，图4(a)为六面体半径是0.0625的有限元网格，图4(b)为六面体半径是0.03125的有限元网格，图4(c)为六面体半径是0.015625的有限元网格；

图5为对大形变处理示意图，其中，在相同条件下，图5(a)为不处理大形变时在第9帧的结果，图5(b)为使用旋转应变坐标方法处理大形变在第9帧的结果，图5(c)为使用共旋方法处理大形变在第9帧的结果，图5(d)为使用旋转应变坐标方法处理大形变在第23帧的结果，图5(e)为使用共旋方法处理大形变在第23帧的结果；

图6为关键帧位移与目标姿态对比的示意图；

图7为对armadillo模型的两个关键帧姿态间恢复连续形变的示意图；

图8为对beam模型恢复连续形变的示意图(从左到右分别是初始模型、目标姿态、关键帧位移、中间帧恢复结果)。

具体实施方式

图1给出了基于逆有限元优化分析的连续形变恢复方法的处理流程，下面结合其他附图及具体实施方式进一步说明本发明。

本发明提供一种基于逆有限元优化分析的连续形变恢复方法，主要步骤介绍如下：

1、生成输入目标姿态对应的符号距离场序列

该方法根据输入的多边形网格序列或是体数据序列生成对应的符号距离场序列。三维符号距离场定义在一个三维空间点的集合之上，并以空间中所有场值为0的空间点的集合表示其内嵌的闭合曲面。符号距离场中某点的场值就是这个场中该点到这个闭合曲面的最近的距离。可由公式表示为，

式中，dt(p)表示着从空间三维点p到有向距离场内嵌的闭合曲面上的最短的欧几里得距离。

如图2所示是通过将心脏的体数据序列生成对应的符号距离场后，使用marchingcubes算法进行等值面提取并显示的结果，marchingcubes算法是三维离散数据场中提取等值面的经典算法。

图3所示是内嵌闭合曲面为一球面的符号距离场中，每一空间点上的梯度方向示意图，展现出符号距离场的重要属性，即符号距离场在空间中某一点的梯度方向指向它定义的闭合曲面上离此空间点最近的一点。从图3可以看到，对于一个球体的符号距离场，空间中的每一点上的梯度方向正好指向它的球心。

2、有限元仿真模型的建立

使用有限元方法对物体进行建模仿真，首先要将连续体离散化为一个个小单元，单元之间通过结点连接，且在边界上连续。本发明采用八结点六面体单元模型对物体进行离散化，以表示物体初始姿态的多边形网格文件为输入，将多边形网格离散化为有限个六面体单元作为近似，其中近似程度可由单元稀疏层次决定，如图4所示，六面体单元越细小，越能更好地近似原多边形网格的形状，然而相应的计算复杂度也会增大。

此后，为每个六面体单元赋予材质信息并建立运动方程，并将所有单元的运动方程集成全局运动方程，

式中代表着相对于静止姿态的物体姿态的离散化表示中的所有n个结点的位移，分别是速度和加速度向量，分别是质量矩阵，刚度矩阵，以及rayleigh阻尼矩阵，并且fext(t)是在时刻t下的控制外力。

由于上述方程自由度3n通常较大，对于大型矩阵的求解运算的计算复杂度相对较高，将会导致后续的时空优化问题求解过程中需要进行多次有限元仿真时计算复杂度巨大。因此本发明采用模态分析技术进行模型化简。

模态分析通过求解如下一般特征化问题得到基底矩阵

kφ＝λmφ(3)式中λ和φ分别是特征值和特征向量。通过将最小的r个特征值集成到一个对角矩阵λ＝diag(λ1,λ2,...,λr)中，并将它们对应的特征向量存储在矩阵w＝(φ1,...,φr)中，且具有w^tmw＝i,w^tkw＝λ的性质，若设则公式(2)化简为

式中z为子空间位移向量，是子空间中的阻尼矩阵投影，fr,ext(t)＝w^tfext(t)代表着外力在子空间的投影。

最后将线弹性模型扩展至大形变，将子空间位移向量z先转化为旋转应变坐标向量后，通过泊松重建将旋转应变坐标向量重建为全局位移向量其中p为常矩阵，为由旋转应变坐标向量计算出的形变梯度。如图5所示，(a)为不处理大形变的结果，(b)、(d)为使用旋转应变坐标方法处理大形变的结果，(c)、(e)为使用共旋方法处理大形变的结果。结果表明旋转应变坐标方法在处理大形变时具有较好的鲁棒性。

3、有限元模型关键帧位移的生成

以步骤1中生成的符号距离场序列作为输入，构造需要最小化的目标函数，定义为有限元模型位移表示的模型姿态在输入的符号距离场序列中的场值和的大小，

式中nkey是关键帧的总数目，表示在时刻ti下，六面体有限元模型内嵌的网格形变后的m个采样顶点的3d坐标组成的向量，x(ti)＝aq(ti)+x0是从六面体体素网格位移q到内嵌网格位移x的映射,是常量的映射矩阵，dist(x)是距离场值构成的向量，qi是单位矩阵的比例缩放。

之后通过共轭梯度优化方法求解时空优化问题对其最小化，其中将有限元仿真中的每一时间步下的fr,ext和λ作为控制参数，目标函数相对于控制参数的梯度由伴随方法求解。最后在最优的材质信息和控制施力场下进行有限元仿真可以得出尽可能符合输入姿态序列的有限元模型关键帧位移如图6所示，上排灰色为给定的目标姿态序列，下排为生成的对应的关键帧位移。

4、连续形变空间的恢复方法

以步骤3中生成的有限元模型关键帧位移作为输入，首先将其转化为旋转应变坐标向量并对旋转应变坐标向量进行插值之后通过泊松重建将插值得到的旋转应变坐标向量恢复为连续平滑的连续位移场q0_tar,q1_tar,...,qnum_tar，num为重建出的位移数目。此后构造需要最小化的目标函数，定义为有限元模型位移与输入的关键帧位移场的偏差的加权和和控制施力量的总和，

式中fr,ext是各个时间步下的控制外力fr,ext组合成的长向量，ri是单位矩阵的比例缩放。

并通过共轭梯度优化方法求解时空优化问题对其最小化，其中目标函数相对于控制参数的梯度由伴随方法求解。最后在最优的材质信息和控制施力场下进行有限元仿真即可得出符合输入姿态序列的连续平滑的形变仿真输出。如图7、8所示为在关键帧位移之间恢复出的连续形变。

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

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