76、多边形一些基本操作（自相交、尖刺、保证逆时针、求交）

最新推荐文章于 2024-05-31 18:36:38 发布
chaoyuzhang
最新推荐文章于 2024-05-31 18:36:38 发布
阅读量1.1k
点赞数 3
分类专栏：图形学
本文链接：https://blog.csdn.net/chaoyuzhang/article/details/115918738
版权
图形学专栏收录该内容
3 篇文章 0 订阅
订阅专栏
这里求交集采用的思路是找到内边，沿内边不断延展直到形成多边形
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)

    #并集处理
chaoyuzhang
关注
3
点赞
踩
8

收藏

觉得还不错? 一键收藏
1
评论
76、多边形一些基本操作（自相交、尖刺、保证逆时针、求交）

pro = """1、造一些多边形和一个大的范围2、多边形交集部分归属问题3、多边形未达部分归属问题"""import random as rdimport mathimport matplotlib.pyplot as pltfrom tool.Convex_hull_find import boundary_find#多边形生成def Polygon_generate(x_min, y_min, x_max, y_max, num_point): L = [[int(.
复制链接

扫一扫