谢增广

（华北计算技术研究所，北京100083）

0 引 言

三角剖分是计算几何中领域的重要课题之一。平面点集Delaunay三角剖分与该点集的Voronoi图是对偶图，具有许多优良性质，在图形网格化技术领域有着广泛应用。三角剖分算法理论研究已经相当成熟，已经能够证明算法时间复杂度的上界和下界［1］，然而由于区域拓扑划分问题，算法在实际工程中应用的难易程度有所不同。根据实现过程，Delaunay三角剖分的算法可以分为逐点插入法、三角网生长法，分治算法等，其中分治算法最适合实际工程应用。

论文引入了一种数据结构——双链接边表 （doublyconnected edge list，DCEL），利用该数据结构可以方便地对平面区域进行划分；引入了3个工具算法，简化了Delaunay三角剖分算法分治的实现过程；另外对特殊退化情况进行分类处理，增强了该算法的实用性。

1 算法设计与实现

1.1 数据结构

双链接 边 表 （doubly－connected edge，DCEL）［4－5］是 用来表示平面图拓扑信息的一种数据结构，如图1所示。

在DCEL中，将每条边的两端分别作为一条半边。半边有方向性，由起点出发，指向另一个端点。一条边的两条半边互为孪生半边。

DCEL由三组记录构成，分别存储顶点、面、半边的信息。

（1）在顶点V的记录中，除了存储点的坐标信息，还有一个指向半边的指针，指向以V为起点的某一条半边。

（2）在面f的记录中，只存储了一个指向半边的指针，指向该面边界上的某一条半边。

图1 DCEL数据结构

（3）在半边edge的记录中，存储了半边的起点，一个半边的指针指向其孪生半边，一个面指针指向位于半边左侧的面face。此外还有2个半边指针，分别指向沿着face边界方向的前一条半边和后一条半边。

根据DCEL提供的信息，可以很方便地对平面进行子区域划分。为了简化分治Delauany三角剖分算法实现过程，不必存储面的信息。

1.2 预备工具算法

介绍分治Delauany三角剖分算法之前，先引入3个预备工具算法：判断点是否在指定三点外接圆内，判断夹角不大于π，求不相交的2个凸多边形的上下公切线。

1.2.1 判断点是否在指定三角形的外接圆内

设不共线三点a，b，c构成的三角形记作△ （a，b，c），点的d｛a，b，c｝。判断d是否在△ （a，b，c）外接圆内的方法是［6］：记ret＝，ret＜0，d在△ （a，b，c）外接圆外部；ret＝0，d在△ （a，b，c）外接圆边界上；ret＞0，d在△ （a，b，c）外接圆内部。

定义操作inCircle（a，b，c，d），如果d在a，b，c确定的三角形外接圆内，返回TRUE；否则，返回FALSE。显然地，inCircle（a，b，c，d）的时间复杂度为O （1）。

1.2.2 判断夹角不大于π

设不共线三点为a，b，c，记逆时针方向夹角∠abc为θ，θ与π大小关系判断方法是［6］：记ret＝，假定人沿着方向行走，如果c在其左侧，ret＞0，θ＜π；如果c在其右侧，ret＜0，θ＞π；如果c在方向上，ret＝0，θ＝π。

1.2.3 求不相交的2个凸多边形的上下公切线

设左边的凸多边形为LP，右边的凸多边形为RP，如图2所示。记L为LP边界上的最右边的顶点，R为RP边界上的最左边的顶点；记LP和RP的上公切线是tct（top common tangent），TL在LP的边界上，为tct的左端点，TR在RP的边界上，为tct的右端点；记LP和RP的下公切线是bct（bottom common tangent），BL在LP的边界上，为bct的左端点，BR在RP的边界上，为bct的右端点。为上公切线，当且仅当LP边界上的点和RP边界上的点均不在的左侧；同理，为下公切线，当且仅当LP边界上的点和RP边界上的点均在的左侧。

图2 求2个不相交凸多边形的上下公切线

为了更好地说明该算法，设点V为凸多边形P边界上的一点，定义下面2个操作：①ccw （V，P），V沿着P边界逆时针方向移动到达的第一个点；②cw （V，P），V沿着P边界顺时针方向移动到达的第一个点。

zig－zag算法

输入：左侧凸多边形LP，右侧凸多边形RP。LP和RP不相交。

输出：LP和RP的上下公切线。

