python 實現德洛內三角剖分的操作

Posted on by

我也不知道這玩意主要是幹啥用的

實現如下

我用剖分的三角形的三個頂點到中心點的距離和作為顏色， 結果顯示： 點越密集的地方， 圖片上的顏色越深。

from scipy.spatial import Delaunay
import numpy as np
import matplotlib.pyplot as plt
width = 80
height = 40 
pointNumber = 50
points = np.zeros((pointNumber, 2)) 
points[:, 0] = np.random.randint(0, width, pointNumber) 
points[:, 1] = np.random.randint(0, height, pointNumber)
tri = Delaunay(points)
center = np.sum(points[tri.simplices], axis=1)/3.0 
'''
color = []
for sim in points[tri.simplices]:
    x1, y1 = sim[0][0], sim[0][1]
    x2, y2 = sim[1][0], sim[1][1]
    x3, y3 = sim[2][0], sim[2][1]
    
    s = ((x1-x2)**2+(y1-y2)**2)**0.5 + ((x1-x3)**2+(y1-y3)**2)**0.5 + ((x3-x2)**2+(y3-y2)**2)**0.5
    color.append(s)
color = np.array(color)
'''
color = []
for index, sim in enumerate(points[tri.simplices]):
    cx, cy = center[index][0], center[index][1]
    x1, y1 = sim[0][0], sim[0][1]
    x2, y2 = sim[1][0], sim[1][1]
    x3, y3 = sim[2][0], sim[2][1]
    
    s = ((x1-cx)**2+(y1-cy)**2)**0.5 + ((cx-x3)**2+(cy-y3)**2)**0.5 + ((cx-x2)**2+(cy-y2)**2)**0.5
    color.append(s)
color = np.array(color)
plt.figure(figsize=(20, 10)) 
plt.tripcolor(points[:, 0], points[:, 1], tri.simplices.copy(), facecolors=color, edgecolors='k') 
plt.tick_params(labelbottom='off', labelleft='off', left='off', right='off', bottom='off', top='off') 
ax = plt.gca() 
plt.scatter(points[:,0],points[:,1], color='r')
#plt.grid()
plt.savefig('Delaunay.png', transparent=True, dpi=600)

補充：生長算法實現點集的三角剖分（ Python（Tkinter模塊））

關於三角剖分

假設V是二維實數域上的有限點集，邊e是由點集中的點作為端點構成的封閉線段, E為e的集合。那麼該點集V的一個三角剖分T=(V,E)是一個平面圖G，該平面圖滿足條件：

1.除瞭端點，平面圖中的邊不包含點集中的任何點。

2.沒有相交邊。

3.平面圖中所有的面都是三角面，且所有三角面的合集是散點集V的凸包。

在實際中運用的最多的三角剖分是Delaunay三角剖分，它是一種特殊的三角剖分。

【定義】Delaunay邊：假設E中的一條邊e（兩個端點為a,b），e若滿足下列條件，則稱之為Delaunay邊：

存在一個圓經過a,b兩點，圓內(註意是圓內，圓上最多三點共圓)不含點集V中任何其他的點，這一特性又稱空圓特性。

【定義】Delaunay三角剖分：如果點集V的一個三角剖分T隻包含Delaunay邊，那麼該三角剖分稱為Delaunay三角剖分。

關於Delaunay三角剖分算法可以參考百度百科Delaunay三角剖分算法

我做三角剖分的目的——構建TIN，不規則三角網

不規則三角網（TIN）是DEM的重要形式之一，相較於規則格網，其具有數據冗餘小、細節丟失少的特點。

在分佈不規則的高程點之間構建出三角網，其關鍵技術就是三角剖分

算法步驟

1、首先任選一點，在點集中找出距離改點最近的點連成一條線，以該線為基線。

2、在所有點中尋找能與該基線構成具有空圓性三角形的點，並構成三角形。

3、以新生成的邊為基線，重復第二步，直至點集構網完成。

具體代碼如下

所使用的python版本為python3.6，編輯器為Pycharm2018.3.1

