matlab中绘制 维诺图（Voronoi Diagram）

1.专业术语（相关概念）：
基点Site：具有一些几何意义的点

细胞Cell：这个Cell中的任何一个点到Cell中基点中的距离都是最近的，离其他Site比离内部Site的距离都要远。

Cell的划分：基点Site与其它的n-1个点所对应的那个平分线所确定的那个离它更近的那个半平面。把所有这些半平面公共的交集求出来就是这个cell.

因此每个Cell都是凸的，图中一定会有一些Cell是无界的

 二维平面上的Voronoi图，存储Voronoi图时只需要存储点和边界信息。

 Voronoi图中包括Voronoi点（线段与线段之间的交点）、线段Segment、射线Ray或直线，这些线都被称为Voronoi edge

2.相关特性：

（1）非空性：每个Site都拥有一个非空的Site

（2）空圆性：在任何一个Site所对应的Cell中，任何一个点以这个点为圆心，以到这个Site的距离为半径的圆必然是空的。

这里的空指的是：圆其中它的内部不包含任何的Site

空圆的另一种形式：点x位于Vertex Edge上，x处于几个Cell的公共交上，那么x到这些Cell对应的Site的距离都是相等的，以x为圆心以这个距离为半径做出来的圆也是空圆。这个空圆也是最大的。

 

辨析两个概念：disk是实心的包括内部的整个部分，Circle是一个圆环，只考虑它的边界。

讨论：以一个x为圆心做出的空圆的圆环上最多会有几个Site？（度数为几）

①、如果x在Cell内部，那么它一定是以到该Cell内部的Site为距离作圆，因此只会有1个，度为1。

②、如果x在一条边界上（不在Voronoi Vertex处），那么就至少有两个Site是最近邻，度至少为2

③、如果x在一个Voronoi Vertex上，那么它是多少个Cell的公共交处那么就有多少度。度至少为3

但是多个Site共圆的情况是一种退化情况，本课程不考虑这个，因此我们在这里认为如果是在Voronoi Vertex上的点那么它的度数就为3。

3、复杂性

为了对Voronoi图进行复杂性的讨论，我们将那些射线强行指向一个无穷远的点，那么此时的Voronoi图就是一个连通图。

那么此时Voronoi图顶点的数目大致在O(2n-4)的范围，边大致在O(3n-6)的范围。无论是顶点还是边，都是线性的。

证明：令D是Voronoi图中所有的Cell（即面face）的周长的总和。3f≤D=2e。引入欧拉公式v-e+f-c=1。c是连通域。在Voronoi图里连通域c=1。将这两个公式联立可以得到顶点数目与边的数目都与n个Site成线性关系

 4.代表性

我们可以用计算机高效的表示以Voronoi图为代表的子区域剖分的几何结构。

平面图嵌入：可为图中的每个顶点安排一个具体的位置，并且能保证它们之间的那些连边可以画在这个平面上，且互相 不会再内部相交。

subdivision（子区域剖分）：将平面图所依附的平面这个空间做一个气的剖分，把它剖分成一个个类似于Voronoi图的单元部分。如此划分出来的一个个单元格称之为face。（对于Voronoi图而言face就是Cell）

face与face之间的关系：它们联合在一起能够将整个平面覆盖住；任何两张face在内部都不会有交。

Voronoi图是subdivision的一个特例。

 Fary Theorem（法里定理）：如果能将一幅图用一种曲线形式画在平面上的话，那么也必然能够以这种完全使用直线的方式把它画在平面上。如下图

对于一个Voronoi图，我们不仅要存储每个Vertex，每条edge，每个Cell，还需要存储一些点与边与面之间的关联关系，这些关联关系的incidence需要记录下来。另外，我们可能还需要对图进行某种遍历。为此，我们借助DCEL结构来进行操作。

5、DCEL

下面是DCEL的一系列数据结构

表边

边分解：将subdivision中连接于以任何两个顶点中之间的一条无向边拆解成一对有向边，这对边是彼此共存，又叫做twin edges。其中一半边叫half edge。我们用edge list边表来记录拆分出来的这一系列半边。

