2010年第7期福建电脑
申永源，曹布阳张叶帆
（同济大学软件学院上海201804）
【摘要】：为了加快大规模二维平面点集的泰森多边形生成速度，本文设计并实现了一种并行优化算法。该算法在保证与串行算法具有相同的精准度的条件下，利用串行算法的分治特征对其有效的进行了并行化优化。经过实验证实，该算法在并行计算的环境下有效地提高了计算速度，减少了执行时间，并且获得了较高的计算加速比。 【关键词】：泰森多边形,Delaunay三角形，并行计算，三角剖分
1、引言
在城镇化与村镇建设动态监测的过程中，常常需要借助于泰森多边形[1](Thiessen Polygon，又被称为Dirichlet图、Voronoi 图)进行地理信息分析。对于大规模的数据点，串行的泰森多边形的生成算法速度较慢，因而需要一种并行算法提高生成速度。
生成Delaunay三角形是Thiessen多边形生成的基础，二者可以通过一定规则互相转换[2]。而生成Delaunay三角形具有多种算法，依照参考文献[3]的研究，在串行计算的条件下，Dwyer算法[4]无论对于均匀分布或者非均匀分布点集，在输入点数小于一定数量(216)的情况下，都相对于其它的算法具有较佳的性能。并且其本身的分治特征也比较容易进行并行化，所以本文在此的基础上设计与改进并行化的泰森多边形生成算法。
2、算法流程
一个较为实用的泰森多边形生成程序的算法流程如下所示：
1.输入点集读取，并且去除点集中的重复点，计算MBR(最小外包矩形)
2.在远离点集MBR的四角加入四个额外点
幽灵楼道
3.Delaunay三角剖分
4.将Delaunay三角形转换为泰森多边形
5.将生成的泰森多边形裁切到合适的范围内
6.输出最终结果
钨铜电触头
米脂的婆姨性功能由于输入点集存储于ESRI Shape文件中，其采取了四叉树的存储形式，因而MBR可以直接获得。而去除重复点时如果采取数组作为查存储结构，则不得不多次线性扫描数组出重复点，数据量大时性能十分低下。因此可以建立一个哈希表，以X坐标作为键值，Y坐标的链表作为查存储结构。经过实验证明，对于2万个输入点的问题规模，该改进可以获得约30倍的性能提升。而对泰森多边形的裁切，使用的是Sutherland Hodg-man多边形裁切算法[5]。由于对于每一个泰森多边形的裁切都是独立的，因而可以使用并行计算进行优化。
但是经过实验可知，对于较大规模的问题(输入点数1万以上)，步骤3与步骤4的计算时间占了总时间的90%以上，因而对这两个步骤进行并行化，对于总计算性能的提高最为关键。3、Quad-Edge
考虑到Delaunay三角形与Thiessen三角形的快速转换，本算法使用Quad-Edge[6](四方边缘结构)作为存储数据结构。这是由于每个Quad-Edge单元都储存了两条主边(e0与e2)与两条对偶边(e1与e3)，并且按照逆时针顺序可以顺序访问。即对于每个顶点，包围其的对偶边就是其泰森多边形的边，并且可以直接转换为Shape文件格式所要求的顺时针顺序，无需再进行额外的排序工作。所以这种情况下，Delaunay三角形与泰森多边形的转换是非常方便与快速的，并且由于每个泰森多边形单元的转换独立性，因此该过程也是可并行化的。
Quad-Edge的结构十分简单，仅仅使用点与有向边即可完备的描述空间数据模型的拓扑与几何特性。
对于Quad-Edge结构，Org表示当前有向边的起点，eRot表示当前Quad-Egde中逆时针旋转90度的有向边即对偶边。eONext为以当前Quad-Edge 的基点逆时针旋转的下一条Quad-Edge。除此之外，其定义了终点End，左剖分区域Left与右剖分区域Right[7]。另外Quad-Edge 在此基础上定义了基本的创建(Make)、连接(Connect)、接合与分解(Splice)等基本操作过程。具体可参见参考文献[6]。
图1:Quad-Edge示例
经过实验可知，采用了Quad-Edge结构后，整个转换过程的性能得到极大提升，相对于总时间，步骤4的时间几乎可以忽略。因而一种并行化的Delaunay三角剖分算法，对于整体计算性能的提高更为关键。
3、并行Delaunay三角化
本文选用的算法首先按照
，首先将按照之前计算的MBR，个子区域。然后对于每个区域进行分治法求解，再加以合并。由于无论是区域的划分还是每个区域的求解，都存在独立的计算过程，因而都可以进行并行优化。另外这里值得指出的是，由于输入点集的分布未必均匀分布，每个区域内的点数未必均衡。因此需要一定的措施保证各个线程的计算负载均衡，一个较为简单的办法是产生远多于CPU内核数的计算线程，通过一个线程池的管理，保证CPU的各个核心一直都处于计算状态[8]。
算法对于每个小区域采用分治法求解。对于划分的区域仅有2个或者3个点的情况下，生成操作比较简单，直接连接这2个或者3个点构成Delaunay三角形即可。如果每个区域的点数
本文受"十一五"国家科技支撑计划课题：城镇化与村镇建设动态监测关键技术研究(2006BAJ11B)资助。
1
福建电脑2010年第7期
大于等于4个，则需要对其进行划分，并且对于划分计算后的结果进行合并。
图2:2个或者3个点的Delaunay三角形生成
本算法使用的合并的方法如下所示。首先使用Zig-Zag算法[9]获得两子凸包的上下底边线，分别记为B1B2与T1T2。在B1B2上寻一点B3，如果B3属于左子凸包，那么新的下底边为B2B3，如果B3属于右子凸包，则新的下底边为B1B3。对于B3的选择，尽量选择左子凸包的右侧或者右子凸包左侧Y坐标较低的点，以减少后续的优化步骤。同时保证新的底边与原有左右凸包不能有相交，重复这一步骤直至底边与顶边重合。这时合并后的结果还不是Delaunay三角形，还需要执行LOP局部优化过程[10]。其主要过程是将有共同边的三角形合并成为一个多边形，然后使用最大空圆准则[2]检验其第四个点是否在其它三个点构成的三角形外接圆之内。如果是，则对调对角线，完成局部优化处理。
图3：合并过程
值得注意的是，点集中有大量的点位于同一垂直线上的时候，按照此合并算法，可能无法保证算法的鲁棒性，因为底边可能会变成一条垂直线。为了避免这个问题，可以在划分子区域的时候，每一次迭代时变化划分区域的坐标轴方向，即一次横向划分，一次纵向划分，避免过于狭长的子区域出现。
除此之外，给定圆上的三个顶点，为了求出第四个点是否在其外接圆内，较为常见的做法是求出既有点的两两连线的垂直平分线交点即圆心，进而以第四个点到圆心的距离与半径做比较，判定其是否在园内。但是如果既有点三点共线，垂直平分线两两无交点，则必须考虑这种特殊情况。
分析合并步骤可知，无论是上下公切线的查，还是边连接与优化，其时间复杂度都是O(n)。而由于
分解步骤的时间复杂度为O(1)，合并步骤的时间复杂度为O(n)，设总时间复杂度为T (n)，则T(n)=2T(n/2)+O(n)，因此可以推算出算法的整体时间复杂度为O(nlogn)。
4、实验结果
按照前文描述的算法，我们使用.NET Framework4.0(C#)平台实现了相应的计算程序。同时我们也使用ESRI ArcGIS Desktop9.3中泰森多边形生成工具对于各组数据进行了计算，用于计算速度的对比以及计算结果的正确性检验。其中各组测试数据都是均匀分布的随机生成点集。
本文所使用的测试环境为，CPU：Intel Core2Quad Q9400(2. 66GHz*4,1333Mhz FSB,6M L2Cache)，内存：2G(DDR3 1333MHz)，操作系统为Windows7(32bit)。实验的结果如下表所示，其中时间单位为秒，采取计算5次的平均值。
表1：计算时间对比
表2：计算加速比
5、结束语
由实验结果可知，本文的单线程生成算法在点集少于8万时较ArcGIS性能具有一定的优势，但是点数
继续增加时，算法的速度较低。但是经过多线程并行优化，在多核处理器平台的支持下，算法获得了较高的加速比，并且依据参考文献[11]，可以推测出这种加速性可以随着处理器的个数或核数的增加而能够得到进一步的提高。
参考文献：
[1]Thiessen A H.Precipitation Averages for Large Areas[J].Monthly Weather Review.1911(39):1082～1084.
[2]Delaunay B.Sur la Sphere Vide.Bulletin of the Academy of Sciences of the USSR[J].Classe des Sciences Mathematiques et Naturelles,1934(8):793～800.
[3]Peter Su,Robert L.Scot Drysdale.A Comparison of Sequential De-launay Triangulation Algorithms[J].Computational Geometry7(1997) 361-385.
[4]Rex A.Dwyer.A Faster Divide-and-Conquer Algorithm for Con-structing Delaunay Triangulations[J].Algorithmica(1987)2:137-151. [5]Ivan E.Sutherland,Gary W.Hodgman.Reentrant polygon clipping. Communications of the ACM[J].Volume17,Issue1(January1974):32-42.
[6]Leonidas Guibas,Jorge Stolfi.Primitives for the manipulation of general subdivisions and the comp
utation of Voronoi[J].ACM Transactions on Graphics(TOG).Volume4,Issue2(April1985):74-123.
[7]李翔，王卫安.基于四方边缘结构的实用TIN快速构建[J].测绘工程. 200716(6):29-32.
[8]张明敏，潘志庚，郑文庭等.散乱点集Delaunay三角剖分的分布并行算法[J].计算机辅助设计与图形学学报.200012(7):484-487.
[9]Jonathan Richard Shewchuk.Stabbing Delaunay Tetrahedralizations[J]. Discrete&Computational Geometry.200432(3):339-343.
[10]Lawson,C.L.Software for C1surface interpolation:in Rice,JR(ed.), MathematicalSoftwareIII[M].Academic Press,New York.1977:161-194. [11]GE Blelloch,GL Miller,JC Hardwick,D.Talmor.Design and Imple-mentation of a Practical Parallel Delaunay Algorithm[J].Algorithmica.1999
柳智惠
24(3-4):243-269. 2
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议