CN105243688A - 基于重心Voronoi图的非钝角的重新网格化方法 - Google Patents

Abstract

本发明实施例公开了一种基于重心Voronoi图的非钝角的重新网格化方法，该方法包括以下步骤：步骤S1、获取待处理网格的采样点集；步骤S2、对采样点集进行CVT优化；步骤S3、对CVT优化后得到的网格的度数进行优化；步骤S4、基于经过度数优化后得到的网格，进行特征敏感的CVT优化，并引入对偶惩罚项，惩罚较短的Voronoi边，以移除经过度数优化后得到的网格中的钝角三角形。本发明实施例通过避免坏三角形(＜30°或者＞90°)的产生，得到了高品质的重新网格化结果，非常适用于物理模拟以及几何建模等应用领域。

Description

基于重心Voronoi图的非钝角的重新网格化方法

技术领域

本发明实施例涉及计算机图形处理技术领域，尤其是涉及一种基于重心Voronoi图的非钝角的重新网格化方法。

背景技术

多边形几何处理提供了很多有用的技术，例如网格平滑、网格变形、参数化、表面近似和网格分割等，可用来产生、编辑和制造三维模型。

网格是三维几何表示的基本形式，网格生成技术是连接计算机图形学和工业界的重要桥梁。借助激光扫描仪得到网格数据是常用的网格获取方法，但是通过该方法得到的模型具有数据量大、网格质量差等缺点，不能直接用于计算机动画、有限元模拟等高级应用。此外，曲面网格生成方法则研究如何从一个已有曲面出发，通过计算得到高质量的输出网格来近似表示输入曲面。因其得到的网格质量好而被广泛使用和研究。如果输入曲面用三角形网格表示，则这一过程称之为曲面的重新网格化(remeshing)，重新网格化是获得高质量网格模型的重要手段。

很多的重新网格化方法基于求解偏微分方程(PDE)。如果网格中含有小角度(＜30°)的三角形，将会影响Laplacian的条件数，继而影响网格化的品质。因此有很多工作都在侧重于如何消除小角度的三角形。Mullen2011(P.Mullen,P.Memari,F.deGoes,andM.Desbrun,“Hot:Hodge-optimizedtriangulations,”ACMTrans.onGraphics(Proc.SIGGRAPH),vol.30,no.4,pp.103:1–103:12,July2011.)提出使用powerdiagram和regulartriangulation来优化网格质量。deGoes2014(F.deGoes,P.Memari,P.Mullen,andM.Desbrun,“Weightedtriangulationsforgeometryprocessing,”ACMTrans.onGraphics,vol.33,no.3,pp.28:1–28:13,2014.)使用带权三角化来提高网格品质。但是，这些工作并不能消除网格中的钝角三角形。在求解偏微分方程过程中，标准的做法是使用余切(cotan)Laplacian来建立Laplacian矩阵。如果含有钝角三角形，则权重可能产生负值，影响方程求解。另外，非钝角的三角网格也有很好的品质。例如，对于一些有限元方法，非钝角网格具有更多好的数值属性；非钝角性可以确保通过离散Harmonic映射得到的平面网格嵌入的有效性等等。关于各向同性的表面重新网格化，Li2006(J.LiandH.Zhang,“Nonobtuseremeshinganddecimation,”inProc.ofEurographicsSymposiumonGeometryProcessing,2006.)使用基于二次误差的网格简化算法来产生非钝角网格。但是，这种方法不能去除小角度的三角形。

有鉴于此，特提出本发明。

发明内容

本发明实施例的主要目的是提供一种基于重心Voronoi图的非钝角的重新网格化方法，以至少解决如何使得曲面的重新网格化的结果中不含有钝角三角形和小角度三角形的技术问题。

为实现上述目的，本发明实施例提供一种基于重心Voronoi图的非钝角的重新网格化方法，该方法包括以下步骤：

步骤S1、获取待处理网格的采样点集；