Begin

步骤1 L为LP边界上的最左点，R为RP边界上的最右点。L1＝L，R1＝R，L2＝L，R2＝R。

步骤2 若toLeft（L1，R1，cw （R1，RP））为真，R1＝cw （R1，RP）；若toLeft（L1，R1，ccw （L1，LP））为真，L1＝ccw （L1，LP）。迭代这一过程，直至LP边界上的点和RP边界上的点均不在的左侧。

步骤3 若toLeft（L2，R2，ccw （R2，RP））为假，R2＝ccw （R2，RP）；若toLeft （L2，R2，cw （L2，LP））为假，L2＝cw （L2，LP）；迭代这一过程，直至LP边界上的点和RP边界上的点均在的左侧。

End

zig－zag算法执行的次数只与LP的右半边界上的点和RP的左半边界上的点规模有关。设点集规模为N，zigzags算法时间复杂度不超过O （N）。

定义操作zig－zag（LP，RP），得到LP和RP的上公切线为tct，下公切线为bct。

1.3 分治Delauany三角剖分算法

分治Delauany三角剖分算法的基本步骤如下：

（1）点集初始化：将点集以横坐标为主，纵坐标为辅按升序排序。

（2）将点集划分为近似相等的两个子点集。

（3）分别对两个子点集进行Delauany三角剖分。

（4）合并两个子点集的Delauany三角网。

（5）递归执行步骤 （2）到步骤 （4），直至点集的所有点都参加了Delauany三角网的构造。

由排序算法时间复杂度可知，点集初始化算法时间复杂度为O （N＊logN）。点集初始化之后，将点集合划分，直至每个子点集合的规模不超过3个；对这些子点集三角化，得到点、线段或者三角形，如图3所示。

图3 点集初始化和点集合划分

点集初始化后，设点集存储为数组P［］。定义以下2个操作：

（1）Merge （LDT，RDT）：LDT 为 左 侧 子 点 集Delauany三角剖分后得到的结果，RDT为右侧子点集Delauany三角剖分后得到的结果，Merge是合并操作。该操作返回合并LDT和RDT之后得到的Delauany三角剖分后的结果DT。合并是分治Delauany三角剖分算法的核心，会在下一节合详细阐述。

（2）Delauany＿Triangulate（start，end）：start为点集的起始点编号，end为点集的终结点编号，Delauany＿Triangulate是分治操作。该操作将点集划分为若干子点集，使得每个子点集的数量不超过3个，并对子点集进行三角化；对相邻的点集构成的简单凸多边形迭代地进行Merge操作合并，直至得到点集的Delauany三角剖分的结果。该操作返回点集的Delauany三角剖分的结果，设为DT。

分治算法伪代码如下：

Delauany＿Triangulate（start，end）

｛

if（1＝＝end－start＋1）//只有1个点

｛

将Pstart放入DT，return DT；

｝

else if（2＝＝end－start＋1）//只有2个点

｛

｝

esle if（3＝＝end－start＋1）//只有3个点 （3点共线为退化情况，暂不考虑）

｛

｝

else //多于3个点，进行分治

｛

mid＝ （start＋end）/2；

//LDT表示左边点子集Delauany＿Triangulate后得到的结果

LDT＝ Delauany＿Triangulate（start，mid）；

//RDT表示右边点子集Delauany＿Triangulate后得到的结果

RDT＝ Delauany＿Triangulate（mid＋1，end）；

//合并LDT和RDT，得到结果DT

DT＝ Merge（LDT，RDT）；

｝

return DT；

｝

1.3.1 合并算法

设DT是合并LDT和RDT的结果。如图4所示，DT中的边分为3类：

（1）LL－edge：在Delauany三角剖分左子点集结束得到LDT时，LL－edge已经被创建，且LL－edge∈LDT；显然的，LL－edge的两个端点 ∈LDT。

（2）RR－edge：在Delauany三角剖分右子点集结束得到RDT时，RR－edge已经被创建，且RR－edge∈RDT；显然的，RR－edge的两个端点∈RDT。

（3）LR－edge：LL－edge是在合并LDT和 RDT过程中被新创建的，且LR－edge∈DT；在DT中，LR－edge是Delauany边；在创建LR－edge之前，LR－edge的左端点 ∈LDT，LR－edge的右端点∈RDT。

