CN103942836A - 三维网格模型四面体化方法 - Google Patents

Abstract

本发明公开了一种三维网格模型四面体化方法，包括以下步骤：(1)模型初始位置预处理；(2)构建体心立方结构；(3)计算节点符号及两端节点符号相反的四面体边的切点；(4)切点移动；(5)边界重四面体化。本发明公开的三维模型四面体化方法，基于体心立方的网格四面体化算法，增加了三维模型主成分分析的预处理，改进了模型边界切点的移动方式转向模型特征点移动，并对最终的四面体网格使用了密度能量误差函数来优化节点的位置。三维模型运用主成分分析后，提高了初始四面体单元的质量；改进切点向模型特征点移动保证了四面体化后模型的局部特征；密度能量误差函数优化了最终四面体网格的质量。

Description

三维网格模型四面体化方法

技术领域

本发明涉及一种三维网格模型四面体化方法，特别是一种适用于有限元分析、模型逼近、空间网格剖分等应用的三维网格模型四面体化方法。 

背景技术

有限元分析方法是一种能够有效地求解各种工程计算问题的通用数值分析方法。其中网格生成是有限元方法的关键步骤，应用领域十分广泛，因此有限元网格生成算法成为众多国内外学者研究的重点。其中，四面体网格单元简单灵活，对复杂边界有较强的适应能力，是最常用的空间单元表示结构之一，可用于有限元分析、模型逼近、空间网格剖分等。因此，三维模型四面体化已成为一种重要的有限元网格生成技术。 

20世纪80年代中后期，国内外学者将二维平面网格生成算法的研究转向三维空间，对四面体网格生成算法进行了广泛的研究，相关技术可分为网格模板法、拓扑分解法、映射法和几何分解法等。网格模板法适用于自适应网格生成，完全自动且速度快，但对边界单元的质量无法保证，对模型位置方向也特别敏感。拓扑分解法过于依赖于几何体的拓扑结构，使生成的四面体网格质量不理想，甚至很差。映射法生成四面体网格单元很方便，但对模型形状要求过高，过于复杂就很难处理。几何分解法所生成的网格单元形状和分布均较好但是自动化程度低，不利于复杂模型的网格生成。 

目前，出现了以下几种四面体化方式：(1)基于四叉树法的思想，提出八叉树法，用于实现三维空间的网格划分，但是得到的仅是初始四面体，生成的四面体有待优化。(2)基于体心立方(body-centered cubic,BCC)结构，并结合自适应八叉树的思想，首先得到初始的四面体，再对边界四面体重四面体化得到最终的四面体逼近结果。该方法操作简单，易于实现，且具有较好的四面体化效果，但对模型位置方向较为敏感，且四面体单元质量还可进一步优化。(3)把拓扑分解法延伸到三维空间的网格划分，该方法虽然理论比较简单，实现起来也比较容易，但是在几何因素方面考虑欠缺，生成的网格形状不理想。(4)将几何分解法应用到三维空间的网格生成，首先找出一个分割面，然后用该面将待划分区域用递归的方法将其分割成两个子体，循环操作，直至所有的子体都变成四面体，但是该方法操作繁琐。 

发明内容

本发明公开了一种三维模型四面体化方法，提高了初始四面体单元的质量，保证了四面体化后模型的局部特征，优化了最终四面体网格的质量。 

为解决上述技术问题，本发明提供一种三维模型四面体化方法，包括以下步骤： 

(1)、模型初始位置预处理 

根据三维网格模型数据中面的法向，运用表面法向主成分分析法(NPCA)提取三维网格模型的主要成分，可有利于提高四面体单元的质量。 

(2)、构建体心立方结构 

根据三维网格模型表面曲率，自适应构建基于欧式距离变换的细分八叉树结构的体心立方，得到初始四面体。 

(3)、计算节点符号及两端节点符号相反的四面体边的切点 

计算体心立方的节点符号及切点，其中，节点符号计算方法为：根据角度权伪法矢量判断节点在三维网格模型内部或外部，定义位于三维网格模型内部的节点的符号为正，位于三维网格模型外部的节点的符号为负，位于三维网格模型表面的节点的符号为零；或者，定义位于三维网格模型外部的节点的符号为正，位于三维网格模型内部的节点的符号为负，位于三维网格模型表面的节点的符号为零；切点计算方法为：模型中节点符号不同的边所在的三角面组成一个三角面集合，分别计算每条节点符号不同的边与这些三角面所在的平面的交点，切点为在三角面内的交点。 

(4)、切点移动 