步骤S2、对所述采样点集进行CVT优化；

步骤S3、对所述CVT优化后得到的网格的度数进行优化；

步骤S4、基于经过度数优化后得到的网格，进行特征敏感的CVT优化，并引入对偶惩罚项，惩罚较短的Voronoi边，以移除所述经过度数优化后得到的网格中的钝角三角形。

本发明实施例采用网格几何处理的技术，通过上述技术方案所采用的优化手段使得新的网格可以很好地近似原始网格，避免了在网格化的结果中出现小角度(＜30°)或者钝角(＞90°)的三角形。在物理模拟以及几何建模等领域有很好的应用价值。

附图说明

图1是本发明基于重心Voronoi图的非钝角的重新网格化方法的流程图；

图2A为根据一示例性实施例示出的输入三角网格的示意图；

图2B为对图2A中的输入三角网格进行CVT优化后的网格的示意图；

图3A为根据一示例性实施例示出的对输入网格进行CVT优化后的网格示意图；

图3B为对CVT优化后的网格中的顶点加入微小扰动后的网格示意图；

图3C是对图3B得到的网格进行CVT优化后的网格示意图；

图4为钝角三角形的最长边和短的Voronoi边的对偶关系示意图；

图5为本发明实施例的方法中特征敏感CVT优化中不同的s的重新网格化结果示意图；

图6为本发明实施例的方法中不同的λ的重新网格化结果示意图；

图7为原始的输入网格、原始的CVT方法的结果、特征敏感的CVT方法的结果和本发明实施例的方法的结果示意图；

图8是本发明实施例的方法与deGoes2014(F.deGoes,P.Memari,P.Mullen,andM.Desbrun,“Weightedtriangulationsforgeometryprocessing,”ACMTrans.onGraphics,vol.33,no.3,pp.28:1–28:13,2014.)的加权三角形的方法对于输入网格进行重新网格化的结果对比图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实例，并参照附图，对本发明进一步详细说明。

需要说明的是，在没有明确限定或不冲突的情况下，本申请的实施例及实施例中的技术特征可以相互组合。下面将参考附图并结合实施例对本发明进行进一步详细地说明。显然，所描述的实施例仅仅是本申请的一部分实施例，而不是全部的实施例。基于本申请的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都应当属于本申请的保护范围。

本发明实施例的核心思想是在传统的基于CVT(CentroidalVoronoiTessellation)的重新网格化方法中引入一个惩罚项，利用能量最小化方法来优化点集分布，继而生成高品质三角网格的方法。

图1是本发明基于重心Voronoi图的非钝角的重新网格化方法的流程图，如图1所示，该方法可以包括步骤S1至步骤S4。

步骤S1、获取待处理网格的采样点集。

在一个可选的实施例中，步骤S1具体包括：由用户指定采样点数，根据定义在网格上的密度函数进行采样。

在步骤S1中，本发明实施例的输入是一个三角网格M以及定义在M上的密度函数ρ(x)(x∈M)，首先由用户指定采样点数N(N>0)，然后根据密度函数ρ(x)随机产生采样点集。如果密度函数为一个常数，则对应均匀采样；对于非均匀采样的情况，密度函数使用局部特征大小(即localfesturesize—简称lfs)，lfs小的地方采样点比较密集，反之则稀疏，即ρ(x)＝1/lfs(x)²。

步骤S2、对所述采样点集进行CVT优化。

在步骤S2中，使用CVT优化(没有惩罚项和特征敏感约束)，对采样点集进行优化，一般需要30～50次迭代。定义在R³(三维欧式空间)上的N个不同的采样点(或者种子点)的重心Voronoi(沃罗诺伊)图是一类特殊的Voronoi图，通过最小化下面的能量函数可以得到：

$F_{CVT}\left(X\right)=\underset{i=1}{\overset{n}{\Σ}}{\∫}_{{\mathrm{\Ω}}_{i|M}}\mathrm{\ρ}\left(x\right)\left|\right|x-x_i|{|}^{2}d\mathrm{\σ};$

