CN107752228B - 基于位置动力学的三维脚模型与鞋楦模型的匹配方法 - Google Patents

Abstract

本发明提供了一种基于位置动力学的三维脚模型与鞋楦模型的匹配方法，涉及虚拟试鞋技术领域，该方法以人脚和鞋楦的三维网格模型为基础，使用基于位置动力学的物理仿真算法，模拟人脚穿鞋的最终形态与受力状况，该方法能够实现自动的脚楦匹配，脚模型的变形符合物理规律，能比较真实的模拟人脚穿鞋的最终形态和受力状况，方法运行效率高、稳定性强。

Description

基于位置动力学的三维脚模型与鞋楦模型的匹配方法

技术领域

本发明涉及虚拟试鞋技术领域，尤其涉及一种基于位置动力学的三维脚模型与鞋楦模型的匹配方法。

背景技术

目前人们在网购平台上选购鞋子时往往是根据自己的脚的长度码数和宽度码数，但是人脚是复杂的三维形状，由于每款鞋子的鞋型不尽相同，因此这样的选择方式往往不能保证鞋子尺寸合脚。

鞋楦是对人脚的抽象，是鞋的母体和制鞋的基础和重要模具，一双鞋穿着是否舒适关键要看鞋楦是否合适。所以通过三维鞋楦模型和脚模型的匹配模拟消费者穿鞋的最终形态，再对匹配的结果进行受力分析，得出合脚性评价，从而帮助人们在网络上选购合适的鞋子应该是一个可行的方案。

现有的脚楦匹配方案有很多种，例如由Li等人提出了一种基于有向包围盒(OBB)和轴向包围盒(AABB)的脚楦模型配准方法。Leon等人提出了用迭代最近点算法(IterativeClosest Point，ICP)处理脚楦匹配问题。冯粮城等人提出了基于轴变形的脚楦匹配方法，将脚楦匹配分为空间位置对准和脚型变形两个主要步骤。李萌坚人提出的脚楦匹配过程综合了以上方法，第一步使用包围盒方法做脚楦模型的初始对齐，第二步使用ICP算法做脚楦迭代配准，第三步使用变形算法对脚做变形处理，针对普通跟高的鞋使用轴变形算法，，针对高跟鞋使用分段变形融合算法达到脚楦匹配的目的。Jia Yu等人利用有限元分析软件ABAQUS，通过对人脚三维有限元模型的不同部位施加不同的力，对人脚穿入高跟鞋进行了物理仿真。以上这几种脚模型与鞋楦模型的匹配方法还存在一些不足，例如基于包围盒技术和最近迭代点算法的脚楦匹配还需人工调整才能达到比较好的匹配效果。基于轴变形和分段融合变形的脚楦匹配方法中，脚模型的变形基于几何规律，变形状态不够真实。而基于有限元的动态仿真效率低，仿真不稳定等。

发明内容

本发明要解决的技术问题是：现有技术中几种脚模型与鞋楦模型的匹配方法存在的不足，例如：需人工调整才能达到比较好的匹配效果、脚楦匹配最终形态中脚模型的变形不够真实、动态仿真效率低，本发明提供了一种基于位置动力学的三维脚模型与鞋楦模型的匹配方法来解决上述问题。

本发明解决其技术问题所采用的技术方案是：一种基于位置动力学(PositionBased Dynamics，PBD)的三维脚模型与鞋楦模型的匹配方法，包括以下步骤：

S1、将三角网格数据格式的脚模型进行四面体有限元网格划分；

S2、将三角网格数据格式的鞋楦模型与所述脚模型放入同一空间坐标系中并进行几何对齐；

S3、对脚模型和鞋楦模型进行初始匹配仿真：

S301、设置初始匹配仿真的算法参数，包括杨氏模量E，泊松比υ，约束方程的求解迭代次数，时间步长△t，约束方程刚度系数α为1；

S302、初始化脚模型顶点p_i的位置x_i，速度v_i和质量m_i

x_i＝X_i,v_i＝0,m_i＝1

其中i表示顶点序号，X_i表示第i个顶点的初始位置坐标；

S303、在每个时间步长△t更新顶点的速度：

v'_i＝v_i+△tw_if_ext (p_i)

v'_i表示顶点p_i新的速度，f_ext(p_i)表示顶点p_i受到的合外力，w_i＝1/m_i表示顶点p_i质量的倒数；

求取顶点p_i的预测位置：

x'_i＝x_i+Δtv'_i

