Convex hull

hull.py

import random

# Determinant predicate (line sidedness test)
def det3x3_pt(p, q, r):
    a = (q[0] - p[0], q[1] - p[1])
    b = (r[0] - p[0], r[1] - p[1])
    return a[0]*b[1] - a[1]*b[0]

def convex_hull(points):
    # sorts points by x then y which works for us (we only need the sort by x part here, but it doesn't hurt)
    # NOTE: I'm not handling two points with equal x here. This needs care (which I'm not using right now). :)
    sorted_points = sorted(points)
    upper_contour = []
    lower_contour = []

    for p in sorted_points:
        # reject previous points in the upper contour that produce a concavity with the new point
        # added, until we're convex again
        while len(upper_contour) >= 2 and det3x3_pt(upper_contour[-2], upper_contour[-1], p) > 0:
            upper_contour.pop()

        upper_contour.append(p)

        # analogous for lower contour, which is the same with signs reversed
        while len(lower_contour) >= 2 and det3x3_pt(lower_contour[-2], lower_contour[-1], p) < 0:
            lower_contour.pop()

        # add the new point
        lower_contour.append(p)

    # concatenate the upper contour with the inner portion of the reversed lower contour
    # (since the min/max x points are on both contours)
    return (upper_contour + lower_contour[-2:0:-1])

# Generate a bunch of random points in the plane, at integer coords for simplicity
rnd = random.seed(10)
points = []
count = 20

for i in range(count):
    x = random.randint(0, 99)
    y = random.randint(0, 99)
    points.append((x,y))

hull = convex_hull(points)

with open('points.dat', 'w') as f:
    for p in points:
        print(p[0], p[1], file=f)

with open('hull.dat', 'w') as f:
    for p in hull + [hull[0]]:
        print(p[0], p[1], file=f)

#print(points)
print(convex_hull(points))

hull2.py

# Original version preserved for posterity, but this one tries to be a bit more robust and tries
# to avoid emitting degenerate hull segments. It also starts out with a nicer point
# distribution is more circular so everything doesn't end up looking like a misshapen rectangle.
import random

# Determinant predicate (line sidedness test)
def det3x3_pt(p, q, r):
    a = (q[0] - p[0], q[1] - p[1])
    b = (r[0] - p[0], r[1] - p[1])
    return a[0]*b[1] - a[1]*b[0]

def convex_hull(points):
    # sorts points by x then y which works for us (we only need the sort by x part here, but it doesn't hurt)
    sorted_points = sorted(points)
    upper_contour = []
    lower_contour = []

    for p in sorted_points:
        # Reject duplicate points
        if len(upper_contour) > 0 and upper_contour[-1] == p:
            continue

        # reject previous points in the upper contour that produce a concavity with the new point
        # added, until we're convex again
        while len(upper_contour) >= 2 and det3x3_pt(upper_contour[-2], upper_contour[-1], p) >= 0:
            upper_contour.pop()

        upper_contour.append(p)

        # analogous for lower contour, which is the same with signs reversed
        while len(lower_contour) >= 2 and det3x3_pt(lower_contour[-2], lower_contour[-1], p) <= 0:
            lower_contour.pop()

        # add the new point
        lower_contour.append(p)

    # concatenate the upper contour with the inner portion of the reversed lower contour
    # (since the min/max x points are on both contours)
    return (upper_contour + lower_contour[-2:0:-1])

# Generate a bunch of random points in the plane, at integer coords for simplicity
rnd = random.seed(1)
points = []
count = 50
radius = 20

for i in range(count):
    while True:
        x = random.randint(-radius, radius)
        y = random.randint(-radius, radius)
        if x*x + y*y <= radius*radius:
            break
    points.append((x,y))

hull = convex_hull(points)

with open('points.dat', 'w') as f:
    for p in points:
        print(p[0], p[1], file=f)

with open('hull.dat', 'w') as f:
    for p in hull + [hull[0]]:
        print(p[0], p[1], file=f)

print(hull)