其中， ${\mathrm{\Ω}}_i=\{x\∈R^3|\left|\right|x-x_i\left|\right|\≤\left|\right|x-x_j\left|\right|,\∀j\≠i\}$ 是点x_i的Voronoi单元；Ω_i|M是定义在三角网格M上的限制Voronoi图；ρ(x)是定义在M上的密度函数，ρ(x)≥0。

通过最小化能量函数F_CVT(Χ)，可以使原始的采样点集与每个采样点的Voronoi单元的重心重合。

图2A为根据一示例性实施例示出的输入三角网格的示意图；图2B为对图2A中的输入三角网格进行CVT优化后的网格的示意图。从图2B中可以发现，经过CVT优化后，结果中存在钝角三角形(如矩形框所示)。但是该方法可以有效地消除小角度的三角形。

步骤S3、对所述CVT优化后得到的网格的度数进行优化。

因为CVT能量函数是非凸的，存在很多局部最小值，有时这些局部最小值会导致一些顶点的度数小于5或者大于7(如图3A所示)。将这些顶点称之为坏度数顶点。

为此，在步骤S3中，采用的处理方法是：首先对于这些顶点加入微小的扰动(如图3B所示)。在一个可选的实施例中，对坏度数顶点引入扰动的具体做法可以是：在坏度数顶点的1-ring三角形中随机采样产生新的顶点，将新的顶点映射到输入网格上，并将其插入Delaunay三角形中。同时，从Delaunay三角形中移除相应的坏度数顶点。然后对扰动后的顶点进行局部的CVT优化(如图3C所示)。这种方法可以有效地移除坏度数(<5或者>7)的顶点而不需要改变其他顶点的位置，有利于优化的全局收敛。

CVT优化可以很好地避免小角度三角形的产生。通过分析CVT的结构可以发现，钝角三角形的最长边总是对偶于一个短的Voronoi边(如图4所示)。在图4中，对偶于钝角三角形最长边的短的Voronoi边用白色表示。

步骤S4、基于经过度数优化后得到的网格，进行特征敏感的CVT优化，并引入对偶惩罚项，惩罚较短的Voronoi边，以移除经过度数优化后得到的网格中的钝角三角形。

在步骤S4中，为了在网格上计算CVT，我们需要把Voronoi图限制在网格表面上，叫做restrictedvoronoidiagram—简称RVD。它被定义为三维Voronoi图Ω＝{Ω_i}与网格表面的交集，即Ω_|M＝Ω∩M＝{Ω_i|M}。RVD的对偶是三角网格，叫做restricteddelaunaytriangulation—简称RDT。

在获得好的分布的采样点集的基础上，通过在新网格(restrictedDelaunaytriangulation—简称RDT)的CVT能量函数中引入对偶惩罚项λR(Χ)来惩罚较短的Voronoi边，并进行特征敏感的CVT优化(引入法向各向异性A_s(Ν^f)到原始的CVT能量函数中，对特征敏感的CVT能量函数进行最小化)，最后对全局能量函数进行最小化，这样做的主要目的是自动地对特征敏感的地方(如，尖锐的边缘等)进行很好的采样，使得重新网格化后的结果可以很好地保持原始网格的重要特征。其中，对应的公式为：

$\begin{array}{c}{F}_{{\mathrm{CVT}}_{fs}}\left(X\right)=\underset{i=1}{\overset{n}{\&Sigma;}}\underset{f\&Element;{\mathrm{\&Omega;}}_{i|M}}{\&Sigma;}{\&Integral;}_f\mathrm{\&rho;}\left(x\right)\left|\right|A_s\left(N^f\right){(x-x_i)}^2\left|\right|d\mathrm{\&sigma;};\\ F\left(X\right)=F_{{\mathrm{CVT}}_{fs}}\left(X\right)+\mathrm{\&lambda;}R\left(X\right).\end{array}$