S304、遍历脚模型的所有四面体，第i个四面体四个顶点的预测位置分别为x'_i1，x'_i2，x'_i3，x'_i4；迭代求解约束方程，计算四个顶点的位置改变量Δx_i1，Δx_i2，Δx_i3，Δx_i4，并更新四面体顶点的位置：

x″_ij＝x′_ij+Δx_ij (25)

其中x″_ij表示第i个四面体的第j个顶点更新后的位置；

S305、使用更新后的位置x″_ij更新第i个四面体每个顶点的速度：

v″_ij＝(x″_ij-x_ij)/Δt (26)

其中v″_ij表示第i个四面体第j个顶点的更新速度。

S306、进行连续碰撞检测和响应；

S307、重复步骤S303～S306，直到脚模型的应变能不再变化；

S4、对脚模型和鞋楦模型进行精确匹配仿真：

S401、设置初始匹配仿真的算法参数，包括杨氏模量E，泊松比υ，约束方程的求解迭代次数，时间步长△t，约束方程刚度系数α为0.02；

步骤S402～S407与S302～S307相同；

S5、脚模型和鞋楦模型精确匹配仿真后，通过此时脚模型每个顶点的位置x_i和对应的初始位置X_i，计算得到每个四面体单元的应变张量ε和应力张量σ。

本发明的有益效果是，这种基于位置动力学的三维脚模型与鞋楦模型的匹配方法，该方法以人脚和鞋楦的三维网格模型为基础，使用基于位置动力学的物理仿真算法，模拟人脚穿鞋的最终形态与受力状况，该方法能够实现自动的脚楦匹配，脚模型的变形符合物理规律，能比较真实的模拟人脚穿鞋的最终形态和受力状况，方法运行效率高、稳定性强。

附图说明

下面结合附图和实施例对本发明进一步说明。

图1是本发明的基于位置动力学的脚楦匹配方法的最优实施例的结构示意图。

图2：脚模型的三维示意图；

图3：鞋楦模型的面片显示图和线框显示图；

图4：四面体网格化后的脚模型的剖面图。

图5：鞋楦模型和脚模型的初始状态示意图。

图6：从三个视角观察的经过空间几何对齐的鞋楦模型与脚模型的示意图。

图7：从三个视角观察的脚模型与鞋楦模型初始匹配仿真后的状态示意图。

图8：脚模型与鞋楦模型精确匹配仿真后的状态示意图。

图9：脚模型的应力状态图。

具体实施方式

下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，仅用于解释本发明，而不能理解为对本发明的限制。

在本发明的描述中，需要理解的是，术语“中心”、“纵向”、“横向”、“长度”、“宽度”、“厚度”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”“内”、“外”、“轴向”、“径向”、“周向”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。

此外，术语“第一”、“第二”等仅用于描述目的，而不能理解为指示或暗示相对重要性。在本发明的描述中，需要说明的是，除非另有明确的规定和限定，术语“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。此外，在本发明的描述中，除非另有说明，“多个”的含义是两个或两个以上。

流程图中或在此以其他方式描述的任何过程或方法描述可以被理解为，表示包括一个或更多个用于实现特定逻辑功能或过程的步骤的可执行指令的代码的模块、片段或部分，并且本发明的优选实施方式的范围包括另外的实现，其中可以不按所示出或讨论的顺序，包括根据所涉及的功能按基本同时的方式或按相反的顺序，来执行功能，这应被本发明的实施例所属技术领域的技术人员所理解。

如图1所示，本发明提供了一种基于位置动力学的三维脚模型与鞋楦模型的匹配方法，包括以下步骤：

S1、将三角网格数据格式的脚模型进行四面体有限元网格划分(预处理)。如果输入计算机中的脚模型和鞋楦模型的原始数据不是三角网格格式的，则需要将其它格式的模型数据转换成三角网格数据格式，三角网格数据格式脚模型和鞋楦模型如图2、图3所示。例如可以通过均匀采样得到模型的三维点云信息，再使用泊松曲面重建算法得到三角网格形式的三维模型数据。三角网格形式数据的表示形式是模型表面顶点集合P(p₁,p₂,...,p_i,...,p_a)(其中p_i表示模型的第i个顶点，a表示模型的顶点数)和三角网格集合F(f₁,f₂,...,f_i,...,f_b)(其中f_i表示模型的第i个三角网格，b表示模型的网格数量)，每一个三角网格f包含三个索引值i₁,i₂,i₃，分别表示模型的第i₁,i₂,i₃个顶点。此实例中使用Delaunay四面体网格生成算法对脚模型进行有限元网格划分。脚模型进行有限元网格划分后的模型剖面图如图4所示。

