简单多边形的三角剖分

1 折线、闭折线、简单闭折线、简单多边形

折线：将平面上一些点的序列连接起来，可以得到一条折线。

如：有点列$A_1,A_2,...,A_n$，则将其顺次连接，得到的一条折线。

点列中的点，称为折线的控制点。

闭折线：如果一条折线的起点和终点相同，则称该折线为闭折线。以连续两控制点为端点的线段，称为折线的边。有一个共同端点的边称为相邻边。

如：顺次连接点列$A_1,A_2,...,A_n,A_1$，得到的一条闭折线。

简单闭折线：一条闭折线，如果其控制点和边满足如下四点，则该闭折线称为简单闭折线。

1、闭折线包含至少三个不同的控制点
2、控制点两两不同
3、连续的三个控制点不共线
4、不相邻的两条边，不相交

简单闭折线的控制点，也称为其顶点。

简单多边形：由简单闭折线围成的图形称为简单多边形。简单多边形按照边的数目，给其命名为简单n边形。
简单多边形，是单连通的。

下列图形都是简单多边形：

图一

说明：若无特殊交待，下述多边形均指简单多边形。

2 多边形的三角剖分

多边形的三角剖分：一个多边形，如果可以表示为一些内部两两不交的三角形的并，并且这些三角形满足如下条件

1、三角形的顶点，是多边形的顶点
2、若两三角形相交，则其要么交于一点，要么交于一边（即恰有一条公共边）。

则称这样的三角形划分方法为该多边形的一个三角剖分。

例如

图二

多边形的分类：多边形可以分为凸多边形和凹多边形。若连接多边形内任意两点的线段也在多边形内，就称这样的多边形为凸多边形；不是凸多边形的多边形称为凹多边形。

剩余图形：超过三个点的多边形，可以擦去一个点，形成一个新的多边形。比如，多边形$ABCDA$，擦去$B$点，形成新的图形$ACDA$，擦去$A$点则形成新图形$BCDB$。如果记原图形$ABCDA$为$\mathbf{G}$，则擦去一点$B$后的图形记为${\mathbf{G}}_{\mathbf{B}}$。

图三

拆出三角形：擦去的一个点，必然能和其相邻的两个点构成一个三角形，这个三角形称为该擦去点的拆出三角形。上图中的$\Delta ABC$就是$B$点的拆出三角形。

凸顶点、凹顶点：多边形$\mathbf{G}$的一个顶点$V$，如果不包含在擦除该顶点后的剩余图形中，即$V\notin {\mathbf{G}}_{\mathbf{V}}$，则称$V$为$\mathbf{G}$的一个凸顶点；否则称为凹顶点。

命题1：若一个多边形的每个顶点都是凸顶点，则该多边形为凸多边形。
三角形是最简单的凸多边形。

命题2：一个凹多边形，至少有一个凹顶点，至少有三个凸顶点。

命题3：一个凸$n$边形（$n>3$），擦去一个点后的得到的剩余图形，为凸$n-1$边形。

可划分顶点：多边形的一个顶点，如果擦去这个顶点得到的剩余图形为多边形，且这个顶点的拆出三角形与其剩余图形内部不交，则称这个顶点为可划分顶点。

命题4：凸多边形的每个顶点都是可划分顶点。凹多边形的凹顶点一定是不可划分顶点。可划分顶点一定是凸顶点。
图二，就是一个凹四边形，其中A为凹顶点，B,C,D为凸顶点，其中B,D为可划分顶点。

  显然，要实现对一个多边形$\mathbf{G}$的三角剖分，只需找到其一个可划分顶点$P$，得到$P$的拆出三角形并记录，对其剩余图形${\mathbf{G}}_{\mathbf{P}}$重复上述过程，就可以得到多边形的一个三角剖分。由于凸多边形的良好性质，即每个顶点都是可划分顶点，剩余图形也是凸多边形，因此对凸多边形的三角剖分是非常简单的。对凹多边形，关键就是要寻找其可划分顶点。

3 可划分顶点的计算机判定¹

  用有序顶点给出一个多边形$\mathbf{G}=A_1,A_2,...,A_n,A_1$，给定一个顶点$A_j$，首先要判定这个顶点是否是凸顶点。
$A_1,A_2,...,A_n,A_1$要么是顺时针的，要么是逆时针的。定义一个方向因子$\delta$，如果是逆时针方向因子为$1$，是顺时针方向因子为$-1$。

from shapely.geometry import LinearRing
LinearRing(顶点序列).is_ccw 
可以判断一个顶点序列是顺时针的还是逆时针的，逆时针返回True，顺时针返回False

定义二维向量的叉积：$(x_1,y_1)\times (x_2,y_2)=x_1{y}_2-y_1{x}_2$

