使用scipy.spatial.Delaunay 三角网的构建

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 plt
import matplotlib.tri as tri
import numpy as np


plt.triplot(points[:,0], points[:,1], delaunay, linewidth=1.5)
plt.show()