一种软组织局部有限元模型的快速生成方法与流程

本发明涉及一种人体组织有限元仿真技术，特别涉及一种软组织局部有限元模型的快速生成方法领域。

背景技术：

在传统的手术培训中，新医生完全熟练掌握一项手术技能的时间长达几年，而这对手术的成功率具有较大的影响。直到虚拟手术培训系统的诞生，上述弊端才能够得以解决，用户通过操作虚拟手术培训系统，进行手术训练，大大减少了手术培训周期和培训成本。在虚拟手术培训系统的构建中存在两个基本要求，即真实性和实时性。另一方面，人体软组织具有复杂的力学特性，利用有限元法能够精确地预测人体组织的变形，但在其求解过程中，往往对包括载荷区域稍远的整个模型进行迭代求解，导致其计算规模庞大，难以满足虚拟手术培训系统实时性的要求。

对现有技术的文献检索发现，闫桢南在上海交通大学的硕士学位论文《软组织非线性形变快速模型与应用》中，基于有限元法将软组织模型分为操作区域和非操作区域，并只对操作区域的形变进行计算的方法来减少软组织模型计算规模，从而达到加速仿真的效果。而这种方法需要用户在模拟手术演练前，根据自己以往的经验，对软组织的模型进行手动划分操作区域，因此缺乏灵活性，尤其是对手术路径较长的情况，其加速效果并不理想。

综上所述，在软组织实时仿真领域中，生成局部的形变模型是提高计算速度的有效手段之一，并也是一项有待于克服的难题。

技术实现要素：

本发明的目的在于为了克服现有技术中的不足，提供一种软组织局部有限元模型的快速生成方法，能够同时实现虚拟手术培训系统的真实性和实时性。

本发明为解决上述技术问题采取的技术方案是：

一种软组织局部有限元模型的快速生成方法，其特征在于包括以下步骤：

步骤1）：根据软组织的网格模型，创建一个map容器，map容器的关键字定义为int64型变量A，map容器关键字对应的值定义为结构体B，结构体B中包括两个int型变量，分别为B1和B2。将64位的变量A分为三个部分，分别为0位到20位的A1部分、21位到41位的A2部分以及42位到63位的A3部分。根据软组织网格模型的每个三角形单元所组成的节点序号信息，在每个map容器的迭代器中，int64型变量A的不同部分A1、A2、A3中分别存放每个三角单元的三个节点序号的二进制数值，以由小到大的顺序，A1中存放的节点序号值为最小，A3中存放的节点序号值最大，实现通过一个int64型变量能够表达一个三角形单元。而在每个map容器的迭代器中，结构体B包含的int型变量B1和B2分别存放共用变量A所对应三角形单元的两个相邻四面体单元序号，如果变量A所对应的三角形单元为软组织模型的边界单元，则只在B1中存放包含变量A所对应三角形单元的四面体单元，并将变量B2赋值为-1。

步骤2）：筛选参与计算的单元的集合，以受力点所在单元的内心为基准点，给出如下所示的筛选条件的关系式：

(1)

式中，l为当前四面体单元的内心与受力点所在单元的内心之间的距离，a为筛选系数，R为受力单元的外接圆半径；

判断当前的四面体单元是否满足筛选条件关系式（1），如果是，将当前四面体单元列入参与计算单元的集合，判断首先从受力点所在四面体单元的相邻四面体单元开始进行，然后再对每个相邻四面体单元相邻的四面体单元进行判断，参与计算单元集合的判断以受力点所在四面体单元为中心，由里到外进行，直到不满足筛选条件关系式（1）为止。筛选参与计算的单元时，根据筛选判断的先后顺序，将参与计算的单元集合的四面体单元序号和节点序号进行重新排序。当受力点所在四面体单元个数大于1时，依次进行受力点所在四面体单元个数次的循环进行判断，并统一封装到参与计算单元的集合。每个四面体单元的相邻四面体单元，根据每个四面体单元所组成的节点，生成四个三角形面的节点组合，节点序号由小到大的顺序排列，然后通过步骤1）中生成的map容器查找。

