构造voronoi图

出处http://blog.csdn.net/shunan/archive/2007/05/28/1628238.aspx

Typedef struct Node

{

  Int x;

  Int y;

}MyNode;

Typedef struct Nodel

{

  Int start;

  Int num;

  Edge *firstEdge;

}MyNodel;

Typedef struct Edge

{

Int mainNum;

Int adhereNum;

Double angle;

Struct Edge *pre;

Struct Edge *next;

Struct Edge *twin;

}edge;

Typedef struct MPIEDGE

{

 Int mainNum;

 Int adhereNum;

 Double angle;

}MyEdge;

 

 

1，MyNode结构，节点数组保存的是平面上点集的横坐标和纵坐标

2，edge结构，边节点，保存与序号为mainNum的点相邻的点信息：adhereNum为此相邻点的序号；angle为序号mainNum的点到此相邻点的向量到序号mainNum的点到序号为start(在结构MyNodel中)的点的向量的夹角的弧度(0<=,,,<=2* )；pre，next域分别为双向链表中的前一个和后一个边节点；twin域为与此节点对应的一个边节点(即序号为mainNum的点成为序号为adhereNum的一个相邻点)

3，MyNodel结构:start为第一个与此点相邻的点的序号；num为当前与此点相邻的点的个数；firstEdge域指向与此点相邻的边节点中angle最小的那个节点(序号为start的那个相邻边节点的angle为0)

4, MyEdge结构，只是提取了edge结构的前三个元素组成一个节点类型，主要是为了满足MPI函数在通信时(自定义通信元素的类型)确定边的类型用的

 

如果只需实现串行的voronoi图算法，则可以把MyNode和MyNodel结构合并，也用不到MyEdge结构了

 

输入：坐标平面内的n>1个无重合的点组成的点集P

输出：构造好的Voronoi图V(P)

 

1                对这n个点按照x坐标值为主关键字，y坐标值为次关键字进行排序

2                开始构造：

2.1        if  n==2  then 构造两个点的Voronoi图，并返回

2.2        if  n==3  then 构造三个点的Voronoi图，并返回

2.3        把n个点平均分成两部分Pl,Pr

2.4        把点集Pl构造成Voronoi图 V(Pl)

2.5        把点集Pr构造成Voronoi图 V(Pr)

2.6        合并V(Pl)，V(Pr)为V(P)

2.7        返回V(P)

输入：已生成的V(Pl),V(Pr)

输出：合并后的V(P),P=PlUPr

 

1        分别求Pl,Pr的凸包CH(Pl)，CH(Pr)

2        BCT ß bottomConnected_tangent(CH(Pl)，CH(Pr))

3        UCT ß upConnected_tangent(CH(Pl)，CH(Pr))

4        L ß left point of BCT 

Rß right point of BCT

5        while (not (BCT =UCT)) do

5.1        isRight ß isLeft ß false

5.2        initEdge(edgeL，edgeR，L，R)

5.3        insertEdge(L，edgeL)

insertEdge(R，edgeR)

