在Python中计算两个旋转矩形的相交区域面积

Question

我有两个2D旋转矩形，由一个（中心x，中心y，高度，宽度）和旋转角度（0-360°）定义。如何计算这两个旋转矩形的交集面积。

- user824624

---

将顶点投影到轴上并检查这些点的交点。 - Netwave

可能是旋转矩形的两个交集的面积的重复问题。 - custom_user

这并不是一个软件工程问题，而是一个数学问题。 - Philipp

作为数学问题，对于数学家来说很正常。作为一名软件工程师，我在碰撞检测、OCR、计算机视觉和其他需要功能代码的系统中遇到过这种类型的问题，而不是数学理论。 - Bug

我们可以假设两条边不共线吗？ - Jared Goguen

这里有一个关于多边形面积计算的很好的解释：https://www.youtube.com/watch?v=0KjG8Pg6LGk。结合链接和提示，以及相关问题（区域：交集=多边形1+多边形2-聚合多边形），缺失的部分是将您的数据转换为顶点的（x，y）坐标。即一些简单的正弦和余弦计算。 - VPfB

Answer 1

这种任务可以使用计算几何包来解决，例如Shapely：

import shapely.geometry
import shapely.affinity

class RotatedRect:
    def __init__(self, cx, cy, w, h, angle):
        self.cx = cx
        self.cy = cy
        self.w = w
        self.h = h
        self.angle = angle

    def get_contour(self):
        w = self.w
        h = self.h
        c = shapely.geometry.box(-w/2.0, -h/2.0, w/2.0, h/2.0)
        rc = shapely.affinity.rotate(c, self.angle)
        return shapely.affinity.translate(rc, self.cx, self.cy)

    def intersection(self, other):
        return self.get_contour().intersection(other.get_contour())


r1 = RotatedRect(10, 15, 15, 10, 30)
r2 = RotatedRect(15, 15, 20, 10, 0)

from matplotlib import pyplot
from descartes import PolygonPatch

fig = pyplot.figure(1, figsize=(10, 4))
ax = fig.add_subplot(121)
ax.set_xlim(0, 30)
ax.set_ylim(0, 30)

ax.add_patch(PolygonPatch(r1.get_contour(), fc='#990000', alpha=0.7))
ax.add_patch(PolygonPatch(r2.get_contour(), fc='#000099', alpha=0.7))
ax.add_patch(PolygonPatch(r1.intersection(r2), fc='#009900', alpha=1))

pyplot.show()

- Leon

---

嗨，您如何计算与此交集的面积？ - Bug

@Bug 一旦你获得了交叉区域的形状 r1.intersection(r2) （它是一个 Shapely 对象），你只需要访问它的 area 属性：r1.intersection(r2).area。 - Leon

shapely真是太棒了！ - Rick

我已经通过pip install shapely安装了shapely，它显示Requirement already satisfied: shapely in /usr/local/lib/python2.7/dist-packages，但是当我运行代码时，它总是告诉我们ImportError: No module named affinity，有什么线索吗？ - user824624

如果您有每个矩形的4个坐标，如何做到这一点？ - Majid Azimi

Answer 2

这里有一个解决方案，不使用Python标准库之外的任何库。确定两个矩形相交的区域可以分为两个子问题：

- 找到相交的多边形（如果有）； - 确定相交多边形的面积。

当您使用矩形的顶点（角落）时，这两个问题都相对容易。因此，首先必须确定这些顶点。假设坐标原点位于矩形的中心，顶点从左下角开始按逆时针方向排列： (-w/2, -h/2)、(w/2, -h/2)、(w/2, h/2) 和 (-w/2, h/2)。 将其旋转角度a并将它们平移到矩形中心的正确位置，这些变成： (cx + (-w/2)cos(a) - (-h/2)sin(a), cy + (-w/2)sin(a) + (-h/2)cos(a))，其他角落点也是类似的。

一种简单的确定交集多边形的方法如下： 将一个矩形作为候选交集多边形。 然后应用顺序切割过程（如此处所述）。 简而言之：依次处理第二个矩形的每条边， 并删除在该边定义的“外”半平面上的所有部分， （向两个方向延伸）从候选交集多边形中移除。 对于所有边都这样做，留下候选交集多边形 仅包含第二个矩形内部或边界上的部分。