分别计算切点和两端节点之间的距离与边长的比值，比值小于阈值的切点进行移动，切点移动方向为三维网格模型表面的特征点，且该特征点对应的虚拟原点与切点之间的距离最小，其中，切点移动方向的确定方法为：计算切点所在三维网格三角面的特征点，使用拉普拉斯算法计算各特征点对应的拉普拉斯坐标虚拟原点，计算切点和各虚拟原点之间的距离，移动方向为与切点距离最小的虚拟原点对应的特征点，所述拉普拉斯算法具体为：三维网格上某一点的拉普拉斯坐标定义为该点指向其相邻点的所有向量之和，拉普拉斯坐标具有平移不变性，使用网格上相邻顶点的线性组合来表示顶点的网格坐标，描述网格的细节特性和局部特征，具体的拉普拉斯坐标公式表示如下： 

${\mathrm{\σ}}_i=\underset{j\∈N\left(i\right)}{\mathrm{\Σ}}{w}_{\mathrm{ij}}(v_i-v_j)$

其中N(i)＝{j|{i,j}∈E}是与顶点v_i的相邻的顶点的集合，w_ij表示顶点v_i和v_j之间的权重，满足等式

(5)、基于网格质量的密度能量误差函数对四面体的边界重新四面体化，构成逼近模型的四面体网格；其中，基于网格质量的密度能量误差函数如下： 

$E_{\mathrm{ODT}}=\frac{1}{4}\underset{i}{\mathrm{\Σ}}{x_i}^2\left({\mathrm{\Σ}}_{T_j\∈{\mathrm{\Ω}}_i}{q}_j\right|T_j\left|\right)-{\∫}_M{x}^2\mathrm{dx}$

其中|T_j|是T_j的面积，q_j是网格质量，通过求解上式的极小值可得优化后节点的位置为： 

${x_i}^{*}=\frac{1}{{\mathrm{\Σ}}_{T_j\∈{\mathrm{\Ω}}_i}{q}_{j\left|T_j\right|}}\underset{T_j\∈{\mathrm{\Ω}}_i}{\mathrm{\Σ}}{q}_j\left|T_j\right|c_j$

其中x_i ^*是x_i的优化后的位置，c_j是T_j的外接圆圆心； 

在四面体网格中，网格单元的质量函数可设为： 

$q=\sqrt{24}\frac{r}{l_{\max }}$

其中r为四面体的内接球半径，l_max是四面体的最长边。 

本发明公开的三维模型四面体化方法，基于体心立方的网格四面体化算法，增加了三维模型主成分分析的预处理，改进了模型边界切点的移动方式转向模型特征点移动，并对最终的四面体网格使用了密度能量误差函数来优化节点的位置。三维模型运用主成分分析后，提高了初始四面体单元的质量；改进切点向模型特征点移动保证了四面体化后模型的局部特征；密度能量误差函数优化了最终四面体网格的质量。 

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。 

图1为二维图形的初始位置。 

图2为图1经NPCA主成分分析法预处理后的示意图。 

图3为图2的初始四面体示意图。 

图4为图3中初始四面体的节点和切点示意图。 

图5为图4中需要移动的切点的示意图。 

图6为图5中需要移动的切点移动之后的示意图。 

图7为图6边界重四面体化后的示意图。 

图8为兔子模型的初始位置。 

图9为图8经基于顶点位置主成分分析法(CPCA)后的示意图。 

图10为图8经表面法向主成分分析法(NPCA)后的效果图。 

图11为点到三角面的最近点在三角面中间的示意图。 

图12为点到三角面的最近点在三角面边上的示意图。 

图13为点到三角面的最近点在三角面顶点的示意图。 

图14为切点与两端节点距离的示意图。 

图15为点在拉普拉斯坐标下的虚拟原点。 

图16为切点移动方式的示意图。 

图17为一个网格节点星形结构的示意图。 

图18为节点优化参数示意图。 

图19为兔子模型的初始模型。 

图20为图19经NPCA主成分分析法预处理后的示意图。 

图21为图20的初始四面体示意图。 

图22为图21中初始四面体的节点和切点示意图。 

图23为图22经切点移动后的示意图。 

图24为图23边界重四面体化后的剖面示意图。 

图25为兔子模型基于现有四面体化方法的内部剖面效果图。 

图26为兔子模型基于本发明四面体化方法的内部剖面效果图。 

图27为犰狳模型基于现有四面体化方法的内部剖面效果图。 

图28为犰狳模型基于另一现有四面体化方法的内部剖面效果图。 

图29为犰狳模型基于本发明四面体化方法的内部剖面效果图。 

图30为龙模型基于现有四面体化方法的内部剖面效果图。 

图31为龙模型基于另一现有四面体化方法的内部剖面效果图。 

图32为龙模型基于本发明四面体化方法的内部剖面效果图。 

图33为恐龙模型基于现有四面体化方法的内部剖面效果图。 

图34为恐龙模型基于另一现有四面体化方法的内部剖面效果图。 

图35为恐龙模型基于又一现有四面体化方法的内部剖面效果图。 

图36为恐龙模型基于本发明四面体化方法的内部剖面效果图。 

图37为牙齿模型基于现有四面体化方法的内部剖面效果图。 

图38为牙齿模型基于另一现有四面体化方法的内部剖面效果图。 

