一种面向任意分布大规模点云数据的快速Delaunay构网方法

摘要

本发明提供一种面向任意分布大规模点云数据的快速Delaunay构网方法，包括以下步骤：划分多重网格，创建初始三角网；逐点插入；三角网更新；删除辅助三角形完成Delaunay三角网的构建。本发明通过多重网格划分和Hilbert曲线顺序遍历网格的排序方法，减少了点定位过程的搜索步长以及构网过程中产生的最终要删除的狭长三角形数量，并且通过对点定位过程中，采用新的算法，避免了求三角形重心和相交边的过程，减少了计算量。大大提高了点云在数据量较大以及分布不均匀的情况下构建Delaunay三角网的效率。

说明书

技术领域

本发明涉及信息处理技术领域，尤其涉及一种面向任意分布大规模点云数据的快速Delaunay构网方法。

背景技术

三角剖分是计算几何领域的重要研究内容之一，而其中Delaunay三角剖分由于具有许多优良性质，在地理信息系统、数字高程插值、有限元分析，几何重构等领域有着广泛的应用。

Delaunay三角网的构建方法主要可以分为逐点插入法、分治算法和三角网生长法三种。其中三角网生长法由于要频繁遍历数据点，时间效率低且算法复杂，目前已很少采用。分治算法的优点是时间效率高，但占用空间大，算法复杂。虽然分治算法可以通过并行计算的方法使算法效率得到很大提高，但是并行计算对硬件要求较高，不具有普遍性。逐点插入法由于其算法简单、占用空间小等特点成为了目前最流行的一种算法。近年来，又有许多学者在传统逐点插入法的基础上进行了改进。徐明海等利用存储三角形的邻接关系的方法在Bowyer和Watson的基础上对逐点插入法进行改进，使插入点进行外接圆判断的次数大大减少，算法效率有了很大提高，也使算法效率更大程度上取决于点插入三角网的顺序。徐道柱等提出了一种对点集划分均匀的网格，按网格加点的方法，可以使时间进一步降低。Liu和Snoeyink、Zhou和Jones、Hornus和Boissonnat、Buchin在划分均匀网格的基础上提出了依照不同的空间填充曲线(如Hilbert曲线、Peano曲线等)顺序遍历插入网格中点的方法，进一步提高了划分均匀网格加点的效率，对均匀分布点集的Delaunay构网有较好的效果。

现有方法虽然比较成熟，但仍普遍存在许多不足之处，如：

第一：处理对象针对少量离散点云对算法的处理效率要求不高，但当点云的数目较大时，算法处理效率普遍较低，对构网影响很大；

第二：由于点云数据的分布并不均匀，规则网格方法在网格划分过程中会形成很多点集分布过密或者过疏的网格，使网格划分程度无法统一，影响了后续的处理效率；

第三：目前的点定位过程中计算重心坐标求相交边的过程大大增加了计算量，降低了对海量点云数据Delaunay构网的速度。

发明内容

针对现有技术中存在的上述的问题，本发明提出一种面向任意分布大规模点云数据的快速Delaunay构网方法，对现有算法进行了改进优化，大大提高了构网效率，其具体的技术方案如下：

一种面向任意分布大规模点云数据的快速Delaunay构网方法，其特征在于：包括以下步骤：

步骤一：划分多重网格，并以Hilbert曲线遍历网格的方式，对插入点的顺序进行排序：

A：建立控制点数组和非控制点数组，用于存储点排序数据；

B：根据数据点的个数建立一级网格，并将点分配到网格中；

C：对划分好的一级网格按照Hilbert曲线从左下角到右下角的顺序进行排序；

D：对排序后的一级网格，按顺序对每个网格进行控制点和非控制点的提取；

E：设定网格中的点数的上限阈值，如果当前网格中的点数超过该阈值，则按照步骤B至步骤D网格迭代划分，在迭代划分过程中，对网格开口方向进行调整，但不进行控制点的提取。

步骤二：创建初始三角网：

通过选择把数据点坐标最大值加一个较小的数，最小值减一个较小的数，找到能够包围所有数据点的初始矩形凸壳连接初始矩形凸壳的一条对角线，将其分为两个三角形，作为初始三角网，并为初始三角网添加邻接关系；

步骤三：逐点插入：按照先插入控制点数组，插入顺序是原存储顺序的逆序、后插入非控制点数组，插入顺序与原存储顺序相同，逐点进行如下操作：