S2、鞋楦模型和脚模型都有自己的坐标系，将三角网格数据格式的鞋楦模型与所述脚模型放入同一空间坐标系中时，两个模型相距较远，它们的初始状态如图5所示，此实例中，使用最近点迭代算法(ICP)对两个模型进行几何对齐。几何对齐状态如图6所示。

S3、对脚模型和鞋楦模型进行初始匹配仿真：

基于位置动力学仿真中，物体的表示形式包括一个顶点集合，脚模型的顶点集合包括模型表面的顶点和四面体有限元划分后插入的顶点，第i个顶点p_i包含的属性值有：质量m_i、位置坐标x_i和速度矢量v_i。物体的表示形式还包括一个基于顶点位置的约束方程C_j(x₁,x₂,...,x_n)，j为约束方程的序号，n表示约束方程的顶点个数。基于应变能建立约束方程的过程如下：

在连续介质力学中，物体的形变是通过连续的位移场u来描述的。当位移用拉格朗日描述法时，物体变形前顶点的坐标X＝(X₁,X₂,X₃)^T为基本未知量，物体变形后的位置为x＝(x₁,x₂,x₃)^T。每个顶点变形前后坐标之间的关系可以用形变函数φ:R³→R³来表示，则有

x＝φ(X)＝X+u (1)

φ(X)的雅克比式为

F称为变形梯度，我们用F来确定非线性格林应变张量

其中I表示3阶单位矩阵。

由广义胡克定律给出的各向同性材料的应力-应变关系，则应力张量

σ＝2με+λtr(ε)I (4)

其中I表示三阶单位矩阵，tr(ε)表示应变张量ε的迹，μ，λ表示拉梅常数，其计算公式为

其中E表示杨氏模量，υ表示泊松比。应变能密度函数函数

则因变形存储在物体内的能量就是在物体空间Ω内对ψ_s积分，计算公式为

E_s＝∫_Ωψ_sdx (7)

为了进行柔性体的数值仿真，我们选用四面体网格和线性拉格朗日形函数将物体离散化。四面体单元的变形梯度

其中

D_s＝(x₁-x₄,x₂-x₄,x₃-x₄) (9)

D_m＝(X₁-X₄,X₂-X₄,X₃-X₄) (10)

由公式(5)可计算存储在四面体单元内的应变能

V表示未变形前的四面体的体积。这样就可以将应变能函数E_s作为基于位置动力学的约束方程。定义应变能约束

C(x)＝E_s(x)＝0 (12)

由公式(2)(3)还需要计算应变能梯度计算公式为

其中x₁，x₂，x₃，x₄为四面体单元四个顶点的位置，V表示四面体单元的体积。以上就是应变能约束方程的建立过程。

S301、设置初始匹配仿真的算法参数，包括杨氏模量E为1Mpa，泊松比υ为0.4，约束方程的求解迭代次数为3，时间步长△t为1.5s，约束方程刚度系数α为1。

S302、初始化脚模型顶点p_i的位置x_i，速度v_i和质量m_i

x_i＝X_i,v_i＝0,m_i＝1

其中i表示顶点序号，X_i表示第i个顶点的初始位置坐标；

S303、在每个时间步长△t更新顶点的速度：

v'_i＝v_i+△tw_if_ext (p_i)

v'_i表示顶点p_i新的速度，f_ext(p_i)表示顶点p_i受到的合外力，w_i＝1/m_i表示顶点p_i质量的倒数；

求取顶点p_i的预测位置：

x'_i＝x_i+△tv'_i

这些预测位置会被修改到重新满足约束方程的位置。基于位置的求解器目标就是计算顶点的预测位置修正量△x，使之满足C(x+△x)＝0。Δx通过求解下述代数方程组得到：

为应变能约束函数的梯度。如果将△x的方向限制在E(x)的梯度方向上，并考虑每个顶点的质量，则顶点p_i的位置修正量

其中λ₂为拉格朗日乘子，w_i＝1/m_i表示第i个顶点质量的倒数。将式(16)带入(15)式可得

其中j表示约束方程中顶点的序号。

S304、遍历脚模型的所有四面体，第i个四面体四个顶点的预测位置分别为x'_i1，x'_i2，x'_i3，x'_i4；迭代求解约束方程，在每次迭代时利用公式(13)(14)和公式(15)(16)(17)计算四个顶点的位置改变量△x_i1，△x_i2，△x_i3，△x_i4，并更新四面体顶点的位置：

x″_ij＝x′_ij+△x_ij

其中x″_ij表示第i个四面体的第j个顶点更新后的位置；

S305、使用更新后的位置x″_ij更新第i个四面体每个顶点的速度：

