Python实现Delaunay三角剖分与TIN构建详解

晚来风

发布时间：2026-10-15| 3人浏览过

原创

本文详细介绍如何使用Python实现Delaunay三角剖分及不规则三角网（TIN）的构建。文章首先解析Delaunay三角剖分的定义与空圆特性，随后通过生长算法和逐点插入算法对比，展示基于Tkinter的三角剖分实现过程。内容涵盖随机点生成、坐标处理、界面展示及算法优缺点分析，为DEM数据处理提供参考。

关于该技术的实际应用价值，此处暂不展开讨论。

具体实现过程如下。

在可视化结果中，我采用三角形三个顶点到中心点的距离之和作为颜色映射依据。观察结果显示，点分布越密集的区域，对应图片上的颜色越深。


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为所有边的集合。那么该点集V的一个三角剖分T=(V,E)是一个平面图G，该平面图需满足特定条件。

首先，除了端点外，平面图中的边不包含点集中的任何其他点。

其次，图中不存在相交的边。

最后，平面图中所有的面均为三角面，且所有三角面的并集构成散点集V的凸包。

在实际应用中，应用最为广泛的三角剖分是Delaunay三角剖分，它属于一种特殊的三角剖分形式。

【定义】Delaunay边：假设E中的一条边e（其两个端点为a和b），若e满足下列条件，则称其为Delaunay边。

存在一个圆经过a、b两点，且该圆内部（注意是圆内，圆上最多允许三点共圆）不包含点集V中的任何其他点，这一特性被称为空圆特性。

【定义】Delaunay三角剖分：如果点集V的一个三角剖分T仅包含Delaunay边，那么该剖分即称为Delaunay三角剖分。

关于Delaunay三角剖分的具体算法细节，可参考百度百科中“Delaunay三角剖分算法”的相关介绍。

进行三角剖分的主要目的——构建TIN（不规则三角网）。

不规则三角网（TIN）是数字高程模型（DEM）的重要表现形式之一。与规则格网相比，TIN具有数据冗余小、细节丢失少的显著特点。

在分布不规则的高程点之间构建三角网，其核心技术关键在于三角剖分算法。

算法执行步骤。

第一步，任意选取一点，在点集中找出距离该点最近的点，将这两点连成一条线，以此线作为基线。

第二步，在所有剩余点中寻找能与该基线构成具有空圆性三角形的点，并据此构成三角形。

第三步，以新生成的边为基线，重复执行第二步，直至整个点集构网完成。

具体实现代码如下。

本代码使用的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)

运行上述程序后，我们得到一组坐标数据如下。

将上述生成的点在界面中进行展示。