步骤3）：创建数组C，将步骤2）中生成的参与计算单元集合的新四面体单元序号和原模型中的四面体单元序号相连。创建数组D，将步骤2）中生成的参与计算单元集合的新节点序号和原模型中的节点序号相连。

步骤4）：基于步骤2）所生成的参与计算的单元集合，并结合步骤3）中所生成的数组C和数组D，封装整体刚度矩阵，在施加的外加载荷和边界条件下，只对软组织模型中所筛选出的参与计算单元的集合进行形变计算。

上述的一种软组织局部有限元模型的快速生成方法，所述步骤3)中的数组C为int型的一维数组，数组C的下标表示原模型的四面体单元序号，每个参与计算单元的新四面体单元序号分别存储到对应的原模型的四面体单元序号下标的数组C的元素中，而原模型的四面体单元序号中，非参与计算单元下标的元素值赋为-1。数组D为int型的一维数组，数组D的下标表示原模型的节点序号，每个参与计算单元的新节点序号分别存储到对应的原模型节点序号下标的数组D的元素中，而原模型节点序号中，非参与计算节点下标的对应元素值赋为-1。

上述的一种软组织有限元模型的加速计算方法，所述步骤4）中软组织形变计算的过程中，软组织的模型采用指数多项式混合形式的超弹性模型，模型的表达式如下：

(2)

式中，W为应变能密度函数，C₁、β为模型的材料常数，I₁和I₂为应变张量的主不变量；

基于软组织的指数多项式混合形式的超弹性模型建立平衡方程，然后施加边界条件和外加载荷，可进行形变计算。在进行形变计算的过程中，平衡方程采用完全拉格朗日法，进行分步载荷求解，每一载荷步的非线性方程组采用修正的牛顿迭代法进行求解。在软组织有限元计算的过程中，整体刚度矩阵采用GPU计算，非线性方程组采用CPU进行求解，来进行进一步地提高计算效率。

本发明的有益效果是：

本发明基于软组织的指数多项式混合形式的超弹性模型，并结合有限元计算方法保证了软组织仿真的真实性。同时通过算法改进，建立软组织的局部有限元模型，减少了不必要的计算，从而提高了计算的速度。用户通过合理地选取筛选系数a的值，能够有效地减少系统的计算规模，同时能够真实地仿真出人体软组织的生物力学特性。

在建立软组织局部有限元模型的过程中，用map容器和int64位变量的结合，实现了搜索相邻四面体单元的功能，相比于数组的存储方法进一步地加快了搜索相邻四面体单元的速度，并同时提高了软组织局部有限元模型的生成速度。

附图说明

图1是本发明流程框图，即软组织加速计算流程图；

图2是筛选参与计算单元集合的判断示意图；

图3是本发明应用的平面示意图。

具体实施方式

具体实施方式一：如图1所示，本实施方式所述的一种软组织局部有限元模型的快速生成方法，包括以下步骤：

步骤1）：根据软组织的网格模型，创建一个map容器，map容器的关键字定义为int64型变量A，map容器关键字对应的值定义为结构体B，结构体B中包括两个int型变量，分别为B1和B2。将64位的变量A分为三个部分，分别为0位到20位的A1部分、21位到41位的A2部分以及42位到63位的A3部分。根据软组织网格模型的每个三角形单元所组成的节点序号信息，在每个map容器的迭代器中，int64型变量A的不同部分A1、A2、A3中分别存放每个三角单元的三个节点序号的二进制数值，以由小到大的顺序，A1中存放的节点序号值为最小，A3中存放的节点序号值最大，实现通过一个int64型变量能够表达一个三角形单元。而在每个map容器的迭代器中，结构体B包含的int型变量B1和B2分别存放共用变量A所对应三角形单元的两个相邻四面体单元序号，如果变量A所对应的三角形单元为软组织模型的边界单元，则只在B1中存放包含变量A所对应三角形单元的四面体单元，并将变量B2赋值为-1。本步骤在离线状态下实现，因此对软组织计算速度无影响。

步骤2）：筛选参与计算的单元的集合，以受力点所在单元的内心为基准点，给出如下所示的筛选条件的关系式：