v″_ij＝(x″_ij-x_ij)/△t

其中v″_ij表示第i个四面体第j个顶点的更新速度。

S306、进行连续碰撞检测和响应：首先判断脚模型中在鞋楦外的点；接着修改这些点的速度。判断一个顶点是否在鞋楦外使用k近邻算法，具体方法是对脚模型上的点p_i(i为顶点的序号)，找到鞋楦上距离该点最近的k个点，按距离由近到远依次排列为p′₁,p′₂,...,p′_k,p′₁的点法线为n，方向指向鞋楦外面，计算鞋楦顶点的法线。顶点p_i和p′₁构成的向量

a＝p_i-p′₁ (18)

如果

a·n>0 (19)

则认为点p_i在鞋楦外面，在修改点p_i的速度为

其中min(a,b)表示两个向量模长中较小的向量。

如果

v·n≤0 (21)

则认为点v_i在鞋楦里面，将点v_i的速度改为0。

S307、重复步骤S303～S306，直到脚模型的应变能不再变化。初始匹配仿真后的的结果如图7所示。

S4、对脚模型和鞋楦模型进行精确匹配仿真：

S401、设置初始匹配仿真的算法参数，包括杨氏模量E为1Mpa，泊松比υ为0.4，约束方程的求解迭代次数位3，时间步长△t为1.5s，约束方程刚度系数α为0.02；

步骤S402～S407与S302～S307相同。精确匹配仿真后的的结果如图8所示。

S5、脚模型和鞋楦模型精确匹配完成后，通过此时脚模型每个顶点的位置x_i和对应的初始位置X_i，根据公式(3)(4)(8)(9)(10)计算得到每个四面体单元的应变张量ε和应力张量σ。使用openGL对应力状态可视化如图9所示。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对所述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

以上述依据本发明的理想实施例为启示，通过上述的说明内容，相关工作人员完全可以在不偏离本项发明技术思想的范围内，进行多样的变更以及修改。本项发明的技术性范围并不局限于说明书上的内容，必须要根据权利要求范围来确定其技术性范围。

Claims (1)

1.一种基于位置动力学的三维脚模型与鞋楦模型的匹配方法，其特征在于，包括以下步骤：

S1、将三角网格数据格式的脚模型进行四面体有限元网格划分；

S2、将三角网格数据格式的鞋楦模型与所述脚模型放入同一空间坐标系中并进行几何对齐；

S3、对脚模型和鞋楦模型进行初始匹配仿真：

S301、设置初始匹配仿真的算法参数，包括杨氏模量E，泊松比υ，约束方程的求解迭代次数，时间步长Δt，约束方程刚度系数α为1；

S302、初始化脚模型顶点p_i的位置x_i，速度v_i和质量m_i

x_i＝X_i,v_i＝0,m_i＝1

其中i表示顶点序号，X_i表示第i个顶点的初始位置坐标；

S303、在每个时间步长Δt更新顶点的速度：

v'_i＝v_i+Δtw_if_ext(p_i)

v'_i表示顶点p_i新的速度，f_ext(p_i)表示顶点p_i受到的合外力，w_i＝1/m_i表示顶点p_i质量的倒数；

求取顶点p_i的预测位置：

x'_i＝x_i+Δtv'_i

S304、遍历脚模型的所有四面体，第i个四面体四个顶点的预测位置分别为x'_i1，x'_i2，x'_i3，x'_i4；迭代求解约束方程，计算四个顶点的位置改变量Δx_i1，Δx_i2，Δx_i3，Δx_i4，并更新四面体顶点的位置：

x″_ij＝x'_ij+Δx_ij (25)

其中x″_ij表示第i个四面体的第j个顶点更新后的位置；

S305、使用更新后的位置x″_ij更新第i个四面体每个顶点的速度：

v″_ij＝(x″_ij-x_ij)/Δt (26)

其中v″_ij表示第i个四面体第j个顶点的更新速度；

S306、进行连续碰撞检测和响应；

S307、重复步骤S303～S306，直到脚模型的应变能不再变化；

S4、对脚模型和鞋楦模型进行精确匹配仿真：

S401、设置初始匹配仿真的算法参数，包括杨氏模量E，泊松比υ，约束方程的求解迭代次数，时间步长Δt，约束方程刚度系数α为0.02；

步骤S402～S407与S302～S307相同；

S5、脚模型和鞋楦模型精确匹配仿真后，通过此时脚模型每个顶点的位置x_i和对应的初始位置X_i，计算得到每个四面体单元的应变张量ε和应力张量σ。