1. 前言

　　众所周知，Delaunay三角剖分算法在测绘工程中有着重要的应用。如果你是使用C#，Java之流的编程语言，不好意思你得自己去实现这个算法。好在python有着非常强大的开源库，python+numpy+scipy+matplotlib可谓科学计算”黄金搭档“，当然诸如pandas之流的高性能数据分析库，掌握他们你就不必重复造轮子了。

2. 进入正题

　　这篇随笔主要介绍如何使用Python+scipy+numpy+matplotlib来演示Delaunay三角网的构建并最终以图形显示出来。

• 　　首先先贴一下主要的数据结构代码：

class Vector2d():    def __init__(self, x, y, name=""):        self.name = name        self.x = x        self.y = y        self.rnx_path = None    # 重载运算符    def __eq__(self, other):        if isinstance(other, types.NoneType) or not isinstance(other, Vector2d):            raise ValueError('Expected type of: %s' % (type(Vector2d)))        else:            return (other.x == self.x) and (other.y == self.y)    def __add__(self, other):        if isinstance(other, types.NoneType) or not isinstance(other, Vector2d):            raise ValueError('Expected type of: %s' % (type(Vector2d)))        else:            return Vector2d(self.x + other.x, self.y + other.y)    def __sub__(self, other):        if isinstance(other, types.NoneType) or not isinstance(other, Vector2d):            raise ValueError('Expected type of: %s' % (type(Vector2d)))        else:            return Vector2d(self.x - other.x, self.y - other.y)    # 将乘法重载为叉积    def __mul__(self, other):        if isinstance(other, types.NoneType) or not isinstance(other, Vector2d):            raise ValueError('Expected type of: %s' % (type(Vector2d)))        else:            return other.x * self.y - self.x * other.y    @property    def length(self):        return math.sqrt(self.x ** 2 + self.y ** 2) def __str__(self):        return '(%s,%s,%s)' % (self.name, self.x, self.y)

　　这个类的主要作用就是来表示二维向量，你也可以理解为二维点。在实现主要是重载了几个操作符。

class Triangle():    '''    形参point的类型均为 Vector2d    '''    def __init__(self, point1, point2, point3):        self.p1 = point1        self.p2 = point2        self.p3 = point3    def __str__(self):        return '%s->%s->%s' % (self.p1.name, self.p2.name, self.p3.name)    def is_in_triangle(self, point):        if isinstance(point, types.NoneType) or not isinstance(point, Vector2d):            raise ValueError('Expected type of: %s' % (type(Vector2d)))        else:            o2a = self.p1 - point            o2b = self.p2 - point            o2c = self.p3 - point            return ((o2a * o2b > 0 and o2b * o2c > 0 and o2c * o2a > 0) or (o2a * o2b < 0 and o2b * o2c < 0 and o2c * o2a < 0))

　　这个类主要是用来表示三角单元的，不用多说。

• 　　下面就是关键的生成三角网部分的代码，首先你必须先导入Delaunay

from scipy.spatial import Delaunay

def triangulate(vertex):    '''    Get delauney triangles from points    :param vertex: list of Vector2d    :return: Delauney triangles    '''    triangles = []    delauney = Delaunay([[pt.x, pt.y] for pt in vertex]).simplices.copy()    for tri in delauney:        triangles.append(Triangle(vertex[tri[0]], vertex[tri[1]], vertex[tri[2]]))    return triangles，delauney

　　triangulate方法的形参为Vector2d类的列表类型。 方法返回了所有的三角单元和delaunay对象。delaunay在使用matplotlib绘图的时候需要用到。

• 　　绘图

points = np.array([[pt.x, pt.y] for pt in vertex])

import matplotlib.pyplot as pltimport matplotlib.tri as triimport numpy as npplt.triplot(points[:,0], points[:,1], delaunay, linewidth=1.5) plt.show()