寻找给定点的最小矩形算法，以计算主轴和次轴长度。

Question

我有一组点（地理坐标值中的黑点），这些点来自多边形（红色）的凸包（蓝色）。请参见图：

[(560023.44957588764,6362057.3904932579), 
 (560023.44957588764,6362060.3904932579), 
 (560024.44957588764,6362063.3904932579), 
 (560026.94957588764,6362068.3904932579), 
 (560028.44957588764,6362069.8904932579), 
 (560034.94957588764,6362071.8904932579), 
 (560036.44957588764,6362071.8904932579), 
 (560037.44957588764,6362070.3904932579), 
 (560037.44957588764,6362064.8904932579), 
 (560036.44957588764,6362063.3904932579), 
 (560034.94957588764,6362061.3904932579), 
 (560026.94957588764,6362057.8904932579), 
 (560025.44957588764,6362057.3904932579), 
 (560023.44957588764,6362057.3904932579)]

我需要按照以下步骤计算主轴和次轴长度（来自这个帖子，用R-project和Java编写，或者遵循这个示例过程)：

1. 计算云的凸包。
2. 对于凸包的每条边： 2a. 计算边缘方向， 2b. 使用该方向旋转凸包，以便使用旋转后的凸包的x/y最小/最大轴易于计算边界矩形面积， 2c. 存储对应于找到的最小面积的方向，
3. 返回与找到的最小面积相对应的矩形。

之后我们知道角度Theta（表示边界矩形相对于图像y轴的方向）。所有边界点上的a和b的最小值和最大值均已找到：

• a(xi，yi)=xi*cos Theta + yi sin Theta
• b(xi，yi)=xi*sin Theta + yi cos Theta

(a_max - a_min)和(b_max - b_min)的值定义了边界矩形在Theta方向上的长度和宽度。

- Gianni Spear

---

你已经找到了算法 - 你有什么问题？ - Eric

所以你的问题是：“是否存在一个已经能够实现这个功能的模块？” - Eric

@scravy：不，这只是工作的一小部分。我的目标是通过景观竞赛中的每个多边形提取几个指标（例如：凸度、固实度、圆度等）。之后我需要训练机器学习来对每个多边形进行分类。主要目标是了解我们是否正在失去土地以及如何保护自然。 - Gianni Spear

在http://gis.stackexchange.com/questions/22895/how-to-find-the-minimum-area-rectangle-for-given-points上发布了交叉帖。 - whuber

我已经为这个问题创建并测试了一个程序包。文章链接。仓库链接。 - William Rusnack

Answer 1

我刚刚自己实现了这个功能，所以我想在这里分享我的版本供其他人查看：

import numpy as np
from scipy.spatial import ConvexHull

def minimum_bounding_rectangle(points):
    """
    Find the smallest bounding rectangle for a set of points.
    Returns a set of points representing the corners of the bounding box.

    :param points: an nx2 matrix of coordinates
    :rval: an nx2 matrix of coordinates
    """
    from scipy.ndimage.interpolation import rotate
    pi2 = np.pi/2.

    # get the convex hull for the points
    hull_points = points[ConvexHull(points).vertices]

    # calculate edge angles
    edges = np.zeros((len(hull_points)-1, 2))
    edges = hull_points[1:] - hull_points[:-1]

    angles = np.zeros((len(edges)))
    angles = np.arctan2(edges[:, 1], edges[:, 0])

    angles = np.abs(np.mod(angles, pi2))
    angles = np.unique(angles)

    # find rotation matrices
    # XXX both work
    rotations = np.vstack([
        np.cos(angles),
        np.cos(angles-pi2),
        np.cos(angles+pi2),
        np.cos(angles)]).T
