python - 找到给定点的最小面积矩形以计算长轴和短轴长度的算法

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)]

我需要按照这些步骤计算长轴和短轴长度(形成这个 post 写在 R-project 和 Java 中)或按照 this example procedure

1. 计算云的凸包。
2. 对于凸包的每条边:2a.计算边缘方向，2b。使用该方向旋转凸包，以便轻松计算具有旋转凸包的 x/y 的 min/max 的边界矩形区域，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 的边界矩形。

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 个随机点并找到了边界框。

(由@heltonbiker 编辑)一个简单的绘图代码:

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

Link to the same answer over on gis.stackexchange供我自己引用。