三维重建技术是根据物体的二维图像信息恢复其三维模型, 多年来一直是计算机图形学领域的研究热点^[1-2].三角剖分是三维重建过程中最重要的基础之一, 通过将散点连接形成许多互不相交的三角形或四面体, 从而使模型面化或者体化^[3].经过多年探索, 研究者们对于二维三角剖分的研究已经取得了很多成果, 提出了多种三角网格的优化准则和构造方法.其中, Delaunay三角剖分由于具有良好的几何特性-可以最大化最小角^[4], 即能使得到的三角网格最平均, 使用最广泛^[5].然而, 随着计算机图形学技术的发展, 扫描设备的精度越来越高, 用于三维重建的点集对象的规模也越来越大.因此, 在保证三角网格的整体质量的提前下, 有必要寻找一种更高效的Delaunay三角剖分, 从而满足三维重建的实时性.

Delaunay三角剖分算法有很多种类, 主要分为三大类:逐点插入算法^[6]、分治算法^[7]、三角网生长算法.其中, 逐点插入法由于实现相对简单且占用空间相对较小^[8], 被大量研究者所采用; 分治算法在执行时间方面可以达到最优, 但占用内存较大, 不适用于一般的计算机平台; 三角网生长算法相对于前两种算法效率较低, 不适用于大规模点集, 在实际应用中很少被采用.逐点插入法的主要流程如下: 1)构造一个包含所有插入点的超级三角形; 2)将存储在链表中的点按序插入, 遍历三角形链表找出包含插入点的三角形(即目标三角形), 利用插入点将目标三角形拆分生成新的三角形, 在三角形链表中删除目标三角形, 完成一个点的插入; 3)利用空圆特性对新生成的三角形进行Delaunay规则判断; 4)循环执行2)和3), 直至完成对所有插入点的处理.经典的逐点插入法有Bowyer算法^[9]、Watson算法^[10]以及Lawson算法^[11]等.

虽然逐点插入法易于实现, 但在处理大规模的散点集数据时, 实时性不是很好, 为了提高逐点插入法的效率, 研究者们进行了大量的研究.一部分研究者认为可以对插入点进行预处理, 通过改进插入点的序列来提高三角剖分的效率. Liu等^[12]提出了一种基于广度优先搜索的确定性插入序列, 使用k-d树来构造Delaunay三角关系, 但是构造k-d较为复杂, 需耗费大量时间.在多重网格插入法^[13]的基础上, Su等^[14]提出了使用Hilbert曲线遍历插入点, 减少了三角剖分过程中狭长三角形出现的数量, 从而降低局部优化时间. Zalik等^[15]提出一种两级均匀细分加速技术, 将目标三角形问题转化为最近点问题, 在他们的方案中, 从最近的Delaunay点开始, 然后从第二个最近的Delaunay点开始, 以这种递归的方式直到找到目标三角形.在此基础上, Zadravec等^[16]提出结合哈希表与跳跃表来寻找插入点的最近点, 从而快速定位目标三角形.一些研究人员提出基于重心方向来定位目标三角形, 但是, 当出现分界点时, 这种方法的搜索路径可能不唯一.针对此问题, Xi等^[17]提出可以通过移动重心来避免截点, 从而使目标三角形可以连续搜索.随着计算机技术的快速发展, 许多研究者们提出了并行Delaunay三角剖分, 将点集划分成许多独立的分区, 这些分区可以同时进行Delaunay三角剖分, Kohout等^[18]、Rong等^[19]、Cuong等^[20]、Lo ^[21]均对此进行了研究及改进.此外, 杨昊禹等^[22]提出了并行动态Delaunay三角剖分算法, 可以解决新增点位于原来三角网之外的情况.另外李国庆等^[23]提出了基于凸多边形的Delaunay三角剖分算法, 使用生成的凸多边形代替传统算法中的超级三角形.常见的生成凸包算法有Graham扫描法^[24]、分区算法^[25]等.刘斌等^[26]提出将主成分分析法(PCA)与二分法结合, 通过快速确定凸包边缘点来计算平面点集的凸包.