A：对插入点进行定位；

B：构建包含当前插入点所有影响域的空腔；把当前插入点所在三角形的三条边作为空腔边，分别对空腔边的邻接三角形做外接圆检测；若当前插入点位于某三角形外接圆内部，则将该三角形对应的空腔边删除，将该三角形的另两边加入到空腔边中，继续对空腔边进行判断，直到所有空腔边的邻接三角形的外接圆都不包含当前插入点；

C：三角网更新：将当前插入点与得到的空腔节点相连，形成新的三角形，并更新由空腔形成的新三角形及其相邻三角形的邻接关系；

步骤四：删除辅助三角形：全部点集插入后，将顶点中含有矩形凸壳四个顶点的三角形删除，即完成了Delaunay三角网的构建。

进一步地，步骤一中的步骤B中建立一级网格，并将点分配到网格中按如下方法进行：

A：设定平均每个网格大致的点数，用总点数除以平均点数即为网格的大致个数；

B：将网格的个数调整至4的整数次幂的形式，并保持行与列的网格数相同；

C：以点云数据中x,y坐标的最大最小值为边界，划分网格，并根据每个数据点的位置将点分配到一级网格中去。

进一步地，步骤一中的步骤D中对于控制点和非控制点的提取按如下方法进行：

A：如果当前网格有点，则将第一个点放入控制点数组；

B：设定网格中的点数的上限阈值，如果当前网格中的点数小于等于所设定的阈值，则将除控制点外的其他点放入非控制点数组。

进一步地，步骤一中的步骤E中对网格的迭代划分按如下方法进行：

A：若当前网格中的点数超过阈值，则对目前操作的网格区域，按照与一级网格相同的划分规则，进行下级网格的迭代划分；

B：在迭代划分过程中，不再进行控制点的提取，已排序的网格中的点全部放入非控制点数组；

C：在对当前网格进行Hilbert排序时，通过Hilbert曲线的开口方向调整，保证对于总体的遍历顺序，把所有网格按照放入非控制点集的顺序连接起来，所形成的Hilbert曲线没有交叉。

进一步地，步骤三的步骤A中，对插入点进行定位按如下方法进行：

A：.找到最新更新的三角形，即对上一个插入点进行处理时三角网中最后一个新生成的三角形。设该三角形顶点为A、B、C，沿逆时针排列；

B：设当前插入点坐标为P，P与三角形三边AB、BC、CA构成的三角形面积分别为$S_A=\frac{1}{2}\overrightarrow{\mathrm{AB}}\times \overrightarrow{\mathrm{AP}},S_B=\frac{1}{2}\overrightarrow{\mathrm{BC}}\times \overrightarrow{\mathrm{BP}},S_C=\frac{1}{2}\overrightarrow{\mathrm{CA}}\times \overrightarrow{\mathrm{CP}},$计算S_A，S_B，S_C，判断其中小于0的个数。若S_A，S_B，S_C均大于等于0，则ΔABC即为当前待插点P所在三角形，点定位过程结束；当S_A，S_B，S_C有一个小于0时，进行步骤C，当S_A，S_B，S_C有两个小于0时，进行步骤D；

C：找到相应小于0的S，此时P应位于该小于0的S所对应的边的外侧，把该边的临边三角形作为ΔABC，返回步骤B继续判断；

D：找到相应小于0的两个S中最小的S，此时把最小的S所对应的边的临边三角形作为ΔABC，返回步骤B继续判断。

本发明的有益效果体现在：

第一：在本发明中采用Hilbert曲线遍历网格的方式，且在遍历过程中调整曲线开口方向，减少了点定位过程的搜索步长；采用添加控制点的方式，避免一次性插入网格中所有点时插入点与初始矩形凸壳顶点以及已完成Delaunay构网的区域边界形成的大量狭长三角形，降低了构建和删除这些三角形占用的额外时间，从而提高构网效率；

第二：通过对点定位过程中，采用新的算法，避免了求三角形重心和相交边的过程，减少了计算量。实现了大规模点云数据的Delaunay三角网快速构建，大大提高了点云在数据量较大以及分布不均匀的情况下构建Delaunay三角网的速度，本发明的算法在构网效率上，具有明显优势。

附图说明

图1为本发明Delaunay整体构网过程流程图；

图2为本发明点排序过程流程图；

图3为本发明点定位过程流程图；

