CN109949420A - 适用于GPU的Delaunay三角剖分网格细化方法、GPU及系统 - Google Patents

Abstract

本公开提供了适用于GPU的Delaunay三角剖分网格细化方法、GPU及系统。其中，该细化方法包括步骤(1)：计算给定二维点集的Delaunay三角剖分网格，并按照顺序标记当前网格中各个三角形的序号；步骤(2)：判断Delaunay三角剖分网格中是否存在坏三角形，若存在，则标记Delaunay三角剖分网格中所有的坏三角形，进入下一步；否则，输出Delaunay三角剖分网格；步骤(3)：并行计算所有坏三角形的外接圆圆心，将这些圆心记为Steiner点；步骤(4)：在当前网格中并行插入所有Steiner点；步骤(5)：并行做翻转边操作；在做翻转边操作过程中，若包含两个Steiner点的多边形区域有重叠，则标记并删除其中一个冗余点，持续做翻转边操作，直至所有冗余点被删除，且三角剖分网格满足Delaunay属性，返回步骤(2)。

Description

适用于GPU的Delaunay三角剖分网格细化方法、GPU及系统

技术领域

本公开属于数据处理领域，尤其涉及一种适用于GPU的Delaunay三角剖分网格细化方法、GPU及系统。

背景技术

本部分的陈述仅仅是提供了与本公开相关的背景技术信息，不必然构成在先技术。

在生产和研究的过程中人们得到的空间数据往往杂乱且无规则，因此必须采用一定方法对空间散乱数据进行网格化处理，在此基础进行插值，才能够客观的反映出数据在空间内的演化规律。如图1所示的Delaunay三角剖分网格具有结构良好、易于更新、可适应各种分布密度的数据等特点，因此，Delaunay三角剖分网格成为重要的网格之一，可广泛应用于数值分析(例如有限元分析)、几何建模、图形学等领域。在实际应用中，人往往还需要生成的Delaunay三角剖分网格满足可以满足人们对于网格中三角形的数量、形状等的要求，如图2所示的狭长的三角形(通常需要被消除，此时人们需要对Delaunay三角剖分网格进行进一步的细化处理，即采用Delaunay三角剖分网格的细化算法以生成高质量的三角剖分网格。但是该算法实现比较复杂，使得在处理大规模的离散数据时，速度比较慢。

越来越多的学者将拥有强大浮点计算能力的图形处理器GPU(GraphicsProcessing Unit)列为加速算法的解决方案。为了能够在图形计算之外的更多领域发挥GPU强大的计算功能，人们开始研究如何能够利用GPU完成通常意义上的数据运算，这被称之为GPGPU(General-Purpose Computing on Graphics Processing Units，基于GPU的通用计算)。在程序开发上，GPU一般采用NIDIA推出的CUDA(Compute Unified DeviceArchitecture)，作为通用并行计算架构，该架构使GPU能够解决复杂的计算问题。

发明人发现如何在普通计算机上利用GPU编程实现面向海量数据的高质量的Delaunay三角网格的高效构建，是实现对离散数据进行研究与分析的关键问题。上述问题的解决方案仍然以下存在空白与缺陷：

1)传统的基于CPU的网格细化算法如circumcenter-insertion、sink-insertion和off-center-insertion，并不能直接应用于在GPU体系架构上。因为直接将串行算法应用于GPU上，很难保证算法在并行执行时的收敛性，使得算法很大程度上无法正确结束并得到正确的结果。并且计算几何的算法往往非常复杂，循环分支跳转极为频繁，且数据间的关联性很强，这一切使得当将经典的CPU算法并行地实施到GPU上时，往往并不能够得到高效且正确的算法结果。

2)传统的网格细化算法的并行算法，例如分而治之算法，也不能直接用于GPU体系。因为在并行算法直接套用于GPU算法时，在“分”的部分，很难使得每个核都负载均衡；在“治”的部分每个核面临的任务都非常繁复；而更困难的则是如何将分而治之的结构合并起来，且合并的过程往往需要大量的内存。此外，算法中需要进行大量的嵌套递归，而嵌套递归对于缺少cache的GPU是非常不利的。因此，目前也还没有该算法在GPU上成功实现的案例。

