百度百科
泰森多邊形又叫馮洛諾伊圖（Voronoi diagram），得名於Georgy Voronoi，是由一組由連接兩鄰點線段的垂直平分線組成的連續多邊形組成。
泰森多邊形是對空間平面的一種剖分，其特點是多邊形內的任何位置離該多邊形的樣點（如居民點）的距離最近，離相鄰多邊形內樣點的距離遠，且每個多邊形內含且僅包含一個樣點。由於泰森多邊形在空間剖分上的等分性特征，因此可用於解決最近點、最小封閉圓等問題，以及許多空間分析問題，如鄰接、接近度和可達性分析等。
特征：
- 每個泰森多邊形內僅含有一個離散點數據；
- 泰森多邊形內的點到相應離散點的距離最近；
- 位於泰森多邊形邊上的點到其兩邊的離散點的距離相等。
泰森多邊形圖例：
詳細介紹見百度百科泰森多邊形
算法實現
算法一：
算法二：
算法二是基於算法一的優化。本文將着重介紹算法二。
1）算法二同樣采用特征點以均勻的速度向外擴張的方式進行。既然我們速度一定，那我們不妨設置為1，那么諸多特征點同時同步地向外擴張，說明在相遇時，相遇的特征點是本着相同的速度，以相同的時間到達的相遇地點（在這里是相遇的單元格），那么根據大眾熟知的速度路程公式s=v*t，我們知道這兩個相遇的特征點距離該相遇點的路程是相等的，也就是距離一樣，說明該相遇點是這兩個特征點兩線的中點。這就符合了泰森多邊形的定義（是由一組由連接兩鄰點線段的垂直平分線組成的連續多邊形組成）
2）擴張的速度我們假定為1，基准的核心點，我們設定為特征點所在像素單元的幾何中心。我們規定，在以該特征點為圓心的輻射區域內的其他像素單元格的幾何中心距離該特征點所在的幾何中心的直線距離（根據勾股定理計算即可）小於或者等於某時間點該特征點以速度為1外擴的半徑長度（即路程）時（也就是該像素單元格的大部分面積都包含在該時間點特征點輻射半徑所划的圓中的時候），將該像素單元格賦值為特征點外擴運動中此時的時間點信息。
計算中心點：（其中紅色點為核心特征點，黑色點為其外擴過程中的一個示例點）
模擬泰森多邊形由核心特征點外擴：（單元格中數值代表外擴到該單元格所用的時間）
3）以步驟2）中介紹的方式，我們對所有特征點進行相應外擴，假設每個特征點都輻射滿整張矩形區域，共輻射出與特征點數量等量的矩形表數量。（其實可以在代碼中把這些特征點輻射的不同數值放在一個數組集合中）。
以其中一個特征點輻射滿全圖為例：（圖二是用數值填充滿整個矩陣區域后，輸出excel，在excel中使用查找替換將相同的數值單元附成相同顏色的背景色）
（注意輸出excel時，要輸出為“.csv”格式而不是“.xsl”和“.xslx”格式，因為后兩種格式有輸出行列不超過256的限制，如果你使用矩陣行列小於256，那么就可以使用后兩種格式了）
（圖一，選中單元即為核心單元點）
（圖二，相關查找替換后的模擬圖表）
4）依據步驟3）中的操作，我們將所有特征值的輻射全圖都計算出后，現在我們依據步驟1）中的理論，當所有的輻射全圖中的某個相同單元格中，所有的輻射全圖中該單元格賦值有兩個並且是最小值相同時，那我們就把它確定為泰森多邊形的邊界線。
為什么是最小呢？
因為如果兩個最近的特征點外延先相遇了，那么后期有人再在此處相遇，對於我們而言就沒有什么意義了。
5）現在我們詳細說一下步驟2）的某一個特征點擴充到全圖的具體操作，回到前面看算法一中的描述，我們要在某特征值的周圍進行無窮無盡的試探，以及對於邊界單元兩邊的各特征點外擴部隊彼此試探，互通信息困難的問題，我們對算法一進行了優化，也就是算法二，為了提高效率，也本着一次遍歷運算完所有單元的原則，我們使用步驟2）中的理論也就是大眾熟知的“勾股定理”，對整個矩陣的單元格進行遍歷（事先我們不是已經知道特征點的坐標了么，也就是矩陣中的行列號）。對於計算出來的當前單元格中心到特征點中心的直線距離，我們進行一下處理，如果是整數，我們就直接將此數附與該單元格；如果有小數部分，例如3.4，也就是說時間為3的時候，輻射區域還未到此單元格，時間為4的時候已經超越了該單元格中心點，那么我們就為該單元格附上4，也就是對於小數取整加一。
代碼實例（python）：
小編拿到的初始特征點坐標文件為“xsl”文件，如下：
共29個特征點。
我們用到的矩陣大小為533*401。
我們的目標是對這些特征點進行構造泰森多邊形，輸出該泰森多邊形的excel表和圖像文件。
我們的大體流程就是：
- 先將這些特征值分別計算出自身輻射全圖，分別為每個矩陣單元附上自身輻射到此處的時間
- 對比某個特定單元格中的29個值，如果其中有兩個值相等並且是在這29個值中是最小的，那么我們就設定該單元格為邊界區域，附上特定值
- 將上步結果附在一個新矩陣中，用此矩陣進行excel文件和圖像的輸出
#! /usr/bin/env python #coding=utf-8 #************************************輸出圖像和excel表 # pyexcel_xls 以 OrderedDict 結構處理數據 from collections import OrderedDict from pyexcel_xls import get_data from pyexcel_xls import save_data import cmath import numpy as np from PIL import Image import math def read_xls_file(): xls_data = get_data(r"C:\Users\123\Desktop\Voronoi\result.xlsx") # print("Get data type:", type(xls_data)) #for sheet_n in xls_data.keys(): #print(sheet_n, ":", xls_data[sheet_n]) # print(len(xls_data["Sheet1"])) # print(xls_data["Sheet1"][0]) # print(len(xls_data["Sheet1"][0])) # print(type(xls_data["Sheet1"][0][0])) # print(xls_data["Sheet1"][0][0]) # 獲取特征點行列號 VList=[] for i in range(len(xls_data["Sheet1"])): if i==0: continue; list1=[] for j in range(len(xls_data["Sheet1"][i])): if j==0: list1.append(xls_data["Sheet1"][i][j]) if j==1: list1.append(xls_data["Sheet1"][i][j]) VList.append(list1) break; # 打印獲取的行列號 # for i in range(len(VList)): # for j in range(len(VList[i])): # print(VList[i][j],end="\t") # print() # 建立柵格矩陣 VMatrix=[None]*533 for i in range(len(VMatrix)): VMatrix[i] = [0]*401 for j in range(401): VMatrix[i][j]=[0]*30 print(type(VMatrix[1][1])) # print(VMatrix) # for i in range(533): # for j in range(401): # VMatrix[i][j].append(0); # print(VMatrix[i][j],end="\t") # print() # 為柵格矩陣附上特征點 sum=0; for i in range(len(VList)): VMatrix[VList[i][0]][VList[i][1]][0]=66666 sum=sum+1 print(sum) #打印一下，共29個點 # VMatrix[col][row+v][1]=v;#右 # VMatrix[col][row-v][1]=v;#左 # VMatrix[col+v][row][1]=v;#下 # VMatrix[col-v][row][1]=v;#上 # 在該矩形區域內進行遍歷計算 # for i in range(col-1,col+1): # for j in range(row-1,row+1): # Value2=(i-col)*(i-col)+(j-row)*(j-row) # if i==col&j==row: # continue; # elif abs(cmath.sqrt(Value2))<=v: # VMatrix[i][j][1]=v; # 所謂相同速度不同時間，不就是路程么？計算每個單元距離核心單元的具體距離，以速度為1算。 # 若為整數則保留，若有余數，則舍掉，加一 for k in range(len(VList)): col=VList[k][0] row=VList[k][1] for i in range(len(VMatrix)): for j in range(len(VMatrix[i])): if i==col&j==row: continue; Value2=(i-col)*(i-col)+(j-row)*(j-row) distance=cmath.sqrt(abs(Value2)).real if isinstance(distance,int): VMatrix[i][j][k+1]=distance; else: VMatrix[i][j][k+1]=math.ceil(distance); # 單獨提取三維列表中的第三維的30個變量中的一個變量 # 該30個變量中，第一個表示特征點，余下29個為29個特征點以自己為核心勻速向外擴展的層次矩陣 # 在某一個單元內，29個特征點流過的路程值中，如果有最小值不止一個，就說明此處為邊界上的一個點 SingleVM=[None]*533 for i in range(len(SingleVM)): SingleVM[i] = [255]*401 for j in range(len(SingleVM)): for k1 in range(len(SingleVM[j])): count=0 #29個數中最小值的計數 minV=999999 for k2 in range(len(VMatrix[j][k1])): if k2==0: continue;#要把第一個元素跳過去，因為第一個元素存儲的是特征值 # print(VMatrix[j][k1][k3],end="\t") if VMatrix[j][k1][k2]<minV: minV=VMatrix[j][k1][k2] count=1 elif VMatrix[j][k1][k2]==minV: count=count+1 if count>=2: SingleVM[j][k1]=0;#邊界矩陣 # print() for i in range(len(VList)): SingleVM[VList[i][0]][VList[i][1]]=0 # 保存到csv文件 data = OrderedDict() # 添加sheet表 data.update({u"Sheet1": SingleVM}) #打印矩陣 # 保存成xls文件 有256行的限制 # 故存成csv文件 addr=r"C:\Users\123\Desktop\Voronoi\ResultBorder2.csv" save_data(addr, data) #保存成圖像輸出 data_list=np.asarray(SingleVM) pic=Image.fromarray(data_list.astype(np.uint8)) pic.save(r"C:\Users\123\Desktop\Voronoi\Voronoi.tif") if __name__ == '__main__': read_xls_file()
結果：
1）excel文件：（由於數值過多，此處只附上將邊界單元格替換成紅色背景后的一部分表格）
2)圖像文件：
驗證代碼的正確性：
下面代碼，大家可以替換掉開始的VList數組進行運行生成新的泰森多邊形。
#! /usr/bin/env python #coding=utf-8 # pyexcel_xls 以 OrderedDict 結構處理數據 # 隨便點幾點試試 from collections import OrderedDict from pyexcel_xls import get_data from pyexcel_xls import save_data import cmath import numpy as np from PIL import Image import math def read_xls_file(): VList=[[300,200],[100,200],[350,90],[60,80],[150,400]] # 打印獲取的行列號 # for i in range(len(VList)): # for j in range(len(VList[i])): # print(VList[i][j],end="\t") # print() # 建立柵格矩陣 VMatrix=[None]*533 for i in range(len(VMatrix)): VMatrix[i] = [0]*401 for j in range(401): VMatrix[i][j]=[0]*6 print(type(VMatrix[1][1])) # print(VMatrix) # for i in range(533): # for j in range(401): # VMatrix[i][j].append(0); # print(VMatrix[i][j],end="\t") # print() # 為柵格矩陣附上特征點 sum=0; for i in range(len(VList)): VMatrix[VList[i][0]][VList[i][1]][0]=0 sum=sum+1 print(sum) #打印一下，共29個點 # VMatrix[col][row+v][1]=v;#右 # VMatrix[col][row-v][1]=v;#左 # VMatrix[col+v][row][1]=v;#下 # VMatrix[col-v][row][1]=v;#上 # 在該矩形區域內進行遍歷計算 # for i in range(col-1,col+1): # for j in range(row-1,row+1): # Value2=(i-col)*(i-col)+(j-row)*(j-row) # if i==col&j==row: # continue; # elif abs(cmath.sqrt(Value2))<=v: # VMatrix[i][j][1]=v; # 所謂相同速度不同時間，不就是路程么？計算每個單元距離核心單元的具體距離，以速度為1算。 # 若為整數則保留，若有余數，則舍掉，加一 for k in range(len(VList)): col=VList[k][0] row=VList[k][1] for i in range(len(VMatrix)): for j in range(len(VMatrix[i])): if i==col&j==row: continue; Value2=(i-col)*(i-col)+(j-row)*(j-row) distance=cmath.sqrt(abs(Value2)).real if isinstance(distance,int): VMatrix[i][j][k+1]=distance; else: VMatrix[i][j][k+1]=math.ceil(distance); # 單獨提取三維列表中的第三維的30個變量中的一個變量 # 該30個變量中，第一個表示特征點，余下29個為29個特征點以自己為核心勻速向外擴展的層次矩陣 # 在某一個單元內，29個特征點流過的路程值中，如果有最小值不止一個，就說明此處為邊界上的一個點 SingleVM=[None]*533 for i in range(len(SingleVM)): SingleVM[i] = [255]*401 for i in range(len(VList)): SingleVM[VList[i][0]][VList[i][1]]=0 for j in range(len(SingleVM)): for k1 in range(len(SingleVM[j])): count=0 #29個數中最小值的計數 minV=999999 for k2 in range(len(VMatrix[j][k1])): if k2==0: continue;#要把第一個元素跳過去，因為第一個元素存儲的是特征值 # print(VMatrix[j][k1][k3],end="\t") if VMatrix[j][k1][k2]<minV: minV=VMatrix[j][k1][k2] count=1 elif VMatrix[j][k1][k2]==minV: count=count+1 if count>=2: SingleVM[j][k1]=0;#邊界矩陣 # print() # 保存到csv文件 data = OrderedDict() # 添加sheet表 data.update({u"Sheet1": SingleVM}) #打印矩陣 # 保存成xls文件 有256行的限制 # 故存成csv文件 addr=r"C:\Users\123\Desktop\Voronoi\zzshishiBorder2.csv" save_data(addr, data) data_list=np.asarray(SingleVM) pic=Image.fromarray(data_list.astype(np.uint8)) pic.save(r"C:\Users\123\Desktop\Voronoi\zzshishiVoronoi.tif") if __name__ == '__main__': read_xls_file()
本測試代碼生成結果：
