[IDL]自动构建泰森多边形（沃罗诺多边形）
发表于 2014-02-14
首先，祝大家元宵节快乐，情人节快乐！
说一下为什么关注到Voronoi Polygon，是因为在ENVI 5.1版本中提供的无缝镶嵌工具可以自动生成接边线，而使用的方法就类似于Voronoi。感兴趣的话可以去下面的博文学习和感受一下此工具的魅力：http://blog.sina.com.cn/s/blog_764b1e9d0101cbf0.html
1. 什么是泰森多边形
下面来科普下泰森多边形（Voronoi），来自百度百科。
左边的图就是泰森多边形，看着是不是有点眼熟？对了，水立方便是基于此原理设计的。
图：左为泰森多边形，右为水立方
long long ago，美国气候学家A·H·Thiessen提出了一种根据离散分布的气象站的降雨量来计算平均降雨量的方法，即将所有相邻气象站连成三角形，作这些三角形各边的垂直平分线，于是每个气象站周围的若干垂直平分线便围成一个多边形。用这个多边形内所包含的一个唯一气象站的降雨强度来表示这个多边形区域内的降雨强度，并称这个多边形为泰森多边形。
泰森多边形的特性是：
1、每个泰森多边形内仅含有一个离散点数据；
2、泰森多边形内的点到相应离散点的距离最近；
3、位于泰森多边形边上的点到其两边的离散点的距离相等。
泰森多边形可用于定性分析、统计分析、邻近分析等。例如，可以用离散点的性质来描述泰森多边形区域的性质；可用离散点的数据来计算泰森多边形区域的数据；判断一个离散点与其它哪些离散点相邻时，可根据泰森多边形直接得出，且若泰森多边形是n边形，则就与n个离散点相邻；当某一数据点落入某一泰森多边形中时，它与相应的离散点最邻近，无需计算距离。
在泰森多边形的构建中，首先要将离散点构成三角网。这种三角网称为Delaunay三角网。
泰森多边形又叫冯洛诺伊图（Voronoi diagram），得名于Georgy Voronoi。
2. 泰森多边形的建立步骤
建立泰森多边形算法的关键是对离散数据点合理地连成三角网，即构建Delaunay三角网。建立泰森多边形的步骤为：
1) 离散点自动构建三角网，即构建Delaunay三角网。对离散点和形成的三角形编号，记录每个三角形是由哪三个离散点构成的。
2) 找出与每个离散点相邻的所有三角形的编号，并记录下来。这只要在已构建的三角网中找出具有一个相同顶点的所有三角形即可。
3) 对与每个离散点相邻的三角形按顺时针或逆时针方向排序，以便下一步连接生成泰森多边形。设离散点为o。找出以o为顶点的一个三角形，设为A；取三角形A除o以外的另一顶点，设为a，则另一个顶点也可找出，即为f；则下一个三角形必然是以of为边的，即为三角形F；三角形F的另一顶点为e，则下一三角形是以oe为边的；如此重复进行，直到回到oa边。
4) 计算每个三角形的外接圆圆心，并记录之。
5) 根据每个离散点的相邻三角形，连接这些相邻三角形的外接圆圆心，即得到泰森多边形。对于三角网边缘的泰森多边形，可作垂直平分线与图廓相交，与图廓一起构成泰森多边形。
3. IDL中建立泰森多边形的方法
从上边可以看到，需要先提供离散点，然后构建Delaunay三角网，最后根据规则构建泰森多边形，用到的IDL函数为TRIANGULATE和VORONOI。分几种情况进行介绍。
注：下图中黑点为离散点位置，可以更好理解泰森多边形的构建过程。
3.1 随机离散点
通过如下代码生成随机离散点，效果图如下：
N = 36
X = RANDOMN(seed, N)
Y = RANDOMN(seed, N) 
图：随机离线点自动构建的泰森多边形
3.2 规则离散点
坐标为正规网格，例如[0, 0]、[0, 1]、[0, 2]、[1, 0]、[1, 1]、[1, 2]… …
使用如下代码生成规则离散点：
X = REBIN(INDGEN(6)+10, 36,1)
Y = REFORM(REBIN(INDGEN(6)+10, 6,6),36,1)
X为：
10 10 10 10 10 10
11 11 11 11 11 11
12 12 12 12 12 12
13 13 13 13 13 13
14 14 14 14 14 14
15 15 15 15 15 15
Y为：
10 11 12 13 14 15
10 11 12 13 14 15
10 11 12 13 14 15
10 11 12 13 14 15
10 11 12 13 14 15
10 11 12 13 14 15 生成的泰森多边形就是这样的规则网格了。
图：规则离散点生成的泰森多边形
3.3 较规则离散点
其实就是给规则离散点加入一些随机噪声，利用Randomu或Randomn可以生成。代码如下：
X = REBIN(INDGEN(6)+10, 36,1)+ RANDOMU(seed, 36,1)
Y = REFORM(REBIN(INDGEN(6)+10, 6,6),36,1)+ RANDOMU(seed, 36,1)
图：较规则离散点构建的泰森多边形（类似水立方）
最后来一张叠加三角网显示效果图，更有利于了解泰森多边形的构建原理：
图：叠加三角网显示效果
附IDL源代码：
PRO testVoronoi
idx = 0
; 创建离散点
CASE idx OF
; 随机离散点
0: BEGIN
N = 36
X = RANDOMN(seed, N)
Y = RANDOMN(seed, N)
END
; 规则离散点
1: BEGIN
X = REBIN(INDGEN(6)+10, 36,1)
Y = REFORM(REBIN(INDGEN(6)+10, 6,6),36,1)
N = N_ELEMENTS(X)
END
; 较规则离散点
2: BEGIN
X = REBIN(INDGEN(6)+10, 36,1)+ RANDOMU(seed, 36,1)
Y = REFORM(REBIN(INDGEN(6)+10, 6,6),36,1)+ RANDOMU(seed, 36,1)
N = N_ELEMENTS(X)
END
ELSE:
ENDCASE
; 构建Delaunay三角网
TRIANGULATE, X, Y, tr, CONN=C
; 绘制离散点
Points = PLOT(x, y, $
LINESTYLE=6, $
SYMBOL='o', SYM_COLOR='black', $
SYM_SIZE=0.5, $
/SYM_FILLED, SYM_FILL_COLOR='black', $
AXIS_STYLE=4, $
MARGIN=[0,0,0,0], $
WINDOW_TITLE='泰森多边形')
FOR I=0, N-1 DO BEGIN
; 获取第i个泰森多边形:
VORONOI, X, Y, I, C, Xp, Yp
; 快速可视化绘制
; 创建Polygon函数中的CONNECTIVITY关键字，指定连通性
n=N_ELEMENTS(Xp)
con = [n,INDGEN(n)]
; 绘制泰森多边形
poly = POLYGON(Xp,Yp,/data,CONNECTIVITY=con, /current, $
FILL_BACKGROUND=1, FILL_COLOR=I*7+5, RGB_TABLE=25, $
COLOR='white')
ENDFOR
; 绘制三角网，如果不想绘制，把这几行注释即可
n_Tr = (SIZE(tr, /DIMENSIONS))[1]
FOR i = 0,n_Tr-1 DO BEGIN
tri = POLYGON(X[tr[*,i]], Y[tr[*,i]], CONNECTIVITY=[3,0,1,2], $
/current, color = 'black', /data, FILL_BACKGROUND=0, $
LINESTYLE=4)
ENDFOR
; 将离散点置顶显示
Points.ORDER, /BRING_TO_FRONT
END