Half-edge（半边缘）的数据结构组成：id，twin边，半边的起始顶点ori，这条边参与围成的面inc（半边的左侧为内侧，记录这个内侧边），半边的前驱pred和后继succ。

Half-edge的数据结构组成：id，twin边，半边的起始顶点ori，这条边参与围成的面inc（半边的左侧为内侧，记录这个内侧边），半边的前驱pred和后继succ。

这里的前驱和后继是顺着在这个面里的，不能够一下子跳到另一个面去

• 点表

顶点的数据结构：id，坐标x,y，第一条向外发出的关联的边inc。

• 面表

face的数据结构：id，一条incident half-edge(inc)（这个inc就是指一个面做遍历的时候的起始边）

这里介绍了更简单的绘制方法

如图所示，我们给定一个face f后，就可以利用其中的信息找到它的一条边，然后通过前驱后继即可遍历这个面围成的边。如果要得到旁边的面，那么通过twin的信息可以翻到另一个面去，它在另一个面也会有相应的前驱后继。

6、Hardness

对于一般意义上的Voronoi图：一维的Voronoi图构造的时间复杂度通过sorting问题来进行reduction。

sorting的输入是一系列的数，可以转化为一维Voronoi图中一系列的点在数轴上。然后我们要将Voronoi图一系列的点进行排序，并把它的输出转化为sorting的输出。输入是数轴上的蓝色的点，然后要将这些点对应的Voronoi图构造出来。而构造出来的Voronoi图的点就是图中一系列红色的点。这些红色的点是两个蓝色点的中点。这样我们得到了Voronoi图的一维的答案，然后要将这些答案转化成排序问题的输出。我们可以通过一次线性扫描，得到第一个0号点的位置，因为知道了红1是蓝0和蓝1的中点，那么必然知道这个中点的距离，那么可以通过蓝0和红1知道蓝1的位置，如此类推，就会得到之前蓝色输入的点的顺序，也就规约到了sorting上。

因为Voronoi图构造可以规约到sorting，那么这个问题的下界就是O(nlogn)

对于二维的Voronoi图，也可以规约到sorting上。

如图所示，我们输入了n个未排序的数，我们的任务是将这n个数映射到一个单位元的圆周上。然后这n个点就是二维Voronoi图在平面上的一个输入点集。再用二维Voronoi图的算法把输出构造出来。这个Voronoi图只有唯一一个vertex就位于这个圆的中心x，有多少个点，就有多少条边，每一条Voronoi边都是在圆周上相邻某一对点的平分线。对于这个Voronoi图的输出，我们只要知道其中一个点，就可以确定它的Cell，然后通过数据结构里存储的信息可以将这个图遍历出来，也就得到了Sorting的排序。

因此二维的Voronoi图的构造的下界也是O(nlogn)

然而上图是一个退化的Voronoi图，对于非退化一般性的Voronoi图而言，也可以设计这样一个reduction。

如图所示，我们用一个抛物线来代替圆弧。将输入的点映射到一个抛物线上。然后可以通过简单的操作在线性时间内找到起点，通过一次一次的flip逐一将所有点枚举出来，那么就将输出映射回sorting的输出。

因此确凿明确的知道2维Voronoi图的问题复杂度下界不会低于O(nlogn)

7、Sorted Sets（分类设置）

讨论：如果在给定一些附加条件后，Voronoi图的难度是否会降低？例如在某些时候我们的输入子集是蕴含着一些排序的信息，像凸包那样的问题就可以降低维度，那么Voronoi图是否也可以降低难度？

结论：即使输入的点是按照高度的次序由低到高给出来的，我们依然不能使Voronoi图这个问题从下界意义上讲的效率有丝毫的改进。

解释：凸包和Voronoi图的结构其实是有本质区别的。凸包的在本质上是一种Combinatorial structures组合结构，而Voronoi图本质上是一种geometric structures几何结构。这二者的区别凸包在仿射变换下具有不变性，即虽然形状发生了变化，但它的极点和非极点以及极点之间的连接关系是不会发生丝毫的变化的，它的拓扑结构是保持的。

而Voronoi图在经过仿射变换后经常会发生拓扑结构的改变。

如图所示是第一个图先压缩再拉伸，Voronoi图的拓扑结构发生了本质上的变化。