图39为牙齿模型基于又一现有四面体化方法的内部剖面效果图。 

图40为牙齿模型基于本发明四面体化方法的内部剖面效果图。 

图41为牛模型基于现有四面体化方法的外部剖面效果图。 

图42为牛模型基于本发明四面体化方法的外部剖面效果图。 

图43为人体模型基于本发明四面体化方法的内四面体化效果图。 

图44为人体模型基于本发明四面体化方法的外部四面体化效果图。 

图45为兔子模型基于现有方法二的网格二面角度分布图。 

图46为兔子模型基于本发明方法的网格二面角度分布图。 

图47为犰狳模型基于现有方法一的四面体网格二面角度分布图。 

图48为犰狳模型基于现有方法二的四面体网格二面角度分布图。 

图49为犰狳模型基于本发明方法的四面体网格二面角度分布图。 

图50为龙模型基于现有方法一的四面体网格二面角度分布图。 

图51为龙模型基于现有方法二的四面体网格二面角度分布图。 

图52为龙模型基于本发明方法的四面体网格二面角度分布图。 

图53为恐龙模型基于现有方法一的四面体网格二面角度分布图。 

图54为恐龙模型基于现有方法二的四面体网格二面角度分布图。 

图55为恐龙模型基于本发明方法的四面体网格二面角度分布图。 

图56为牙齿模型基于现有方法一的四面体网格二面角度分布图。 

图57为牙齿模型基于现有方法二的四面体网格二面角度分布图。 

图58为牙齿模型基于本发明方法的四面体网格二面角度分布图。 

具体实施方式

下面结合实施例对本发明做进一步的详细说明，以下实施例是对本发明的解释而本发明并不局限于以下实施例。 

本发明公开的三维网格模型四面体化方法，是基于体心立方(body-centeredcubic,BCC)结构来构建的。体心立方结构，是材料科学技术组织结构中的一种立方晶系结构，这种结构具有很多优良的特性，在实际的工程中得到了广泛的 应用。 

实施例一 

如图1-7所示，图1为初始模型，采用本发明公开的三维网格四面体化方法进行处理，具体算法步骤如下： 

(1)、模型初始位置预处理 

三维模型主成分分析(principal component analysis,PCA)是模型处理中常用的一种方法，因为它能够有效地提取出模型的主要成分和结构，一定程度上去除噪音和干扰，因此成为了对3D模型施加其它后续操作的基本变换。通常，在对模型施加PCA操作时，模型的3D自由度需要满足以下几个基本要求：平移不变性、旋转不变性、缩放不变性。在进行PCA分析时，平移和缩放不变性比较容易实现，因而各种方法的差别主要就集中在旋转不变性的实现上。 

由参考文献《Efficient3D shape matching and retrieval using a concreteradialized spherical projection representation》可知，三维模型主成分分析可分为基于顶点位置主成分分析法(CPCA)和表面法向主成分分析法(NPCA)两种。如图8所示，图8为模型原始位置，图9为CPCA方法的分析结果，图10为NPCA方法的分析结果，可以看出NPCA方法的旋转不变性较好。 

另外，由于生成的四面体的二面角度和模型表面法向有关，因此本发明根据三维模型的数据中面的法向，运用NPCA主成分分析法，对模型进行旋转变换，能使初始四面体化的二面角度得到优化。 

因此，在本步骤中采用NPCA主成分分析法对图1所示的初始位置模型进行预处理，预处理后得到图2所示的结果，经本步骤后有利于提高步骤(2)中四面体单元的质量。 

(2)、构建体心立方结构 

由参考文献《Feature-sensitive tetrahedral mesh generation with guaranteedquality》可知，根据三维网格模型表面曲率的不同，可自适应构建基于欧式距离变换的细分八叉树结构的体心立方结构。如图3所示，经本步骤后得到了初始四面体。 

(3)、计算节点符号及两端节点符号相反的四面体边的切点 

构建体心立方结构得到初始四面体后，需要计算体心立方中节点的符号，若节点位于模型的内部则符号为正，若位于模型外部则符号为负，若位于模型表面则符号为零；或者若节点位于模型的内部则符号为负，若位于模型外部则符号为正，若位于模型表面则符号为零。若体心立方中某条边两端节点的符号 不同，则该边一定和模型的表面相交，这些交点称为切点。切点的具体计算方法为：模型中节点符号不同的边所在的三角面组成一个三角面集合，分别计算每条节点符号不同的边与这些三角面所在的平面的交点，若交点在三角面内，该交点即为切点。 

