Gianni Spear - 2 months ago 20

Python Question

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

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.

Answer

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
```