计算几何实验（一）：简单多边形的三角剖分算法

njw21377128

866人浏览 · 2023-10-11 19:34:17

njw21377128 · 2023-10-11 19:34:17 发布

预备知识

概念主要参考了Discrete and Computational Geometry 这本书。每一个概念我都会给出原书中的陈述以及我自己翻译的版本。

简单多边形：多边形 $P$ 是平面上的封闭区域，以形成不与自身相交的封闭曲线的有限线段集为边界。

simple polygon: A polygon P is the closed region of the plane bounded by a finite collection of line segments forming a closed curve that does not intersect itself.

多边形的边界：多边形的顶点和边的集合是多边形的边界，记作 $\partial P$

boundary of polygon： The set of vertices and edges of P is called the boundary of the polygon, denoted as ∂P.

多边形约当曲线定理：一个多边形的边界把平面分成两部分。特别地, $R^2\backslash \partial P$ 的两个部分是有界的内部和无界的外部。

Polygonal Jordan Curve: The boundary ∂P of a polygon P partitions the plane into two parts. In particular, the two components of $R^2\backslash \partial P$ are the bounded interior and the unbounded exterior.

对角线：多边形的对角线是一条连接多边形 $P$ 两个顶点且完全在 $P$ 的内部，不与多边形的边界 $\partial P$ 相交除了在它的两个端点处的线段。

diagonal： A diagonal of a polygon is a line segment connecting two vertices of P and lying in the interior of P, not touching ∂P except at its endpoints.

三角剖分：多边形 $P$ 的三角剖分是把 $P$ 分解成三角形，要求分解时添加的对角线不相交且数量最多。

triangulation: A triangulation of a polygon P is a decomposition of P into triangles by a maximal set of noncrossing diagonals.

多边形的耳朵：三个连续的顶点 $a,b,c$ 形成一个多边形的耳朵，如果 $ac$ 是该多边形的对角线。

Ears: Three consecutive vertices a, b, c form an ear of a polygon if ac is a diagonal of the polygon.

双耳定理： 每一个具有超过三个顶点的多边形有至少两个耳朵。

Every polygon with more than three vertices has at least two ears.

算法思路分析

考虑到多边形的耳朵的特殊性，我们可以通过对多边形一直切耳朵来实现三角剖分。

多边形切去一个耳朵仍然是多边形，所以该算法可以写成递归算法。

除了三角形外每个多边形都有至少两个耳朵，所以递归算法的出口就是三角形。

所以这里的重点是如何识别多边形的耳朵，按照定义，我们只要枚举所有连续的三个顶点 $a,b,c$ 判断两个不相邻的顶点连接而成的线段 $ac$ 是不是对角线就可以了。

两顶点连线 $ac$ 有四种可能的情况（以下情况都省略掉 $ac$ 的端点）：

完全在多边形内部，是对角线
完全在多边形外部，不是对角线
部分在多边形内部，部分在多边形外部，不是对角线（包括与多边形的边界相交于顶点处的情况）
与多边形的边界完全重合（这部分是新增的，对应于 $b$ 是平点，即在该点处内角为180度）

第3种情况，可以通过判断 $ac$ 是否与多边形的边界相交（不包括 $ac$ 的端点）来判别。

如果 $ac$ 不与多边形的边界相交，那么可以随意取 $ac$ 上一点（不包括端点，下同），通过判断该点是否在多边形的内部来判别第1种和第2种情况，如果在多边形内部，则整条线段都在多边形内部，如果在多边形外部，则整条线段都在多边形外部。

第4种情况，只需要单独判断是否 $b$ 为平点即可，只需要判断 $ab$ 和 $bc$ 是否平行。注意，当多边形为带平点的四边形时，需要额外考虑算法实现中的细节。例如：

如果只判断 $ab$ 和 $bc$ 是否平行是不足以判断 $ac$ 是否是对角线的，还需要判断 $ad$ 和 $dc$ 是否平行。

所以现在需要实现的算法主要有

线段的相交判定算法（包括端点）
一个点是否在多边形的内部的判定算法（不需要判断该点是否在多边形的边界，因为 $ac$ 不与多边形的边界相交）
针对第4种情况的线段平行判定算法

这里算法实现主要参考了《算法设计与分析》第2版李春葆中的计算几何部分。

