预备知识
概念主要参考了Discrete and Computational Geometry 这本书。每一个概念我都会给出原书中的陈述以及我自己翻译的版本。
简单多边形:多边形是平面上的封闭区域,以形成不与自身相交的封闭曲线的有限线段集为边界。
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.
多边形的边界:多边形的顶点和边的集合是多边形的边界,记作
boundary of polygon: The set of vertices and edges of P is called the boundary of the polygon, denoted as ∂P.
多边形约当曲线定理:一个多边形的边界把平面分成两部分。特别地,的两个部分是有界的内部和无界的外部。
Polygonal Jordan Curve: The boundary ∂P of a polygon P partitions the plane into two parts. In particular, the two components of are the bounded interior and the unbounded exterior.
对角线:多边形的对角线是一条连接多边形两个顶点且完全在
的内部,不与多边形的边界
相交除了在它的两个端点处的线段。
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.
三角剖分:多边形的三角剖分是把
分解成三角形,要求分解时添加的对角线不相交且数量最多。
triangulation: A triangulation of a polygon P is a decomposition of P into triangles by a maximal set of noncrossing diagonals.
多边形的耳朵:三个连续的顶点形成一个多边形的耳朵,如果
是该多边形的对角线。
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.
算法思路分析
考虑到多边形的耳朵的特殊性,我们可以通过对多边形一直切耳朵来实现三角剖分。
多边形切去一个耳朵仍然是多边形,所以该算法可以写成递归算法。
除了三角形外每个多边形都有至少两个耳朵,所以递归算法的出口就是三角形。
所以这里的重点是如何识别多边形的耳朵,按照定义,我们只要枚举所有连续的三个顶点判断两个不相邻的顶点连接而成的线段
是不是对角线就可以了。
两顶点连线有四种可能的情况(以下情况都省略掉
的端点):
-
完全在多边形内部,是对角线
-
完全在多边形外部,不是对角线
-
部分在多边形内部,部分在多边形外部,不是对角线(包括与多边形的边界相交于顶点处的情况)
-
与多边形的边界完全重合(这部分是新增的,对应于
是平点,即在该点处内角为180度)
第3种情况,可以通过判断是否与多边形的边界相交(不包括
的端点)来判别。
如果不与多边形的边界相交,那么可以随意取
上一点(不包括端点,下同),通过判断该点是否在多边形的内部来判别第1种和第2种情况,如果在多边形内部,则整条线段都在多边形内部,如果在多边形外部,则整条线段都在多边形外部。
第4种情况,只需要单独判断是否为平点即可,只需要判断
和
是否平行。注意,当多边形为带平点的四边形时,需要额外考虑算法实现中的细节。例如:
![](https://i-blog.csdnimg.cn/blog_migrate/9e6ab492320bc1c7aaa79dde5c99737b.png)
如果只判断和
是否平行是不足以判断
是否是对角线的,还需要判断
和
是否平行。
所以现在需要实现的算法主要有
-
线段的相交判定算法(包括端点)
-
一个点是否在多边形的内部的判定算法(不需要判断该点是否在多边形的边界,因为
不与多边形的边界相交)
-
针对第4种情况的线段平行判定算法
这里算法实现主要参考了《算法设计与分析》第2版 李春葆 中的计算几何部分。
判定点是否在多边形内部的算法大致思路:
过该点作一条射线,如果射线与多边形相交的点的个数是奇数,那它在多边形内部,如果是偶数,则在多边形外部。为了方便,我们一般使用水平向右的射线。
这里需要考虑一些算法实现方面的细节。
在判断射线与多边形的相交点的个数时,我们肯定会对多边形的每一条边都进行判断。
如果射线与多边形交于顶点,如何判断交点的个数?
有两种情况:
![](https://i-blog.csdnimg.cn/blog_migrate/b0877b2bfd296bb73c768fb8c66a748d.png)
![](https://i-blog.csdnimg.cn/blog_migrate/bf1949b5f50e3d853e44c4abf79c8cfc.png)
对于第1种情况,两条邻接边记成有偶数个交点即可
对于第2种情况,两条邻接边记成有奇数个交点即可
还有可能射线与多边形的水平边相交,如何判断?
![](https://i-blog.csdnimg.cn/blog_migrate/5e2daca04c5e79f07ebf0b4122d0e192.png)
此时射线必然与另外两条边的端点相交,忽视掉水平边(即不计入点的个数),然后转化成之前的情况即可。
《算法设计与分析》中的判断点是否在多边形内部的算法实现了同样的效果。
算法代码
# 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仓库计算几何实验代码上,欢迎大家收藏,我会更新许多课程的实验代码,希望能对读者有用
实现效果:
![](https://i-blog.csdnimg.cn/blog_migrate/8924ee84fb73e1b3aa2df3ec923789dc.png)
![](https://i-blog.csdnimg.cn/blog_migrate/fce6e65b1b4f09a71f1e64f88c50b98b.png)
![](https://i-blog.csdnimg.cn/blog_migrate/4f19b966e49bf890f6ddd1b417535beb.png)
![](https://i-blog.csdnimg.cn/blog_migrate/86b819020fda34e11e74ca7940cadddb.png)
最后如果作者有误,欢迎大家交流指正。