命题5 凸顶点判定的有向面积法：设$A_j$的前驱节点为$A_i$，后继节点为$A_k$，设向量$\vec{x}=A_j-A_i,\vec{y}=A_k-A_j$，那么$\delta (\vec{x}\times \vec{y})>0$则$A_j$为凸顶点；$\delta (\vec{x}\times \vec{y})<0$则$A_j$为凹顶点。

如果一个顶点是凸顶点，那么再根据可划分顶点的定义去判定该顶点是否为可划分顶点。可以利用shapely，构造多边形，然后根据条件进行判定。

完整的python代码实现如下：

import matplotlib.pyplot as plt
import numpy as np
import random
from shapely.geometry import Polygon, LinearRing #多边形模型，和线性环模型

def IsSimplePoly(poly):
    """判断多边形poly是否为简单多边形"""
    poly_ring = poly.boundary
    if poly_ring.is_ring and list(poly.interiors) == []:
        return True
    else:
        return False
def GetPolyVex(poly):
    """得到poly的顶点序列，以numpy数组的形式返回
       :首尾点重合，该数组为n*2的数组
    """
    return np.asarray(poly.exterior)
def VexCCW(poly):
    """判断poly的顶点给出的是顺时针顺序还是逆时针顺序
       :若给出的顶点为逆时针排列则返回1，为顺时针旋转则返回-1
    """
    return 1 if LinearRing(poly.exterior).is_ccw else -1
def GetDivideVexIdx(poly):
    """得到poly中的可划分顶点的下标序列
         :返回1，无重复顶点np数组，顺序与poly.exterior中的顺序相同
         :返回2，其中可划分顶点的下标序列
         :返回3，可划分顶点的多在角的弧度序列
    """
    dividevex_idx_li = [] #存储可划分顶点的下标
    dividevex_arg_li = [] #存储可划分顶点所对应角的弧度值
    vex_arr = GetPolyVex(poly) #顶点序列
    vex_arr = vex_arr[:-1,:] #去掉最后一个回环
    nums = vex_arr.shape[0] #顶点序列的个数
    if nums <= 3: #三角形则不用再处理
        return vex_arr, dividevex_idx_li, dividevex_arg_li
    
    pm = VexCCW(poly) #poly的顺逆时针状态
    for i in range(nums):
        v = vex_arr[i,:]#当前顶点
        l = vex_arr[i-1,:]#前驱顶点
        r = vex_arr[(i+1)%nums,:]#后继顶点
        fir_vector = v - l #用有向面积法计算是否为凸顶点
        sec_vector = r - v
        A = np.array([fir_vector,sec_vector]) #判断矩阵
        if pm*np.linalg.det(A) > 0:#此时的顶点为凸顶点，在此基础上判断其是否为可划分顶点
            remainvex_arr = np.concatenate([vex_arr[:i,:],vex_arr[i+1:,:]],axis=0)
            remain_poly = Polygon(remainvex_arr)
            tri = Polygon([l,v,r])
            if (remain_poly.is_valid
                and remain_poly.intersection(tri).area < 1e-8 #为一个可调整系数
                and poly.equals(remain_poly.union(tri))):#判断一个凸顶点是否为可划分顶点的依据
                
                dividevex_idx_li.append(i) #将可划分的顶点下标压入序列
                #下面计算对应的弧度
                arc = np.arccos(-np.dot(fir_vector,sec_vector)/np.linalg.norm(fir_vector)/np.linalg.norm(sec_vector))
                dividevex_arg_li.append(arc)
    return vex_arr, dividevex_idx_li, dividevex_arg_li
def GetDivTri(poly, tris = []):
    """递归的将多边形，进行三角剖分，每次都以角度最小的可划分顶点为依据"""
    vex_arr, dv_idx_li, dv_arc_li = GetDivideVexIdx(poly)
    nums = vex_arr.shape[0]
    if nums <= 3: #三角形，则直接处理
        tris.append(poly)
        return tris
    idx = dv_idx_li[np.argmin(np.array(dv_arc_li))]#取出最小的一个可划分顶点的下标
    #idx = dv_idx_li[np.random.randint(len(dv_idx_li))]#随机取出一个下标
    v = vex_arr[idx, :]
    l = vex_arr[idx-1, :]
    r = vex_arr[(idx+1)%nums, :]
    tri = Polygon([l,v,r]) #划分出来的三角形
    tris.append(tri) #将这个处理好的三角形压入序列
    #下面为得到新序列，并转化为图形，用于递归
    remain_vex_arr = np.concatenate([vex_arr[:idx,:],vex_arr[idx+1:,:]],axis=0)
    remain_poly = Polygon(remain_vex_arr)
    GetDivTri(remain_poly,tris)
    return tris