在合并过程中，由于创建LR－edge，可能会删除若干LL－edge和 RR－edge，且已经删除的 LL－edge或 RR－edge不会再生成；合并过程中，不会创建新的LL－edge和RR－edge。

对图3中划分的点集合进行合并，得到图4。由于划分的各个子点集三角化后，得到最简单的凸多边形 （点，线段，三角形），此次合并并没有删除LL－edge和RR－edge。

图4 合并后的边分为LL－edge，RR－edge，LR－Edge

为了详细阐述合并算法且不失一般性，下面对图4的左子点集Delauany三角剖分和右子点集Delauany三角剖分进行合并。

合并算法初始化时记DT＝LDT∪RDT。Delauany三角剖分得到的闭合边界就是点集构成的凸包，因此LDT和RDT的上下公切线在DT中。

记下公切线形成的LR－edge为Base LR－edge，上公切线形成的LR－edge为Roof LR－edge，如图5所示。

图5 合并初始化，Base LR－Edge和Roof LR－edge

合并算法的核心是从Base LR－edge向上搜索LR－edge，将其加入到DT中，直至到达Roof LR－edge为止。在这一过程 中，可能使 LL－edge或 RR－edge不 再 是 DT 中 的Delauany边，因此这样的LL－edge或RR－edge需要从DT中删除。

为了阐述合并算法，记DT为平面点集的Delauany三角剖分，定义以下5个操作：

根据DCEL数据结构的特性，以上5个操作可以在O（1）时间完成。

记New LR－edge的左端点为L，右端点为R。初始化时，L为下公切线的左端点，R为下公切线的右端点。New LR－edge从Base LR－edge开始向上寻找下一个LR－edge，记为Next LR－edge。Next LR－edge的一个端点是L或R，另一个端点记为V。因此根据LR－edge定义，只需要确定端点 V （VNew LR－edge），就可以确定Next LR－edge。

可以证明，V是L或R的邻点。所谓L的邻点，也称作LDT关于L的候选点，是指该点和L为端点的线段是LL－edge；同理，所谓R的邻点，也称作RDT关于R的候选点，是指该点和R为端点的线段是RR－edge。否则，如果V不是L或R的邻点，则New LR－edge会与LL－edge，或RR－edge，或其它的LR－edge相交，这不符合三角剖分的基本性质。

为了进一步确定V的搜索范围，V只能在2个点之间选择，一个是L的某个邻点，也称为LDT关于L的最终候选点 （final candidate）；另一个是R的某个邻点，也称为RDT关于R的最终候选点。

1.3.2 选取候选点

为了详细阐述寻找V的过程且不失一般性，先寻找RDT中关于R的最终候选点。记RDT中关于R的第一潜在候选点 （first potential candidate）R1，R1是R邻点Ri中与线段以R为轴点顺时针方向构成夹角∠LRRi最小的点；记RDT中关于R的次潜在候选点 （next potential candidate）R2，R2是R邻点Ri中与线段以R为轴点顺时针方向构成夹角∠LRRi次小的点。如图6所示。

第一潜在候选点是否是最终候选点，需要经过以下2条判定规则检验：

（1）向上搜索判定：第一潜在候选点R1与边顺时针方向夹角∠LRRi＜π。向上搜索判定保证寻找LR－edge的过程是从下公切线逐点向上，直至到达上公切线为止。

图6 RDT关于R的最终候选点的选取

（2）外接圆判定：由第一潜在候选点R1，New LR－edge的左右端点L，R三点确定的外接圆不包含次级潜在候选点R2。外接圆判定保证Delaunay剖分的正确性。

如果 （1）和 （2）都满足，则R1是RDT关于R的最终候选点；如果 （1）不满足，则搜索LR－edge过程已经达到了Roof LR－edge的右端点，RDT中关于R的候选点不必考虑，即没有可以选择的候选点；如果 （1）满足，而 （2）不满足，则从DT中删除线段所表示的 RR－edge，R1＝R2，即次潜在候选点作为第一候选点，迭代这一处理过程直至找到最终候选点或者没有可以选择的候选点结束，如图6所示。寻找RDT候选点的算法描述如下：

Begin

步骤1 R1＝cw＿rotate（R，L，RDT），若toLeft（L，R，R1）为真，执行步骤2；否则执行步骤3。