由前述可知，计算节点符号的问题可归结于判断空间中给定点是否在三维模型内部，而给定点是否在三维模型内部，与给定点到三维模型表面的距离有关。在求给定点P到三角面P₀P₁P₂的距离时，关键是要找到三角面P₀P₁P₂中距点P最近的点，由参考文献《Algorithm for fast calculating the nearest distancebetween space point and arbitrary polyhedron》公开的动态球搜索技术可以快速计算获得点到三角网格的最近点。如图11示，点P到三角面P₀P₁P₂的最近点P′在三角面中间，点P到三角面P₀P₁P₂的距离为PP′，三角面P₀P₁P₂的法向矢量，可以用来判断点P在模型内部还是外部。如图12所示，点P到三角面P₀P₁P₂的最近点P"在三角面的边上，点P到三角面P₀P₁P₂的距离为PP"，公共边为P₁P₂的两个三角面法向矢量的加权，可以用来判断点P在模型内部还是外部。如图13所示，点P到三角面P₀P₁P₂的最近点为三角面的顶点P1，点P到三角面P₀P₁P₂的距离为PP₁，参考文献《Computing vertex normals from polygonal facets》公开的角度权伪法矢量，可以用来判断点P在模型内部还是外部。 

经本步骤后计算得到切点，如图4所示，“+”表示计算得到的内部节点，“-”表示计算得到的外部节点，“。”表示计算得到的切点。 

(4)、切点移动 

现有技术中的移动方向为向距离切点较近的节点移动，但是这样容易破坏模型的特征结构。因此，本发明考虑运用拉普拉斯坐标对部分距离节点相对较近的切点进行移动，并将切点移动方向改为模型表面的特征点，从而保持四面体网格模型的局部特征，避免四面体二面角度过小或过大，在优化四面体二面角度的同时保持了模型的特征。 

由参考文献《Feature-sensitive tetrahedral mesh generation with guaranteedquality》可知需要移动的切点的确定方法为：分别计算切点和两端节点之间的距离与边长的比值，比值小于阈值的切点进行移动。如图14所示，若p为切点，v₁为外部节点，v₂为内部节点，分别计算切点p和两端节点v₁、v₂之间的距离pv₁、pv₂与边长v₁v₂的比值λ₁、λ₂，若比值λ₁小于阈值λ₀或者比值λ₂小于阈值λ₀(阈值λ₀可定为0.2)，则切点p进行移动，图5中“▲”表示需要移动的切点。 

切点移动方向的确定方法为：计算切点所在三维网格三角面的特征点，切点所在三维网格三角面的三个顶点就是该切点的特征点，使用拉普拉斯算法计算各特征点对应的拉普拉斯坐标虚拟原点，计算切点和各虚拟原点之间的距离，移动方向为与切点距离最小的虚拟原点对应的特征点。所述拉普拉斯算法具体为：三维网格上某一点的拉普拉斯坐标定义为该点指向其相邻点的所有向量之和，拉普拉斯坐标具有平移不变性，使用网格上相邻顶点的线性组合来表示顶点的网格坐标，描述网格的细节特性和局部特征，具体的拉普拉斯坐标公式表示如下： 

${\mathrm{\σ}}_i=\underset{j\∈N\left(i\right)}{\mathrm{\Σ}}{w}_{\mathrm{ij}}(v_i-v_j)$

其中N(i)＝{j|{i,j}∈E}是与顶点v_i的相邻的顶点的集合，w_ij表示顶点v_i和v_j之间的权重，满足等式

如图15所示，vo_i为点v_i的拉普拉斯坐标下的虚拟原点。切点移动方式如图16所示，假设v_c为模型表面的一个切点，该切点v_c所在三维网格三角面的特征点为v_i1、v_i2和v_i3，对应的拉普拉斯坐标虚拟原点为vo_i1、vo_i2和vo_i3，计算切点v_c和各虚拟原点之间的距离，切点v_c与特征点v_i3对应的拉普拉斯坐标虚拟原点vo_i3的距离最近，则切点v_c的移动方向为特征点v_i3，即切点v_c的移动方向为v_c→v_i3。 

经本步骤后，图5中“▲”表示的切点移动到了图6所示的位置，完成切点移动。 

(5)、边界重四面体化 

切点移动后，再对四面体的边界重新四面体化得到最终的四面体集合，构成逼近模型的四面体网格。现有技术中，有基于Delaunay三角形网格的ODT光顺算法对每个星形结构的核心节点位置进行优化，网格中节点的星形结构如图17所示，星形结构的能量误差函数为： 

$E_{\mathrm{ODT}}=\frac{1}{4}\underset{i}{\mathrm{\Σ}}{x_i}^2\left|{\mathrm{\Ω}}_i\right|-{\∫}_M{x}^2\mathrm{dx}$

其中Ω_i是x_i所对应的星形结构，|Ω_i|是Ω_i的面积； 

考虑到优化的节点位置和网格单元的角度有关，而单元的角度又是衡量单元质量的一个重要的因素，可将上面的能量误差函数改为基于网格质量的密度能量误差函数，因此改进后的密度能量误差函数如下： 