3)CPU拥有很强的通用性来处理各种不同类型的数据，同时又可兼顾逻辑判断、大量的分支跳转以及终端处理，而GPU的核相比于CPU的核是非常轻量级的，虽然拥有强大于CPU数倍的浮点计算能力，但只有非常简单的控制逻辑。只有面对的是类型高度统一，相互无依赖的大规模数据和不需要被中断的纯净的计算环境时，GPU才能发挥高效的运行速度。Delaunay三角剖分网格的细化算法，在算法上非常复杂，步骤很多，且元素之间彼此关联，因此要设计相关的高效的GPU算法，难度很大。

发明内容

根据本公开的一个或多个实施例的一个方面，提供一种适用于GPU的Delaunay三角剖分网格细化方法，其根据GPU独特的体系结构，不断并行插入Steiner点的方式减少坏三角形的个数，直至网格中所有的三角形都是Delaunay且不是坏三角形为止，可以大大加快细化网格的生成速度。

本公开的一种适用于GPU的Delaunay三角剖分网格细化方法的技术方案为：

一种适用于GPU的Delaunay三角剖分网格细化方法，包括：

步骤(1)：计算给定二维点集的Delaunay三角剖分网格，并按照顺序标记当前网格中各个三角形的序号；

步骤(2)：判断Delaunay三角剖分网格中是否存在坏三角形，若存在，则标记Delaunay三角剖分网格中所有的坏三角形，进入下一步；否则，输出Delaunay三角剖分网格；所述坏三角形满足三角形的半径与最短边比大于

步骤(3)：并行计算所有坏三角形的外接圆圆心，将这些圆心记为Steiner点；

步骤(4)：在当前网格中并行插入所有Steiner点；

步骤(5)：并行做翻转边操作；在做翻转边操作过程中，若包含两个Steiner点的多边形区域有重叠，则标记并删除其中一个冗余点，持续做翻转边操作，直至所有冗余点被删除，且三角剖分网格满足Delaunay属性，返回步骤(2)。

其中，Delaunay属性，是指的每个三角形的外接圆内都不包含其他数据点。

进一步地，所述步骤(2)中判断Delaunay三角剖分网格中是否存在坏三角形的过程为：

步骤(2-1)：设公共变量t＝0，一个初始长度为网格中三角形的个数n的整形数组isBad[]；

步骤(2-2)：对于每个三角形i，都分配一个GPU的线程thread，每个thread负责一个三角形，计算该三角形的半径r与最短边d之比ratio，如果则该三角形为坏三角形，则置t＝1,且isBad[i]＝1；

步骤(2-3)：若t＝1，说明存在坏三角形，进入步骤3，否则，输出Delaunay三角剖分网格。

进一步地，所述步骤(3)中并行计算所有坏三角形的外接圆圆心，将这些圆心记为Steiner点的过程为：

步骤(3-1)：设置初始长度为网格中三角形的个数n的point类型的数组circumcenter[]；

步骤(3-2)：判断isBad[i]＝1是否成立，若是，则该三角形是坏三角形，计算其外接圆的圆心，并存放在circumcenter[i]中。

进一步地，所述步骤(4)中在当前网格中并行插入所有Steiner点的具体过程为：

步骤(4-1)：定位Steiner点j；

步骤(4-2)：插入Steiner点j；

若点j位于某个三角形内，则将该三角形分成三个三角形；

若点j位于某条三角形的边上，且该边不是网格的边界边，则将彼此相邻的两个三角形分成四个三角形；

若点j位于某条三角形的边上且该边是网格的边界边，则将该三角形分成两个三角形；

若点j位于网格外部，则将边界三角形分成两个三角形；

步骤(4-3)：所有新增的三角形，均被标记为新的三角形，并行更新所有新的三角形的邻居信息。