判定点是否在多边形内部的算法大致思路：

过该点作一条射线，如果射线与多边形相交的点的个数是奇数，那它在多边形内部，如果是偶数，则在多边形外部。为了方便，我们一般使用水平向右的射线。

这里需要考虑一些算法实现方面的细节。

在判断射线与多边形的相交点的个数时，我们肯定会对多边形的每一条边都进行判断。

如果射线与多边形交于顶点，如何判断交点的个数？

有两种情况：

对于第1种情况，两条邻接边记成有偶数个交点即可

对于第2种情况，两条邻接边记成有奇数个交点即可

还有可能射线与多边形的水平边相交，如何判断？

此时射线必然与另外两条边的端点相交，忽视掉水平边（即不计入点的个数），然后转化成之前的情况即可。

《算法设计与分析》中的判断点是否在多边形内部的算法实现了同样的效果。

算法代码

# author:njw

import random
import time
import matplotlib.pyplot as plt


# 点类
class Point:
    def __init__(self,x,y):
        self.x=x
        self.y=y
# 线段类
class LineSegment:
    def __init__(self,point1,point2):
        self.point1=point1
        self.point2=point2
#   求两个向量的点积
def DotProduct(point1,point2):
    return point1.x*point2.x+point1.y*point2.y

# 求两个向量的叉积
def CrossProduct(point1,point2):
    return point1.x*point2.y-point1.y*point2.x

# 判断两条线段是否平行(p1p2,p2p3)
def IsParallel(point1,point2,point3):
    return CrossProduct(Point(point2.x-point1.x,point2.y-point1.y),Point(point3.x-point2.x,point3.y-point2.y))==0

# 判断两线段的相对方向
def Direction(point1,point2,point3):
    d = CrossProduct(Point(point3.x-point1.x,point3.y-point1.y),Point(point2.x-point1.x,point2.y-point1.y))
    if d>0:
        return 1
    elif d<0:
        return -1
    else:
        return 0

# 判断一个点是否在线段上
def OnSegment(point1,point2,point3):
    return CrossProduct(Point(point2.x-point1.x,point2.y-point1.y),Point(point3.x-point1.x,point3.y-point1.y))==0 and \
    DotProduct(Point(point2.x-point1.x,point2.y-point1.y),Point(point3.x-point1.x,point3.y-point1.y))<=0

# TODO:判断两线段是否相交
# 判断两线段是否相交
def IsIntersect(lineSegment1,lineSegment2):
    d1 = Direction(lineSegment2.point1,lineSegment2.point2,lineSegment1.point1)
    d2 = Direction(lineSegment2.point1,lineSegment2.point2,lineSegment1.point2)
    d3 = Direction(lineSegment1.point1,lineSegment1.point2,lineSegment2.point1)
    d4 = Direction(lineSegment1.point1,lineSegment1.point2,lineSegment2.point2)
    if d1*d2<0 and d3*d4<0:
        return True
    elif d1==0 and OnSegment(lineSegment1.point1,lineSegment2.point1,lineSegment2.point2):
        return True
    elif d2==0 and OnSegment(lineSegment1.point2,lineSegment2.point1,lineSegment2.point2):
        return True
    elif d3==0 and OnSegment(lineSegment2.point1,lineSegment1.point1,lineSegment1.point2):
        return True
    elif d4==0 and OnSegment(lineSegment2.point2,lineSegment1.point1,lineSegment1.point2):
        return True
    else:
        return False


# 初始化多边形时，点的列表不包括初始点，连接点的顺序按照自然顺序
# 如果自然顺序连接不成多边形，则改变点的顺序
# 多边形类
class Polygon:
    def __init__(self,points):
        self.points=points
        self.n = len(points)
        self.edges=[]
        for i in range(self.n-1):
            self.edges.append(LineSegment(points[i],points[i+1]))
        self.edges.append(LineSegment(points[-1],points[0]))
    #     画多边形
    def plot(self):
        x=[]
        y=[]
        for point in self.points:
            x.append(point.x)
            y.append(point.y)
        x.append(self.points[0].x)
        y.append(self.points[0].y)
        plt.plot(x,y)
        # plt.plot([1,2],[1,2])
        plt.show()
    #     算多边形的面积
    def area(self):
        area=0
        for i in range(self.n):
            area+=(self.points[i].x*self.points[(i+1)%self.n].y-self.points[(i+1)%self.n].x*self.points[i].y)/2
        return abs(area)