发生这个的原因是因为Voronoi图是依赖于距离，如果是发生与距离有关的改变，那么就会对它的拓扑结构造成影响。

8、VD_sorted

考虑一个特例：Voronoi图给出的输入是在某种意义下有序。

ε-closeness问题：每次输入n个实数以及一个ε的正数，问在所有这些数中，是否存在两个它们的间距是不超过ε的。这个问题的复杂度是O(nlogn)

接下来要做的就是VD_sorted这个问题规约到ε-closeness问题上去。

Lifting Transform：对于输入的一系列的点，在ε的高度区间内， 按照输入的次序把每个点拔高到对应的区间上。如下图

这个对应的公式是：

接下来要针对一组任意的输入调用任何的包括优化的Voronoi图算法，去构造出它对应的Voronoi图。（构造是指要得到例如表示为一个DCEL结构的Voronoi图）。这种情况下我们可以依次枚举出与x轴相交的所有的那些cell。接下来将这些Site的点在投影回x轴，其实就是一开始ε-closeness的输入。

下面分两种情况讨论，无论哪种情况我们都可以得到ε-closeness的答案。区分的准则是这里的每一个site在经过投影之后回到当初那个点是否就位于提升后的这个site所对应的Voronoi图中的那个cell中。

如图所示，1号点和2号点是满足要求的，其他则不满足。

Case A：Yes 所有的点都无一例外的落在自己所对应的Cell中。

那么就可以得到所有的Cell的次序和它对应的site的次序。这种情况下我们就可以得到所有输入点排序的结果。知道了这些点的排序那么也就很容易知道相邻点之间的间隙，而所有这些相邻点间隙中的最小者就是MinGap，有了它就可以回答ε-closeness这个问题了。

Case B：其他情况，有点不落在自己对应的Cell中。

如图所示，I点投影下来的i点就并不落在Cell I中。那么连接i点和它所落在的Cell J里，再加上J的投影点j构成了一个直角三角形。我们可以快速得到ij<iJ<Ii<ε，那么在整个的输入的那些实数中，至少存在i和j两个实数它们的间距是小于ε的。也就是说这种情情况下ε-closeness问题的答案是Yes。

综上所述，我们只要将ε-closeness问题转化为VD_sorted这个问题，那么无论是哪种情况都可以反过来根据VD_sorted问题的输出顺利的得到ε-closeness答案，从而完成归约。这就意味着VD_sorted这个问题的难度依然是O(nlogn)。

9、Naive_Construction

一个最直接最简单的方法：一个Cell一个Cell去构造。对于单个Cell构造就是按照定义找出这个Site与其他Site的平分线并且求交，得到的包含Site的封闭的平面就是一个Cell。