可以从顶点的坐标计算出结果多边形的面积， 通过对每条边的顶点进行叉积求和（再次按逆时针顺序）， 并将其除以二。例如，请参见www.mathopenref.com/coordpolygonarea.html 足够的理论和解释。以下是代码：

from math import pi, cos, sin


class Vector:
    def __init__(self, x, y):
        self.x = x
        self.y = y

    def __add__(self, v):
        if not isinstance(v, Vector):
            return NotImplemented
        return Vector(self.x + v.x, self.y + v.y)

    def __sub__(self, v):
        if not isinstance(v, Vector):
            return NotImplemented
        return Vector(self.x - v.x, self.y - v.y)

    def cross(self, v):
        if not isinstance(v, Vector):
            return NotImplemented
        return self.x*v.y - self.y*v.x


class Line:
    # ax + by + c = 0
    def __init__(self, v1, v2):
        self.a = v2.y - v1.y
        self.b = v1.x - v2.x
        self.c = v2.cross(v1)

    def __call__(self, p):
        return self.a*p.x + self.b*p.y + self.c

    def intersection(self, other):
        # See e.g.     https://en.wikipedia.org/wiki/Line%E2%80%93line_intersection#Using_homogeneous_coordinates
        if not isinstance(other, Line):
            return NotImplemented
        w = self.a*other.b - self.b*other.a
        return Vector(
            (self.b*other.c - self.c*other.b)/w,
            (self.c*other.a - self.a*other.c)/w
        )


def rectangle_vertices(cx, cy, w, h, r):
    angle = pi*r/180
    dx = w/2
    dy = h/2
    dxcos = dx*cos(angle)
    dxsin = dx*sin(angle)
    dycos = dy*cos(angle)
    dysin = dy*sin(angle)
    return (
        Vector(cx, cy) + Vector(-dxcos - -dysin, -dxsin + -dycos),
        Vector(cx, cy) + Vector( dxcos - -dysin,  dxsin + -dycos),
        Vector(cx, cy) + Vector( dxcos -  dysin,  dxsin +  dycos),
        Vector(cx, cy) + Vector(-dxcos -  dysin, -dxsin +  dycos)
    )

def intersection_area(r1, r2):
    # r1 and r2 are in (center, width, height, rotation) representation
    # First convert these into a sequence of vertices

    rect1 = rectangle_vertices(*r1)
    rect2 = rectangle_vertices(*r2)

    # Use the vertices of the first rectangle as
    # starting vertices of the intersection polygon.
    intersection = rect1

    # Loop over the edges of the second rectangle
    for p, q in zip(rect2, rect2[1:] + rect2[:1]):
        if len(intersection) <= 2:
            break # No intersection

        line = Line(p, q)

        # Any point p with line(p) <= 0 is on the "inside" (or on the boundary),
        # any point p with line(p) > 0 is on the "outside".

        # Loop over the edges of the intersection polygon,
        # and determine which part is inside and which is outside.
        new_intersection = []
        line_values = [line(t) for t in intersection]
        for s, t, s_value, t_value in zip(
            intersection, intersection[1:] + intersection[:1],
            line_values, line_values[1:] + line_values[:1]):
            if s_value <= 0:
                new_intersection.append(s)
            if s_value * t_value < 0:
                # Points are on opposite sides.
                # Add the intersection of the lines to new_intersection.
                intersection_point = line.intersection(Line(s, t))
                new_intersection.append(intersection_point)

        intersection = new_intersection

    # Calculate area
    if len(intersection) <= 2:
        return 0

    return 0.5 * sum(p.x*q.y - p.y*q.x for p, q in
                     zip(intersection, intersection[1:] + intersection[:1]))


if __name__ == '__main__':
    r1 = (10, 15, 15, 10, 30)
    r2 = (15, 15, 20, 10, 0)
    print(intersection_area(r1, r2))

- Ruud de Jong

---

