计算几何实验(一):简单多边形的三角剖分算法

预备知识

概念主要参考了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为平点即可,只需要判断abbc是否平行。注意,当多边形为带平点的四边形时,需要额外考虑算法实现中的细节。例如:

特殊四边形

如果只判断abbc是否平行是不足以判断ac是否是对角线的,还需要判断addc是否平行。

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

  • 线段的相交判定算法(包括端点)

  • 一个点是否在多边形的内部的判定算法(不需要判断该点是否在多边形的边界,因为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仓库计算几何实验代码上,欢迎大家收藏,我会更新许多课程的实验代码,希望能对读者有用

实现效果:

i.18

i.snake
带平点的多边形
递归切耳后会出现平点

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

  • 3
    点赞
  • 7
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值