Returns the convex hull (separated into upper and lower chains of vertices) and the diameter (farthest pair of points), given input consisting of a list of 2d points represented as pairs (x,y). The convex hull algorithm is Graham's scan, using a coordinate-based sorted order rather than the more commonly seen radial sorted order. A rotating calipers algorithm generates candidate pairs of vertices for the diameter calculation. Care was taken handling tricky cases such as pairs of points with the same x-coordinate and colinear triples of points.
Download
Copy to clipboard1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 | # convex hull (Graham scan by x-coordinate) and diameter of a set of points
# David Eppstein, UC Irvine, 7 Mar 2002
from __future__ import generators
def orientation(p,q,r):
'''Return positive if p-q-r are clockwise, neg if ccw, zero if colinear.'''
return (q[1]-p[1])*(r[0]-p[0]) - (q[0]-p[0])*(r[1]-p[1])
def hulls(Points):
'''Graham scan to find upper and lower convex hulls of a set of 2d points.'''
U = []
L = []
Points.sort()
for p in Points:
while len(U) > 1 and orientation(U[-2],U[-1],p) <= 0: U.pop()
while len(L) > 1 and orientation(L[-2],L[-1],p) >= 0: L.pop()
U.append(p)
L.append(p)
return U,L
def rotatingCalipers(Points):
'''Given a list of 2d points, finds all ways of sandwiching the points
between two parallel lines that touch one point each, and yields the sequence
of pairs of points touched by each pair of lines.'''
U,L = hulls(Points)
i = 0
j = len(L) - 1
while i < len(U) - 1 or j > 0:
yield U[i],L[j]
# if all the way through one side of hull, advance the other side
if i == len(U) - 1: j -= 1
elif j == 0: i += 1
# still points left on both lists, compare slopes of next hull edges
# being careful to avoid divide-by-zero in slope calculation
elif (U[i+1][1]-U[i][1])*(L[j][0]-L[j-1][0]) > \
(L[j][1]-L[j-1][1])*(U[i+1][0]-U[i][0]):
i += 1
else: j -= 1
def diameter(Points):
'''Given a list of 2d points, returns the pair that's farthest apart.'''
diam,pair = max([((p[0]-q[0])**2 + (p[1]-q[1])**2, (p,q))
for p,q in rotatingCalipers(Points)])
return pair
|
The convex hull part of the recipe is essentially the same as the one in http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/66527 although (unlike that one) this can handle multiple copies of the same point without getting an assertion failure. The orientation test has also been somewhat simplified.
This revision replaces a loop in the diameter routine with a list comprehension.
I really like this use of 'yield' I find this code to be especially readable and elegant with its use of 'yield'.
Complex arithmetic. Some of the formulas become a little simpler if we represent points by complex numbers as pairs of reals:
Of course, under the hood, the complex version is doing more work: it's finding the real as well as imaginary components in the first and second formula, and doing an unnecessary square root in the third one.
Note. If you need a result as a single sequence of points, you can join the lists (the first of L is the last of U):
U[:-1] + L
Not Graham's algorithm at all. This is Andrew's monotone chain algorithm, from the paper "Another Efficient Algorithm for Convex Hulls in Two Dimensions". The lexicographic sort and generation of upper and lower portions of the hull are features of Andrew's algorithm that are not present in Graham's.