这里求交集采用的思路是找到内边,沿内边不断延展直到形成多边形
pro = """
1、造一些多边形和一个大的范围
2、多边形交集部分归属问题
3、多边形未达部分归属问题
"""
import random as rd
import math
import matplotlib.pyplot as plt
from tool.Convex_hull_find import boundary_find
#多边形生成
def Polygon_generate(x_min, y_min, x_max, y_max, num_point):
L = [[int(x_min+(x_max-x_min)*rd.random()), int(y_min+(y_max-y_min)*rd.random())] for i in range(num_point)]
return L + [L[0]]
#单个多边形画图
def Polygon_show(P):
x, y = [c[0] for c in P], [c[1] for c in P]
plt.plot(x, y)
plt.show()
#多个多边形画图
def MultiPolygon_show(LP):
for P in LP:
x, y = [c[0] for c in P], [c[1] for c in P]
plt.plot(x, y)
plt.show()
#直线方程,然后计算x,y在直线的哪边
def abc(x1, y1, x2, y2, x, y):
return (y2-y1)*x-(x2-x1)*y+y1*(x2-x1)-x1*(y2-y1)
#判断一个多边形是否自相交
def is_self_intersect(P):
for i in range(len(P)-2):
for j in range(i+1, len(P)-1):
x1, y1, x2, y2 = P[i][0], P[i][1], P[i+1][0], P[i+1][1]
x3, y3, x4, y4 = P[j][0], P[j][1], P[j + 1][0], P[j + 1][1]
#p1,p2在p3p4两边,p3,p4在p1p2两边
if abc(x1, y1, x2, y2, x3, y3)*abc(x1, y1, x2, y2, x4, y4) < 0 and abc(x3, y3, x4, y4, x1, y1)*abc(x3, y3, x4, y4, x2, y2) < 0:
return True
return False
#多边形去重,只有邻点一样才去
def polygon_deduplication(P):
marks = [False] + [True if P[i-1][0] == P[i][0] and P[i-1][1] == P[i][1] else False for i in range(1, len(P))]
P1 = [P[i] for i in range(len(P)) if not marks[i]]
return P1
def polygon_fill(P):
flag = True
while flag:
flag = False
for i in range(len(P) - 2):
for j in range(i + 1, len(P) - 1):
x1, y1, x2, y2 = P[i][0], P[i][1], P[i + 1][0], P[i + 1][1]
x3, y3, x4, y4 = P[j][0], P[j][1], P[j + 1][0], P[j + 1][1]
if abc(x1, y1, x2, y2, x3, y3) * abc(x1, y1, x2, y2, x4, y4) < 0 and abc(x3, y3, x4, y4, x1, y1) * abc(x3, y3, x4, y4, x2, y2) < 0:
x_new = ((x4-x3)*(x2-x1)*(y3-y1)-(y4-y3)*(x2-x1)*x3+(y2-y1)*(x4-x3)*x1)/((y2-y1)*(x4-x3)-(y4-y3)*(x2-x1))
y_new = ((y4-y3)*(y2-y1)*(x3-x1)-(x4-x3)*(y2-y1)*y3+(x2-x1)*(y4-y3)*y1)/((x2-x1)*(y4-y3)-(x4-x3)*(y2-y1))
P = P[:i+1]+[[x_new,y_new]]+P[i+1:j+1]+[[x_new,y_new]]+P[j+1:]
flag = True
break
if flag:
break
return P
def polygon_split(P):
LP = []
i = 1
while i < len(P):
if P[i] in P[:i]:
idex = P[:i].index(P[i])
LP.append(P[idex:i+1])
P = P[:idex+1] + P[i+1:]
i = idex
i += 1
return LP
#三个点组成的三角形面积
def area_triangle(a, b, c):
ab = [b[0]-a[0], b[1]-a[1]]
ac = [c[0]-a[0], c[1]-a[1]]
return 0.5*abs(ab[0]*ac[1]-ab[1]*ac[0])
#多边形面积,拆成三角形的和
def area_polygon(P):
return sum([area_triangle(P[0], P[i], P[i+1]) for i in range(1, len(P)-2)])
#计算每个多边形面积并取最大的一个多边形留下来
def polygon_choose(LP):
LS = [[area_polygon(P), P] for P in LP]
LS.sort(key=lambda x:x[0],reverse=True)
return LS[0][1]
#处理自相交
def process_self_intersect(P):
#相邻的点如果一样,则去重
P = polygon_deduplication(P)
#把交点补上
P = polygon_fill(P)
#按交点拆成多个多边形
LP = polygon_split(P)
#计算每个多边形面积并取最大的一个多边形留下来
P = polygon_choose(LP)
return P
#计算角b
def calcInsideCornor(a, b, c):
ba = [a[0]-b[0], a[1]-b[1]]
bc = [c[0]-b[0], c[1]-b[1]]
degree = 0
if (math.sqrt(ba[0]*ba[0]+ba[1]*ba[1])*math.sqrt(bc[0]*bc[0]+bc[1]*bc[1])) > 0:
degree = math.acos((ba[0]*bc[0]+ba[1]*bc[1])/(math.sqrt(ba[0]*ba[0]+ba[1]*ba[1])*math.sqrt(bc[0]*bc[0]+bc[1]*bc[1])))
degree = degree*180/math.pi
return degree
#去尖刺
def smoothPolygon(polygon, e=0.05):
#首尾点一致
polygon = polygon[:-1]
#循环砍点
loop = 0
while len(polygon) >= 3:
loop += 1
nums0 = len(polygon)
degrees = [calcInsideCornor(polygon[-1], polygon[0], polygon[1])] + \
[calcInsideCornor(polygon[i], polygon[i + 1], polygon[i + 2]) for i in range(len(polygon) - 2)] + \
[calcInsideCornor(polygon[-2], polygon[-1], polygon[0])]
print("\nloop: ", loop, "\ndegrees: ", degrees)
marks = []
for i in range(len(polygon)):
if degrees[i] <= e:
marks.append(i)
polygon = [polygon[i] for i in range(len(polygon)) if i not in marks]
if len(polygon) == nums0:
break
if len(polygon) < 3:
polygon = []
else:
polygon.append(polygon[0])
return polygon
#保证逆时针
def positive_constraint(P):
#求凸包
boundary = boundary_find(P)[:-1]
#根据凸包的走向决定是否调整方向
idex = 0
for i in range(len(P)-1):
if P[i][0] == boundary[0][0] and P[i][1] == boundary[0][1]:
idex = i
break
P1 = P[:-1][idex:] + P[:-1][:idex]
j = 0
for i in range(len(P1)):
if P1[i][0] == boundary[j][0] and P1[i][1] == boundary[j][1]:
j += 1
if j == len(boundary):
break
if j == len(boundary):
return P
else:
return [P[i] for i in range(len(P)-1, -1, -1)]
def P_cases():
cases = []
#一个自相交
cases.append([[0,0],[1,0],[0,1],[1,1],[0,0]])
#两个自相交
cases.append([[0, 0], [1, 0], [0 ,1], [1, 2], [0, 2], [1, 1], [0, 0]])
#角小于0.05
cases.append([[0, 0], [1250, 0], [0, 1], [-10, 1], [0, 0]])
return cases
def Polygon_generate_and_clean(x_min, y_min, x_max, y_max, num_point, show_polygon=True):
#1-构造多边形
P = Polygon_generate(x_min, y_min, x_max, y_max, num_point)
# P = P_cases()[2]
print(P)
if show_polygon:
Polygon_show(P)
#2-自相交判断及处理(这里相交时取大)
P = process_self_intersect(P)
if show_polygon:
Polygon_show(P)
#3-尖刺处理
P = smoothPolygon(P)
if show_polygon:
Polygon_show(P)
#4-保证逆时针
P = positive_constraint(P)
return P
def LP_cases():
LP_cases = []
#简单交集
LP_cases.append([[[0, 0], [2, 0], [2, 2], [0, 2], [0, 0]], [[1, 1], [3, 1], [3, 3], [1, 3], [1, 1]]])
#一般交集
LP_cases.append([[[91, 66], [72, 16], [20, 17], [69, 47], [91, 66]], [[91, 47], [97, 38], [71, 15], [35, 61], [91, 47]]])
# #角小于0.05
# LP_cases.append([[0, 0], [1250, 0], [0, 1], [-10, 1], [0, 0]])
return LP_cases
def polygon_fill_2(P1, P2):
flag = True
while flag:
flag = False
for i in range(len(P1)-1):
for j in range(len(P2)-1):
x1, y1, x2, y2 = P1[i][0], P1[i][1], P1[i + 1][0], P1[i + 1][1]
x3, y3, x4, y4 = P2[j][0], P2[j][1], P2[j + 1][0], P2[j + 1][1]
#如果非端点处相交,则都加点
if abc(x1, y1, x2, y2, x3, y3) * abc(x1, y1, x2, y2, x4, y4) < 0 and abc(x3, y3, x4, y4, x1, y1) * abc(x3, y3, x4, y4, x2, y2) < 0:
x_new = ((x4-x3)*(x2-x1)*(y3-y1)-(y4-y3)*(x2-x1)*x3+(y2-y1)*(x4-x3)*x1)/((y2-y1)*(x4-x3)-(y4-y3)*(x2-x1))
y_new = ((y4-y3)*(y2-y1)*(x3-x1)-(x4-x3)*(y2-y1)*y3+(x2-x1)*(y4-y3)*y1)/((x2-x1)*(y4-y3)-(x4-x3)*(y2-y1))
P1 = P1[:i + 1] + [[x_new, y_new, 1]] + P1[i+1:]
P2 = P2[:j + 1] + [[x_new, y_new, 1]] + P2[j + 1:]
flag = True
break
#3在12上
if abc(x3, y3, x4, y4, x1, y1) * abc(x3, y3, x4, y4, x2, y2) < 0 and abc(x1, y1, x2, y2, x3, y3) == 0:
P1 = P1[:i + 1] + [[x3, y3, 1]] + P1[i+1:]
flag = True
break
#4在12上
if abc(x3, y3, x4, y4, x1, y1) * abc(x3, y3, x4, y4, x2, y2) < 0 and abc(x1, y1, x2, y2, x4, y4) == 0:
P1 = P1[:i + 1] + [[x4, y4, 1]] + P1[i+1:]
flag = True
break
#1在34上
if abc(x1, y1, x2, y2, x3, y3) * abc(x1, y1, x2, y2, x4, y4) < 0 and abc(x3, y3, x4, y4, x1, y1) == 0:
P2 = P2[:j + 1] + [[x1, y1, 1]] + P2[j+1:]
flag = True
break
#2在34上
if abc(x1, y1, x2, y2, x3, y3) * abc(x1, y1, x2, y2, x4, y4) < 0 and abc(x3, y3, x4, y4, x2, y2) == 0:
P2 = P2[:j + 1] + [[x2, y2, 1]] + P2[j+1:]
flag = True
break
#其他情况:不相交或有端点重叠
if flag:
break
return P1, P2
#判断点是否在多边形上
def is_on_edge_point(point, P):
x, y = point[0], point[1]
for i in range(len(P)-1):
x1, y1, x2, y2 = P[i][0], P[i][1], P[i+1][0], P[i+1][1]
if (y-y1)*(x2-x1) == (y2-y1)*(x-x1):
return True
return False
#判断点是否在多边形内部
def is_interior_point(point, P):
if is_on_edge_point(point, P):
return False
vertx = [p[0] for p in P]
verty = [p[1] for p in P]
x, y = point[0], point[1]
res = False
for i in range(1, len(P)):
j = i-1
#点的纵坐标在中间(一边大于y,一边小于等于y)
if (((verty[i] > y) != (verty[j] > y)) and (x < (vertx[j] - vertx[i]) * (y - verty[i]) / (verty[j] - verty[i]) + vertx[i])):
res = not res
return res
#判断边p1p2是否在P内,由于这一步之前会补交点,所以不存在部分在里面,部分在外面的情况
def is_interior_edge(p1, p2, P):
#p1为内点,边即为内边
if is_interior_point(p1, P):
return True
#p1在边上时,p2在内部或边上为内边
elif is_on_edge_point(p1, P) and (is_interior_point(p2, P) or is_on_edge_point(p2, P)):
return True
return False
def intersect(P1, P2):
S = set()
D = {}
for i in range(len(P1)-1):
if is_interior_edge(P1[i], P1[i+1], P2):
S.add((P1[i][0], P1[i][1], P1[i + 1][0], P1[i + 1][1]))
if (P1[i][0], P1[i][1]) not in D:
D[P1[i][0], P1[i][1]] = []
D[P1[i][0], P1[i][1]].append([P1[i + 1][0], P1[i + 1][1]])
for i in range(len(P2)-1):
if is_interior_edge(P2[i], P2[i+1], P1):
S.add((P2[i][0], P2[i][1], P2[i + 1][0], P2[i + 1][1]))
if (P2[i][0], P2[i][1]) not in D:
D[P2[i][0], P2[i][1]] = []
D[P2[i][0], P2[i][1]].append([P2[i + 1][0], P2[i + 1][1]])
L = list(S)
LP = []
while len(L) > 0:
e = L.pop(0)
D[e[0], e[1]].remove([e[2], e[3]])
P = [[e[0], e[1]], [e[2], e[3]]]
start = (e[2], e[3])
#回到起点或者无法构成多边形,结束
while start[0] != P[0][0] or start[1] != P[0][1]:
if start not in D or len(D[start]) == 0:
break
end = D[start].pop(0)
P.append(end)
L.remove((start[0], start[1], end[0], end[1]))
start = (end[0], end[1])
#回到起点便将多边形加入
if start[0] == P[0][0] and start[1] == P[0][1]:
LP.append(P)
print(LP)
return LP
def union(P1, P2):
pass
if __name__ == '__main__':
#1-确定区间范围及多边形点数
x_min, y_min, x_max, y_max, num_point = 0,0,100,100,4
# #2-单个多边形生成并展示
# P = Polygon_generate_and_clean(x_min, y_min, x_max, y_max, num_point, show_polygon=True)
#3-多边形list构造
num_polygon = 2
LP = [Polygon_generate_and_clean(x_min, y_min, x_max, y_max, num_point, show_polygon=False) for i in range(num_polygon)]
LP = LP_cases()[0]
print(LP)
# MultiPolygon_show(LP)
#3.1-交集处理
P1, P2 = LP[0], LP[1]
#补交点,若为交点则点长度为3且最后一位数为1
P1, P2 = polygon_fill_2(P1, P2)
#算交集
I = intersect(P1, P2)
#并集处理