#     rotations = np.vstack([
#         np.cos(angles),
#         -np.sin(angles),
#         np.sin(angles),
#         np.cos(angles)]).T
    rotations = rotations.reshape((-1, 2, 2))

    # apply rotations to the hull
    rot_points = np.dot(rotations, hull_points.T)

    # find the bounding points
    min_x = np.nanmin(rot_points[:, 0], axis=1)
    max_x = np.nanmax(rot_points[:, 0], axis=1)
    min_y = np.nanmin(rot_points[:, 1], axis=1)
    max_y = np.nanmax(rot_points[:, 1], axis=1)

    # find the box with the best area
    areas = (max_x - min_x) * (max_y - min_y)
    best_idx = np.argmin(areas)

    # return the best box
    x1 = max_x[best_idx]
    x2 = min_x[best_idx]
    y1 = max_y[best_idx]
    y2 = min_y[best_idx]
    r = rotations[best_idx]

    rval = np.zeros((4, 2))
    rval[0] = np.dot([x1, y2], r)
    rval[1] = np.dot([x2, y2], r)
    rval[2] = np.dot([x2, y1], r)
    rval[3] = np.dot([x1, y1], r)

    return rval

以下是四个不同的示例，每个示例我生成了4个随机点并找到了边界框。

简单的绘图代码：

import matplotlib.pyplot as plt
for n in range(10):
    points = np.random.rand(4,2)
    plt.scatter(points[:,0], points[:,1])
    bbox = minimum_bounding_rectangle(points)
    plt.fill(bbox[:,0], bbox[:,1], alpha=0.2)
    plt.axis('equal')
    plt.show()

这些样本在4个点上的处理速度相对较快：

>>> %timeit minimum_bounding_rectangle(a)
1000 loops, best of 3: 245 µs per loop

在gis.stackexchange上查看相同答案的链接，以供参考。

- JesseBuesking

---

了不起的解决方案，功能非常准确！对我来说，需要进行一些微小的源代码调整才能让它顺利运行 - 将 hull 改为 hull_points ，并加入 rval = np.zeros((4,2))。非常感谢你的分享。 - Serj Zaharchenko

@SerjZaharchenko 对不起！在我最后一次编辑中，我在代码中使用了两个不同的变量名。我刚刚根据你提到的修复更新了它，并更新了示例图片。谢谢你的提醒！ - JesseBuesking

Answer 2

在Github上已经有一个可以实现此功能的模块了。 https://github.com/BebeSparkelSparkel/MinimumBoundingBox

你只需要将你的点云插入其中即可。

from MinimumBoundingBox import minimum_bounding_box
points = ( (1,2), (5,4), (-1,-3) )
bounding_box = minimum_bounding_box(points)  # returns namedtuple

你可以通过以下方式获得椭圆的主轴和次轴长度：

minor = min(bounding_box.length_parallel, bounding_box.length_orthogonal)
major = max(bounding_box.length_parallel, bounding_box.length_orthogonal)

它还会返回面积、矩形中心、矩形角度和角点。

- William Rusnack

---

我认为这个答案应该得到更多的赞。它对我来说开箱即用，效果很好。 - biomiker

谢谢，我也这么想！ - William Rusnack

非常不错。你的.py文件省下了我很多时间。直接用就可以了。50,000个坐标点时速度还挺快。 - EA304GT

请注意，需要检查corner_points的顺序，以确定角点的顺时针或逆时针方向。 - undefined

Answer 3

给定凸包中的n个点的顺时针排序列表，查找最小面积包围矩形是一个O(n)操作。(对于凸包查找，可以在O(n log n)时间内查找，请参见activestate.com recipe 66527或查看非常紧凑的tixxit.net上Graham扫描代码。)

下面的Python程序使用类似于计算凸多边形最大直径的通常O(n)算法的技术。也就是说，它维护三个索引(iL、iP、iR)，相对于给定基线，分别指向最左侧、对面和最右侧的点。每个索引最多经过n个点。程序的示例输出如下所示(附加标题):

 i iL iP iR    Area
 0  6  8  0   203.000
 1  6  8  0   211.875
 2  6  8  0   205.800
 3  6 10  0   206.250
 4  7 12  0   190.362
 5  8  0  1   203.000
 6 10  0  4   201.385
 7  0  1  6   203.000
 8  0  3  6   205.827
 9  0  3  6   205.640