$E_{\mathrm{ODT}}=\frac{1}{4}\underset{i}{\mathrm{\Σ}}{x_i}^2\left({\mathrm{\Σ}}_{T_j\∈{\mathrm{\Ω}}_i}{q}_j\right|T_j\left|\right)-{\∫}_M{x}^2\mathrm{dx}$

其中|T_j|是T_j的面积，q_j是网格质量，通过求解上式的极小值可得优化后节点的位置为： 

${x_i}^{*}=\frac{1}{{\mathrm{\Σ}}_{T_j\∈{\mathrm{\Ω}}_i}{q}_{j\left|T_j\right|}}\underset{T_j\∈{\mathrm{\Ω}}_i}{\mathrm{\Σ}}{q}_j\left|T_j\right|c_j$

其中x_i ^*是x_i的优化后的位置，c_j是T_j的外接圆圆心。 

在四面体网格中，网格单元的质量函数可设为： 

$q=\sqrt{24}\frac{r}{l_{\max }}$

其中r为四面体的内接球半径，l_max是四面体的最长边。 

类似于二维情况，节点优化相关参数如图17所示，C_i为需优化的节点，C_j1为外接圆圆心，Cg_j1为内接圆圆心，r_j1为内接圆半径。 

经本步骤后，得到如图7所示的节点优化后的最终四面体化结果。 

实施例二 

如图19-24所示，图19为兔子初始模型，采用本发明公开的三维网格四面体化方法进行处理： 

(1)、模型初始位置预处理 

三维模型主成分分析(principal component analysis,PCA)是模型处理中常用的一种方法，因为它能够有效地提取出模型的主要成分和结构，一定程度上去除噪音和干扰，因此成为了对3D模型施加其它后续操作的基本变换。通常，在对模型施加PCA操作时，模型的3D自由度需要满足以下几个基本要求：平移不变性、旋转不变性、缩放不变性。在进行PCA分析时，平移和缩放不变性比较容易实现，因而各种方法的差别主要就集中在旋转不变性的实现上。 

由参考文献《Efficient3D shape matching and retrieval using a concreteradialized spherical projection representation》可知，三维模型主成分分析可分为基于顶点位置主成分分析法(CPCA)和表面法向主成分分析法(NPCA)两种。如图8至10所示，图8为模型原始位置，图9为CPCA方法的分析结果，图10为NPCA方法的分析结果，可以看出NPCA方法的旋转不变性较好。 

另外，由于生成的四面体的二面角度和模型表面法向有关，因此本发明根据三维模型的数据中面的法向，运用NPCA主成分分析法，对模型进行旋转变换，能使初始四面体化的二面角度得到优化。 

因此，在本步骤中采用NPCA主成分分析法对图19所示的初始位置模型进行预处理，预处理后得到图20所示的结果，经本步骤后有利于提高步骤(2)中四面体单元的质量。 

(2)、构建体心立方结构 

由参考文献《Feature-sensitive tetrahedral mesh generation with guaranteedquality》可知，根据三维网格模型表面曲率的不同，可自适应构建基于欧式距离变换的细分八叉树结构的体心立方结构。如图21所示，经本步骤后得到了初始四面体。 

(3)、计算节点符号及两端节点符号相反的四面体边的切点 

构建体心立方结构得到初始四面体后，需要计算体心立方中节点的符号，若节点位于模型的内部则符号为正，若位于模型外部则符号为负，若位于模型表面则符号为零；或者若节点位于模型的内部则符号为负，若位于模型外部则符号为正，若位于模型表面则符号为零。若体心立方中某条边两端节点的符号不同，则该边一定和模型的表面相交，这些交点称为切点。切点的具体计算方法为：模型中节点符号不同的边所在的三角面组成一个三角面集合，分别计算每条节点符号不同的边与这些三角面所在的平面的交点，若交点在三角面内，该交点即为切点。 

由前述可知，计算节点符号的问题可归结于判断空间中给定点是否在三维模型内部，而给定点是否在三维模型内部，与给定点到三维模型表面的距离有关。在求给定点P到三角面P₀P₁P₂的距离时，关键是要找到三角面P₀P₁P₂中距点P最近的点，由参考文献《Algorithm for fast calculating the nearest distancebetween space point and arbitrary polyhedron》公开的动态球搜索技术可以快速计算获得点到三角网格的最近点。如图11所示，点P到三角面P₀P₁P₂的最近点P′在三角面中间，点P到三角面P₀P₁P₂的距离为PP′，三角面P₀P₁P₂的法向矢量，可以用来判断点P在模型内部还是外部。如图12所示，点P到三角面P₀P₁P₂的最近点P"在三角面的边上，点P到三角面P₀P₁P₂的距离为PP"，公共边为P₁P₂的两个三角面法向矢量的加权，可以用来判断点P在模型内部还是外部。如图13所示，点P到三角面P₀P₁P₂的最近点为三角面的顶点P₁，点P到三角面P₀P₁P₂的距离为PP₁，参考文献《Computing vertex normals from polygonal facets》公开的角度权伪法矢量，可以用来判断点P在模型内部还是外部。 