其中， ${\mathrm{\&Omega;}}_i=\{x\&Element;R^3|\left|\right|x-x_i\left|\right|\&le;\left|\right|x-x_j\left|\right|,\&ForAll;j\&NotEqual;i\}$ 是点x_i的Voronoi单元，M表示三角网格，Ω_i|M是定义在M上的限制Voronoi图；f表示三角形积分区域；ρ(x)(x∈M)表示密度函数；N表示三角形的法向量；Ν^f是用于积分的三角形f的法向量；Α_s(Ν)＝(s-1)ΝΝ^t+Ι_3×3，Α_s(Ν)是一个3x3的矩阵，N^t是N的转置；参数s控制法向各向异性的权重，s≥1。优选地，s＝5；I表示3x3的单位矩阵；表示特征敏感的CVT能量函数；F(Χ)表示全局能量函数；λR(Χ)为对偶惩罚项，调节λ可以改变惩罚项的权重，λ∈[0.1,1]， $R\left(X\right)=\underset{i=1}{\overset{n}{\&Sigma;}}\left|\right|\underset{\&lsqb;x_i,x_j\&rsqb;\&Element;RDT}{\&Sigma;}{w}_{i,j}(x_i-x_j)|{|}^2,w_{i,j}=\frac{\left|\right|x_i,x_j\left|\right|}{\left|\right|dual(x_i,x_j)\left|\right|+\mathrm{\&epsiv;}},$ 其中ε是一个很小的数，以避免被除数是0。在这个公式中，钝角三角形的最长的边将会贡献更多的能量，并收缩相应的长边，使得三角形更加地均匀。dual(x_i,x_j)表示限制的Voronoi边，其可能不是一个线段，即它可能是网格表面上的相邻的限制的Voronoi单元的线段序列。

图5比较了特征敏感CVT优化中不同的s的重新网格化结果。最左边的是CVT的优化结果(s＝1,λ＝0)；其余四个从左到右依次为：s＝1,5,8,10,λ＝1。从图中可以看出，更大的s会引入更多的坏形状的三角形，但是可以更好地近似原始网格。在本发明实施例中，优选地，s＝5。

图6是本发明实施例的方法与其他基于CVT优化的重新网格化的结果对比图。如图6所示，比较了不同的λ的重新网格化结果。从左到右：λ＝0,1,10。λ＝0对应的是CVT。从图中可以看出λ＝10时会平滑掉一些细节(如矩形框所示)，而λ＝1可以消除钝角三角形，同时可以很好地近似输入三角形。可见，该方法可以移除钝角三角形同时保证顶点的好的全局分布。

图7是本发明实施例的方法与deGoes2014(F.deGoes,P.Memari,P.Mullen,andM.Desbrun,“Weightedtriangulationsforgeometryprocessing,”ACMTrans.onGraphics,vol.33,no.3,pp.28:1–28:13,2014.)的加权三角形的方法对于输入网格进行重新网格化的结果对比图。图7比较了基于CVT的三种重新网格化方法，从左到右分别示出了：原始的输入网格、原始的CVT方法(CVT)的结果、特征敏感的CVT方法(CVT_fs)的结果和本发明实施例的方法(CVT_nob)的结果，其中第一行是Elk模型的均匀重新网格化结果；第二行是Bunny模型的自适应重新网格化结果。浅灰色对应包含小角度(＜30°)的三角形，黑色对应钝角三角形。从图7中可以看出CVT_fs方法会引入更多的坏三角形(＜30°或者＞90°)，而本发明实施例的方法产生的重新网格化结果不含有任何的坏三角形。

图8是本发明实施例的方法与deGoes(2014)的加权三角形的方法对于输入网格进行重新网格化的结果对比图。从左到右分别示出了：输入的原始网格、加权三角形方法的结果、本发明的结果。从该图中可以看出，加权三角形的结果中含有少量的坏三角形(浅灰色对应小角度三角形，黑色对应钝角三角形)，相比之下，本发明实施例的方法可以得到正确的结果。