5.4        /*对Pr中与R相邻的点进行测试(以顺时针的方向)

5.4.1        pre ß R中edgeR顺时针方向下一个半边节点

5.4.2        if pre != NULL and pre关联的点在LR向量的左侧 then

5.4.2.1  pre1 ß R中pre顺时针方向下一个半边节点

5.4.2.2  while pre1 != NULL and pre1关联的点在L，R，pre组成的三角形外接圆内  do

5.4.2.2.1        deleteEdge(pre)

5.4.2.2.2        pre ß pre1

5.4.2.2.3        pre1 ß R中pre顺时针方向下一个半边节点

5.4.3        else isRight ß true

     

5.5        /*对Pl中与L相邻的点进行测试(以逆时针方向)

5.5.1        next ß L中edgeL逆时针方向下一个半边节点

5.5.2        if next!=NULL and next关联的点在RL向量的右侧

5.5.2.1  next1 ß L中next逆时针方向下一个半边节点

5.5.2.2  while next1 !=NULL and next1关联的点在L，R，next组成的三角形外接圆内 do

5.5.2.2.1        deleteEdge(next)

5.5.2.2.2        next ß next1

5.5.2.2.3        next1ß L中next逆时针方向下一个半边节点

5.5.3        else isLeft ß true

5.6        if isRight then L ß next关联的点序号

5.7        else if isLeft then R ß pre关联的点序号

5.8        else if next关联的点在L，R，pre组成的三角形外接圆内 then

          L ß next关联的点序号

5.9  else R ß pre关联的点序号

6  initEdge(edgeL，edgeR，L，R)

7  insertEdge(L，edgeL)

   insertEdge(R，edgeR)

initEdge(edgeL,edgeR,L,R)函数初始化两个半边节点edgeL和edgeR,其中edgeL将来插入到L的边链表中，使R成为L的相邻点(edgeL->mainNum ß L,edgeL->adhereNum ß R)；edgeR则刚好相反。

insertEdge(number,edge)函数把半边节点edge插入到序号为number的点的边链表中，由于每个点维护的一个循环边链表是按照相邻点的逆时针方向存储的，为了在每次插入后仍然保持这个特性，必须在插入前，先计算要插入的点到第一个点(MyNodel结构中的start序号所示，在第一次插入时赋值，以后就算此节点被删除，也不改变，以作参考)的夹角，按照夹角的大小来判断插入的地方(夹角有序排列即意味着那些相邻点按照逆时针方向排列)。

deleteEdge(edge)函数就是把半边节点edge和它的对应半边(twin域指向的半边节点)删除。

由于篇幅过长，对于算法的说明在此省略(如求凸包，BCT，UCT等)

 

VORONOI图分治算法并行化

目前已经存在的构造Voronoi图的并行算法有两种方式：一种是使用与平面点集同数量级的p=O(n)个处理器个数，对基于MPI的并行算法，能达到O(lgn*lgn*lgn)的复杂度，比较适用于对时间要求非常高的场合，显然这样的处理环境成本也是非常的高，而当平面点集数量较少时，也体现不出此类算法的优越性，因为通信会占用掉大部分的时间，举个bitonic sort的例子，应该可以说明此类情况，由于我的并行机拥有32个处理器，在运行由MPI实现的bitonic sort算法时，对64个数据进行排序(每个处理器存放两个数据)时，时间就已经是毫秒量级的了，而对64个数据进行快速排序时间是可以忽略的，在微秒量级，但是从理论上分析bitonic sort，它的时间复杂度只有O(lgn*lgn)，是一种非常适合并行化的排序算法，可见，由于它使用的处理器多，相对于数据处理的时间，处理器之间的通信时间占了整个算法执行的非常大块的时间。另一种是使用常数p个处理器来对平面点集中的n个点进行处理，此类算法易于实现，而且算法提速明显，非常适用于一般场合的计算。

构造Voronoi图的分治算法的并行化，理论上也可以参考上述的快速排序并行化的方法，但是涉及到图的构造，并行处理时，各处理器之间需要交换的信息非常得多，不像并行快速排序那样，只需要把本身的数据分成两部分，可以整块的发送给对应的处理器，那么方便。

 

Delaunay三角剖分

最近点意义下的Voronoi图的直线对偶图就是Delaunay三角剖分，因此可以把上上篇描述的构造Voronoi图的分治算法用于Delaunay三角剖分。

增量算法：易于实现，使用广泛，适合小规模点集的三角化。具体过程如下：

1 遍历所有散点，生成一个包含所有散点的大三角形(顶点不在点集中)

2 有未处理过的点p，插入之，否则结束算法，退出

3 在已剖分好的三角网中找出包含点p的三角形t，把p与t的三个顶点相连，构成三个三角形

4 根据优化准则对局部生成的三角形进行优化(互换对角线等)

5 返回第2步

局部变换法：根据Delaunay三角剖分的性质2，首先构造一个不满足Delaunay三角剖分条件的三角网络，然后对两个共边三角形构成的凸四边形迭代换边使之满足Delaunay三角剖分的条件(主要是交换对角线的方法)