(1)

式中，l为当前四面体单元的内心与受力点所在单元的内心之间的距离，a为筛选系数，R为受力单元的外接圆半径。

判断当前的四面体单元是否满足筛选条件关系式（1），如果是，将当前四面体单元列入参与计算单元的集合，判断首先从受力点所在四面体单元的相邻四面体单元开始进行，然后再对每个相邻四面体单元相邻的四面体单元进行判断。如图2所示，受力点所在的四面体单元包含i、j、k和p四个点，并与相邻的四面体单元共用由i、k和p组成的三角形单元，点n1到n2的距离即为公式(1)中的，并且满足时，可判断该相邻四面体单元为参与计算的单元。参与计算单元集合的判断以受力点所在四面体单元为中心，由里到外进行，直到不再满足筛选条件关系式（1）为止。筛选参与计算的单元时，根据筛选判断的先后顺序，将参与计算的单元集合的四面体单元序号和节点序号进行重新排序。当受力点所在四面体单元个数大于1时，依次进行受力点所在四面体单元个数次的循环进行判断，并统一封装到参与计算单元的集合。每个四面体单元的相邻四面体单元，根据每个四面体单元所组成的节点，生成四个三角形面的节点组合，节点序号由小到大的顺序排列，然后通过步骤1）中生成的map容器查找。

步骤3）：创建数组C，将步骤2）中生成的参与计算单元集合的新四面体单元序号和原模型中的四面体单元序号相连，创建数组D，将步骤2）中生成的参与计算单元集合的新节点序号和原模型中的节点序号相连，能够实现快速查询筛选出的参与计算的四面体单元新序号和节点新序号与原模型对应的原序号，并能够保证实时更新发生形变的节点坐标。

步骤4）：基于步骤2）所生成的参与计算的单元集合，并结合步骤3）中所生成的数组C和数组D，封装整体刚度矩阵，实现根据施加的外加载荷和边界条件，只对软组织模型中所筛选出的参与计算单元的集合进行形变计算。

具体实施方式二：本实施方式所述的一种软组织局部有限元模型的快速生成方法，所述步骤3)中的数组C为int型的一维数组，数组C的下标表示原模型的四面体单元序号，每个参与计算单元的新四面体单元序号分别存储到对应的原模型四面体单元序号下标的数组C的元素中，而原模型四面体单元序号中，非参与计算的单元下标的对应元素值赋为-1。例如，参与计算单元集合中序号为2的四面体单元在原模型中的序号为7（假设四面体单元序号从0开始），则C[7]赋值为2，如果原模型中序号为7的四面体单元不是参与计算单元集合中的四面体单元，则C[7]赋值为0。数组D为int型的一维数组，数组D的下标表示原模型的节点序号，每个参与计算单元的新节点序号分别存储到对应的原模型节点序号下标的数组D的元素中，而原模型节点序号中，非参与计算节点下标的元素值赋为0。例如，参与计算单元集合中序号为2的节点在原模型中的序号为7（假设节点序号从0开始），则D[7]赋值为2，如果原模型中序号为7的节点不在参与计算单元集合中的四面体单元上，则D[7]赋值为-1。

具体实施方式三：本实施方式所述的一种软组织局部有限元模型的快速生成方法，所述步骤4）中，如图3所示，只对软组织的有限元模型中参与计算单元集合，即图3中的灰色区域进行计算。软组织的模型采用指数多项式混合形式的超弹性模型，模型的表达式如下：

(2)

式中，W为应变能密度函数，C₁、β为模型的材料常数，I₁和I₂为应变张量的主不变量；

基于软组织的指数多项式混合形式的超弹性模型建立平衡方程，然后施加边界条件和外加载荷，可进行形变计算。在进行形变计算的过程中，平衡方程采用完全拉格朗日法，进行分步载荷求解，每一载荷步的非线性方程组采用修正的牛顿迭代法进行求解，而在求解非线性方程组的过程中产生的线性方程组可采用CG法进行求解。在软组织有限元计算的过程中，整体刚度矩阵采用GPU计算，非线性方程组采用CPU进行求解，来进行进一步地提高计算效率。