问题描述和算法背景

算法概述

实验难点

数据分析

软件用法

实例展示

参考文献

概述

在计算几何中，离散点构造Voronoi图是一个非常基础且应用广泛的问题。N个离散点按照最邻近原则划分区域，每个点与它的最近邻区域相关联。本实验重点研究Voronoi图的三维情况，运用分块的方法动态实现了单一Voronoi细胞和三维Voronoi图的构造。同时，本实验对均匀分布的随机点产生的性质进行了一系列研究和分析。

问题描述

在二维（平面）情况下，Voronoi图是由一组连接两邻近点直线的垂直平分线组成的连续多边形组成。N个在平面上有区别的点，按照最邻近原则划分平面。每个离散点则各自拥有一个细胞区域，区域内部的点到相对应的离散点距离最近。在高维情况下，连接两邻近点直线的垂直平分不再是直线，形成Voronoi图的细胞则会从平面上的多边形转变为高维的多面体。
Voronoi图的平面情况和三维情况的区别不光体现在细胞的构成上。平面情况下，由于可以应用欧拉公式，我们能知道构成Voronoi图的顶点和边数和离散点数为相同数量级。但涉及到三维情况：首先，每个细胞区域构成多了面的概念；其次，构成细胞的顶点数和边数在数量级上也有了质变，可以达到离散点数数量级的平方。这为构造三维Vonoroi图带来了困难，导致了时间复杂度和空间复杂度的增加。三维Voronoi图主要有Quickhull[11]和分块动态实现[12]两种算法，但它们的时间复杂性最差情况下都可以达到θ(n^2)。
实验中，我们首先动态实现了单一Voronoi细胞的构造。事实上，在许多物理应用中，人们更关注的仅是几个单一Voronoi细胞，而不是整个Voronoi图的构造。运用单一Voronoi细胞的构造，结合分块方法，我们动态的构造三维Vonoroi图。尽管在最差情况下，时间复杂性可以达到θ(n^2)，但如果考虑离散点是随机均匀分布的，实验结果显示平均复杂度随着离散点数的增加线性增长。

相关工作

自从1907年Georgy Voronoi提出了Voronoi图的概念[1]后，Voronoi图这种几何结构就被广泛的应用到了各个领域，比如形状建模[2]、运动规划[3]、科学可视化[4]、碰撞侦测[5]、几何学[6]以及化学[7]。针对三维情况的N个离散点，每个点的Voronoi细胞是一个三维的多面体。
三维Voronoi图有许多变种问题被人考虑过。如果不给区域限制，那么离散点构成的凸包边界上的离散点的细胞区域是无界的。但在许多应用中，我们考虑的是对某个封闭的区域进行划分，构造Voronoi图。这就涉及到了封闭区域的形状性质——一个很关键的问题，区域是否是凸的。有许多篇文章[8]对区域问题进行了分析。
另一个被广泛考虑的问题是离散点的分布情况，这在一定程度上影响了算法的时间复杂度和空间复杂度。一个熟知的结论是如果离散点是独立均匀分布于空间中的，那么有平均复杂度仅为O(n)的算法能构造出三维Voronoi图。更复杂的，若考虑离散点符合某种条件的高斯分布，则也能在一定程度上降低算法复杂度。[9][10]两篇文章针对不同条件的高斯分布进行了分析和算法设计。

算法部分的代码是用C++实现的，主要围着两个大类展开。Voronoicell类包含了如何构造一个单一Voronoi细胞。该类包含初始化、构造和输出功能，用Voronoi细胞（三维多面体）的顶点以及它们连成的边来表示该Voronoi细胞。另一个大类是container类，主要解决了三维Voronoi图的构造和表示。该类能调用Voronoicell类，并在这基础上包含一系列计算Voronoi图的函数。为完善Voronoi图的构造以及每个细胞邻居的计算，我们需要另外两个类，分别是voronoicell_neighbor和container_poly。

单一Voronoi细胞