难怪你是计算几何书的合著者之一 Mark de Berg、Marc van Kreveld 和 Mark Overmars 的同胞。我有一种感觉，你可能亲自认识他们中的一些人。 - Leon

顺便提一下，在你的回答中应该注意到它适用于计算任何多边形与任何凸多边形的交集。 - Leon

谢谢您的回答，非常清晰和有用！ Line.__call__() 中发生的操作是否有数学名称？它似乎与这个有关？https://math.stackexchange.com/questions/274712/calculate-on-which-side-of-a-straight-line-is-a-given-point-located - dgoldman

Answer 3

intersection, pnt = contourIntersection(rect1, rect2)

在查看了这个问题可能的重复页面后，我没有找到一个完整的Python解决方案，所以这里是我的解决方案，使用遮罩。这个函数可以处理任何角度上的复杂形状，而不仅仅是矩形。

您将旋转矩形的2个轮廓作为参数传入，如果没有交集，则返回“None”，否则返回相交区域的图像以及该图像相对于原始图像轮廓位置的左/上位置。

使用Python、CV2和NumPy。

import cv2
import math
import numpy as np


def contourIntersection(con1, con2, showContours=False):

    # skip if no bounding rect intersection
    leftmost1 = tuple(con1[con1[:, :, 0].argmin()][0])
    topmost1 = tuple(con1[con1[:, :, 1].argmin()][0])
    leftmost2 = tuple(con2[con2[:, :, 0].argmin()][0])
    topmost2 = tuple(con2[con2[:, :, 1].argmin()][0])

    rightmost1 = tuple(con1[con1[:, :, 0].argmax()][0])
    bottommost1 = tuple(con1[con1[:, :, 1].argmax()][0])
    rightmost2 = tuple(con2[con2[:, :, 0].argmax()][0])
    bottommost2 = tuple(con2[con2[:, :, 1].argmax()][0])

    if rightmost1[0] < leftmost2[0] or rightmost2[0] < leftmost1[0] or bottommost1[1] < topmost2[1] or bottommost2[1] < topmost1[1]:
        return None, None

    # reset top / left to 0
    left = leftmost1[0] if leftmost1[0] < leftmost2[0] else leftmost2[0]
    top = topmost1[1] if topmost1[1] < topmost2[1] else topmost2[1]

    newCon1 = []
    for pnt in con1:

        newLeft = pnt[0][0] - left
        newTop = pnt[0][1] - top

        newCon1.append([newLeft, newTop])
    # next
    con1_new = np.array([newCon1], dtype=np.int32)

    newCon2 = []
    for pnt in con2:

        newLeft = pnt[0][0] - left
        newTop = pnt[0][1] - top

        newCon2.append([newLeft, newTop])
    # next
    con2_new = np.array([newCon2], dtype=np.int32)

    # width / height
    right1 = rightmost1[0] - left
    bottom1 = bottommost1[1] - top
    right2 = rightmost2[0] - left
    bottom2 = bottommost2[1] - top

    width = right1 if right1 > right2 else right2
    height = bottom1 if bottom1 > bottom2 else bottom2

    # create images
    img1 = np.zeros([height, width], np.uint8)
    cv2.drawContours(img1, con1_new, -1, (255, 255, 255), -1)

    img2 = np.zeros([height, width], np.uint8)
    cv2.drawContours(img2, con2_new, -1, (255, 255, 255), -1)

    # mask images together using AND
    imgIntersection = cv2.bitwise_and(img1, img2)

    if showContours:
        img1[img1 > 254] = 128
        img2[img2 > 254] = 100

        imgAll = cv2.bitwise_or(img1, img2)
        cv2.imshow('Merged Images', imgAll)

    # end if

    if not imgIntersection.sum():
        return None, None

    # trim
    while not imgIntersection[0].sum():
        imgIntersection = np.delete(imgIntersection, (0), axis=0)
        top += 1

    while not imgIntersection[-1].sum():
        imgIntersection = np.delete(imgIntersection, (-1), axis=0)

    while not imgIntersection[:, 0].sum():
        imgIntersection = np.delete(imgIntersection, (0), axis=1)
        left += 1

    while not imgIntersection[:, -1].sum():
        imgIntersection = np.delete(imgIntersection, (-1), axis=1)

    return imgIntersection, (left, top)