步骤2 R2＝cw＿rotate（R，R1，RDT），若inCircle（L，R，R1，R2）为真为RR－edge，需要从DT中删除，delete＿RR （R，R1，DT），R1＝R2，R2＝cw＿rotate（R，R1，RDT）；迭代这一过程，直至inCircle （L，R，R1，R2）为假。R1是RDT最终候选点。

步骤3 RDT中没有候选点

End

同理，镜像地可得到LDT中关于L的最终候选点，如图7所示。

图7 找到LDT关于L的最终候选点

寻找LDT候选点的算法描述如下：

Begin

步骤1 L1＝ccw＿rotate（L，R，LDT），若toLeft（L，R，L1）为真，执行步骤2；否则执行步骤3。

步骤2 L2＝ccw＿rotate（L，L1，LDT），若inCircle（L，R，L1，L2）为真，为LL－edge，需要从DT中删除，delete＿LL （L，L1，DT），L1＝L2，L2＝ccw＿rotate（L，L1，LDT）；迭代这一过程，直至inCircle（L，R，L1，L2）为假。L1是LDT最终候选点。

步骤3 LDT中没有候选点

End

1.3.3 确定LR－edge

根据LDT关于L的最终候选点和RDT中关于R的最终候选点的存在性，分以下4种情况确定LR－edge：

（1）如果RDT中关于R没有可以选择的候选点，且LDT中关于L没有可以选择的候选点，则线段表示的LR－edge为Roof LR－edge，合并过程结束。

（2）如果RDT中关于R的最终候选点是R1，且LDT中关于L没有可以选择的候选点，则V＝R1，边是LR－edge，需要加入到DT中。

（3）如果RDT中关于R没有可以选择的候选点，且LDT中关于L的最终候选点是L1，则V＝L1，边是LR－edge，需要加入到DT中。

（4）如果RDT中关于R的最终候选点是R1，且LDT中关于L的最终候选点是L1，则需要对L，R，L1，R1这4个点进行Delaunay三角形测试，如图8所示。如果L1在△ （L，R，R1）外接圆的外部，则 V＝L1，边是LR－edge，需要加入到DT中；否则 V＝R1，边是LR－edge，需要加入到DT中。

由此LR－edge向上寻找下一个LR－edge，迭代上述合并处理流程，在到达Roof LR－edge时合并结束，得到点集的Delauany三角剖分，如图9所示。

1.3.4 合并算法伪代码及算法复杂度分析

综上所述，合并算法伪代码如下：

Merge（LDT，RDT）

｛

DT＝LDT∪RDT；//初始化操作

zig－zag （LDT，RDT）得上公切线Roof＿LR－edge，下公切线Base＿LR－edge；

点L＝Base＿LR－edge的左端点，点R＝Base＿LR－edge的右端点；

insert＿LR （L，R，DT）；//将Base＿LR－edge加入到DT中

while（不是 Roof＿LR－edge）//向上搜索LR－edge，在到达Roof＿LR－edge时终止

｛

分别寻找LDT和RDT候选点L1，R1；

根据左右候选点L1，R1的存在性确定新的LR－edge并删除无效的LR－edge，搜索向上进一步；

｝//end while

Merge结束，得到得到点集的Delauany三角剖分，return DT；

｝

Merge算法的执行分为以下3种操作：

（1）Merge初始化。用zig－zag算法寻找上下公切线，算法复杂度不超过O （N）。

（2）删除LL－edge或RR－edge。在搜索LR－edge过程中删除的LL－edge或RR－edge的数量为常数。

（3）加入LR－edge。Merge生成 LR－edge的过程与zigzag类似。记点集U＝左Delauany三角剖分右边界上的点∪右Delauany三角剖分的左边界上的点。U的规模不超过N，而建立LR－edge的次数与点集U的规模是线性关系，所以加入LR－edge的次数不超过N。

综上所述，Merge算法的时间复杂度为O （N），分治的复杂度为O（1）。设分治Delauany三角剖分算法的总的时间复杂度为T （N），则T （N）＝2T （N/2）＋O （N）。根据主定理［9］，分治Delaunay三角剖分算法的时间复杂度为O （N＊LogN）。

1.3.5 退化情况处理

一般地，退化情况在边界值判断时候出现。

（1）对于toLeft（A，B，C）操作，A，B，C三点共线是退化情况。为此，增加一个判断共线的操作coLine（A，B，C）操作。在点集合划分子点集时，如果子点集的规模是3，且这三点共线，则此子点集三角化应该得到2条边。在此基础上，三点共线的退化情况不会影响算法的正确性。