单一Voronoi细胞是三维的一个凸多面体，主要用voronoicell类实现，用一系列相邻的顶点表示。该类包含一系列函数，能计算和输出对应某一中心点的单一Voronoi细胞。我们使用正方体作为该Voronoi细胞的初始架构。
算法是动态的，每次新加入一个离散点相当于引入了一张切平面。主要分成四个步骤：
1）判断切平面是否与目前的Voronoi细胞相交；
2）如果相交，找出原来Voronoi细胞上某条与切平面相交的边（特殊情况下，可能是点）；
3）由该边出发，有方向性的找出原Voronoi细胞的面，重构Voronoi细胞的边和顶点结构；
4）删去原Voronoi细胞多余的边和顶点。
单一Voronoi细胞主要由如下形式表示：
p：public变量，记录Voronoi细胞上的顶点个数
pts[]：记录Voronoi细胞上顶点的坐标
nu[]：记录Voronoi细胞上顶点的度数
ed[][]：记录边，针对第i个顶点，ed[i]有2n_i+1个元素。前n_i个元素表示第i个顶点连出去的边，即ed[i][j]=e(j,i)，这里的e(j,i)表示对第i个顶点来说的第j个相邻点。接着的n_i个元素表示第i个顶点相对于它相邻的顶点是第几个顶点，即e(ed[i][ n_i+j],e(j,i))=i，这是为了方便删除时的维护工作（如果要删除第i个顶点，它的相邻点能迅速定位需要删除的边）。最后一个元素指向自身，作为分配内存使用的一个返回标志。

3D Voronoi图

我们仅考虑正方体区域内离散点的Voronoi图结构的构造，主要用container类来实现。核心方法是compute_cell()，该方法调用voronoicell类并结合分块的方法来构造某一特定的Voronoi细胞。
算法的核心思想是新加入一个离散点v后，不断的寻找邻近的离散点。每个邻近的离散点对应着一个切平面，不断的重新计算新加入离散点的Voronoi细胞，直到确定其它离散点离新加入的离散点v足够远，不会对目前构造出的Voronoi细胞产生影响为止。
可以看出，算法中一个重要的问题是如何判断什么时候可以终止寻找邻近离散点的过程，我们主要有以下两种方法判断：
1）半径测试：如果目前构造出的Voronoi细胞距离细胞中心（新加入的离散点v）的最远距离为R，则距离离散点v超过2R的离散点不会影响到Voronoi细胞的构成。半径测试非常简单，但也有缺陷——没有方向性。假设Voronoi细胞在某一个方向特别长，那么其它方向的离散点也还需要进行搜索检测。
2）区域测试（代码见附录2）：我们可以测试一个特定区域内是否有点可能影响到离散点v的Voronoi细胞的架构。方法是运用一系列平面测试对该特定区域内最接近离散点v（Voronoi细胞中心点）的点进行测试。这个过程比较慢，但具有方向性，即若该区域没有点有足够影响力，则相对于该区域更远的区域也不会有点有足够影响力。
算法结合了这两种测试：每测试几块邻近的区域就会做一次半径测试，如果是离散点独立均匀分布，经常做前几次半径测试就能停止。如果半径测试没有通过，那么就进行区域测试，以防出现极端情况使得程序需要检测非常多区域的点。
核心方法compute_cell()可以大致分为如下五个要点：
1） 初始化。确定维护区域，即一个正方体。将正方体划分为64个小正方体，并为每个小正方体计算出一张工作链表，用来记录针对特定小正方体需要考虑的邻近小正方体；
2） 新加入离散点v时，定位该离散点所位于的小正方体，对该小正方体内部的其它离散点运用voronoicell类构造单一Voronoi细胞；
3） 用工作链表测试相邻的区域，每过几步做一次半径测试看是否能终止计算；
4） 如果工作链表测试到了底部，将邻近区域加入新的测试链表；
5） 每次对链表的第一个区域执行区域测试。如果该区域需要被检测，则检测其内部所有的离散点，并将它在工作链表中的邻近区域加入到链表的末尾。不断重复这些过程直到确定可以终止。

