求给定点最小面积矩形以计算长轴和短轴长度的算法

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中填写此post）或在this example procedure之后计算长轴和短轴长度

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

然后我们知道角度θ（表示边框相对于图像y轴的方向）。所有边界点上a和b的最小值和最大值为 找到：

• a（xi，yi）=xi*cosθ+yi sinθ
• b（xi，yi）=xi*sinθ+yi cosθ

值（a_max-a_min）和（b_max-b_min）分别定义了长度和宽度， θ方向的边界矩形。

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供我参考。