通过对Delaunay三角剖分逐点插入法的研究, 可知在构建三角网过程中最耗时的部分为寻找插入点的目标三角形.本文在传统算法的基础上, 将边指针与区域划分相结合, 很大程度上提高了构建Delaunay三角网的效率.

1.   基于边指针的搜索

在传统的逐点插入法中, 对每个插入点的目标三角形的搜索总是从三角形链表的表头开始.这意味着, 目标三角形若位于链表的表头, 则可以很快被找到; 若位于链表的尾部, 则需要遍历整个链表才能被找到.随着点的数量大幅度增加, 三角形链表的长度也必然会大幅度增加, 从而导致寻找目标三角形会越来越耗时, 降低三角剖分的效率.针对此问题, 本文对三角形的数据结构进行优化, 提出了基于边指针的数据结构.如图 1所示, 在$ \triangle ABC $中, $ A $为起始顶点且$ A $、$ B $、$ C $逆时针分布, $ AC $为边1, $ AB $为边2, $ BC $为边3.在三条边上分别存在边指针$ P1 $、$ P2 $、$ P3 $, 分别指向边指针所在的公共边的相邻三角形, 如$ P1 $指向与边$ AC $相邻的三角形, 对于$ P2 $、$ P3 $情况类似.

图 1  三角形数据结构

Fig. 1  Data structure of triangle

下载: 全尺寸图片 幻灯片

1.1   点与三角形的位置判断

点与三角形之间的位置关系可以根据其坐标的相对关系来判断.对于任意三个点$ A(x_A, y_A ) $, $ B(x_B, y_B ) $, $ C(x_C, y_C ) $, 其中

若$ Flag(A, B, C)>0 $, 则$ A $、$ B $、$ C $按逆时针方向分布; 若$ Flag(A, B, C)<0 $, 则$ A $、$ B $、$ C $按顺时针方向分布; 若$ Flag(A, B, C) = 0 $, 则$ A $、$ B $、$ C $三点共线.如图 2 (a)所示, 若$ Flag(C, A, M) $、$ Flag(A, B, M) $、$ Flag(B, C, M) $三者的值都大于零, 则点$ M $位于$ \triangle ABC $的内部.若三者中一个值为零且其余二者皆为正值, 则位于某条边上, 如图 2 (b)中$ Flag(A, B, M) $为零, $ Flag(C, A, M) $、$ Flag(B, C, M) $为正值, 故点$ M $位于边$ AB $上.若三者中至少有一个为负值, 则点位于三角形的外部.如图 2 (c)所示, $ Flag(A, B, M) $为负值, 故点$ M $位于边$ AC $的外侧.当然, 点位于三角形外侧还会出现图 2 (d)与图 2 (e)的情况.在图 2 (d)中, $ Flag(A, B, M) $、$ Flag(B, C, M) $皆为负值, 在此情况下, 由于$ AB $为边2, 而$ BC $为边3, 本文中优先取编号在前的边, 故认为$ M $位于边$ AB $的外侧.在图 2 (e)中, 点$ M $虽位于三角形外部且与边$ AC $共线, 但是$ Flag(B, C, M) $为负值, 故本文认为点$ M $位于边$ BC $的外侧.

图 2  点与三角形位置关系判断

Fig. 2  Judgment of positional relationship between point and triangle

下载: 全尺寸图片 幻灯片

1.2   基于边指针搜索的路径确定

在搜索目标三角形的过程中, 搜索路径通过判断插入点与当前被搜索三角形的位置关系来确定的.如图 3所示, 点$ V $为插入点, 假设$ \triangle AIJ $存储在三角形链表的表头, 搜索过程如下所述. 1)检测点$ V $与$ \triangle AIJ $的位置关系, 首先确定了点$ V $位于$ \triangle AIJ $的外部, 然后确定点$ V $位于$ \triangle AIJ $的哪条边的外侧, 根据式(1)可知点$ V $位于边$ IJ $外侧, 故下一个被搜索的三角形为$ \triangle AIJ $的边$ IJ $上的边指针$ p_1 $所指的$ \triangle IEJ $. 2)用同样的方法确定接下来的被搜索三角形分别是$ \triangle EFJ $、$ \triangle EHF $, 然后根据式(1)可知点$ V $位于$ \triangle EHF $的内部, 故$ \triangle EHF $为点$ V $的目标三角形.在整个搜索过程中, 共经过4次搜索找到点$ V $的目标三角形, 搜索路径为: $ \triangle AIJ(p_1) $→$ \triangle IEJ(p_2) $→$ \triangle EFJ(p_3) $→$ \triangle EHF $.边指针搜索的方式, 由于在搜索目标三角形的过程中通过边指针的方向性可以快速地定位目标三角形, 使搜索路径大大缩短, 从而大大降低了三角剖分的耗时.