鲁棒性

A．通常情况下，Voronoi细胞顶点度数为3（事实上经过测试，即使加入100000个随机点，顶点最高度数还是3）。但不能排除特殊情况下，新加入离散点产生的切平面使得顶点度数超过3。甚至，当切平面离顶点距离太近时，会导致绘图的困难。我们的解决方案是当顶点和切平面过近时，容忍一个小误差，当做顶点位于面上，并允许顶点的度数超过3。

Figure 1：单一Voronoi细胞顶点度数超过3

Figure 2：3D VD顶点度数超过3
B．当新加入离散点时，构造Voronoi细胞后可能会产生部分顶点的度数为1或者2，这是不允许的。我们的解决方案是最后循环的测试顶点度数：如果为1，则直接删去；如果为2，则将两个相邻点连边，并删去该顶点。这里要注意只遍历一次顶点可能不够，因为删点可能会导致产生新的少于3度的顶点。
C．新加入离散点时，计算切平面时，有可能会产生重复边。针对这种情况，只需最后做一步check_duplicate()即可。

参数选取

构造三维Voronoi图时，我们需要分块，那么分多少块最合适呢？如果分块过多，那工作链表要测试的块过多。但如果分块太少的话，每个分块里有过多的离散点需要计算。事实上，如果离散点数为N，则分块数量为O(√n)比较合适。但输入点数并不是已知的，经测试，64个分块能处理大部分情况。

这一部分，我们主要针对均匀分布的随机点构成单一Voronoi细胞和三维Voronoi图进行初步的分析。我们的数据分析主要分为运行时间分析和特征分析两部分。其中，对顶点、边、面等的特征分析体现了均匀分布对三维Voronoi图构造的影响和特殊性。

运行时间分析

我们实验的离散点数范围从1000到100000（间隔为1000），针对单一Voronoi细胞和三维Voronoi图进行了时间分析。可以看出基本上两个部分运行时间都随着离散点数的增加线性增长。特别是三维Voronoi图，尽管最差情况是θ(n^2)，但平均时间仅为线性。

特征分析

A．单一Voronoi细胞：下一图说明随着离散点数的增加，单一Voronoi细胞的顶点、边和面都会增加。从1到80个点，这三个特征开始时近似于线性增加，在80个离散点时迎来一个拐点，之后增长速度逐渐减慢至趋于平稳。当离散点数增加到足够多时，顶点数趋近150，边数趋近于230，面数趋近于75。下二图说明随着离散点数的增加，单一Voronoi细胞的体积会迅速的下降到趋近于0，在10个离散点前近乎于线性下降，之后下降速度逐渐变慢至趋于平稳。

B．三维Voronoi图：下图是随着离散点数增加，针对最大Voronoi细胞的体积进行研究。结论类似于单一Voronoi细胞，在10个离散点前近乎于线性下降，之后下降速度逐渐变慢至趋于平稳，只是绝对体积要比单一Voronoi细胞里相同离散点数情况大。

可执行文件 和 报告下载

软件运行环境的搭建

算法核心由C++编写(Visual Studio 2013)，3d图形的展示用mathematica编写，窗口的界面用JAVA编写。
程序运行需安装mathematica(9.0.1.0)和java runtime environment(1.7-17)。

程序的运行方法

搭建好程序运行的环境之后，用java -jar 3DvoroDiagram.jar运行打包后的3DvoroDiagram.jar程序

用户可以通过两种方式加入顶点

1.通过设置顶点数目并点击“随机生成点”按钮，随机产生指定个数的顶点。
2.通过点击“增加点”按钮，逐个加入顶点，并可以指定每个顶点具体的坐标。

用户可以通过两种方式删除顶点

1.通过点击“重置”按钮，所有顶点都被清空的顶点。
2.通过在点坐标列表框里选中待删除顶点，并点击“删除点”按钮，逐个，或者通过多选的方式删除被选择