10  0  4  7   187.451
11  0  4  7   189.750
12  1  6  8   203.000

例如，i=10的条目表示相对于从点10到11的基准线，点0最左，点4相对，点7最右，产生一个面积为187.451单位。

请注意，代码使用mostfar()来推进每个索引。 mostfar()的mx，my参数告诉它要测试哪个极端值；例如，使用mx，my = -1,0，mostfar()将尝试最大化-rx（其中rx是点的旋转x），从而找到最左边的点。 请注意，在不精确计算时，if mx * rx + my * ry >= best应该使用epsilon允许：当外壳有众多点时，四舍五入误差可能会成为问题，并导致该方法不正确地不推进索引。

下面显示代码。 外壳数据取自上面的问题，省略了不相关的大偏移和相同的小数位数。

#!/usr/bin/python
import math

hull = [(23.45, 57.39), (23.45, 60.39), (24.45, 63.39),
        (26.95, 68.39), (28.45, 69.89), (34.95, 71.89),
        (36.45, 71.89), (37.45, 70.39), (37.45, 64.89),
        (36.45, 63.39), (34.95, 61.39), (26.95, 57.89),
        (25.45, 57.39), (23.45, 57.39)]

def mostfar(j, n, s, c, mx, my): # advance j to extreme point
    xn, yn = hull[j][0], hull[j][1]
    rx, ry = xn*c - yn*s, xn*s + yn*c
    best = mx*rx + my*ry
    while True:
        x, y = rx, ry
        xn, yn = hull[(j+1)%n][0], hull[(j+1)%n][1]
        rx, ry = xn*c - yn*s, xn*s + yn*c
        if mx*rx + my*ry >= best:
            j = (j+1)%n
            best = mx*rx + my*ry
        else:
            return (x, y, j)

n = len(hull)
iL = iR = iP = 1                # indexes left, right, opposite
pi = 4*math.atan(1)
for i in range(n-1):
    dx = hull[i+1][0] - hull[i][0]
    dy = hull[i+1][1] - hull[i][1]
    theta = pi-math.atan2(dy, dx)
    s, c = math.sin(theta), math.cos(theta)
    yC = hull[i][0]*s + hull[i][1]*c

    xP, yP, iP = mostfar(iP, n, s, c, 0, 1)
    if i==0: iR = iP
    xR, yR, iR = mostfar(iR, n, s, c,  1, 0)
    xL, yL, iL = mostfar(iL, n, s, c, -1, 0)
    area = (yP-yC)*(xR-xL)

    print '    {:2d} {:2d} {:2d} {:2d} {:9.3f}'.format(i, iL, iP, iR, area)

注意：要获取最小面积包围矩形的长度和宽度，请按以下方式修改上述代码。这将生成一个输出行，如下所示。

Min rectangle:  187.451   18.037   10.393   10    0    4    7

其中第二个和第三个数字表示矩形的长度和宽度，四个整数给出了位于其边上的点的索引号。

# add after pi = ... line:
minRect = (1e33, 0, 0, 0, 0, 0, 0) # area, dx, dy, i, iL, iP, iR

# add after area = ... line:
    if area < minRect[0]:
        minRect = (area, xR-xL, yP-yC, i, iL, iP, iR)

# add after print ... line:
print 'Min rectangle:', minRect
# or instead of that print, add:
print 'Min rectangle: ',
for x in ['{:3d} '.format(x) if isinstance(x, int) else '{:7.3f} '.format(x) for x in minRect]:
    print x,
print

- James Waldby - jwpat7

---

