You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
97 lines
3.0 KiB
Python
97 lines
3.0 KiB
Python
11 years ago
|
# Example from http://wiki.scipy.org/Cookbook/Finding_Convex_Hull
|
||
|
|
|
||
|
|
import numpy as n
|
||
|
|
import pylab as p
|
||
|
|
import time
|
||
|
|
|
||
|
|
|
||
|
|
def _angle_to_point(point, centre):
|
||
|
|
'''calculate angle in 2-D between points and x axis'''
|
||
|
|
delta = point - centre
|
||
|
|
res = n.arctan(delta[1] / delta[0])
|
||
|
|
if delta[0] < 0:
|
||
|
|
res += n.pi
|
||
|
|
return res
|
||
|
|
|
||
|
|
|
||
|
|
def _draw_triangle(p1, p2, p3, **kwargs):
|
||
|
|
tmp = n.vstack((p1, p2, p3))
|
||
|
|
x, y = [x[0] for x in zip(tmp.transpose())]
|
||
|
|
p.fill(x, y, **kwargs)
|
||
|
|
# XXX time.sleep(0.2)
|
||
|
|
|
||
|
|
|
||
|
|
def area_of_triangle(p1, p2, p3):
|
||
|
|
'''calculate area of any triangle given co-ordinates of the corners'''
|
||
|
|
return n.linalg.norm(n.cross((p2 - p1), (p3 - p1)))/2.
|
||
|
|
|
||
|
|
|
||
|
|
def convex_hull(points, graphic=True, smidgen=0.0075):
|
||
|
|
'''Calculate subset of points that make a convex hull around points
|
||
|
|
|
||
|
|
Recursively eliminates points that lie inside two neighbouring points until
|
||
|
|
only convex hull is remaining.
|
||
|
|
|
||
|
|
:Parameters:
|
||
|
|
points : ndarray (2 x m)
|
||
|
|
array of points for which to find hull
|
||
|
|
graphic : bool
|
||
|
|
use pylab to show progress?
|
||
|
|
smidgen : float
|
||
|
|
offset for graphic number labels - useful values depend on your data
|
||
|
|
range
|
||
|
|
|
||
|
|
:Returns:
|
||
|
|
hull_points : ndarray (2 x n)
|
||
|
|
convex hull surrounding points
|
||
|
|
'''
|
||
|
|
if graphic:
|
||
|
|
p.clf()
|
||
|
|
p.plot(points[0], points[1], 'ro')
|
||
|
|
n_pts = points.shape[1]
|
||
|
|
assert(n_pts > 5)
|
||
|
|
centre = points.mean(1)
|
||
|
|
if graphic:
|
||
|
|
p.plot((centre[0], ), (centre[1], ), 'bo')
|
||
|
|
angles = n.apply_along_axis(_angle_to_point, 0, points, centre)
|
||
|
|
pts_ord = points[:, angles.argsort()]
|
||
|
|
if graphic:
|
||
|
|
for i in xrange(n_pts):
|
||
|
|
p.text(pts_ord[0, i] + smidgen, pts_ord[1, i] + smidgen,
|
||
|
|
'%d' % i)
|
||
|
|
pts = [x[0] for x in zip(pts_ord.transpose())]
|
||
|
|
prev_pts = len(pts) + 1
|
||
|
|
k = 0
|
||
|
|
while prev_pts > n_pts:
|
||
|
|
prev_pts = n_pts
|
||
|
|
n_pts = len(pts)
|
||
|
|
if graphic:
|
||
|
|
p.gca().patches = []
|
||
|
|
i = -2
|
||
|
|
while i < (n_pts - 2):
|
||
|
|
Aij = area_of_triangle(centre, pts[i], pts[(i + 1) % n_pts])
|
||
|
|
Ajk = area_of_triangle(centre, pts[(i + 1) % n_pts],
|
||
|
|
pts[(i + 2) % n_pts])
|
||
|
|
Aik = area_of_triangle(centre, pts[i], pts[(i + 2) % n_pts])
|
||
|
|
if graphic:
|
||
|
|
_draw_triangle(centre, pts[i], pts[(i + 1) % n_pts],
|
||
|
|
facecolor='blue', alpha=0.2)
|
||
|
|
_draw_triangle(centre, pts[(i + 1) % n_pts],
|
||
|
|
pts[(i + 2) % n_pts],
|
||
|
|
facecolor='green', alpha=0.2)
|
||
|
|
_draw_triangle(centre, pts[i], pts[(i + 2) % n_pts],
|
||
|
|
facecolor='red', alpha=0.2)
|
||
|
|
if Aij + Ajk < Aik:
|
||
|
|
if graphic:
|
||
|
|
p.plot((pts[i + 1][0], ), (pts[i + 1][1], ), 'go')
|
||
|
|
del pts[i+1]
|
||
|
|
i += 1
|
||
|
|
n_pts = len(pts)
|
||
|
|
k += 1
|
||
|
|
return n.asarray(pts)
|
||
|
|
|
||
|
|
if __name__ == "__main__":
|
||
|
|
points = n.random.random_sample((2, 40))
|
||
|
|
hull_pts = convex_hull(points)
|
||
|
|
p.show()
|