基于最小角消除的三角形网格曲面的重新网格化方法

摘要

本发明涉及一种基于最小角消除的三角形网格曲面的重新网格化方法，包括：步骤S1，输入三角形网格模型M，设定最小角度阈值θ；构建初始输出网格模型M'；步骤S2，将M'中三角形存入动态优先队列Q；步骤S3，剔除动态优先队列Q中最小角大于θ的三角形，判断队列Q是否为空；如果是，执行步骤S6；如果否，执行步骤S4；步骤S4，选取队列顶部的三角形t，判断三角形t的最小角对应的边是否可以被折叠；如果是，执行边折叠操作；如果否，重定位顶点v；步骤S5，优化顶点v所在的位置，返回步骤S3；步骤S6，输出重新网格化后的新输出网格模型M”。本发明在保证近似误差的前提下，可以大大提高三角网格的质量。

说明书

技术领域

本发明属于计算机图形处理技术领域，具体涉及一种基于最小角度消除并且自动保持特征的重新网格化方法。

背景技术

在计算机图形学领域，三角网格是三维几何表示中最基本和使用最广泛的形式。目前三角网格模型的获取手段主要包括三维扫描重建，等值面提取(比如Marching Cube)，连续曲面的离散化，或者计算机视觉方法重建(比如Visual Hull)等。但是，这些网格往往具有几何采样不合理、存在大量狭长或者退化的三角形等缺点，不能直接用于比如有限元模拟、计算机动画等高级应用。重新网格化(Remeshing)，或者叫网格重剖，是在输入网格的基础上,在保持拓扑信息不变的前提下对网格的顶点采样、连接关系、面片形状等进行重构,从而获得满足用户或应用需求的结果。

一个高质量的网格模型，应该满足以下标准：(1)网格多边形质量高(内角差异较小)；(2)网格与原输入曲面的几何误差较小(例如Hausdorff距离小)；(3)方向保持，即网格边需要沿着特征、曲率或者用户指定的方向进行分布；(4)特征对齐，即模型上的特征需要被正确地表示为连续的网格边。我们从网格质量、逼近误差、特征保持等几个方面调研并叙述重新网格化的研究现状。

1、质量驱动的网格优化

目前，高质量的网格生成方法主要基于采样和CVT优化[Du99]。早期的方法主要在参数域中使用二维CVT进行优化，然后将结果映射到三维空间中[Alliez03,Alliez08,Fuhrmann10]；[Vorsatz03]则在参数域中采用粒子系统方法来取代CVT优化。一般来说，基于参数化的方法在生成高亏格的曲面时，通常会面临扭曲以及需要将参数曲面进行拼接等问题，导致拼接边界处网格质量较低。[Yan09,Yan15b]通过直接计算网格表面的三维CVT来避免参数化过程。此外，基于蓝噪声采样，他们使用自适应最大泊松盘采样[Yan13]和最远点优化进一步提高网格的质量[Yan14]。但是，这些方法同样也有其局限性，比如，几何保真度和最小三角形角度都不能进行显式控制。同时，为了实现特征保持，上述方法均需要提前检测或者手工指定特征边和特征点。

为了避免参数域方法引起的缝合问题，一些学者研究直接在网格表面进行操作[Frey00,Narain12]。其中，[Botsch04]提出了一种基于局部操作算子的网格优化方法。该方法以边界长度作为输入，通过不断地分离长边、合并短边、均衡点的度数以及调整点的位置，直到所有边都近似满足设定的长度。[Dunyach13]使用基于局部曲率确定的自适应的边长度来取代固定的目标边长度，将该方法扩展到自适应的版本。这类方法既不需要表面参数化过程又不需要密度函数，具有易于实现、适于高亏格网格以及高效率的优点。

2、误差驱动的网格优化