更进一步地，在所述步骤(4-1)的Steiner点j定位过程中，若点j位于网格外部，且点j与网格的其它点分立于某个网格边界三角形的边ab的两侧；则用ab中点替代点j的坐标。

更进一步地，在所述步骤(4-2)中，采用1-3flip操作将该三角形分成三个三角形。

更进一步地，在所述步骤(4-2)中，采用2-4flip操作将彼此相邻的两个三角形分成四个三角形。

更进一步地，在所述步骤(4-2)中，采用1-2flip操作将该三角形分成两个三角形。

进一步地，所述步骤(5)中并行做翻转边操作的具体过程为：

步骤(5-1)：检查所有的新三角形，对于每个新三角形i，都分配一个GPU的线程thread，每个thread负责一个三角形，依次检查三角形i的三个邻居，如果其邻居与三角形i构成non-Delaunay，则标记该对三角形为需要做翻转边的三角形对；其中，non-Delaunay为三角形的外接圆内包含其他数据点；

步骤(5-2)：采用翻转边操作对所有标记的需要做翻转边的三角形对进行并行翻转；

步骤(5-3)：若包含两个Steiner点的多边形区域有重叠，则标记并删除其中一个冗余点；

步骤(5-4)：并行更新所有新的三角形的邻居信息。

进一步地，在所述步骤(5-2)中，一个三角形与一个或多个邻居形成翻转边的三角形对，依据邻居们的三角形序号，选择一个最小的进翻转边操作，其余的保留，等待下一轮迭代。

给每个翻转边三角形对分配一个thread，每个thread负责一个翻转边三角形对。对于每翻转边三角形对中的三角形abc和三角形cbd，如果点a和点d都是该轮插入时插入的Steiner点，说明它们其中之一是冗余点，应该被删除，选择其中序号较小的一个，将其标记为冗余点；否则将该两个三角形翻转为三角形adc和三角形dab。

所述步骤(5)中，所有步骤是迭代进行的，即4个步骤组成一个迭代iteration，迭代持续进行，直至再没有可以进行翻转边操作的翻转边三角形对。

进一步地，在所述步骤(5)中采用flip-flop操作删除冗余点，其具体过程为：

找到与当前冗余点相连的所有边，并沿着其中任意一条边开始遍历所有的以当前冗余点为顶点的三角形；

通过翻转边操作，将以当前冗余点为顶点的三角形个数逐渐减少为三个；

通过1-3flip的逆向操作，将与当前冗余点相连的三条边去除，也即将当前冗余点移除出网格。

根据本公开的一个或多个实施例的另一个方面，提供一种GPU，其利用GPU自身独特的体系结构，不断并行插入Steiner点的方式减少坏三角形的个数，直至网格中所有的三角形都是Delaunay且不是坏三角形为止，可以大大加快细化网格的生成速度。

本公开的一种GPU的技术方案为：

一种GPU，其上存储有计算机程序，该程序被GPU执行时实现上述所述的适用于GPU的Delaunay三角剖分网格细化方法中的步骤。

根据本公开的一个或多个实施例的另一个方面，提供一种Delaunay三角剖分网格细化系统，其利用GPU独特的体系结构，不断并行插入Steiner点的方式减少坏三角形的个数，直至网格中所有的三角形都是Delaunay且不是坏三角形为止，可以大大加快细化网格的生成速度。

本公开的一种Delaunay三角剖分网格细化系统的技术方案为：

一种Delaunay三角剖分网格细化系统，其包括上述所述的GPU。

本公开的有益效果是：