如图22所示，经本步骤后计算得到兔子模型的切点。 

(3)、切点移动 

现有技术中的移动方向为向距离切点较近的节点移动，但是这样容易破坏模型的特征结构。因此，本发明考虑运用拉普拉斯坐标对部分距离节点相对较近的切点进行移动，并将切点移动方向改为模型表面的特征点，从而保持四面体网格模型的局部特征，避免四面体二面角度过小或过大，在优化四面体二面角度的同时保持了模型的特征。 

由参考文献《Feature-sensitive tetrahedral mesh generation with guaranteedquality》可知需要移动的切点的确定方法为：分别计算切点和两端节点之间的距离与边长的比值，比值小于阈值的切点进行移动。如图14所示，若p为切点，v₁为外部节点，v₂为内部节点，分别计算切点p和两端节点v₁、v₂之间的距离pv₁、pv₂与边长v₁v₂的比值λ₁、λ₂，若比值λ₁小于阈值λ₀或者比值λ₂小于阈值λ₀(阈值λ₀可定为0.2)，则切点p进行移动。 

切点移动方向的确定方法为：计算切点所在三维网格三角面的特征点，切点所在三维网格三角面的三个顶点就是该切点的特征点，使用拉普拉斯算法计算各特征点对应的拉普拉斯坐标虚拟原点，计算切点和各虚拟原点之间的距离，移动方向为与切点距离最小的虚拟原点对应的特征点。所述拉普拉斯算法具体为：三维网格上某一点的拉普拉斯坐标定义为该点指向其相邻点的所有向量之和，拉普拉斯坐标具有平移不变性，使用网格上相邻顶点的线性组合来表示顶点的网格坐标，描述网格的细节特性和局部特征，具体的拉普拉斯坐标公式表示如下： 

${\mathrm{\σ}}_i=\underset{j\∈N\left(i\right)}{\mathrm{\Σ}}{w}_{\mathrm{ij}}(v_i-v_j)$

其中N(i)＝{j|{i,j}∈E}是与顶点v_i的相邻的顶点的集合，w_ij表示顶点v_i和v_j之间的权重，满足等式

如图15所示，vo_i为v_i的拉普拉斯坐标下的虚拟原点。切点移动方式如图16所示，v_c为模型表面的一个切点，该切点v_c所在三维网格三角面的特征点为v_i1、v_i2和v_i3，对应的拉普拉斯坐标虚拟原点为vo_i1、vo_i2和vo_i3，计算切点v_c和各虚拟原点之间的距离，切点v_c与特征点v_i3对应的拉普拉斯坐标虚拟原点vo_i3的距离最近，则切点v_c的移动方向为特征点v_i3，即切点v_c的移动方向为v_c→v_i3。 

经本步骤后，完成切点移动后的兔子模型如图23所示。 

(5)、边界重四面体化 

切点移动后，再对四面体的边界重新四面体化得到最终的四面体集合，构 成逼近模型的四面体网格。现有技术中，有基于Delaunay三角形网格的ODT光顺算法对每个星形结构的核心节点位置进行优化，网格中节点的星形结构如图17所示，星形结构的能量误差函数为： 

$E_{\mathrm{ODT}}=\frac{1}{4}\underset{i}{\mathrm{\Σ}}{x_i}^2\left|{\mathrm{\Ω}}_i\right|-{\∫}_M{x}^2\mathrm{dx}$

其中Ω_i是x_i所对应的星形结构，|Ω_i|是Ω_i的面积； 

考虑到优化的节点位置和网格单元的角度有关，而单元的角度又是衡量单元质量的一个重要的因素，可将上面的能量误差函数改为基于网格质量的密度能量误差函数，因此改进后的密度能量误差函数如下： 

$E_{\mathrm{ODT}}=\frac{1}{4}\underset{i}{\mathrm{\Σ}}{x_i}^2\left({\mathrm{\Σ}}_{T_j\∈{\mathrm{\Ω}}_i}{q}_j\right|T_j\left|\right)-{\∫}_M{x}^2\mathrm{dx}$

其中|T_j|是T_j的面积，q_j是网格质量，通过求解上式的极小值可得优化后节点的位置为： 

${x_i}^{*}=\frac{1}{{\mathrm{\Σ}}_{T_j\∈{\mathrm{\Ω}}_i}{q}_{j\left|T_j\right|}}\underset{T_j\∈{\mathrm{\Ω}}_i}{\mathrm{\Σ}}{q}_j\left|T_j\right|c_j$

其中x_i ^*是x_i的优化后的位置，c_j是T_j的外接圆圆心。 

在四面体网格中，网格单元的质量函数可设为： 

$q=\sqrt{24}\frac{r}{l_{\max }}$

其中r为四面体的内接球半径，l_max是四面体的最长边。 