图4为本发明实施例中点云数据局部形成的最终网格遍历顺序图；

图5为本发明点定位步骤中点P位于当前判断的三角形内的示意图；

图6为本发明点定位步骤中点P位于当前判断的三角形一条边外侧的示意图；

图7为本发明点定位步骤中点P位于当前判断的三角形两条边外侧的示意图；

图8为本发明实施例中局部点云数据点定位过程示意图；

图9为本发明实施例中点云数据最终Delaunay构网结果示意图(抽稀10倍)。

具体实施方式

下面结合附图及实施例对本发明的一种面向任意分布大规模点云数据的快速Delaunay构网方法作进一步说明。

本实施案例面向分布较为均匀、27013896个点的点云数据，按照本发明图1的Delaunay整体构网过程流程图中的步骤进行构网:

第一步：通过划分多重网格和以Hilbert曲线遍历网格的方式，对插入点的顺序进行排序，以提高后续过程的处理速度。

其中Hilbert曲线是德国数学家DavidHilbert发现了一种曲线：首先把一个正方形等分成四个小正方形，依次从西南角的正方形中心出发往北到西北正方形中心，再往东到东北角的正方形中心，再往南到东南角正方形中心，这是一次迭代，如果对四个小正方形继续上述过程，往下划分，反复进行，最终就得到一条可以填满整个正方形的曲线。

根据本案例的点数及数据点的分布位置，具体排序过程如下：

(1)建立控制点数组和非控制点数组，用于存储点排序数据。

(2)根据数据点的个数建立一级网格，并将点分配到网格中，具体如下：

a.求出数据点中横坐标、纵坐标的最大最小值，得到X_min＝3685.13，X_max＝4030.06，Y_min＝4441.82，Y_max＝4803.90，并把以此为边界划分一级网格。

b.设置平均每个网格中大致的点数为10。则要划分的一级网格的个数大致为：将网格的个数调整至4的整数次幂，即4194304个网格，每行和每列2048个网格。

c.根据网格数求得行宽和列宽为：

$\mathrm{rowwidth}=\frac{Y_{\max }-Y_{\min }}{2048}=0.1768,\mathrm{colwidth}=\frac{X_{\max }-X_{\min }}{2048}=0.1684,$

进而可以根据每个数据点的坐标(x_i,y_i)得到其所在的网格为横排第，竖排第个网格，将每个数据点都分配到其对应的网格中去。

(3)对划分好的一级网格按照Hilbert曲线从左下角到右下角的顺序进行排序。

(4)对排序后的网格，按顺序对每个网格进行如下操作：

a.如果当前网格有点，则将第一个点放入控制点数组；例如附图4中的一级网格G0、G1、G2、G3中均有点，则这些网格中的第一个点都要放到控制点数组中。

b.如果当前网格中的点数小于等于100，则将除控制点外的其他点放入非控制点数组。例如附图4中的一级网格G0、G1中点数小于100，则将G0、G1中除控制点之外的所有点放入非控制点数组。

c.如果当前网格中的点数超过100，则对目前操作的网格区域，按照步骤(2)至(4)的过程进行下一级网格的迭代划分。在迭代划分过程中，不再进行控制点的提取，且在对当前网格进行Hilbert排序时，要通过Hilbert曲线的开口方向调整，保证对于总体的遍历顺序，把所有网格按照放入非控制点集的顺序连接起来，所形成的Hilbert曲线没有交叉。

例如附图4中的一级网格、G3中点数大于等于100，则在G2所在区域继续迭代划分下级网格，按照划分规则下一级网格有16个，，调整Hilbert曲线的起止方向使其从左下角到左上角以使整体Hilbert曲线没有交叉，按照顺序遍历这16个网格，在遍历过程中同时判断出这16个网格中没有点数仍然超过100的，则不必继续划分网格，将16个网格中的点依次放入非控制点集，此时对一级网格G2处理完毕。之后按照同样的划分规则处理G3。以此类推，直到所有一级网格处理完毕。

Hilbert曲线的开口方向调整过程参照下表进行，其中：表中“前一网格与当前网格的位置关系”的“-”代表当前网格是上级网格区域第一个遍历的网格，“当前网格与后一网格的位置关系”的“-”代表当前网格是上级网格区域最后一个遍历的网格。

表1：Hilbert曲线的开口方向调整对照表