该类方法试图在几何保真度和网格复杂度之间取得最佳平衡，即在满足误差要求的前提下，使得网格复杂度尽量降低。[Cohen04]提出了一种误差驱动的聚类方法来简化网格，他们把网格简化问题转换为变分几何分割问题，通过Lloyd迭代方法来优化一系列的平面，以获得预先定义的近似误差。[Garland97]通过迭代地收缩顶点对来简化模型，并通过二次型函数度量来维持表面的近似误差。[Borouchaki05]定义了Hausdorff信封来度量保真度，简化和优化处于包围体内部的参考网格。虽然该方法同时考虑了几何保真度和元素质量，但是没有显式地提高最小网格化角度。本项目拟在保证最小网格化角度和给定误差界限的前提下，尽可能地降低网格的复杂度，进而达到三者的最佳平衡。基于容纳体的概念，[Mandad15]提出了一种同胚形体近似方法。使用该算法生成的三角形网格都处于给定的容纳体内部。容纳体的误差约束类似于重新网格化中的误差界限。但是，该方法主要解决的是网格的同胚逼近问题，并没有考虑元素的质量问题。

3、网格优化中的特征保持

特征保持是重新网格化算法需要重点解决的问题之一。但是三维网格模型的特征不仅取决于局部区域的几何形状，而且涉及到局部区域的上下文以及模型的语义信息。这使得特征检测成为一个极难解决的问题。虽然已有文献[Jiao02,Kolomenkin09]试图解决该问题，但目前尚没有一种方法对所有网格都有效。因此，很多重新网格化方法[Yan09,Zhong13]假定特征已被显式地检测或者手工定义出来。少量方法[Kobbelt03,Lévy10]视图隐式地保持特征，例如，[Vorsatz01]首先在参数域中优化包括特征点在内的网格顶点，然后重新将特征点恢复到原来位置以实现特征保持。因为网格优化和特征保持过程分开进行单独处理，靠近尖锐特征的网格形状质量会变差。[Valette08]为了消除这个缺陷，在CVT优化过程中嵌入了二次型误差度量。但是，该方法的效果过度依赖于输入网格的质量。此外，当用户需要稀疏网格时，尖锐特征可能不能很好地保持。

有关文献参见如下：

[Alliez03]Alliez,P.,De Verdire,E.C.,Devillers,O.,&Isenburg,M.(2003,May).Isotropic surface remeshing.In Shape Modeling International,2003(pp.49-58).IEEE.

[Alliez08]Alliez,P.,Ucelli,G.,Gotsman,C.,and Attene,M.(2008).Recent advances in remeshing of surfaces.In Shape Analysis and Structuring,pages 53–82

[Borouchaki05]H.Borouchaki and P.Frey,“Simplification of surface mesh using Hausdorff envelope,”Computer methods in applied mechanics and engineering,vol.194,no.48,pp.4864–4884,2005

[Botsch04]M.Botsch and L.Kobbelt,“A remeshing approach to multiresolution modeling,”in Symposium on Geometry Processing,2004,pp.185–192

[Cohen04]D.Cohen-Steiner,P.Alliez,and M.Desbrun,“Variational shape approximation,”ACM Trans.Graph.,vol.23,no.3,pp.905–914,2004

[Du99]Q.Du,V.Faber,M.Gunzburger,Centroidal Voronoi tessellations:applications and algorithms,SIAM Review 41(1999)637–676.

[Dunyach13]M.Dunyach,D.Vanderhaeghe,L.Barthe,and M.Botsch,“Adaptive remeshing for real-time mesh deformation,”Eurographics short papers.Eurographics Association,Girona,Spain,pp.29–32,2013

[Frey00]P.Frey and P.George,Mesh generation,application to finite elements.Hermès Science Publ.,Paris,Oxford,814 pages,2000.

[Fuhrmann10]S.Fuhrmann,J.Ackermann,T.Kalbe,and M.Goesele,“Direct resampling for isotropic surface remeshing,”in VMV,2010,pp.9–16.pp.20–30.

[Garland97]M.Garland and P.S.Heckbert,“Surface simplification using quadric error metrics,”in SIGGRAPH,1997,pp.209–216.

[Jiao02]X.Jiao and M.T.Heath,“Feature detection for surface meshes,”in Proceedings of 8th international conference on numerical grid generation in computational field simulations,2002,pp.705–714.

[Kobbelt03]L.Kobbelt and M.Botsch,“Feature sensitive mesh processing,”in Proceedings of the 19th Spring Conference on Computer Graphics,ser.SCCG’03,2003,pp.17–22.