#     TODO：判断点是否在多边形内部
#    判断点是否在多边形内部，不包括边界
#    对该点作水平向右的射线，判断射线与多边形的交点个数，如果是奇数个，则在多边形内部，否则在多边形外部
    def IsInside(self,point):
        cnt=0
        for i in range(self.n):
            if self.points[i%self.n].y == self.points[(i+1)%self.n].y:
                continue
            if point.y < min(self.points[i%self.n].y,self.points[(i+1)%self.n].y):
                continue
            if point.y >= self.points[i%self.n].y and point.y >= self.points[(i+1)%self.n].y:
                continue
            x = (point.y - self.points[i%self.n].y)*(self.points[(i+1)%self.n].x - self.points[i%self.n].x)/(self.points[(i+1)%self.n].y - self.points[i%self.n].y) + self.points[i%self.n].x
            if x > point.x:
                cnt+=1

        # # special用来记录射线与多边形的交点是顶点以及水平边的特殊情况，在special中的边在最简单的情况中不需要再判断
        # special=[]
        # # 考虑到射线与多边形的交点可能是顶点，该情况比较特殊，单独提前判断
        # # 第一种是判断水平边的情况
        # for i in range(self.n):
        #       # 如果射线与多边形的交点是顶点，则通过判断这个顶点的前后两个边是否位于射线的两侧来确定cnt加多少，同侧加2，两侧加1
        #       # 对于水平边，忽视，考虑与水平边相连的两条边是否在射线的两侧，同侧加2，两侧加1
        #
        #       if point.y == self.points[i%self.n].y and point.y == self.points[(i+1)%self.n].y and point.x < self.points[i%self.n].x and point.x < self.points[(i+1)%self.n].x:
        #             if (self.points[(i-1)%self.n].y - point.y)*(self.points[(i+2)%self.n].y - point.y) > 0:
        #                 cnt+=2
        #             else:
        #                 cnt+=1
        #             #     记录下有关顶点和边，避免重复计算
        #             special.append(i%self.n)
        #             special.append((i-1)%self.n)
        #             special.append((i+1)%self.n)
        # # 第二种是射线与多边形交点是顶点的情况，但排除了水平边的情况
        # for i in range(self.n):
        #     if i not in special:
        #         if point.y == self.points[i%self.n].y and point.x < self.points[i%self.n].x:
        #             if (self.points[(i-1)%self.n].y - point.y)*(self.points[(i+1)%self.n].y - point.y) > 0:
        #                 cnt+=2
        #             else:
        #                 cnt+=1
        #             special.append(i%self.n)
        #             special.append((i-1)%self.n)
        # # 第三种是最简单的情况，射线与多边形的交点不是顶点
        # for i in range(self.n):
        #     if i not in special:
        #         if self.points[i%self.n].y == self.points[(i+1)%self.n].y:
        #             continue
        #         if point.y < min(self.points[i%self.n].y,self.points[(i+1)%self.n].y):
        #             continue
        #         if point.y > max(self.points[i%self.n].y,self.points[(i+1)%self.n].y):
        #             continue
        #         x = (point.y - self.points[i%self.n].y)*(self.points[(i+1)%self.n].x - self.points[i%self.n].x)/(self.points[(i+1)%self.n].y - self.points[i%self.n].y) + self.points[i%self.n].x
        #         if x > point.x:
        #             cnt+=1
        return cnt%2==1

#    TODO：判断两顶点连线是不是对角线
#  判断两顶点连线是不是对角线（判断连续的三个点形成的线段是不是对角线）
    def IsDiagonal(self,point1,point2):
        # 两顶点连线有三种情况：
        # 1.完全在多边形内部，是对角线
        # 2.完全在多边形外部
        # 3.部分在多边形内部，部分在多边形外部。
        # 第3种情况，可以通过判断连线与多边形的边的相交情况来确定
        # 如果与多边形的边相交大于4次，有可能在多边形外部，此时在连线上任取一点判断是否在多边形内部即可。
        # 正常的对角线与多边形的边相交4次（如果我们对多边形的每一条边都单独判断的话）
        cnt=0
        for edge in self.edges:
            if IsIntersect(edge,LineSegment(point1,point2)):
                cnt+=1
                if cnt>4:
                    return False
        return self.IsInside(Point((point1.x+point2.x)/2,(point1.y+point2.y)/2))