图 3  点$V$的目标三角形的搜索路径

Fig. 3  Search path for target triangle of point $V$

下载: 全尺寸图片 幻灯片

2.   基于边指针搜索及区域划分的三角剖分

2.1   区域分块

基于边指针搜索目标三角形的方法, 虽然可以利用插入点与三角形的位置关系来快速定位目标三角形, 但是在点数量很大的情况下, 目标三角形的搜索路径可能仍然会比较长, 此时依旧较为耗时.因此, 在边指针搜索的基础上, 提出了区域划分的方案, 以缩小目标三角形的搜索范围.

设区域分块为$ N\times N $, 构造超级三角形的时候, 首先遍历所有散点找到最大的纵横坐标值$ X_{\max} $、$ Y_{\max} $, 取两者中较大者得到$ max $, 然后根据$ max $构造超级三角形, 其各点坐标为:

区域划分是指将包含超级三角形的正方形划分成$ N\times N $个大小相同的区域, 并将这些区域存储在二维数组中.每个区域的位置坐标值$ X $和$ Y $对应于其在二维数组中的存储位置, 如位置坐标(1, 2)的区域块存储在数组area[1][2]中.对于插入点$ P(x, y) $, 它的区域坐标由其本身的横纵坐标值确定:

在式(3)中, $ L $指包含超级三角形的正方形的边长.

区域可以划分成不同的大小, 对应的区域数量也会不同.对于某个确定的点集, 如果区域太大, 搜索范围会变大.但是, 如果区域太小, 则过多的区域很难被分配到入口三角形, 从而导致入口三角形的使用率下降, 进一步导致搜索目标三角形的时间增加.所以对于处理不同规模的点集, 最优的分块也会不同.

对于点集$ V = \{V_1, V_2, \cdots, V_n\} $, 两个点组成一条边, 三条边组成一个三角形, 三角形通过公共边上的边指针相互连接组成三角形网格.基于边指针搜索及区域划分的Delaunay三角剖分的数据结构如表 1所示.

表 1  基于边指针及区域划分的算法数据结构

Table 1  Data structure of algorithm based on edge-pointer and region-division

Structure
POINT	$ x $	$ y $	entry_flag
EDGE	$ V_1 $	$ V_2 $
TRIANGLE	$ V_1 $, $ V_2 $, $ V_3 $	*$ p_1 $, *$ p_2 $, *$ p_3 $	entry	$ X, Y $
REGION	valid	*entrytri

下载: 导出CSV 

| 显示表格

数据结构包括4部分:点、边、三角形和区域.每个顶点包含$ x $和$ y $坐标以及区域入口标志$ entry\_flag; $连接顶点的边包含顶点$ V_i $和$ V_j $的编号; 三角形采用空间网状结构, 包含顶点$ V_i $、$ V_j $、$ V_k $的编号, 三角形的区域入口标志$ entry $以及该区域的坐标$ X $、$ Y $, 将位置相邻的三角形连接起来的边指针$ p_1 $、$ p_2 $、$ p_3 $, 区域坐标只有在$ entry = 1 $时才意义, 如果$ entry = 1 $, 意味着三角形是该区域的入口三角形; 区域块包含区域是否分配入口三角形的标志valid, 以及区域的入口三角形指针$ entry\_tri $.对于一个入口三角形, 在其三个顶点中存在一个入口顶点, 其entry_flag值被置为1, 并且入口顶点位于于该三角形所记录的区域中.边指针搜索及区域划分数据结构中不同区域之间的关系如图 4所示, 不同局域会分配一个入口三角形作为该区域的搜索入口, 该入口三角形再通过边指针指向存在该区域内或其它区域内的相邻三角形.