第二步：令$t_x=\frac{X_{\max }-X_{\min }}{100}=3.4493,t_y=\frac{Y_{\max }-Y_{\min }}{100}=3.6208.$把(X_min-t_x,Y_min-t_y)、(X_max+t_x,Y_min-t_y)、(X_max+t_x,Y_max+t_y)、(X_min-t_x,Y_max+t_y)作为初始矩形凸壳的四个顶点。连接(X_max+t_x,Y_min-t_y)与(X_min-t_x,Y_max+t_y)两个顶点，将初始矩形凸壳分为两个三角形，作为初始三角网，并为初始三角网添加邻接关系。

第三步：按照先插入控制点数组，插入顺序是原存储顺序的逆序、后插入非控制点数组，插入顺序与原存储顺序相同，逐点进行如下操作：

(1)点定位

a.找到最新更新的三角形，即对上一个插入点进行处理时三角网中最后一个新生成的三角形。设该三角形顶点为A、B、C，沿逆时针排列。

b.设当前插入点坐标为P，P与三角形三边AB、BC、CA构成的三角形面积分别为$S_A=\frac{1}{2}\overrightarrow{\mathrm{AB}}\times \overrightarrow{\mathrm{AP}},S_B=\frac{1}{2}\overrightarrow{\mathrm{BC}}\times \overrightarrow{\mathrm{BP}},S_C=\frac{1}{2}\overrightarrow{\mathrm{CA}}\times \overrightarrow{\mathrm{CP}},$计算S_A，S_B，S_C，判断其中小于0的个数。若S_A，S_B，S_C均大于等于0，则ΔABC即为当前待插点P所在三角形，点定位过程结束；当S_A，S_B，S_C有一个小于0时，进行步骤c，当S_A，S_B，S_C有两个小于0时，进行步骤d；

c.找到相应小于0的S，此时P应位于该小于0的S所对应的边的外侧。把该边的临边三角形作为ΔABC，返回步骤b继续判断。

d.找到相应小于0的两个S中最小的S，此时把最小的S所对应的边的临边三角形作为ΔABC，返回步骤b继续判断。

例如附图8是由要插入P点时该点云数据形成的局部三角网，T1为第一个判断的三角形。首先求P与T1的三边形成面积的大小，得到下一个判断的三角形为T3，则继续判断P与T3三边形成的三角形面积大小，得到下一个判断的三角形为T6,。以此类推，直到找到P点与T13三边形成的三角形面积均大于等于0，得到P位于T13中，P的点定位过程结束。

(2)空腔建立

构建包含当前插入点所有影响域的空腔。

首先把当前插入点所在三角形的三条边作为空腔边，分别对空腔边的邻接三角形做外接圆检测。若当前插入点位于某三角形外接圆内部，则将该三角形对应的空腔边删除，将该三角形的另两边加入到空腔边中，继续对空腔边进行判断，直到所有空腔边的邻接三角形的外接圆都不包含当前插入点。

(3)三角网更新

将当前插入点与得到的空腔节点相连，形成新的三角形，并更新由空腔形成的新三角形及其相邻三角形的邻接关系。

第四步：全部点集插入后，将顶点中含有矩形凸壳四个顶点的三角形删除，即完成了Delaunay三角网的构建。

下面对本发明的构网效率进行检测。

以上述实例点集抽稀100倍、10倍和原始点集为测试数据集，采用目前已有的、较常用的串行Delaunay构网算法，对点排序时间、点插入时间和运行总时间进行对比，测试环境为便携式笔记本电脑，Intel(R)Core(TM)i7-4700HQ2.40GHz，8GB内存，win8，64位操作系统。检测结果见下表：

表2：不同构网方法时间比较

其中，CGAL为计算几何算法库，它提供了一种目前公认的比较快的Delaunay构网方法；RG是规则网格划分算法，MG是多重网格划分算法，HG是本发明所述算法。从对比结果可以看出，本发明所述算法效率比其他三种算法要高。本发明使用方便，在不依赖并行硬件结构的基础上，提高了大规模点云数据的Delaunay构网效率，经实验测试，速度最高可以达到每秒处理约90万个点，明显高于现有的串行Delaunay构网算法，提高了大规模点云数据的后处理效率。此外，经验证，本发明针对线型分布、圆型分布、交叉型分布、螺旋型分布等大规模、不均匀分布的点集，在Delaunay构网效率上也具有明显的优势，能够满足任意分布的大规模点云数据的高效Delaunay构网要求。

以上所述仅为本发明的实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。