#   TODO：多边形三角剖分算法(递归切耳算法)
#   返回值是对角线的列表
    def triangulation(self):
        diagonal = []
        # 递归出口，三角形的情况
        if self.n == 3:
            return diagonal
        seed = int(time.time())
        permutation = list(range(self.n))
        random.Random(seed).shuffle(permutation)
        # 遍历所有的点，找到一个对角线 或者 采用随机化的方法，随机选取一个点，找到一个对角线
        # for i in range(self.n):
        for i in permutation:
            # 现在考虑一种新的情况，即出现“平点”（在该点处的内角为180度）的情况
            if self.IsDiagonal(self.points[i % self.n], self.points[(i + 2) % self.n]) and IsParallel(self.points[i % self.n], self.points[(i + 1) % self.n], self.points[(i + 2) % self.n]) == False:
                # 还需要考虑四边形有一个平点导致退化成三角形的情况
                if self.n == 4 and IsParallel(self.points[i % self.n], self.points[(i - 1) % self.n], self.points[(i - 2) % self.n]):
                    continue
                diagonal.append(LineSegment(self.points[i % self.n], self.points[(i + 2) % self.n]))
                # 递归求解剩下的多边形
                newpoints = []
                for j in range(self.n):
                    if j != (i + 1) % self.n:
                        newpoints.append(self.points[j])
                newpolygon = Polygon(newpoints)
                diagonal.extend(newpolygon.triangulation())
                # 找到一个对角线就可以了
                break
        return diagonal
    # 把三角剖分的结果画出来
    def plotTriangulation(self):
        diagonal = self.triangulation()
        for edge in diagonal:
            plt.plot([edge.point1.x,edge.point2.x],[edge.point1.y,edge.point2.y])
        self.plot()


if __name__ == '__main__':
    # i.18
    # points = [
    #     Point(0, 0),
    #     Point(10, 7),
    #     Point(12, 3),
    #     Point(20, 8),
    #     Point(13, 17),
    #     Point(10, 12),
    #     Point(12, 14),
    #     Point(14, 9),
    #     Point(8, 10),
    #     Point(6, 14),
    #     Point(10, 15),
    #     Point(7, 18),
    #     Point(0, 16),
    #     Point(1, 13),
    #     Point(3, 15),
    #     Point(5, 8),
    #     Point(-2, 9),
    #     Point(5, 5)
    # ]
    # i.snake
    # points = [
    #     Point(10, 0),
    #     Point(20, 10),
    #     Point(30, 0),
    #     Point(40, 10),
    #     Point(50, 0),
    #     Point(50, 10),
    #     Point(40, 20),
    #     Point(30, 10),
    #     Point(20, 20),
    #     Point(10, 10),
    #     Point(0, 20),
    #     Point(0, 10)
    # ]
    # 正方形
    # points = [
    #     Point(0, 0),
    #     Point(10, 0),
    #     Point(10, 10),
    #     Point(0, 10)
    # ]
    # 带平点的三角形
    # points = [
    #     Point(0, 0),
    #     Point(10, 0),
    #     Point(20, 0),
    #     Point(30, 0),
    #     Point(40, 0),
    #     Point(20, 10),
    # ]
    # 五角星
    # points = [
    #     Point(20.0, 0.0),
    #     Point(8.090169943749475, 5.877852522924732),
    #     Point(6.180339887498949, 19.02113032590307),
    #     Point(-3.0901699437494736, 9.510565162951536),
    #     Point(-16.180339887498945, 11.755705045849465),
    #     Point(-10.0, 1.2246467991473533e-15),
    #     Point(-16.18033988749895, -11.75570504584946),
    #     Point(-3.0901699437494754, -9.510565162951535),
    #     Point(6.180339887498945, -19.021130325903073),
    #     Point(8.090169943749473, -5.877852522924734),
    # ]
    # 六边形
    # points = [
    #     Point(17.320508075688775, 9.999999999999998),
    #     Point(17.32050807568877, -10.000000000000002),
    #     Point(3.673940397442059e-15, -20.0),
    #     Point(-17.320508075688767, -10.000000000000009),
    #     Point(-17.32050807568877, 10.0),
    #     Point(-6.123233995736766e-15, 20.0)
    #
    #
    # ]
    # 特殊的多边形
    points = [
        Point(0, 0),
        Point(50, 0),
        Point(50, 30),
        Point(40, 30),
        Point(40, 20),
        Point(30, 20),
        Point(30, 30),
        Point(20, 30),
        Point(20, 20),
        Point(10, 20),
        Point(10, 30),
        Point(0, 30)
    ]
    #
    # points = [
    #     Point(0, 0),
    #     Point(1, 0.3),
    #     Point(0.7, 0.7),
    #     Point(1, 1),
    #     Point(0.7, 1.3),
    #     Point(0.3, 0.7),
    #     Point(0.3, 1),
    #     Point(0, 1),
    #     Point(-0.3, 0.7),
    #     Point(-0.7, 1),
    #     Point(-1, 1.3),
    #     Point(-1, 0.7),
    # ]
    # 六角星
    # points = [
    #     Point(6, 3),
    #     Point(7, 4),
    #     Point(9, 4),
    #     Point(7.5, 5.5),
    #     Point(9, 7),
    #     Point(7, 7),
    #     Point(6, 8),
    #     Point(5, 7),
    #     Point(3, 7),
    #     Point(4.5, 5.5),
    #     Point(3, 4),
    #     Point(5, 4)
    # ]
    polygon=Polygon(points)
    print(polygon.area())
    # polygon.plot()
    polygon.plotTriangulation()