图 4  结构体关系图

Fig. 4  Structure diagram

下载: 全尺寸图片 幻灯片

2.2   搜索目标三角形

实现Delaunay三角剖分最重要也最耗时的部分就是搜索目标三角形, 本文对目标三角形的搜索过程如算法1所示.在基于边指针搜索的基础上, 区域划分的搜索方式为每个区域分配了入口三角形.这使得对散点的目标三角形的搜索不是从三角形链表的表头开始, 而是从入口三角形开始, 故搜索的起始三角形更接近目标三角形, 目标三角形的搜索路径进一步被简化.

算法1.搜索插入点$ V_i $的目标三角形

S1:确定$ V_i $所在区域;

S2: if该区域已经分配入口三角形, then

S3:   从入口三角形开始搜索;

  else

S4:   从最新生成的三角形开始搜索;

S5: if $ V_i $位于当前搜索三角形内部, then

S6:   找到并返回$ V_i $的目标三角形;

S7: else if $ V_i $位于当前三角形的边上, then

S8:   根据边指针搜索共边三角形并返回这两个三角形;

  else

S9:根据边指针确定下一个要搜索的三角形, 转到S5.

图 5展示了基于边指针搜索及区域划分算法中搜索目标三角形的过程.如图 5 (a)所示, 包含超级三角形的正方形被划分成2 $ \times $ 2个区域.插入点$ V $时, 首先根据点$ V $的坐标确定其所在区域为area[1][0].由于在点$ V $之前, 点$ H $与点$ K $已经被插入该区域, 故该区域一定存在入口三角形, 假设$ \triangle GCH $为入口三角形, 对点$ V $的目标三角形的搜索从$ \triangle GCH $开始.接下来, 按照第1.2节中所述边指针搜索的方式, 确定目标三角形的搜索路径为: $ \triangle GCH(p_1 $)→$ \triangle DGH(p_2 $)→$ \triangle EDH(p_3 $)→$ \triangle EHF $, 目标三角形为$ \triangle EHF $.

图 5  搜索点$V$的目标三角形

Fig. 5  Searching for target triangle of point $V$

下载: 全尺寸图片 幻灯片

如果检测到点$ V $位于三角形的某条边上, 则需要根据边指针搜索另一个包含该点的一个三角形.如图 5 (b)所示, 点$ V $位于边$ HF $上, 当找到$ \triangle EHF $之后, 根据边$ HF $上的边指针$ p_4 $找到共边三角形$ \triangle FHK $.故$ \triangle EHF $与$ \triangle FHK $为点$ V $的目标三角形, 并且点V将会被插入这两个三角形所构成的四边形内.

如果在点$ V $所在区域中没有入口三角形, 则从最新插入的点生成的三角形开始搜索目标三角形, 搜索路径由边指针确定.由于该算法主要用于三维重建, 而三维重建中扫描得到的连续点在位置上可能是相邻的, 因此, 由最新插入点生成的三角形可能更接近当前插入点的目标三角形.

在传统算法中, 假设当前所有三角形按一下顺序存储在三角形链表中: $ \triangle AIJ $、$ \triangle AJF $、$ \triangle IEJ $、 $ \triangle EFJ $、$ \triangle IBL $、$ \triangle ILD $、$ \triangle IDE $、$ \triangle EDH $、 $ \triangle EHF $、$ \triangle FHK $、$ \triangle CFK $、$ \triangle BGL $、$ \triangle DLG $、$ \triangle DGH $、$ \triangle GCH $.插入散点$ V $时, 则需要依次判断链表中在$ \triangle EHF $之前的三角形是否包含点$ V $, 可以发现一共需要经过9次搜索才能找到目标三角形, 而本文的算法只需要经过4次搜索就能找到目标三角形, 搜索深度为传统算法的44.4 %.随着点集的规模不断变大, 传统算法的执行时间将会呈指数增长, 本文算法的优越性会越来越明显.