类似于二维情况，节点优化相关参数如图18所示，C_i为需优化的节点，C_j1为外接圆圆心，Cg_j1为内接圆圆心，r_j1为内接圆半径。 

经本步骤后，得到如图24所示的节点优化后的最终四面体化结果。 

本发明提出的三维网格模型四面体化方法，首先在四面体化前对三维模型进行基于NPCA的主成分分析，调整模型的初始位置；其次使用基于体心立方的四面体化方法，在切点移动时使用拉普拉斯坐标，并将切点移动方向改为模型表面的特征点，在优化四面体二面角度的同时保持了模型的特征；最后使用改进的密度能量误差函数优化四面体网格单元。本发明方法在一定精度下减少了四面体的个数，并能较好的保持模型局部特征，进一步提高了四面体网格单元质量。 

实验结果： 

1、使用VS2008和OpenGL在2.40GHz Intel Core2CPU，2GB内存的PC 机上实现了本发明公开的方法，并对不同模型的四面体化结果进行验证，实验结果如图25-44所示。图25-26为兔子模型基于不同方法的四面体化内部剖面效果图，图27-29为犰狳模型基于不同方法的四面体化内部剖面效果图，图30-32为龙模型基于不同方法的四面体化内部剖面效果图，图33-36为恐龙模型基于不同方法的四面体化内部剖面效果图，图37-40为牙齿模型基于不同方法的四面体化内部剖面效果图，其中，图26、图29、图32、图36、图40均为本发明提供的四面体化方法的效果图；图41-42为牛模型基于不同方法的外部四面体化效果图，图42为本发明提供的四面体化方法的效果图；图43-44为人体模型基于本发明方法的内、外部四面体化效果图。从实验图可以发现本发明公开的三维网格模型四面体化方法效果较好。 

2、四面体二面角度的大小可用来衡量四面体单元质量的好坏，角度过小或过大，则认为单元质量较差。兔子、犰狳、龙、恐龙、牙齿模型在不同方法四面体化后四面体二面角度分布图如图45-58所示，其中图46、图49、图52、图55、图58为本发明方法的分布直方图；犰狳、龙、恐龙、牙齿模型在不同方法四面体化后四面体网格的最大二面角和最小二面角值对比如表1所示，表一中，现有方法一是指文献《A quality tetrahedral mesh generator and three-dimensionaldelaunay triangulator》中公开的方法，现有方法二是指文献《Feature-sensitivetetrahedral mesh generation with guaranteed quality》中公开的方法。由表一可知，本发明公开的方法进一步提高了四面体网格的最小二面角度，有效地优化了网格的质量。 

表1本发明方法和其他现有方法四面体二面角度比较 

本发明的参考文献如下： 

[1]Zhang Su,Shi Fazhong.Implementation of Finite Element Method MeshGeneration from Multiple Trimmed Free Surfaces[J].Journal of Software,2005,16(11):2008-2013(张苏,施法中.多裁剪自由曲面生成有限元网格的实现[J].软件学报,2005,16(11):2008-2013) 

[2]Chen Xin,Xiong Yueshan.Tetrahedral Mesh Generation from MedicalVolume Data[J].Journal of Software,2008,19(z1):78-86(陈欣,熊岳山.基于医学体数据生成四面体网格的方法[J].软件学报,2008,19(z1):78-86) 

[3]Wang J,Yu Z.Feature-sensitive tetrahedral mesh generation withguaranteed quality[J].Computer-Aided Design,2012,44(5):400-412 

[4]Woo T C,Thomasma T.An algorithm for generating solid elements inobjects with holes[J].Computer&Structure,1984,18(2):333-342 

[5]Shephard M S,Georges M K.Automatic three-dimensional meshgeneration by the finite octree technique[J].International Journal for NumericalMethods in Engineering,1991,32(4):709-749 

[6]Yerry M A,Shephard M S.Automatic three-dimensional mesh generationby the modified-octree technique[J].International Journal for Numerical Methodsin Engineering,1984,20(11):1965-1990 

[7]Amresh J.Modern methods for automatic FE mesh generation[C]//ModernMethods for Automating Finite Element Mesh Generation.ASCE,1986:19-28 

[8]Papadakis P,Pratikakis I,Perantonis S,et al.Efficient3D shape matchingand retrieval using a concrete radialized spherical projection representation[J].Pattern Recognition,2007,40(9):2437-2452 

[9]Thürmer G,Wüthrich C.Computing vertex normals from polygonalfacets[J].Journal of Graphics Tools,1998,3(1):43-46 

[10]Fang Xiang,Bao Hujun,Wang Pingan,et al.Algorithm for fastcalculating the nearest distance between space point and arbitrary polyhedron[J].Journal of Computer-Aided Design&Computer Graphics,2001,13(9):788-792(inChinese)(方向,鲍虎军,王平安,等.点到任意多面体距离的快速计算方法[J].计算机辅助设计与图形学学报,2001,13(9):788-792) 