这个算法复杂度在找平分线求交时会有n²的难度，另外还有n个Cell，最后会高达O(n³）

10、Incremental Construction 递增式建筑

Incremental algorithm：算法核心是如何从一个规模为k-1的Voronoi图通过引入一个新的site构造出规模为k的Voronoi图。

 算法描述：如上图，这里原先的Voronoi图是已经按照DCEL的数据结构存储好，并且图中的编号顺序只是为了描述方便不代表输入时是有这个次序。假设引入一个新的site点4，那么首先和0号site一起找到它们的平分线，并且通过0Cell的边界点找到这条平分线的端点，然后在通过DCEL结构中存储的内容翻到隔壁1号Cell中继续找1号site与4号site的平分线，再找到它的端点。就这样依次寻找直到封闭。

如图所示，在引入黄色的点后产生新的平分线的顺序是1234，那么Voronoi边的遍历顺序则是acdefghijb。

算法复杂度：从求交开始每次沿着一个平分线的方向去做一次求交，由于Cell本身已经存为DCEL结构，根据Voronoi图的性质每一个Cell都是凸的。在一个凸多边形中去对边界做一个这样的穿刺的求交是在O(logn)的时间内完成的。用O(1)的时间翻墙，再用O(1)的时间确定平分线的方向继续前行。而每一个新引入的点边界的规模不超过n，整个算法是引入n个点，因此最后算法的复杂度是O(n²logn)。

11、Divided-And-Conquer

分而治之：首先将整个坐标的点集沿着x轴坐标做一个预排序，然后将点集大概平均地分割成更小的两个点集进行递归构造Voronoi图，最后Merge。

伪代码：

DacVoronoi(S, n)
{x-sort all sites into S = {p1, p2, ..., pn};return dacVD(S, 1, n);
}
dacVD(S, i, j)
{return (j - i < 3) ? trivialVD(S, i, j) : merge(dacVD(S, i, (i + j)/2), dacVD(S,  (i + j)/2 + 1, j);
}

处理重叠部分

在Merge两个Voronoi子图的时候经常会产生重叠，对于重叠我们要根据它们之间的平分线的平面来对它们进行分割。平分线对应的左平面切除右边的Cell，反之右平面切除左边的Cell。

如果单纯用这种方法强行处理一系列的重叠的话，复杂度会比较高，下面介绍一种优化的方法。

轮廓线（Contour）：介于左边的子图和右边子图之间的边界，它到左侧site和右侧site的距离是相等的，它具有等距性。因此Contour上的任何一个点相对于最后要合成的Voronoi图来说必然是某一条edge上的点，所以Contour必然是由最终的Voronoi图上的若干条edge前后依次相连而构成的一条裂缝。

Contour的性质：

①具有y方向上的单调性，是y方向上单调的一条折线，即任何一个水平高度上这样的contour line只会至多有一个交点。

②一条contour在自向而下的行进过程中，它的第一段即上面通向无穷的那一段必然是左右两个子集的公切线的平分线。对称的最后一段即通往负无穷深度的奶段也必然是lower common tangent的平分线。

③若左右各有n和m个site，那么contour的长度不会超过n+m

下面进行Merge的讲解

算法描述：首先找到Contour入手的第一条边，即先找左右子图的公切线（与凸包算法中的类似只需线性时间），然后求它的平分线记做b；接着自顶而下的沿着contour链追踪下去，对于任意一对左右site，我们都找出它们之间的平分线，这里是一个迭代的过程；对于每一次找到的平分线b，要对左边的Cell和右边的Cell同时用b做一次裁切（求交找方向）；裁切后要更新下一个Cell，这里的的关键是看这个平分线是最先从哪个cell出来的，如果一个先出来那么另一个就要保持。而先出来的Cell要替换为它同侧的一个后继。我们要用它的后继来构成下一对需要排解冲突的Cell pair。对于新的Cell，它们的平分线又要重新确定。这个过程一直持续到从contour的最后一条边出来。（这里的描述比较乱...）

伪代码描述：

ComputeContourBetween(VD(SL),VD(SR))
{compute the upper tagent of conv(SL) and conv(SR);//let l∈SL and r∈SR be tagent sites and b be the bisector of segment lrtrace the contour from top down;find b ∩ Cell(l) and b ∩ Cell(r);clip and then flip the cell whose boundary intersects b first;update l or r;update b;
}

对于Merge的算法复杂度主要体现在求交上，因为有可能会在某一个Cell上反复的与其他平分线求交，这里的反复计算就对复杂度造成了很大影响。

由于Contour具有凸性，那么我们可以利用这一性质来对Merge算法进行优化。

如图所示，我们可以发现Contour每次拐弯都是同一个方向的，这意味着contour上这样一系列的边参与构成了这个site经过剪裁以后更小的那个cell的边界，它们会构成一条凸的环路。那么具体来说如果是在左侧的某一个cell中，则所有的拐向都是朝右，反之亦然。

而我们在求交过程中，虽然能确定拐的方向，但是如上图所示我们还会得到类似点1那样的无用点，其实我们只需要用更近的那个交点。正是有这些貌似无用点的求交导致复杂度的增加，因为每次都要从起点遍历。由于contour的凸性，我们可以得到这些貌似无用的交点在整个这个cell的边界上的分布是前后单调的。

因此我们的算法优化就在于，假设黄色点1是我们找到的第一个交点，那么接下来的遍历只需要从点1开始即可而不需要再从头开始遍历。

改进后，所有这些交点的求法都可以在前一个交点的基础上继续搜索，它的成本就是连续的两个交点之间的一段弧的长度，总体不会超过这个Cell的轴长，是一个线性的操作。

那么基于这样的一个合并的基本算法，按照分而治之的框架所得到的总体的Voronoi图的构造算法复杂度不超过O(nlogn)。

12、Plane-Sweep

假设用水平线自顶而下扫描平面，再引入了一些抛物线线段如下图。

我们在曲线上面的部分叫做封存部分（Frozen region），曲线和扫描线之间的叫做未冻结部分（Unfrozen region）。Frozen部分的特性从数学上来说是它到扫描线的距离都要大于它到以上已扫描过的点集的距离；反过来，Unfrozen部分都是到这个扫描线的距离更短，至少相对于已经扫过的那些site而言更短。

我们将那条曲线取名叫做Beach Line（海浪线）。数学上的定义是这条线上的点到扫描线的距离等于它扫过的site的最下的距离。而在Frozen内部的点都已经是确定的落在某个Cell，因为它到扫描线的距离都要大于它到已扫过的点集中离它最近的那个site的距离，因此Frozen region是已经确定了的部分。

准备工作：

抛物线的定义：到准线和焦点距离相等的点的集合。性质：准线离焦点越近，抛物线开口越狭窄。当准线与焦点重合时，此时的抛物线是一条直线。

对于我们这里而言，准线就是那条扫描线，每个扫过的site都可以是焦点，可以说扫过多少个site就有多少条抛物线。beach line就是最底下的轮廓，称作包络（envelope）。这些抛物线有的并不贡献Beach line有的贡献好几段，一条抛物线最多贡献n条弧。而在任意时刻beach line上的弧的数目不超过2n-1段。

在beach line中两段弧的交点称作breakpoint，这个交点一定会位于某条Voronoi edge上。随着扫描线的推进，breakpoint会不断前进构成Voronoi edge，直到扫描线碰到另一个site，此时的breakpoint就是Voronoi vertex。

当然实际算法过程中不可能处理无数个点，我们需要针对一些特定的位置进行处理，这些我们也称作Events事件。这个Events指的是整个扫描线上的抛物线的拓扑结构发生变化的时刻。这些变化的时刻分为了两类。

（1）、Circle Events

circle event是动态生成的，随着算法的进行不断的发现circle event。

若在beach line上出现了新的三个相邻的伙伴弧，分别对应于i,j,k三个site，它们所对应的三段抛物线弧首尾相连，构成这样的局部三人组必会定义一个圆，那么这个圆最下面的点就是一个circle event。

由Voronoi图定义可知，若出现这么一个圆，那么扫描线上这个圆一定是空的。若抵达最低点时圆仍是空的，那么此时i的抛物线和k的抛物线的交汇点就是一个新的Voronoi vertex。在这个缝合处会长出一条新的边，这条边的方向就是i和k的平分线的方向。这个过程会不断的演化下去。

在碰到circle event后①将j对应的弧从beach line中删除②生成一个新的vertex③将原来的两条绿色的平分线的终点绑定在新的vertex上④生成一条新的边，方向为i k的平分线方向⑤删除包含j的circle event⑥检测i k的两边是否会有新的三元组产生。

（2）Site Events

扫描点经过适当扫描后碰到一个site，①确定这个site上方对应的弧是哪段弧②将对应的弧p用一段无限小的弧split成两段③随着扫描线往下，焦点s与扫描线生成的抛物线在弧p中间对应的那一小段弧会不断扩大，此时对应的两个breakpoint连接的线称作悬垂线④将原先p和q的三人组对应的circle events删除⑤对于每个新产生的三元组，创造新的circle events。

当最后抵达新产生的circle event c时，之前的悬垂线和以前p和q的平分线就会汇聚到点v中并产生黄色的新的平分线。

%采用matlab自带的函数进行绘制
clear
xdot=gallery('uniformdata',[200 2],5);
%delaunay三角形
figure(1)
DT=delaunayTriangulation(xdot);
triplot(DT,'color','k')
%voronoi三角形
figure(2)
voronoi(xdot(:,1),xdot(:,2));
xlim([0,1])
ylim([0,1])

A GPU Approach to Voronoi Diagrams