目标三角形的定位是根据插入点与当前搜索三角形的位置关系来确定的, 通过判断插入点位于三角形的哪条边的外侧来确定下一个被搜索的三角形.每一次被搜索的三角形都比上一次的三角形更接近插入点, 因此通过这种方式一定能找到目标三角形, 故该算法是收敛的.

2.3   区域入口三角形及边指针的更新

随着点地不断插入, 区域的入口三角形也需要不断更新以确保其不会随着某些三角形的拆分而消失.同时, 与新三角形有关的边指针也需要被更新以确保三角形之间能够互相找到相邻的三角形.找到目标三角形之后, 需要将点插入到目标三角形中, 并使用空圆特性对新生成的三角形进行Delaunay规则判断及调整违背规则的三角形.在此过程中, 需要更新入口三角形及边指针, 入口三角形的更新分为本区域的更新以及他区域的更新.

如图 6所示, 找到目标三角形$ \triangle EHF $之后, 插入点$ V $, 生成了三个新的三角形($ \triangle EHV $, $ \triangle HFV $与$ \triangle FEV $), 并将原目标三角形$ \triangle EHF $从三角形链表中删掉.首先, 对本区域的入口三角形进行更新.如果在此之前, 区域area[1][0]没有分配到入口三角形, 则选择$ \triangle EHV $为该区域的入口; 如果area[1][0]之前已经分配了入口三角形如第2.2节中所提到的$ \triangle GCH $, 则$ \triangle EHV $代替$ \triangle GCH $作为该区域的新入口, 并将点$ V $指定为入口顶点.然后, 根据被删除三角形的标志位及区域坐标信息检测其是否为其他区域的入口, 若是, 则需要指定一个新的三角形来更新此区域的入口.例如, 如果三角形$ \triangle EHF $是区域area[1][1]的入口, 判断$ \triangle EHF $的哪个顶点是入口顶点, 根据点的入口标志位得出点$ F $是入口顶点, 因此, 包含点$ F $的$ \triangle FEV $成为了area[1][1]的新入口.

图 6  入口三角形及边指针的更新

Fig. 6  Update of entry triangle and edge-pointer

下载: 全尺寸图片 幻灯片

此外, 还需要更新和原三角形$ \triangle EHF $有关的边指针.在图 6中, 随着$ \triangle EHF $的删除, 原本指向该三角形的边指针应该指向新的三角形, $ \triangle EDH $的边$ EH $上的边指针$ p_3 $现在应该指向$ \triangle EHV $, 同理, $ \triangle EFJ $的边$ EF $上的边指针$ p_5 $、$ \triangle FHK $的边$ FH $上的边指针$ p_6 $现在分别指向其相邻三角形$ \triangle FEV $、$ \triangle HFV $.

在进行Delaunay规则判断的过程中, 区域入口及边指针会以同样的方式更新, 此处不再赘述.

3.   实验结果及分析

为了验证基于边指针搜索及区域划分的三角剖分算法的效率, 本文对随机生成的不同规模大小的点集进行三角剖分, 并比较了传统的Delaunay三角剖分(TD)、使用CGAL库的三角剖分、基于边指针搜索的三角剖分(EPD)、基于边指针搜索和区域划分相结合的三角剖分(EPRDD)几种算法之间的执行时间, 此处的执行时间不包括读取插入点数据的时间.

在进行区域划分的时候, 包含超级三角形的正方形被划分成了15 $ \times $ 15个区域.本实验的实验环境为Intel(R) Core(TM) CPU i5-8500 @ 3.00 GHz, 16 GB内存, Windows 10, VS2013, 64位操作系统, 得到的结果如表 2所示:

表 2  算法的执行时间(s)

Table 2  Running time of algorithms (s)

散点数量	2万	3.5万	5万	7万	10万
TD	113.515	298.687	756.656	1 943.39	3 972.81
CGAL	1.719	3.074	4.281	6.265	8.594
EPD	0.563	1.547	2.562	3.218	6.172
EPRDD	0.117	0.234	0.325	0.531	0.813

下载: 导出CSV 

| 显示表格