修改3D视角

1.通过点击“修改视角”按钮，填入想修改的视角坐标。

修改显示模式

1.通过复选框，修改3d显示方式为single cell 还是 VD。
2.通过复选框，修改线条和顶点的显示方式是单个点，简单直线，还是球面点，圆柱状直线。

显示算法的不同步骤

通过滑动窗体最下方进度条，控制，查看Voronoi图的生成进度。
特别需要注意的是，调整进度条的数值之后，需点击“重绘”按钮，左边图片才会更新。

备注

1.在程序运行的过程中会在命令行里同步输出mathematica绘图语句的内容。
2.只有在点击“修改视角”和“重绘”按钮时，3d显示区域才会发生重绘。

500个点，显示Voronoi Diagram，500个点已经产生有数万个面,
图像的显示上面用户已经可以开始感觉到很大的延时，
但是算法的底层产生数据的时间依然非常的短，
实际上，核心的c++算法是可以计算更大规模问题的。

500个点，显示single cell

500个点，显示Voronoi Diagram,修改了视角

4个点，显示single cell

4个点，显示Voronoi Diagram

4个点，显示single cell

4个点，显示Voronoi Diagram

4个点，显示single cell

4个点，显示Voronoi Diagram

6个点，显示single cell

6个点，显示Voronoi Diagram

8个点，显示single cell

8个点，显示Voronoi Diagram

27个点，显示single cell

27个点，显示Voronoi Diagram

27个点，显示single cell

27个点，显示single cell，运行到第12步时的结果

27个点，显示Voronoi Diagram

27个点，显示Voronoi Diagram，运行到第12步时的结果

附录：参考文献

[1] Georgy Voronoi, Nouvelles applications des parametres continus a la theorie des formes quadratiques, Journal fur die Reine und Angewandte Mathematik 133, 97–178, 1907.

[2] Pierre Alliez, Giuliana Ucelli, Craig Gotsman, and Marco Attene. Recent advances in remeshing of surfaces. Shape Analysis and Structuring, pages 53-82, 2008.

[3] Avneesh Sud, Erik Andersen, Sean Curtis, Ming C. Lin, and Dinesh Manocha. Real-time path planning in dynamic virtual environments using multiagent nav- igation graphs. IEEE Transactions on Visualization and Computer Graphics, 14(3):526-538, 2008.

[4] Michael Balzer and Oliver Deussen. Voronoi treemaps. In Proceedings of the 2005 ACM Symposium on Software Visualization, pages 165-172, 2005.

[5] Avneesh Sud, Naga K. Govindaraju, Russell Gayle, Ilknur Kabul, and Dinesh Manocha. Fast proximity computation among deformable models using discrete Voronoi diagrams. ACM Transactions on Graphics (Proceedings of ACM SIG-GRAPH 2006), 25(3):1144-1153, 2006.

[6] Christopher M. Gold. What is GIS and what is not? Transactions in GIS, 10(4):505-519, 2006.

[7] Anne Poupon. Voronoi and Voronoi-related tessellations in studies of protein structure and interaction. Current Opinion in Structural Biology, 14(2):233-241, 2004.

[8] Yan D M, Wang W, Lévy B, et al. Efficient computation of 3d clipped voronoi diagram, Advances in Geometric Modeling and Processing. Springer Berlin Heidelberg, 269-282, 2010.

[9] Lucarini V. Three-dimensional random Voronoi tessellations: From cubic crystal lattices to Poisson point processes[J]. Journal of Statistical Physics, 134(1): 185-206, 2006.

[10] Golin M J, Na H S. On the average complexity of 3d-Voronoi diagrams of random points on convex polytopes[J]. Computational Geometry, 25(3): 197-231, 2003.

[11] Qhull code for convex hull, Delaunay triangulation, Voronoi diagram, and halfspace intersection about a point, http://www.qhull.org/.

[12] Rycroft, Chris. "Voro++: a three-dimensional voronoi cell library in c++." (2009).

B组