亲爱的jwpat7，首先感谢您，您真是个好人。我正在编写您的函数（并将在您的名字下发布，因为您是作者），以便拥有一个返回给定点列表的长度、宽度和面积的单个函数：L，W，A = get_minimum_area_rectangle(hull)。我希望分享这种方法，因为在5天内我没有找到这个在Python中实现的函数。 - Gianni Spear

我添加了一个关于最小矩形长度和宽度的注释。将代码修改为函数以删除所有打印并返回最小矩形的长度和宽度非常简单。在UIUC的方法中，连续的包围矩形有几个共同点。UIUC代码分别处理了几种情况，但我不确定这是否必要或可取。请随意分享代码--并接受答案 :)。 - James Waldby - jwpat7

亲爱的jwpat7，这里有一个新的SO帖子，我在其中报告了您是此解决方案的作者。http://stackoverflow.com/questions/13553884/python-implement-a-script-in-a-function-some-suggestions - Gianni Spear

我会写一个 def get_minimum_area_rectangle(hull) 函数 [其中 hull 是指定方向的凸多边形]，并在该函数内定义 mostfar()。不要让 mostfar() 在 get_minimum_area_rectangle() 函数之外可见。请注意，我已经使用另一个 CW 多边形 hull = [(0,1), (0,2), (1,3), (2,3), (3,2), (3,1), (2,0), (1,0), (0,1)] 测试了代码，但没有测试任何 CCW 方向的多边形；我想它会因为 CW 方向的多边形上的符号相反而失败，但这可能是错误的。 - James Waldby - jwpat7

如果最小面积矩形的基地恰好是hull [n]和hull [0]之间的边，则此算法将失败。一个简单的解决方案是用 for i in range(n): 替换 for i in range(n-1):，并用dx = hull[(i+1)%n, 0] - hull[i, 0] (相应地) 替换 dx = ...。至少对于像 hull = np.array([[0,0],[1,2],[4,2.2],[7,2],[8,0]], np.float32) 这样的情况，这种情况是适用的 (而且解决方案有效)。 - Raketenolli

Answer 4

OpenCV 中有这个功能。请参考下面的链接：

http://docs.opencv.org/trunk/dd/d49/tutorial_py_contour_features.html

7.b. 旋转矩形

使用

cv2.minAreaRect(cnt)

您可以获取矩形的长度、宽度和角度。如果需要绘制矩形，您还可以计算其各个角落。

- Trevor

Answer 5

我找到了计算凸包的方法。

如果我们谈论“完整解决方案”（一个函数完成所有工作），我只找到了arcpy，它是ArcGIS程序的一部分。它提供了MinimumBoundingGeometry_management函数，看起来就像你要找的东西。但它不是开源的。不幸的是，缺乏Python开源GIS库。

- Chris Medrela

---

谢谢。我知道arcpy模块，但我更喜欢不使用非开源（Python哲学）的工具。顺便说一句，奇怪的是没有一个模块能够从一组点中计算这些指数。 - Gianni Spear

Answer 6

原帖发布于2013年2月：

上面的代码示例不够健壮。 我使用真实数据（许多点的凸包）进行测试，它产生了接近正确的结果。然而，对于简单的4至6边形无法正常工作。

这里提供了一个自包含的解决方案，用Python编写： https://github.com/dbworth/minimum-area-bounding-rectangle/

使用qhull（QuickHull）算法的2D实现来查找凸包。 解决方案是计算多边形所有边的角度，然后仅在旋转的第一象限（90度）中唯一的那些角度上操作。 在找到最小面积的边界矩形后，它输出所有数据，包括中心点和角点。 提供了一个简单的测试程序。使用Matlab验证了答案。

请注意，该解决方案依赖于一个假设，即边界框将与输入多边形共享至少一条边。这适用于我的应用程序，但我听说有一篇论文作者展示了一些罕见的解决方案，这并不成立。如果这个结果很重要，你应该研究一下！

- David