从表 2可以看出, 当插入点的数量为2万时, 基于边指针搜索的三角剖分耗时为传统算法的1/201, 加入区域划分之后, 算法的执行时间为传统算法1/970.当插入点的数量为10万时, 基于边指针的搜索算法执行时间分别为传统算法的1/643, 而基于边指针及区域划分算法的执行时间为传统算法的1/4 886.此外, 与CGAL算法相比, 只使用边指针的算法优势并不是很明显; 但是, 边指针与区域划分相结合之后, 算法的耗时明显低于CGAL算法.当对10万个点进行三角剖分的时候, CGAL算法需要8.594 s, 而基于边指针及区域划分的算法仅仅需要0.813 s就能完成.

为了比较几种算法的执行时间与散点数量的关系, 绘制了图 7所示折线图.由于传统算法的执行时间远远高于其余几种算法, 故不与其他算法作比较.从图 7可以看出, 随着散点数量的增多, CGAL算法, 基于边指针搜索的算法与基于边指针搜索及区域划分算法的执行时间基本上都是呈线性增长的, CGAL算法与基于边指针搜索的算法的执行时间随散点数量变化的增长速度还是较快, 但是基于边指针及区域划分算法的执行时间的增长速度较为缓慢, 斜率远远小于前两种算法.可以得知, 随着插入点数量的增长, 边指针搜索与区域划分结合的方法在时间上的优势将会越来越大.

图 7  执行时间比较

Fig. 7  Comparison of run time

下载: 全尺寸图片 幻灯片

进一步对几种算法的平均搜索深度进行比较, 结果如表 3所示.搜索深度是指在将散点插入到Delaunay网格中的过程中搜索三角形的次数, 搜索深度越小, 三角剖分的实时性越好.由于CGAL的源代码没有找到, 在此处未对CGAL的搜索深度进行比较.从表中可以看出, 在传统算法中平均搜索深度与散点数目也是线性增长的, 而基于边指针搜索及区域划分的算法的平均搜索深度很小.即使散点数量成倍增加, 平均搜索深度也只是稍微增加, 证明边指针与区域划分结合的算法实时性要高于传统Delaunay三角剖分算法.其原因在于, 区域划分缩小了目标三角形的搜索范围, 边指针优化了搜索方向, 最终使得搜索深度大大减小, 执行时间得到降低.

表 3  算法的平均搜索深度

Table 3  Average search depth of algorithms

散点数量	2万	3.5万	5万	7万	10万
TD	10 003	17 543	24 989	34 893	49 854
EPD	108	229	233	247	278
EPRDD	8	11	13	15	18

下载: 导出CSV 

| 显示表格

4.   结论

为了满足三维重建的实时性, 针对其基础部分三角剖分处理时间过长的问题, 本文提出了一种基于边指针搜索及区域划分的Delaunay三角剖分算法.在对传统Delaunay三角剖分算法进行研究及实现的基础上, 进一步优化了三角形存储的数据结构及搜索过程, 设计了一种通过边指针反映三角形空间相邻关系的三角形数据结构来取代传统算法中的一维三角形链表.由于边指针具有方向性, 可以用于确定距离插入点最近的相邻三角形, 从而快速确定目标三角形的搜索方向.在此基础上将整个超级三角形进行区域划分, 且随着散点的插入, 每个区域将会分配一个区域入口三角形, 作为搜索目标三角形的起始三角形, 且入口三角形在构建Delaunay三角网的过程中会不断更新.这使得目标三角形的搜索范围由以前的整个三角形链表缩小到插入点所在区域, 从而进一步缩短搜索的起始三角形与目标三角形之间的距离, 大大缩小了搜索范围, 降低了搜索目标三角形的耗时, 以满足三维重建的实时性需求.实验结果表明, 该算法的执行时间随散点数量的增长较为缓慢, 且散点数量为10万时, 所耗费的执行时间仅为0.813 s.

本文中的区域划分为固定划分, 当散点分布密度变化大时, 对于密度大的区域搜索深度可能会相应增加.为了进一步扩大算法的适用范围, 后续研究中将引入自适应算法, 寻求一种能够根据散点分布密度或者搜索深度的变化对区域划分方式进行调整的优化算法.