def PolyPretreatment(poly_arr):
    """用于对poly_arr进行归一化处理"""
    temp = poly_arr - np.min(poly_arr,axis=0)
    return temp / np.max(temp)
def MinAngle(tri):
    """计算一个三角形的最小角的弧度[0,pi/2]"""
    point = np.asarray(tri.exterior)
    arc_li = []
    for i in range(3):
        j = (i+1)%3; k=(i+2)%3
        a = np.linalg.norm(point[i,:] - point[j,:])
        b = np.linalg.norm(point[j,:] - point[k,:])
        c = np.linalg.norm(point[k,:] - point[i,:])
        arc = np.arccos((a**2 + b**2 - c**2)/(2*a*b))
        arc_li.append(arc)
    return min(arc_li)
def OptDiv(poly4_vex_arr):
    """对四边形进行优化划分，返回其最优化的两个三角形"""
    tri1 = Polygon(poly4_vex_arr[[0,1,2]])
    tri2 = Polygon(poly4_vex_arr[[0,2,3]])
    arc1 = min([MinAngle(tri1),MinAngle(tri2)])

    tri3 = Polygon(poly4_vex_arr[[0,1,3]])
    tri4 = Polygon(poly4_vex_arr[[1,2,3]])
    arc2 = min([MinAngle(tri3),MinAngle(tri4)])

    if arc1 >= arc2:
        return tri1,tri2
    else:
        return tri3,tri4
def OptAlltris(tris):
    """对已经给出的三角剖分进行进一步的优化，使得最小角最大
        :对剖分出的三角形序列进行优化
        :通常需要运行两次，才能保证充分优化
    """
    random.shuffle(tris)
    nums = len(tris)
    for i in range(nums):
        tri_i = tris[i]
        for j in range(i+1,nums):
            tri_j = tris[j]
            if tri_i.intersection(tri_j).length > 1e-10:
                u = tri_i.union(tri_j)
                vex_arr, dv_vex_li, _=GetDivideVexIdx(u)
                if len(dv_vex_li) == 4:
                    a,b = OptDiv(vex_arr)
                    flag = True
                    for idx in set(range(nums)) - {i,j}:
                        if a.intersection(tris[idx]).area > 0. or b.intersection(tris[idx]).area > 0.:
                            flag = False
                    if flag:
                        tris[i],tris[j] = a,b
    return tris
##-----------------------------------------------##
poly_arr = np.array([(0,0),(1,2),(0,4),(3,6),(2,5),(3,5),(4,3),(0,0)]) #顶点序列
poly = Polygon(PolyPretreatment(poly_arr)) #构造多边形
#运算,绘图脚本
if IsSimplePoly(poly):
    plt.figure(figsize=(16,16))
    tris = []
    tris =  GetDivTri(poly,tris = tris)
    #用mpl画出，原来图形的线框
    plt.subplot(2,2,1)
    plt.plot(*poly.exterior.xy)
    plt.axis("equal")
    #用线框画出剖分
    plt.subplot(2,2,2)
    for tri in tris: #triangulate得到的所有三角形，这是对凸包的一个划分
        plt.plot(*tri.exterior.xy)
    plt.axis("equal")
    #用色块画出剖分
    plt.subplot(2,2,3)
    for tri in tris:
        color = np.random.rand(4)
        color[3] *= 0.6#以下两句用于调节透明度
        color[3] += 0.4
        tri = plt.Polygon(np.asarray(tri.exterior),facecolor = color)
        plt.gca().add_patch(tri)
    #进行优化，并用色块画出新的剖分
    newtris = tris.copy()
    newtris = OptAlltris(newtris) # 进行优化
    newtris = OptAlltris(newtris)
    plt.subplot(2,2,4)
    for tri in newtris:
        color = np.random.rand(4)
        color[3] *= 0.6#以下两句用于调节透明度
        color[3] += 0.4
        tri = plt.Polygon(np.asarray(tri.exterior),facecolor = color)
        plt.gca().add_patch(tri)
        

    plt.show()
else:
     print("输入的多边形，不是定义要求的简单多边形！")

运行效果。最后一个图对三角剖分进行了优化，使得在一个凸四边形中，剖分出的最小角较大的一种剖分方法。

说明：在进行这样的剖分之后可以插入中间节点，然后对每个三角形或者其组成的凸区域进行更进一步的Delaunay三角剖分。

---

【参考条目】

1. https://www.cnblogs.com/lan-yt/p/9200621.html
2. https://www.doc88.com/p-661156855371.html²

---

1. 当然这里不是从最底层开始的，而是利用了Python的shapely包来实现。 ↩︎

2. 这篇文献中的粗三角剖分方法应该是错误的，错误的原因是寻找可划分顶点的方法错误，请读者注意。 ↩︎