Title: Convex hull and diameter of 2d point sets
Submitter: David Eppstein (other recipes)
Last Updated: 2003/10/16
Version no: 1.1
Category: Algorithms

Description:

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.

Source: Text Source

# 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

Discussion:

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.