# end function

为了完成答案，以便您可以使用上述函数与2个旋转矩形的CenterX、CenterY、Width、Height和Angle值，我添加了以下功能。只需要简单地更改代码底部的Rect1和Rect2属性即可适用于您自己的情况。

def pixelsBetweenPoints(xy1, xy2):
    X = abs(xy1[0] - xy2[0])
    Y = abs(xy1[1] - xy2[1])

    return int(math.sqrt((X ** 2) + (Y ** 2)))
# end function


def rotatePoint(angle, centerPoint, dist):
    xRatio = math.cos(math.radians(angle))
    yRatio = math.sin(math.radians(angle))
    xPotted = int(centerPoint[0] + (dist * xRatio))
    yPlotted = int(centerPoint[1] + (dist * yRatio))
    newPoint = [xPotted, yPlotted]

    return newPoint
# end function


def angleBetweenPoints(pnt1, pnt2):
    A_B = pixelsBetweenPoints(pnt1, pnt2)

    pnt3 = (pnt1[0] + A_B, pnt1[1])
    C = pixelsBetweenPoints(pnt2, pnt3)

    angle = math.degrees(math.acos((A_B * A_B + A_B * A_B - C * C) / (2.0 * A_B * A_B)))

    # reverse if above horizon
    if pnt2[1] < pnt1[1]:
        angle = angle * -1
    # end if

    return angle
# end function


def rotateRectContour(xCenter, yCenter, height, width, angle):
    # calc positions
    top = int(yCenter - (height / 2))
    left = int(xCenter - (width / 2))
    right = left + width

    rightTop = (right, top)
    centerPoint = (xCenter, yCenter)

    # new right / top point
    rectAngle = angleBetweenPoints(centerPoint, rightTop)
    angleRightTop = angle + rectAngle
    angleRightBottom = angle + 180 - rectAngle
    angleLeftBottom = angle + 180 + rectAngle
    angleLeftTop = angle - rectAngle

    distance = pixelsBetweenPoints(centerPoint, rightTop)
    rightTop_new = rotatePoint(angleRightTop, centerPoint, distance)
    rightBottom_new = rotatePoint(angleRightBottom, centerPoint, distance)
    leftBottom_new = rotatePoint(angleLeftBottom, centerPoint, distance)
    leftTop_new = rotatePoint(angleLeftTop, centerPoint, distance)

    contourList = [[leftTop_new], [rightTop_new], [rightBottom_new], [leftBottom_new]]
    contour = np.array(contourList, dtype=np.int32)

    return contour
# end function


# rect1
xCenter_1 = 40
yCenter_1 = 20
height_1 = 200
width_1 = 80
angle_1 = 45

rect1 = rotateRectContour(xCenter_1, yCenter_1, height_1, width_1, angle_1)

# rect2
xCenter_2 = 80
yCenter_2 = 25
height_2 = 180
width_2 = 50
angle_2 = 123

rect2 = rotateRectContour(xCenter_2, yCenter_2, height_2, width_2, angle_2)

intersection, pnt = contourIntersection(rect1, rect2, True)

if intersection is None:
    print('No intersection')
else:
    print('Area of intersection = ' + str(int(intersection.sum() / 255)))
    cv2.imshow('Intersection', intersection)
# end if

cv2.waitKey(0)

- Bug

---

这种解决方案的缺点是它通过光栅化计算结果，这会引入不准确性，速度较慢，并且结果可能无法直接使用（如果我们需要将其作为多边形而不是图像）。请参见我的答案以获取更好的选择。 - Leon

正在尝试计算值而不是形状。同时，考虑到轮廓可能是从图像中提取的，因此返回了一个结果，允许将其映射到交叉像素的个别交点。只有一行代码可以将图像转换为多边形，不想使代码混乱。 - Bug

在上述代码中使用time.time()进行运行时间测试并关闭所有输出屏幕，它的运行时间比time.time()能记录的最小时间还要短。如果没有可能发生交集时，还可以避免不必要的处理。 - Bug