（2）对于inCircle（A，B，C，D）操作，A，B，C，D四点共圆是退化情况。在合并过程中，四点共圆形的退化情况不会影响整个算法的正确性。

2 实验与结果分析

针对不同点集规模，每个数量级进行5次运算，算法运行的平均时间见表1，根据实验结果进行曲线拟合如图10所示。

表1 实验结果数据

图10 时间复杂度曲线拟合

顶部实线表示由公式0.0005＊ N＊LOG （N）得到的曲线，而虚线表示实际中测得的平均运行时间。从图中实际数据拟合曲线与标准曲线的接近程度，表明程序运行时间的变化规律符合O （N＊logN）这一趋势，与算法的复杂度为O （N＊logN）的理论结果相一致。

中间的实线是经典分治算法实验者测得的平均运行时间。因为采用了zig－zag算法，简化了求不想交凸多边形上下公切线的过程，并且DCEL数据结构特性使合并算法中的若干操作都是O （1），从而整个算法的效率较经典算法提升了约5%。

3 结束语

本文使用DCEL数据结构表示平面点集Delaunay三角剖分的结果；阐述了分治Delaunay三角剖分算法，并利用zig－zag算法寻找不相交凸包上下公切线简化了算法的实现。并且对退化情况进行分类处理，增强了算法的实用性。

在实际工程中，Delauany三角剖分的模型会更复杂，实现难度大，其中具有代表性的模型有平面点集含特征约束的Delauany三角剖分和三维空间的Delauany三角剖分，值得进一步深入研究。

由于DCEL数据结构不仅能有效地对平面点集进行子区域划分，而且可以存储三维空间的拓扑信息，本论文为复杂的问题模型研究与实际应用奠定了基础。

［1］Mark de Berg，Otfried Cheong，Marc van Kreveld，et al，Computational Geometry：Algorithm and Applications ［M］.3rd ed.DENG Junhui，transl.New York：Springer－Verlag Berlin Heidelberg，2008：197－218 （in Chinese）. Mark de Berg，Otfried Cheong，Marc van Kreveld，等.计算几何—算法与应用 ［M］.3版.邓俊辉，译.北京：清华大学出版社，2009：250－285.

［2］ZHOU Jia－wen，XUE Zhi－xin，WAN Shi.Survey of triangulation methods ［J］.Computer and Modernization，2010，26（7）：75－78 （in Chinese）.［周佳文，薛之昕，万施.三角剖分综述 ［J］.计算机与现代化，2010，26 （7）：75－78.］

［3］Doubly－connected edge list ［EB/OL］. ［2011－04－18］.http：//en.wikipedia.org/wiki/Doubly＿connected＿edge＿list.

［4］Ryan Holmes.The DCEL data structure for 3Dgraphics［EB/OL］.［2009－09－27］.http：//www.holmes3d.net/graphics/dcel/.

［5］Max McGuire.The half－edge data structure ［EB/OL］.［2000－08－07］.http：//www.flipcode.com/archives/The ＿ Half－Edge ＿Data＿Structure.shtml.

［6］Guibas L，Stolfi J.Primitives for the manipulation of general subdivisions and the computation of Voronoi diagrams ［J］.ACM Transactions on Graphics，1985，4 （2）：75－123.

［7］Rex A Dwyer.A faster divide－and conquer algorithm for constructing Delaunay triangulations［J］.Algorithmica，1987，2（2）：137－151.

［8］Lee D T，Schachter B J.Two algorithms for constructing a Delaunay triangulation ［J］.International Journal of Computer and Information Sciences，1980，9 （3）：219－242.

［9］Thomas H Cormen，Charles E Leiserson，Ronald L Rivest，et al.Introduction to algorithms［M］.PAN Jingui，GU Tiecheng，LI Chengfa，et al transl.2nd ed.MIT Press，2001：43－49 （in Chinese）.Thomas H Cormen，Charles E Leiserson，Ronald L Rivest，等.算法导论 ［M］.2版.潘金贵，顾铁成，李成法，等译.北京：机械工业出版社，2011：43－49.

［10］WANG Li.Research on triangulation algorithm ［D］.Shanghai：Shanghai Jiaotong University，2010 （in Chinese）. ［王力.三角剖分算法研究 ［D］.上海：上海交通大学，2010.］