平面中任意方位矩形的重叠面积——python

import numpy as np
import matplotlib.pyplot as plt
from copy import deepcopy

class Rec():
    def __init__(self,center,l,w,theta):
        # 角度为长l与水平边的夹角，逆时针方向，弧度角
        self.center=center
        self.l=l
        self.w=w
        self.theta=theta
        # 计算出四个顶点的坐标
        l_vec=np.array([l/2*np.cos(theta),l/2*np.sin(theta)])
        w_vec=np.array([w/2*np.cos(3/2*np.pi+theta),w/2*np.sin(3/2*np.pi+theta)])
        # 将四个顶点存起来
        self.points=[]
        self.points.append(l_vec+w_vec+np.array(self.center))
        self.points.append(l_vec-w_vec+np.array(self.center))
        self.points.append(-l_vec-w_vec+np.array(self.center))
        self.points.append(-l_vec+w_vec+np.array(self.center))

    def draw(self):
        #画出矩形
        rec_points = deepcopy(self.points)
        rec_points.append(self.points[0])

        rec_points=np.array(rec_points)

        plt.plot(rec_points[:,0],rec_points[:,1],'-')



def point2rec(point,rec):
    #判断一个点是否在一个矩形中
    rec_points = deepcopy(rec.points)
    rec_points.append(rec.points[0])
    rec_points = np.array(rec_points)

    flag=None
    for i in range(4):
        #主要是利用叉乘判断一个点是否在一个矩形内
        tmp=np.cross(rec_points[i+1,:]-rec_points[i,:],point-rec_points[i+1,:])
        if flag is None:
            flag=tmp
        else:
            if flag*tmp<-1e-6:
                return False
    return True




def line2line(line1,line2):
    xa,ya = line1[0,0],line1[0,1]
    xb,yb = line1[1,0],line1[1,1]
    xc,yc = line2[0,0],line2[0,1]
    xd,yd = line2[1,0],line2[1,1]
    #判断两条直线是否相交，矩阵行列式计算
    a = np.matrix(
        [
            [xb-xa,-(xd-xc)],
            [yb-ya,-(yd-yc)]
        ]
    )
    delta = np.linalg.det(a)
    #不相交,返回两线段
    if np.fabs(delta) < 1e-6:
        return []
    #求两个参数lambda和miu
    c = np.matrix(
        [
            [xc-xa,-(xd-xc)],
            [yc-ya,-(yd-yc)]
        ]
    )
    d = np.matrix(
        [
            [xb-xa,xc-xa],
            [yb-ya,yc-ya]
        ]
    )
    lamb = np.linalg.det(c)/delta
    miu = np.linalg.det(d)/delta
    #相交
    if lamb <= 1 and lamb >= 0 and miu >= 0 and miu <= 1:
        x = xc + miu*(xd-xc)
        y = yc + miu*(yd-yc)
        return [x,y]
    #相交在延长线上
    else:
        return []



def line2rec(line,rec):
    #线段与矩形的交点
    rec_points=deepcopy(rec.points)
    rec_points.append(rec.points[0])
    rec_points=np.array(rec_points)

    cross_points=[]

    for i in range(4):
        cross_point=line2line(line,rec_points[i:i+2,:])
        if len(cross_point) !=0:
            cross_points.append([cross_point[0],cross_point[1]])

    if len(cross_points)==1:
        #把端点也作为多边行的顶点
        if point2rec(line[0,:],rec):
            cross_points.append(list(line[0,:]))
        else:
            cross_points.append(list(line[1,:]))

    if len(cross_points)==0 and point2rec(line[0,:],rec) and point2rec(line[1,:],rec):
        #如果一条线与矩形没有交点，需要考虑点在矩形内部
        cross_points.append(list(line[0,:]))
        cross_points.append(list(line[1, :]))


    return cross_points




def rec2rec(rec1,rec2):
    #矩形与矩形的交点，形成多边形
    rec1_points = deepcopy(rec1.points)
    rec1_points.append(rec1.points[0])
    rec1_points=np.array(rec1_points)

    rec2_points = deepcopy(rec2.points)
    rec2_points.append(rec2.points[0])
    rec2_points = np.array(rec2_points)

    cross_points=[]



    for i in range(4):
        #求出矩形1每条边与矩形2的交点
        cross_point=line2rec(rec1_points[i:i+2,:], rec2)
        if len(cross_point)!=0:
                cross_points.extend(cross_point)

    for i in range(4):
        #求出矩形2每条边与矩形1的交点
        cross_point = line2rec(rec2_points[i:i + 2, :], rec1)
        if len(cross_point) != 0:
            cross_points.extend(cross_point)
    return cross_points


def refine_points(points):
    if len(points)==0:
        return [],0
    #将重复的点删除，按照顺序形成多边形
    new_points=np.array([points.pop(0),points.pop(0)])
    while np.sum((new_points[-1,:]-new_points[0,:])**2)>1e-6:
        tmp_point=new_points[-1,:]

        ind=np.nonzero(np.sum((np.array(points)-tmp_point)**2,1)<1e-6)[0]

        if ind[0]%2==0:
            new_points=np.concatenate((new_points,
                                       np.array(points.pop(ind[0]+1)).reshape((1,-1))),0)
            points.pop(ind[0])
        else:
            new_points = np.concatenate((new_points,
                                         np.array(points.pop(ind[0]-1)).reshape((1,-1))), 0)
            points.pop(ind[0]-1)

    S=0#求面积
    for i in range(1,new_points.shape[0]-2):
        a=np.linalg.norm(new_points[0,:]-new_points[i,:])
        b=np.linalg.norm(new_points[i,:]-new_points[i+1,:])
        c=np.linalg.norm(new_points[0,:]-new_points[i+1,:])
        S+=tri_area(a,b,c)



    return new_points,S







def tri_area(a,b,c):
    #三角形面积
    # 计算半周长
    s = (a + b + c) / 2
    # 计算面积
    area = (s * (s - a) * (s - b) * (s - c)) ** 0.5
    return area




#定义两个矩形
Rec1=Rec(center=[1,2],l=5,w=3,theta=0)
Rec2=Rec(center=[2,2],l=3,w=4,theta=np.pi/3)
Rec1.draw()
Rec2.draw()


#求出矩形1与矩形2的交点
cross_points=rec2rec(Rec1,Rec2)
#将重复的交点去掉，并顺序表示
points,S=refine_points(cross_points)

if len(points)!=0:
    #只有当有交点时才画出相交面积
    plt.plot(points[:,0],points[:,1],'ro-')
    plt.fill(points[:,0],points[:,1], facecolor='r',alpha=0.5)

plt.title('S={:.4f}'.format(S))
plt.axis('equal')
plt.show()