使用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()