[Kolomenkin09]M.Kolomenkin,I.Shimshoni,and A.Tal,“On edge detection on surfaces,”in Computer Vision and Pattern Recognition,2009.CVPR 2009.IEEE Conference on.IEEE,2009,pp.2767–2774.

[Lévy10]B.Lévy and Y.Liu,“L p centroidal Voronoi tesselation and its applications,”vol.29,no.4,pp.119:1–11,2010.

[Mandad15]M.Mandad,D.Cohen-Steiner,and P.Alliez,“Isotopic approximation within a tolerance volume,”ACM Trans.Graph.,vol.34,no.4,pp.64:1–64:12,2015

[Valette08]S.Valette,J.-M.Chassery,and R.Prost,“Generic remeshing of 3d triangular meshes with metric-dependent discrete Voronoi diagrams,”IEEE Trans.Vis.Comput.Graph,vol.14,no.2,pp.369–381,2008

[Vorsatz01]J.Vorsatz,C.L.Kobbelt,and H.-P.Seidel,“Feature sensitive remeshing,”Comput.Graph.Forum,vol.20,no.3,pp.393–401,2001

[Vorsatz03]J.Vorsatz,C.and H.-P.Seidel,“Dynamic remeshing and applications,”J.Comput.Inf.Sci.Eng.,vol.3,no.4,pp.338–344,2003.

[Yan09]Yan,D.M.,Lévy,B.,Liu,Y.,Sun,F.,&Wang,W.(2009,July).Isotropic remeshing with fast and exact computation of restricted Voronoi diagram.In Computer graphics forum(Vol.28,No.5,pp.1445-1454).

[Yan13]D.-M.Yan,P.Wonka,Gap Processing for Adaptive Maximal Poisson-Disk Sam-pling,ACM Trans.on Graphics 32(5)(2013)148:1–148:15.

[Yan14]D.-M.Yan,J.Guo,X.Jia,X.Zhang,P.Wonka,Blue-Noise Remeshing with Farthest Point Optimization,Computer Graphics Forum(Proc.SGP)33(5)(2014)167–176.

[Yan15b]D.Yan and P.Wonka,“Non-obtuse remeshing with centroidal voronoi tessellation,”IEEE Transactions on Visualization and Computer Graphics,2015.

[Zhong13]Z.Zhong,X.Guo,W.Wang,B.Lévy,F.Sun,Y.Liu,and W.Mao,“Particle-based anisotropic surface meshing,”ACM Trans.Graph.,vol.32,no.4,p.99,2013.

发明内容

为了解决现有技术中的上述问题，本发明提出一种基于最小角消除的三角形网格曲面的重新网格化方法，在保证近似误差的前提下，可以大大提高三角网格的质量。

为实现上述目的，本发明提的一种基于最小角消除的三角形网格曲面的重新网格化方法，包括以下步骤：

步骤S1，输入三角形网格模型M，设定最小角度阈值θ；构建初始输出网格模型M'；

步骤S2，将M'中三角形存入动态优先队列Q；

步骤S3，剔除动态优先队列Q中最小角大于θ的三角形，判断队列Q是否为空；

如果是，执行步骤S6；

如果否，执行步骤S4；

步骤S4，选取队列顶部的三角形t，判断三角形t的最小角对应的边是否可以被折叠；

如果是，执行边折叠操作；

如果否，重定位顶点v；

步骤S5，优化顶点v所在的位置，返回步骤S3；

步骤S6，输出重新网格化后的新输出网格模型M”。

优选的，所述的初始输出网格模型M'与步骤S1中输入三角形网格模型M一致。

优选的，所述的初始输出网格模型M'的构建方法为：

步骤S11，设定近似误差阈值δ；令M'＝M；

步骤S12，计算M和M'之间的豪斯多夫距离dH(M，M′)，并判断dH(M，M′)是否小于δ，若不小于则执行步骤S13，若小于则执行步骤S14；

步骤S13，对输入三角形网格模型M进行简化，并令M'＝M；执行步骤S12；

步骤S14，输出初始输出网格模型M'。

优选的，M和M'之间的豪斯多夫距离dH(M，M′)的计算方法为：

d(p，M′)＝min_p′∈M′d(p，p′)

dh(M，M′)＝max_p∈M(d(p，M′))