[11]Sorkine O,Cohen-Or D,Lipman Y,et al.Laplacian surfaceediting[C]//Proceedings of the2004Eurographics/ACM SIGGRAPH symposiumon Geometry processing.ACM,2004:175-184 

[12]Alliez P,Cohen-Steiner D,Yvinec M,et al.Variational tetrahedralmeshing[J].ACM Transactions on Graphics(TOG),2005,24(3):617-625 

[13]Liu Yan,Chang Jihai,Guan Zhenqun.Enhanced3D optimal delaunaytriangulation optimization method for tetrahedral mesh quality[J].Journal ofComputer-Aided Design&Computer Graphics,2012,24(7):949-953(in Chinese) 

(刘岩,昌继海,关振群.改进的三维ODT四面体网格质量优化算法[J].计算机辅助设计与图形学学报,2012,24(7):949-953) 

[14]Si H,TetGen A.A quality tetrahedral mesh generator andthree-dimensional delaunay triangulator[J].Weierstrass Institute for AppliedAnalysis and Stochastic,Berlin,Germany,2006 

[15]Joachim S,Hannes G,Robert G.NETGEN-automatic mesh generator[J].Austria:Johannes Kepler University Linz(2004-04-15)E2009-11-161.http://www.hpfem.jku.at/netgen,2012。 

Claims (1)

1.一种三维网格模型四面体化方法，包括以下步骤：

(1)、模型初始位置预处理

根据三维网格模型数据中面的法向，运用表面法向主成分分析法(NPCA)提取三维网格模型的主要成分；

(2)、构建体心立方结构

根据三维网格模型表面曲率，自适应构建基于欧式距离变换的细分八叉树结构的体心立方，得到初始四面体；

(3)、计算节点符号及两端节点符号相反的四面体边的切点

计算体心立方的节点符号及切点，其中，节点符号计算方法为：根据角度权伪法矢量判断节点在三维网格模型内部或外部，定义位于三维网格模型内部的节点的符号为正，位于三维网格模型外部的节点的符号为负，位于三维网格模型表面的节点的符号为零；或者，定义位于三维网格模型外部的节点的符号为正，位于三维网格模型内部的节点的符号为负，位于三维网格模型表面的节点的符号为零；切点计算方法为：模型中节点符号不同的边所在的三角面组成一个三角面集合，分别计算每条节点符号不同的边与这些三角面所在的平面的交点，切点为在三角面内的交点；

(4)、切点移动

分别计算切点和两端节点之间的距离与边长的比值，比值小于阈值的切点进行移动，切点移动方向为三维网格模型表面的特征点，且该特征点对应的虚拟原点与切点之间的距离最小；其中，切点移动方向的确定方法为：计算切点所在三维网格三角面的特征点，使用拉普拉斯算法计算各特征点对应的拉普拉斯坐标虚拟原点，计算切点和各虚拟原点之间的距离，移动方向为与切点距离最小的虚拟原点对应的特征点，所述拉普拉斯算法具体为：三维网格上某一点的拉普拉斯坐标定义为该点指向其相邻点的所有向量之和，拉普拉斯坐标具有平移不变性，使用网格上相邻顶点的线性组合来表示顶点的网格坐标，描述网格的细节特性和局部特征，具体的拉普拉斯坐标公式表示如下：

${\mathrm{\σ}}_i=\underset{j\∈N\left(i\right)}{\mathrm{\Σ}}{w}_{\mathrm{ij}}(v_i-v_j)$

其中N(i)＝{j|{i,j}∈E}是与顶点vi的相邻的顶点的集合，w_ij表示顶点v_i和v_j之间的权重，满足等式

(5)、边界重四面体化

基于网格质量的密度能量误差函数对四面体的边界重新四面体化，构成逼近模型的四面体网格；其中，基于网格质量的密度能量误差函数如下：

$E_{\mathrm{ODT}}=\frac{1}{4}\underset{i}{\mathrm{\Σ}}{x_i}^2\left({\mathrm{\Σ}}_{T_j\∈{\mathrm{\Ω}}_i}{q}_j\right|T_j\left|\right)-{\∫}_M{x}^2\mathrm{dx}$

其中，|T_j|是T_j的面积，q_j是网格质量，通过求解上式的极小值可得优化后节点的位置为：

${x_i}^{*}=\frac{1}{{\mathrm{\Σ}}_{T_j\∈{\mathrm{\Ω}}_i}{q}_{j\left|T_j\right|}}\underset{T_j\∈{\mathrm{\Ω}}_i}{\mathrm{\Σ}}{q}_j\left|T_j\right|c_j$

其中x_i ^*是x_i的优化后的位置，c_j是T_j的外接圆圆心；

在四面体网格中，网格单元的质量函数可设为：

$q=\sqrt{24}\frac{r}{l_{\max }}$

其中r为四面体的内接球半径，l_max是四面体的最长边。