python:帮助实现一个算法来找到给定点的最小面积矩形,以便计算主轴和副轴的长度 [英] python: Help to implement an algorithm to find the minimum-area-rectangle for given points in order to compute the major and minor axis length

查看:216
本文介绍了python:帮助实现一个算法来找到给定点的最小面积矩形,以便计算主轴和副轴的长度的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

我有一组由多边形(红色)的凸包(蓝色)派生的点(地理坐标值中的黑点)。见图:

  [(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)]
code>

我需要计算这些步骤后的主轴和副轴长度(形式为或查看相当紧凑的

下面的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()来推进每个索引。 mx,my 参数指向 mostfar()告诉它要测试什么极端;作为一个例子,使用 mx,my = -1,0 mostfar()将尝试最大化-rx(其中rx是点的旋转x),从而找到最左边的点。请注意,如果在不精确的算术中完成mx * rx + my * ry> = best ,那么在时应该使用epsilon allowance:当船体有很多点时,舍入误差可能是一个问题,并导致该方法错误地不提前索引。



代码如下所示。船体数据取自上面的问题,不相关的大偏移量和相同的小数位被删除。

 #!/ usr / bin / python 
进口数学

hull = [(23.45,57.39),(23.45,60.39),(24.45,63.39),
(26.95,68.39),(28.45, (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
best = mx * rx + my * ry
而真:
x,y = rx,ry
xn,yn =船体[(j + 1)%n] [0],船体[( 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(船体)
iL = iR = iP = 1#指数向左,向右,相反
pi = 4 * math.atan(1)
对于范围(n-1)中的i:
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 = 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-格式(i,iL,iP,iR)*(xR-xL)

print'{:2d} {:2d} {:2d} {:2d} {:9.3f}'。 ,area)

注意:要获取最小值的长度和宽度-area包围矩形,修改上面的代码,如下所示。这将产生一个输出行,如

  Min矩形:187.451 18.037 10.393 10 0 4 7 

其中第二个和第三个数字表示矩形的长度和宽度,并且这四个整数给出了位于

 #在pi = ...行后加上:
minRect =(1e33,0,0,0 ,0,0,0)#area,dx,dy,i,iL,iP,iR

#在area = ...之后添加:line:
如果area < minRect [0]:
minRect =(area,xR-xL,yP-yC,i,iL,iP,iR)

#在打印后添加...行:
打印'Min rectangle:',minRect
#或者代替那个打印,添加:
打印'Min:矩形:',
for ['{:3d}'.format x)如果isinstance(x,int)else'{:7.3f}'.format(x)for x in minRect]:
print x,
print


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

  1. Compute the convex hull of the cloud.
  2. 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,
  3. 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

这篇关于python:帮助实现一个算法来找到给定点的最小面积矩形,以便计算主轴和副轴的长度的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