dH(M，M′)＝max(dh(M，M′)，dh(M′，M))

其中，d(p，p′)为三维空间中两个点p、p′的欧式距离，p∈M，p′∈M′，dh(M，M′)为M和M′之间的单向豪斯多夫距离。

优选的，动态优先队列Q按照其内各三角形的最小角的数值升序排列。

优选的，步骤S4中判断三角形t的最小角对应的边是否可以被折叠的方法为：

步骤S411，判断边折叠操作能否引起h的1-邻域三角形发生翻转；若是则执行步骤S414；

步骤S412，判断M和M'的豪斯多夫距离是否小于δ；若否则执行步骤S414；

步骤S413，计算执行边折叠操作造成的最小角，并判断该最小角是否小于执行操作之前的最小角；若是则执行步骤S414；若否则输出判断结果为三角形t的最小角对应的边可以被折叠；

步骤S414，输出判断结果为三角形t的最小角对应的边不可以被折叠。

优选的，步骤S4中重定位顶点v的方法为：

步骤S421，判断执行重定位顶点操作是否能引起三角形t中的某个顶点vt的1-邻域三角形发生翻转；若是则执行步骤S424；

步骤S422，判断M和M'的豪斯多夫距离是否小于δ；若否则执行步骤S424；

步骤S423，计算执行重定位顶点操作造成的最小角，并判断该最小角是否小于执行操作之前的最小角；若是则执行步骤S424；若否则输出判断结果为将v设置为三角形t的顶点vt对应的重心Voronoi图的重心；

步骤S424，输出判断结果为将v设置为沿着h传播到局部最长边hl的中点。

优选的，步骤S5中优化顶点v所在的位置的方法为：

优化顶点v所在的位置通过最小化能量函数F计算得到，

其中，S_in和S_out分别为v在M和M'的局部点集集合，a_i是S_out中的采样点，是a_i的最近的点，是S_in中的采样点，b_j是的参考点，和是a_i与b_j各自的权重。

优选的，最小化能量函数F的方法为：

步骤S51，计算最优位置v^*，

其中，是a_i所在的三角形(v，v_i1，v_i2)的重心，是所在的三角形(v，v_j1，v_j2)的重心。

步骤S52，移动v到v+μ(v^*-v)，其中权值μ∈(0，1]；

步骤S53，判断顶点v是否收敛到最优的位置，若否则执行步骤S51，若是则顶点v所在的位置。

本发明在保证近似误差的前提下，可以大大提高三角网格的质量。解决了现有方法不能同时既满足近似误差又提高网格质量这两个方面的缺陷。根据用户指定的误差限制以及最小角度，利用一个动态优先队列，不断调整最小角对应的三角形顶点位置或者连接关系，使得网格的整体质量不断提高，最终可以应用于物理模型，变形动画等应用中。

附图说明

图1是本发明最小角消除的三角形网格曲面的重新网格化方法的流程图；

图2A是边折叠操作引起的三角形翻转的示意图；

图2B是顶点重定位操作引起的三角形翻转的示意图；

图3是最小角度阈值θ对重新网格化结果的影响示意图；

图4A是近似误差δ对网格复杂性的影响示意图；

图4B是近似误差δ对运行时间的影响示意图；

图5是本发明的方法与其他的方法的比较示意图。

具体实施方式

下面参照附图来描述本发明的优选实施方式。本领域技术人员应当理解的是，这些实施方式仅仅用于解释本发明的技术原理，并非旨在限制本发明的保护范围。

如图1所示，本发明的基于最小角消除的三角形网格曲面的重新网格化方法，包括以下步骤：

步骤S1，输入三角形网格模型M，设定最小角度阈值θ；构建初始输出网格模型M'；

步骤S2，将M'中三角形存入动态优先队列Q；

步骤S3，剔除动态优先队列Q中最小角大于θ的三角形，判断队列Q是否为空；

如果是，执行步骤S6；

如果否，执行步骤S4；

步骤S4，选取队列顶部的三角形t，判断三角形t的最小角对应的边是否可以被折叠；

如果是，执行边折叠操作；

如果否，重定位顶点v；

步骤S5，优化顶点v所在的位置，返回步骤S3；