#-*- coding:utf-8 -*-
import tkinter
from tkinter import filedialog
import csv
#根據兩點坐標計算距離
def caldis(x1,y1,x2,y2):
    return ((x1-x2)**2+(y1-y2)**2)**0.5
#輸入三角形三個頂點，計算外接圓圓心及半徑
def calcenter(x1,y1,x2,y2,x3,y3):
    y1=-y1  #計算公式是根據平面直角坐標推算的，原點在左下角，但是計算機屏幕坐標原點在右上角，所以計算式y坐標取負
    y2=-y2
    y3=-y3
    if (y1 != y3 and y1 != y2 and y2 != y3): #判斷是否有y坐標相等，即三角形某邊斜率為0的情況，避免出現墳分母為0的錯誤
        if(((x3-x1)/(y3-y1))-((x2-x1)/(y2-y1)))==0:
            x2=x2+1
        x=(((y1+y3)/2)+((x1+x3)/2)*((x3-x1)/(y3-y1))-((y1+y2)/2)-((x1+x2)/2)*((x2-x1)/(y2-y1)))/(((x3-x1)/(y3-y1))-((x2-x1)/(y2-y1)))
        y=-((x3-x1)/(y3-y1))*x+((y1+y3)/2)+(((x1+x3)/2)*((x3-x1)/(y3-y1)))
        return (x, -y, caldis(x, y, x1, y1))
    elif (y1 == y3 and y1 != y2 and y2 != y3):#若存在斜率為0的邊則計算可簡化
        x=(x1+x3)/2
        y=-((x2-x1)/(y2-y1))*x+((y1+y2)/2)+((x2-x1)/(y2-y1))*((x1+x2)/2)
        return (x, -y, caldis(x, y, x1, y1)) #返回值為元組（圓心橫坐標x，圓心縱坐標y，外接圓半徑r），計算出來的y值要返回屏幕坐標所以再次取負
    elif (y1 != y3 and y1 == y2 and y2 != y3):
        x = (x1 + x2) / 2
        y = -((x3 - x1) / (y3 - y1)) * x + ((y1 + y3) / 2) + ((x3 - x1) / (y3 - y1)) * ((x1 + x3) / 2)
        return (x, -y, caldis(x, y, x1, y1))
    elif (y1 != y3 and y1 != y2 and y2 == y3):
        x = (x3 + x2) / 2
        y = -((x3 - x1) / (y3 - y1)) * x + ((y1 + y3) / 2) + ((x3 - x1) / (y3 - y1)) * ((x1 + x3) / 2)
        return (x, -y, caldis(x, y, x1, y1))
    else:
        return None