本发明实施例通过对原始的CVT加入了对偶惩罚项，惩罚较短的Voronoi边，从而有效地去除钝角三角形。本发明实施例可以用于流体仿真等应用领域。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种基于重心Voronoi图的非钝角的重新网格化方法，其特征在于，该方法包括以下步骤：

步骤S1、获取待处理网格的采样点集；

步骤S2、对所述采样点集进行CVT优化；

步骤S3、对所述CVT优化后得到的网格的度数进行优化；

步骤S4、基于经过度数优化后得到的网格，进行特征敏感的CVT优化，并引入对偶惩罚项，惩罚较短的Voronoi边，以移除所述经过度数优化后得到的网格中的钝角三角形。

2.根据权利要求1所述的方法，其特征在于，所述步骤S1具体包括：由用户指定采样点数，根据定义在网格上的密度函数进行采样。

3.根据权利要求1所述的方法，其特征在于，所述步骤S3具体包括：对坏度数顶点引入扰动，对所述引入扰动后的坏度数顶点进行CVT优化，其中所述坏度数顶点为度数小于5或者大于7的顶点。

4.根据权利要求3所述的方法，其特征在于，所述对坏度数顶点引入扰动，具体包括：

在所述坏度数顶点的1-ring三角形中随机采样产生新的顶点，将所述新的顶点映射到所述待处理网格上，并将所述新的顶点插入Delaunay三角形中，同时，从所述Delaunay三角形中移除所述坏度数顶点。

5.根据权利要求1所述的方法，其特征在于，所述步骤S4具体包括：

将经过度数优化后得到的网格的全局能量函数表示为特征敏感的CVT能量函数与对偶惩罚项之和，对所述全局能量函数进行最小化，其中，

所述特征敏感的CVT能量函数表示为：

$F_{{\mathrm{CVT}}_{fs}}\left(X\right)=\underset{i=1}{\overset{n}{\&Sigma;}}\underset{f\&Element;{\mathrm{\&Omega;}}_{i|M}}{\&Sigma;}{\&Integral;}_f\mathrm{\&rho;}\left(x\right)\left|\right|A_s\left(N^f\right)(x-x_i)|{|}^{2}d\mathrm{\&sigma;};$

所述对偶惩罚项表示为：

$R\left(X\right)=\underset{i=1}{\overset{n}{\&Sigma;}}\left|\right|\underset{\&lsqb;x_i,x_j\&rsqb;\&Element;RDT}{\&Sigma;}{w}_{i,j}(x_i-x_j)|{|}^2;$

所述全局能量函数表示为：

$F\left(X\right)=-F_{{\mathrm{CVT}}_{fs}}\left(X\right)+\mathrm{\&lambda;}R\left(X\right);$

其中， ${\mathrm{\&Omega;}}_i=\{x\&Element;R^3|\left|\right|x-x_i\left|\right|\&le;\left|\right|x-x_j\left|\right|,\&ForAll;j\&NotEqual;i\}$ 是点x_i的Voronoi单元，M表示三角网格，Ω_i|M是定义在M上的限制Voronoi图；f表示三角形积分区域；ρ(x)(x∈M)表示密度函数；Α_s(Ν)＝(s-1)ΝΝ^t+Ι_3×3，Α_s(Ν)是一个3x3的矩阵；N表示三角形的法向量；Ν^f是用于积分的三角形f的法向量；N^t是N的转置；参数s控制法向各向异性的权重；I表示3x3的单位矩阵；λ表示惩罚项的权重，λ∈[0.1,1]；ε是一个很小的数，以避免被除数是0；dual(x_i,x_j)表示限制的Voronoi边，表征网格表面上的相邻的限制的Voronoi单元的线段序列。