步骤S6，输出重新网格化后的新输出网格模型M”。

下面对上述每个步骤进行更加详细的说明。

步骤S1，输入三角形网格模型M，设定最小角度阈值θ；构建初始输出网格模型M'。

本实施例中，初始输出网格模型M'与输入三角形网格模型M一致，即初始化输出三角网格模型M'＝M。

如果输入模型M过于稠密，可以通过保持误差的网格简化方法先对M'进行一定简化，既保证dH(M，M′)＜δ的前提下，尽量减少M'的顶点数目；其中，dH(M，M′)是网格模型M和M'之间的豪斯多夫(Hausdorff)距离，δ为设定的近似误差阈值。

据此，在输入模型M过于稠密的情况下，初始输出网格模型M'的构建方法为：

步骤S11，设定近似误差阈值δ；令M'＝M；

步骤S12，计算M和M'之间的豪斯多夫距离dH(M，M′)，并判断dH(M，M′)是否小于δ，若不小于则执行步骤S13，若小于则执行步骤S14；

M和M'之间的豪斯多夫距离dH(M，M′)的计算如公式(1)、(2)、(3)所示：

d(p，M′)＝min_p′∈M′d(p，p′)>

dh(M，M′)＝max_p∈M(d(p，M′))>

dH(M，M′)＝max(dh(M，M′)，dh(M′，M)) (3)

其中，d(p，p′)为三维空间中两个点p、p′的欧式距离，p∈M，p′∈M′，dh(M，M′)为M和M′之间的单向豪斯多夫距离。

步骤S13，对输入三角形网格模型M进行简化，并令M'＝M；执行步骤S12；

步骤S14，输出初始输出网格模型M'。

步骤S2，将M'中三角形存入动态优先队列Q；

动态优先队列Q按照其内各三角形的最小角的数值升序排列。

步骤S3，剔除动态优先队列Q中最小角大于θ的三角形，判断队列Q是否为空；如果是，执行步骤S6；如果否，执行步骤S4；

这一步用来判断是否达到满足用户要求的重新网格化。其中动态优先队列Q更新后仍按照其内各三角形的最小角的数值升序排列。

步骤S4，选取队列顶部的三角形t，判断三角形t的最小角对应的边是否可以被折叠；如果是，执行边折叠操作；如果否，重定位顶点v；

这个步骤中，首先从队列顶部选择三角形t，三角形t具有当前动态优先队列Q内所有三角形中的最小角度θ_min，然后做进一步优化。具体的，针对该三角形的顶点做一些局部操作，并且进一步更新Q，将局部操作所影响的三角形重新插入Q中。

本实施例中判断三角形t的最小角对应的边是否可以被折叠的方法为：

步骤S411，判断边折叠操作能否引起h的1-邻域三角形发生翻转，；若是则执行步骤S414；

如图2A的示例可以看出，黑色边PQ被折叠，黑色三角形STQ发生翻转；

步骤S412，判断M和M'的豪斯多夫距离是否小于δ；若否则执行步骤S414；

步骤S413，计算执行边折叠操作造成的最小角，并判断该最小角是否小于执行操作之前的最小角；若是则执行步骤S414；若否则输出判断结果为三角形t的最小角对应的边可以被折叠；

步骤S414，输出判断结果为三角形t的最小角对应的边不可以被折叠。

本实施例中，重定位顶点v的方法为：

步骤S421，判断执行重定位顶点操作是否能引起三角形t中的某个顶点vt的1-邻域三角形发生翻转；若是则执行步骤S424；

如图2B示例，重定位顶点P，造成黑色三角形STQ发生翻转；

步骤S422，判断M和M'的豪斯多夫距离是否小于δ；若否则执行步骤S424；

步骤S423，计算执行重定位顶点操作造成的最小角，并判断该最小角是否小于执行操作之前的最小角；若是则执行步骤S424；若否则输出判断结果为将v设置为三角形t的顶点vt对应的重心Voronoi图的重心；其中Voronoi图又叫泰森多边形或Dirichlet图，它是由一组由连接两邻点直线的垂直平分线组成的连续多边形组成。

步骤S424，输出判断结果为将v设置为沿着h传播到局部最长边hl的中点。