class getTIN: #定義窗口及操作類
    def __init__(self):
        self.path=str() #坐標文件路徑
        self.pointlist=[] #存放所有點坐標的列表
        self.linelist=[] #存放線的列表，每條線用兩個點號表示連線
        self.tk=tkinter.Tk() #定義主窗口
        self.tk.title('MyTIN')
        self.tk.geometry('1200x720')
        self.shengzhang=tkinter.Button(self.tk,text='生長算法',width=15,command=self.drawTIN_shengzhang)
        self.shengzhang.place(x=1050,y=100)  #定義按鈕，關聯到生長算法計算TIN的的函數
        self.readin=tkinter.Button(self.tk,text='讀入坐標文件',width=15,command=self.getfile)
        self.readin.place(x=1050,y=50)
        self.can=tkinter.Canvas(self.tk,width=950,height=620,bg='white')
        self.can.place(x=50,y=50)
        self.tk.mainloop()
    def getfile(self):  #選擇坐標文件（*.csv），從文件中讀入坐標存入pointlist列表並在繪圖區展示出來
        self.path=filedialog.askopenfilename()
        f=open(self.path,'r')
        fd=csv.reader(f)
        self.pointlist=list(fd)
        for i in range(0,len(self.pointlist)):
            self.can.create_oval(int(self.pointlist[i][0])-2,int(self.pointlist[i][1])-2,int(self.pointlist[i][0])+2,int(self.pointlist[i][1])+2,fill='black')
            self.can.create_text(int(self.pointlist[i][0])+7,int(self.pointlist[i][1])-7,text=str(i))
    def drawTIN_shengzhang(self):
        j = 1
        k = 0
        mindis = ((int(self.pointlist[0][0]) - int(self.pointlist[1][0])) ** 2 + (int(self.pointlist[0][1]) - int(self.pointlist[1][1])) ** 2) ** 0.5
        x = len(self.pointlist)
        for i in range(1, x):
            dis = ((int(self.pointlist[0][0]) - int(self.pointlist[i][0])) ** 2 + (int(self.pointlist[0][1]) - int(self.pointlist[i][1])) ** 2) ** 0.5
            if dis < mindis:
                mindis = dis
                j = i
        self.linelist.append((k,j)) #首先計算出距起始點（點號為0）距離最短的點，以這兩點的連線作為基線開始生長
        self.shengzhangjixian(k,j)
    def drawTIN(self): #根據線文件在繪圖區繪制出TIN
        for i in self.linelist:
            self.can.create_line(self.pointlist[i[0]][0], self.pointlist[i[0]][1], self.pointlist[i[1]][0], self.pointlist[i[1]][1])
    def shengzhangjixian(self,i,j): #根據某一基線開始生長的函數
        x = len(self.pointlist)
        for k in range(0,x): #遍歷沒一個點，判斷是否與基線構成D三角形
            n = 0 #n用於統計外接圓內的點數
            if ((k,i) not in self.linelist) and ((i,k) not in self.linelist) and ((j,k) not in self.linelist) and ((k,j) not in self.linelist):
                for y in range(0,x): #遍歷每一個點，判斷
                    if y==i or y==j or y==k:
                        continue
                    if(calcenter(int(self.pointlist[i][0]),int(self.pointlist[i][1]),int(self.pointlist[j][0]),int(self.pointlist[j][1]),int(self.pointlist[k][0]),int(self.pointlist[k][1]))==None):
                        continue
                    else:
                        xyr=calcenter(int(self.pointlist[i][0]),int(self.pointlist[i][1]),int(self.pointlist[j][0]),int(self.pointlist[j][1]),int(self.pointlist[k][0]),int(self.pointlist[k][1]))
                    if caldis(xyr[0],xyr[1],int(self.pointlist[y][0]),int(self.pointlist[y][1])) < xyr[2]: #判斷點是否在外接圓內
                        n=n+1
                    else:
                        continue
            else:continue
            if n == 0: #判斷是否為D三角形
                self.linelist.append((k,i)) #將新生成的邊的端點號加入線列表
                self.drawTIN() #調用繪制函數繪制TIN
                self.shengzhangjixian(k,i) #以生成的新邊作為基線，迭代計算
                self.linelist.append((k,j))
                self.drawTIN()
                self.shengzhangjixian(k,j)
            else:continue
if __name__ == '__main__':
    MyTIN=getTIN()

通過如下代碼生成一組隨機的點並存入文件

import random
import csv
from tkinter import filedialog
path=filedialog.askopenfilename()
OutAddress=open(path,'a',newline='')
csv_write = csv.writer(OutAddress,dialect='excel')
for i in range(0,20):
    x=random.randint(30,920)
    y=random.randint(30,590)
    out=(x,y)
    print(out)
    csv_write.writerow(out)

通過上面的程序我們得到一組坐標如下

550,432
81,334
517,265
842,408
369,123
502,169
271,425
213,482
588,248
94,295
344,350
500,385
912,527
801,491
838,455
104,479
760,160
706,77
227,314
764,576

將以上的點在界面中展示出來

點擊生長算法運行得到結果

小結

生長算法在三角剖分算法中並不是什麼高效的算法，其特點在於算法簡單易行，但是計算量大，並且對於新插入的點無法更新，必須重新計算。

相比之下，逐點插入算法雖然計算量仍然較大（似乎三角剖分計算量都不小），但是能實現對新插入點的更新而不用重頭計算。

註：文中部分圖片及介紹來自百度百科。

以上為個人經驗，希望能給大傢一個參考，也希望大傢多多支持WalkonNet。如有錯誤或未考慮完全的地方，望不吝賜教。

推薦閱讀：

    None Found