找到给定点的最小面积矩形以计算长轴和短轴长度的算法
问题描述
我有一组点(地理坐标值中的黑点)来自多边形(红色)的凸包(蓝色).见图:
[(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)]我需要按照这些步骤计算长轴和短轴的长度(形成这个 或查看相当紧凑的 格雷厄姆在 tixxit.net 上扫码.)
以下 python 程序使用类似于通常 O(n) 算法的技术来计算凸多边形的最大直径.也就是说，它维护三个索引(iL、iP、iR)到相对于给定基线的最左边、对面和最右边的点.每个索引最多前进 n 个点.程序的示例输出如下所示(添加了标题):
iL iP iR 区域0 6 8 0 203.0001 6 8 0 211.8752 6 8 0 205.8003 6 10 0 206.2504 7 12 0 190.3625 8 0 1 203.0006 10 0 4 201.3857 0 1 6 203.0008 0 3 6 205.8279 0 3 6 205.64010 0 4 7 187.45111 0 4 7 189.75012 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导入数学船体 = [(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): # 推进 j 到极值点xn, yn = 船体[j][0], 船体[j][1]rx, ry = xn*c - yn*s, xn*s + yn*c最佳 = mx*rx + my*ry而真:x, y = rx, ryxn, yn = 船体[(j+1)%n][0], 船体[(j+1)%n][1]rx, ry = xn*c - yn*s, xn*s + yn*c如果 mx*rx + my*ry >= 最好:j = (j+1)%n最佳 = mx*rx + my*ry别的:返回 (x, y, j)n = len(船体)iL = iR = iP = 1 # 索引左、右、相反pi = 4*math.atan(1)对于范围内的 i (n-1):dx = 船体[i+1][0] - 船体[i][0]dy = 船体[i+1][1] - 船体[i][1]theta = pi-math.atan2(dy, dx)s, c = math.sin(theta), math.cos(theta)yC = 船体[i][0]*s + 船体[i][1]*cxP, yP, iP = mostfar(iP, n, s, c, 0, 1)如果 i==0:iR = iPxR, yR, iR = mostfar(iR, n, s, c, 1, 0)xL, yL, iL = mostfar(iL, n, s, c, -1, 0)面积 = (yP-yC)*(xR-xL)print ' {:2d} {:2d} {:2d} {:2d} {:9.3f}'.format(i, iL, iP, iR, area)注意:要得到最小面积包围矩形的长度和宽度，修改上面的代码如下所示.这将产生类似
的输出行最小矩形:187.451 18.037 10.393 10 0 4 7其中第二个和第三个数字表示矩形的长度和宽度，四个整数给出位于矩形边上的点的索引号.
# 在 pi = ... 之后添加:minRect = (1e33, 0, 0, 0, 0, 0, 0) # 面积, dx, dy, i, iL, iP, iR# 在 area = ... 行之后添加:如果面积 <最小矩形[0]:minRect = (面积, xR-xL, yP-yC, i, iL, iP, iR)# 在打印后添加 ... 行:打印'最小矩形:'，minRect# 或者代替打印，添加:打印'最小矩形:'，for x in ['{:3d} '.format(x) if isinstance(x, int) else '{:7.3f} '.format(x) for x in minRect]:打印 x，打印I have a set of points (black dots in geographic coordinate value) derived from the convex hull (blue) of a polygon (red). see Figure:
[(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)]
I need to calculate the the major and minor axis length following these steps (form this post write in R-project and in Java) or following the this example procedure
- Compute the convex hull of the cloud.
- For each edge of the convex hull: 2a. compute the edge orientation, 2b. rotate the convex hull using this orientation in order to compute easily the bounding rectangle area with min/max of x/y of the rotated convex hull, 2c. Store the orientation corresponding to the minimum area found,
- Return the rectangle corresponding to the minimum area found.
After that we know the The angle Theta (represented the orientation of the bounding rectangle relative to the y-axis of the image). The minimum and maximum of a and b over all boundary points are found:
- a(xi,yi) = xi*cos Theta + yi sin Theta
- b(xi,yi) = xi*sin Theta + yi cos Theta
The values (a_max - a_min) and (b_max - b_min) defined the length and width, respectively, of the bounding rectangle for a direction Theta.
解决方案Given a clockwise-ordered list of n points in the convex hull of a set of points, it is an O(n) operation to find the minimum-area enclosing rectangle. (For convex-hull finding, in O(n log n) time, see activestate.com recipe 66527 or see the quite compact Graham scan code at tixxit.net.)
The following python program uses techniques similar to those of the usual O(n) algorithm for computing maximum diameter of a convex polygon. That is, it maintains three indexes (iL, iP, iR) to the leftmost, opposite, and rightmost points relative to a given baseline. Each index advances through at most n points. Sample output from the program is shown next (with an added header):
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
For example, the i=10 entry indicates that relative to the baseline from point 10 to 11, point 0 is leftmost, point 4 is opposite, and point 7 is rightmost, yielding an area of 187.451 units.
Note that the code uses mostfar() to advance each index. The mx, my parameters to mostfar() tell it what extreme to test for; as an example, with mx,my = -1,0, mostfar() will try to maximize -rx (where rx is the rotated x of a point), thus finding the leftmost point. Note that an epsilon allowance probably should be used when if mx*rx + my*ry >= best is done in inexact arithmetic: when a hull has numerous points, rounding error may be a problem and cause the method to incorrectly not advance an index.
Code is shown below. The hull data is taken from the question above, with irrelevant large offsets and identical decimal places elided.
#!/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)
Note: To get the length and width of the minimal-area enclosing rectangle, modify the above code as shown below. This will produce an output line like
Min rectangle: 187.451 18.037 10.393 10 0 4 7
in which the second and third numbers indicate the length and width of the rectangle, and the four integers give index numbers of points that lie upon sides of it.
# 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
相关文章