(1)对于任意满足条件的输入(合成数据或者源自真实世界的GIS数据，本公开的该细化方法均能鲁棒性的计算出正确结果，而且相比于CPU算法，该算法的运方法算速到提高了4-6倍。

(2)本公开的该细化方法可以大大加快细化网格的生成速度，该方法对于需要使用高质量的三角网格作为数据预处理的各个应用领域都具有非常重要的应用价值。

附图说明

构成本公开的一部分的说明书附图用来提供对本公开的进一步理解，本公开的示意性实施例及其说明用于解释本公开，并不构成对本公开的不当限定。

图1(a)为本公开实施例提供的Delaunay三角剖分网格示意图。

图1(b)为本公开实施例提供的Delaunay三角剖分网格中每个三角形符合Delaunay属性示意图。

图2(a)为本公开实施例提供的第一种狭长三角形示意图。

图2(b)为本公开实施例提供的第二种狭长三角形示意图。

图3为本公开实施例提供的随机二维点集示意图。

图4为本公开实施例提供的随机点集的Delaunay三角剖分网格示意图。

图5为本发明的计算坏三角形外接圆圆心的示意图。

图6为本公开实施例提供的1-3flip的示意图。

图7为本公开实施例提供的2-4flip的示意图。

图8为本公开实施例提供的1-2flip的示意图。

图9为本公开实施例提供的更新新的三角形的邻居信息的示意图。

图10为本公开实施例提供的edge-flip操作的示意图。

图11为本公开实施例提供的用flip-flop删除冗余点的示意图。

图12(a)为本公开实施例提供的随机分布示意图。

图12(b)为本公开实施例提供的高斯分布示意图。

图12(c)为本公开实施例提供的网格分布示意图。

图12(d)为本公开实施例提供的圆盘分布示意图。

图12(e)为本公开实施例提供的圆环分布示意图。

图13为本公开实施例提供的在真实世界GIS数据上的运算结果的示意图。

图14为本公开实施例提供的在真实世界图像数据上的运算结果的示意图。

图15(a)为本公开实施例提供的网格细化方法对不同分布的数据进行网格细化的运行时间的比较图。

图15(b)为本公开实施例提供的网格细化方法对不同分布的数据进行网格细化的加速比示意图。

图16为本公开实施例提供的一种适用于GPU的Delaunay三角剖分网格细化方法流程图。

具体实施方式

应该指出，以下详细说明都是例示性的，旨在对本公开提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本公开所属技术领域的普通技术人员通常理解的相同含义。

需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本公开的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。

术语解释部分：

Steiner点：它到其他所有点的距离之和最小。空间可以是二维或更多维，距离可以是欧氏距离，也可以是曼哈顿距离等。

1-3flip：是计算几何中的操作，其用于将一个三角形分成三个三角形。

2-4flip：是计算几何中的操作，其用于将彼此相邻的两个三角形分成四个三角形。

1-2flip：是计算几何中的操作，其用于将一个三角形分成两个三角形。

edge-flip：是计算几何中的操作，翻转边操作，其用于翻转边。

Delaunay属性是指的每个三角形的外接圆内都不包含其他数据点；

non-Delaunay是指三角形的外接圆内包含其他数据点，这里指的是，两个邻接的三角形(共4个点)，其中一个三角形的外接圆内包含第4点。

如图16所示，本实施例提供的一种适用于GPU的Delaunay三角剖分网格细化方法，包括：

步骤(1)：计算给定二维点集的Delaunay三角剖分网格，并按照顺序标记当前网格中各个三角形的序号。

如图1(a)所示，是二维点集的Delaunay三角剖分网格(左图)，图中的每一个三角形都符合Delaunay属性，即每个三角形的外接圆内都不包含其他数据点，如图1(b)所示。

三角形在构建时，都会有顺序的序号，例如第一个三角形是0号，后面的依次是1号、2号、3号……。

如图3所示，是随机生成的包含均匀分布点的二维点集。

二维点集内的二维点都是离散点，具有(x,y)坐标。几何建模数据中的点可以是地形数据的点；图形数据中的点可以是图象中点或者是特征点等。

在所述步骤(1)中，可采用任意算法，例如增量插入法等计算二维点集的Delaunay三角剖分网格。

步骤(2)：判断Delaunay三角剖分网格中是否存在坏三角形，若存在，则标记Delaunay三角剖分网格中所有的坏三角形，进入下一步；否则，输出Delaunay三角剖分网格；所述坏三角形满足三角形的半径与最短边比大于

如图2(a)和图2(b)所示，存在的这两种形式的狭长三角形，这两种三角形在本实施例内均称为坏三角形，需要采用本实施例提出的基于GPU的Delaunay三角剖分网格的细化方法进行消除。

具体地，所述步骤(2)中判断Delaunay三角剖分网格中是否存在坏三角形的过程为：

步骤(2-1)：设公共变量t＝0，一个初始长度为网格中三角形的个数n的整形数组isBad[]；

步骤(2-2)：对于每个三角形i，都分配一个GPU的线程thread，每个thread负责一个三角形，计算该三角形的半径r与最短边d之比ratio，如果则该三角形为坏三角形，则置t＝1,且isBad[i]＝1；

如果则该三角形中存在小于20.7度的角，或者存在大于138.6度角，该三角形为坏三角形。

步骤(2-3)：若t＝1，说明存在坏三角形，进入步骤3，否则，输出Delaunay三角剖分网格。

如图4所示，采用步骤(1)而生成的Delaunay三角剖分网格，此时网格中存在大量坏三角形，需要通过添加Steiner点进而消除。

步骤(3)：并行计算所有坏三角形的外接圆圆心，将这些圆心记为Steiner点。

具体地，所述步骤(3)中并行计算所有坏三角形的外接圆圆心，将这些圆心记为Steiner点的过程为：

步骤(3-1)：设置初始长度为网格中三角形的个数n的point类型的数组circumcenter[]；

步骤(3-2)：判断isBad[i]＝1是否成立，若是，则该三角形是坏三角形，计算其外接圆的圆心，并存放在circumcenter[i]中。

如图5所示，对于一个坏三角形t,计算它的外接圆圆心c。

步骤(4)：在当前网格中并行插入所有Steiner点。

具体地，所述步骤(4)中并行插入所有Steiner点的过程为：

步骤(4-1)：定位Steiner点j；

在所述步骤(4-1)的Steiner点j定位过程中，若点j位于网格外部，且点j与网格的其它点分立于某个网格边界三角形的边ab的两侧；则用ab中点替代点j的坐标。

步骤(4-2)：插入Steiner点j；

若点j位于某个三角形内，则用1-3flip将该三角形分成三个三角形；如图6所示，Steiner点位于某个三角形内，用1-3flip将一个三角形分成三个三角形。

若点j位于某条三角形的边上，且该边不是网格的边界边，则用2-4flip将彼此相邻的两个三角形分成四个三角形；如图7所示，Steiner点刚好位于某条三角形的边上(该边不是网格的边界边)，用2-4flip将彼此相邻的两个三角形分成四个三角形。

若点j位于某条三角形的边上且该边是网格的边界边，则用1-2flip将该三角形分成两个三角形；

若点j位于网格外部，则用1-2flip将边界三角形分成两个三角形；

如图8所示，Steiner点刚好位于某条三角形的边上(该边是网格的边界边)，或者该点落于网格外，此时用边界边的重点替代Steiner点，用1-2flip将边界三角形分成两个三角形。

所有新增的三角形，均被标记为新的三角形；

步骤(4-3)：并行更新所有新的三角形的邻居信息。

如图9所示，插入Steiner点后，需要更新三角形的邻居关系。例如在进行完1-3flip后，一个三角形t1被分成t1、t5和t6三个三角形。在进行1-3flip前，t1的邻居按照逆时针顺序是t2、t3和t4；在1-3flip后，t1的邻居变为t5、t6和t3,t5的邻居变为t4、t6和t1,t6的邻居变成t2、t1和t5。所有发生变化的三角形，例如t1、t2、t3、t4以及新生成的三角形t5和t6，他们的邻居信息都发生了改变，都需要在算法中进行更新。在插入Steiner点时，一共采用了1-3flip、2-4flip和1-2flip，一共三种类型的操作，这三种操作都导致了新的三角形的生成以及旧三角形的改变，只要三角形发生了改变，其邻居信息都要进行更新。

步骤(5)：并行做翻转边操作；在做翻转边操作过程中，若包含两个Steiner点的多边形区域有重叠，则标记并删除其中一个冗余点，持续做翻转边操作，直至所有冗余点被删除，且三角剖分网格满足Delaunay属性，返回步骤(2)。

具体地，所述步骤(5)中并行做翻转边操作的具体过程为：

步骤(5-1)：检查所有的新三角形，对于每个新三角形i，都分配一个GPU的线程thread，每个thread负责一个三角形，依次检查三角形i的三个邻居，如果其邻居与三角形i构成non-Delaunay，则标记该对三角形为需要做翻转边的三角形对；其中，non-Delaunay为三角形的外接圆内包含其他数据点；

如图10所示，一对三角形abc和cbd，其中点a位于三角形cbd的外接圆内，而点d位于三角形abc的外接圆内，所以这对三角形不符合Delaunay属性，需要用edge-flip操作进行边的翻转。在将边bc翻转为ad后，原三角形abc和cbd变成了三角形adc和dab。反转后，其中点b位于三角形adc的外接圆外，而点c位于三角形dab的外接圆外，此时两个三角形都满足Delaunay属性。

步骤(5-2)：采用翻转边操作-edge-flip操作对所有标记的需要做翻转边的三角形对进行并行翻转；

在所述步骤(5-2)中，一个三角形有可能与多于一个邻居形成edge-flip对，依据邻居们的三角形序号，选择一个最小的进行edge-flip，其余的保留，等待下一轮迭代。

给每个edge-flip对分配一个thread，每个thread负责一个edge-flip对。对于每个edge-flip对三角形abc和三角形cbd，如果点a和点d都是该轮插入时插入的Steiner点，说明它们其中之一是冗余点，应该被删除，选择其中序号较小的一个，将其标记为冗余点；否则将该两个三角形翻转为三角形adc和三角形dab；

步骤(5-3)：若包含两个Steiner点的多边形区域有重叠，则标记并删除其中一个冗余点；

在具体实施中，可采用flip-flop操作删除冗余点，其具体过程为：

找到与当前冗余点相连的所有边，并沿着其中任意一条边开始遍历所有的以当前冗余点为顶点的三角形；

通过翻转边操作，将以当前冗余点为顶点的三角形个数逐渐减少为三个；

通过1-3flip的逆向操作，将与当前冗余点相连的三条边去除，也即将当前冗余点移除出网格。

如图11所示，当一个点v被mark为冗余点后，我们可以通过flip-flop的方法将其去除，具体步骤为：找到与点v相连的所有边，并沿着其中任意一条边开始遍历所有的以点v为顶点的三角形。如图所示，开始时有va、vb、vc、vd和ve五条边连接着点v，也即有五个三角形以点v为顶点。通过翻转边ve为ad，可以将以v点为顶点的三角形个数减少为四个，再通过翻转边vc为bd，可以将以v点为顶点的三角形个数减少为三个，之后可以通过1-3flip的逆向操作，将与点v相连的三条边去除，也即将点v移除出网格。

步骤(5-4)：并行更新所有新的三角形的邻居信息。

所述步骤(5)中，所有步骤是迭代进行的，即4个步骤组成一个迭代iteration，迭代持续进行，直至再没有可以进行edge-flip的edge-flip对。

图12(a)-图12(e)为本实施例在合成数据上的运算结果的示意图。为了验证本实施例的该细化方法的有效性，我们随机生成了多种分布的数据，并在这些合成数据上测试了我们的算法。这些点的分布包含：随机分布如图12(a)所示、高斯分布如图12(b)所示、网格分布(如图12(c)所示，所有点均匀分布在512*512大小的网格上)、圆盘分布(如图12(d)所示，所有点均匀分布在一个圆盘上)以及圆环分布(如图12(e)所示,所有点均匀分布在圆环上)。

图13为本实施例在真实世界GIS数据上的运算结果的示意图。为了验证本实施例中提出的该细化方法的有效性，采用了GIS数据，并在这些真实的contour数据上测试了本实施例中的细化方法。

这些数据来自开源数据库http://www.ga.gov.au/。

图14为本实施例在真实世界图像数据上的运算结果的示意图。为了验证本实施例中提出的该细化方法的有效性，提取了图片中的数据，并在这些数据上测试了本实施例中的细化方法。

在图15(a)和图15(b)中，Uniform表示数据符合随机分布；Gaussian表示数据符合高斯分布；Disk表示数据符合圆盘分布；Cirele表示数据符合圆环分布；Grid表示数据符合网格分布。其中，在图15(b)中，speedup表示运算加速比。

图15(a)和图15(b)分别表示采用本实施例提供的网格细化方法对上述不同分布的数据进行网格细化的运行时间的比较图以及运算加速比。

表1为本实施例在真实世界GIS数据上的运算时间与传统算法的运行时间的比较表格。

表1运行时间的比较表格

对于任意满足条件的输入(合成数据或者源自真实世界的GIS数据，本实施例的该细化方法均能鲁棒性的计算出正确结果，而且相比于CPU算法，该算法的运方法算速到提高了4-6倍。

本实施例的该细化方法可以大大加快细化网格的生成速度，该方法对于需要使用高质量的三角网格作为数据预处理的各个应用领域都具有非常重要的应用价值。

在另一实施例中，提供了一种GPU，该GPU上存储有计算机程序，该程序被GPU执行时实现如图16所示的适用于GPU的Delaunay三角剖分网格细化方法中的步骤。

在另一实施例中，提供了一种Delaunay三角剖分网格细化系统，其包括上述所述的GPU。

本实施例的Delaunay三角剖分网格细化系统，利用GPU独特的体系结构，不断并行插入Steiner点的方式减少坏三角形的个数，直至网格中所有的三角形都是Delaunay且不是坏三角形为止，可以大大加快细化网格的生成速度。

本领域内的技术人员应明白，本公开的实施例可提供为方法、系统、或计算机程序产品。因此，本公开可采用硬件实施例、软件实施例、或结合软件和硬件方面的实施例的形式。而且，本公开可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器和光学存储器等)上实施的计算机程序产品的形式。

本公开是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体(Read-Only Memory，ROM)或随机存储记忆体(RandomAccessMemory，RAM)等。

上述虽然结合附图对本公开的具体实施方式进行了描述，但并非对本公开保护范围的限制，所属领域技术人员应该明白，在本公开的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本公开的保护范围以内。

Claims (11)

1.一种适用于GPU的Delaunay三角剖分网格细化方法，其特征在于，包括：

步骤(1)：计算给定二维点集的Delaunay三角剖分网格，并按照顺序标记当前网格中各个三角形的序号；

步骤(2)：判断Delaunay三角剖分网格中是否存在坏三角形，若存在，则标记Delaunay三角剖分网格中所有的坏三角形，进入下一步；否则，输出Delaunay三角剖分网格；所述坏三角形满足三角形的半径与最短边比大于

步骤(3)：并行计算所有坏三角形的外接圆圆心，将这些圆心记为Steiner点；

步骤(4)：在当前网格中并行插入所有Steiner点；

步骤(5)：并行做翻转边操作；在做翻转边操作过程中，若包含两个Steiner点的多边形区域有重叠，则标记并删除其中一个冗余点，持续做翻转边操作，直至所有冗余点被删除，且三角剖分网格满足Delaunay属性，返回步骤(2)。

2.如权利要求1所述的一种适用于GPU的Delaunay三角剖分网格细化方法，其特征在于，所述步骤(2)中判断Delaunay三角剖分网格中是否存在坏三角形的过程为：

步骤(2-1)：设公共变量t＝0，一个初始长度为网格中三角形的个数n的整形数组isBad[]；

步骤(2-2)：对于每个三角形i，都分配一个GPU的线程thread，每个thread负责一个三角形，计算该三角形的半径r与最短边d之比ratio，如果则该三角形为坏三角形，则置t＝1,且isBad[i]＝1；

步骤(2-3)：若t＝1，说明存在坏三角形，进入步骤3，否则，输出Delaunay三角剖分网格。

3.如权利要求1所述的一种适用于GPU的Delaunay三角剖分网格细化方法，其特征在于，所述步骤(3)中并行计算所有坏三角形的外接圆圆心，将这些圆心记为Steiner点的过程为：

步骤(3-1)：设置初始长度为网格中三角形的个数n的point类型的数组circumcenter[]；

步骤(3-2)：判断isBad[i]＝1是否成立，若是，则该三角形是坏三角形，计算其外接圆的圆心，并存放在circumcenter[i]中。

4.如权利要求1所述的一种适用于GPU的Delaunay三角剖分网格细化方法，其特征在于，所述步骤(4)中在当前网格中并行插入所有Steiner点的具体过程为：

步骤(4-1)：定位Steiner点j；

步骤(4-2)：插入Steiner点j；

若点j位于某个三角形内，则将该三角形分成三个三角形；

若点j位于某条三角形的边上，且该边不是网格的边界边，则将彼此相邻的两个三角形分成四个三角形；

若点j位于某条三角形的边上且该边是网格的边界边，则将该三角形分成两个三角形；

若点j位于网格外部，则将边界三角形分成两个三角形；

步骤(4-3)：所有新增的三角形，均被标记为新的三角形，并行更新所有新的三角形的邻居信息。

5.如权利要求4所述的一种适用于GPU的Delaunay三角剖分网格细化方法，其特征在于，在所述步骤(4-1)的Steiner点j定位过程中，若点j位于网格外部，且点j与网格的其它点分立于某个网格边界三角形的边ab的两侧；则用ab中点替代点j的坐标。

6.如权利要求4所述的一种适用于GPU的Delaunay三角剖分网格细化方法，其特征在于，在所述步骤(4-2)中，采用1-3flip操作将该三角形分成三个三角形；

或在所述步骤(4-2)中，采用2-4flip操作将彼此相邻的两个三角形分成四个三角形；

或在所述步骤(4-2)中，采用1-2flip操作将该三角形分成两个三角形。

7.如权利要求1所述的一种适用于GPU的Delaunay三角剖分网格细化方法，其特征在于，所述步骤(5)中并行做翻转边操作的具体过程为：

步骤(5-1)：检查所有的新三角形，对于每个新三角形i，都分配一个GPU的线程thread，每个thread负责一个三角形，依次检查三角形i的三个邻居，如果其邻居与三角形i构成non-Delaunay，则标记该对三角形为需要做翻转边的三角形对；其中，non-Delaunay为三角形的外接圆内包含其他数据点；

步骤(5-2)：采用翻转边操作对所有标记的需要做翻转边的三角形对进行并行翻转；

步骤(5-3)：若包含两个Steiner点的多边形区域有重叠，则标记并删除其中一个冗余点；

步骤(5-4)：并行更新所有新的三角形的邻居信息。

8.权利要求1所述的一种适用于GPU的Delaunay三角剖分网格细化方法，其特征在于，在所述步骤(5-2)中，一个三角形与一个或多个邻居形成翻转边的三角形对，依据邻居们的三角形序号，选择一个最小的进翻转边操作，其余的保留，等待下一轮迭代。

9.权利要求1所述的一种适用于GPU的Delaunay三角剖分网格细化方法，其特征在于，在所述步骤(5)中采用flip-flop操作删除冗余点，其具体过程为：

找到与当前冗余点相连的所有边，并沿着其中任意一条边开始遍历所有的以当前冗余点为顶点的三角形；

通过翻转边操作，将以当前冗余点为顶点的三角形个数逐渐减少为三个；

通过1-3flip的逆向操作，将与当前冗余点相连的三条边去除，也即将当前冗余点移除出网格。

10.GPU，其特征在于，其上存储有计算机程序，该程序被GPU执行时实现如权利要求1-8中任一项所述的适用于GPU的Delaunay三角剖分网格细化方法中的步骤。

11.Delaunay三角剖分网格细化系统，其特征在于，其包括如权利要求9所述的GPU。