一些注释：

判定点是否在多边形的内部的算法，作者之前针对特殊情况采用了另一种方法，后来发现书上的原算法就可以实现同样的效果
i.18和i.snake都是点的数据集
对于特殊的多边形（某个顶点处的内角为平角），可能不能实现三角剖分（已改进）
多边形初始化时，输入的点集需要保证其自然的顺序（列表顺序）可以连接成为一个多边形，不需要在点集末尾再包括起点
导入了random和time模块，用于三角剖分算法初始化，可以改变随机种子来找到不同的三角剖分方案
代码已同步发布在作者的github仓库计算几何实验代码上，欢迎大家收藏，我会更新许多课程的实验代码，希望能对读者有用

实现效果：

最后如果作者有误，欢迎大家交流指正。

DAMO开发者矩阵

DAMO开发者矩阵，由阿里巴巴达摩院和中国互联网协会联合发起，致力于探讨最前沿的技术趋势与应用成果，搭建高质量的交流与分享平台，推动技术创新与产业应用链接，围绕“人工智能与新型计算”构建开放共享的开发者生态。

更多推荐

第 1 章：为什么嵌入式工程师要学 ROS2？

DAMO开发者矩阵

[具身智能-612]：RAW / NV12 / JPG 变换链路、转换关系与工程用途（适配 RDK X5 MIPI+AI 检测链路）

✅本质：摄像头感光二极管直接输出的原始采样数据，只有单色感光点不是可视图像。表格格式压缩属性上游来源下游输出核心用途能否直接给 AI 推理RAW(Bayer)无压缩仅 ISP底层调试、标定❌ 不能NV12无压缩ISP 输出（RAW 处理后）VPS/BPU/H264/JPG 编码器AI 推理、视频编码、本地预览✅ 推荐首选JPG有损压缩NV12 编码生成网络传输、图片存储网页预览、截图存档❌ 不推荐

DAMO开发者矩阵

[具身智能-611]：JPG / RAW / NV12 格式深度解析（适配 RDK X5 MIPI 相机场景）

才能生成可视图像：去马赛克 (Demosaic) → 白平衡 → 降噪 → 伽马校正 → 输出 YUV (NV12)：mipi_cam 底层 VIO 接收层；一旦经过 ISP，输出就不再是 RAW。网络只能接收 YUV/RGB。基于 YUV420 进行离散余弦有损压缩，带文件头、元信息，通用图片格式。有损压缩，损失图像细节，不建议送入 AI 目标检测网络；采样比例 4:2:0：色度分辨率宽高各减