传播过程具体如下：h₁是构成最小角θ_min的两条边中的较长者，迭代地更新h₁：如果h₁是边界边或者h₁相邻的两个三角形中的最长边，则更新过程终止；否则，h₁被更新为它相邻的两个三角形中的最长边。这一更新过程继续进行，直到上面所述的两个终止条件的其中一个被满足。

步骤S5，优化顶点v所在的位置，返回步骤S3；

优化顶点v所在的位置通过最小化能量函数F计算得到，能量函数F如公式(4)所示：

其中，S_in和S_out分别为v在M和M'的局部点集集合，a_i是S_out中的采样点，是a_i的最近的点，是S_in中的采样点，b_j是的参考点，和是a_i与b_j各自的权重。

上面的能量函数F通过两步迭代方式进行优化：

步骤S51，计算最优位置v^*，如公式(5)、(6)所示：

其中，是a_i所在的三角形(v，v_i1，v_i2)的重心，是所在的三角形(v，v_j1，v_j2)的重心。

步骤S52，移动v到v+μ(v^*-v)，其中权值μ∈(0，1]；

步骤S53，判断顶点v是否收敛到最优的位置，若否则执行步骤S51，若是则顶点v所在的位置。

步骤S6，输出重新网格化后的新输出网格模型M”。

另外，本发明还使用一个3.4GHZ CPU和16G内存的台式机上进行了算法的测试和分析；

最小角度阈值θ对结果的影响：θ越大，元素品质越好。θ从0°到40°，Fandisk模型的重新网格化结果如图3所示。最左边是输入模型，剩下的依次为θ＝0°，10°，20°，30°，35°，40°。图3下部的直方图统计了角度分布(0°～60°)，曲线是近似误差。当θ≤35°时，网格的复杂性和程序的运行时间平滑的增加。但是，当θ＝40°时，二者显著增大。同时，随着θ的增大，近似误差不断地减小。

近似误差δ对结果的影响：如图4A所示，展示了近似误差δ对重新网格化结果的网格复杂性的影响：横坐标表示近似误差δ，纵坐标表示网格复杂性(顶点数目)。图4B展示了近似误差δ对运行时间的影响：横坐标表示近似误差δ，纵坐标表示运行时间。结果显示，给定θ，近似误差δ没有明显地影响网格的复杂性和运行时间。但是，有两个有趣的现象被发现：

(1)当θ≤35°时，θ越大，网格的复杂性越低。但是，当θ＝4θ°，网格复杂性和运行时间显著地增加。这是因为当θ很小时，边折叠操作被优先地执行，增加最小角的同时降低了网格的复杂性。然而，当θ较大时，更多的边剖分操作被执行，改变了局部拓扑。

(2)给定一个θ，当δ增加时，网格的复杂性轻微地减小。这是因为，近似误差越松弛，更多的边折叠操作被执行。

如图5所示，比较了RAR[Dunyach13]，MMGS[Frey00]，CVT[Yan09]，MPS[Yan13]和本算法的重新网格化结果。最左边的是输入的Triceratops模型，然后依次是：RAR，MMGS，CVT，MPS和本发明方法的结果。从比较结果可以看出，本方法的结果在网格复杂性和三角形的品质上都是最优的。

本发明的方法的特色和创新在于，针对现有方法不能同时既满足近似误差又提高网格质量这两个方面的缺陷，提出了一种基于最小角消除的重新网格化方法。根据用户指定的误差限制以及最小角度，利用一个动态优先队列，不断调整最小角对应的三角形顶点位置或者连接关系，在保证近似误差的前提下，使得网格的整体质量不断提高。可以用于数值仿真。

本领域技术人员应该能够意识到，结合本文中所公开的实施例描述的各示例的模块、单元及方法步骤，能够以电子硬件、计算机软件或者二者的结合来实现，为了清楚地说明电子硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以电子硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专本领域技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。

至此，已经结合附图所示的优选实施方式描述了本发明的技术方案，但是，本领域技术人员容易理解的是，本发明的保护范围显然不局限于这些具体实施方式。在不偏离本发明的原理的前提下，本领域技术人员可以对相关技术特征作出等同的更改或替换，这些更改或替换之后的技术方案都将落入本发明的保护范围之内。