From 6ca9e985148ab8f017f82b9b3b851225161bbab4 Mon Sep 17 00:00:00 2001 From: Johann Dreo Date: Thu, 3 Apr 2014 22:17:08 +0200 Subject: [PATCH] Adds the Chan's convex hull algorithm --- hull.py | 231 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 231 insertions(+) create mode 100644 hull.py diff --git a/hull.py b/hull.py new file mode 100644 index 0000000..28ba974 --- /dev/null +++ b/hull.py @@ -0,0 +1,231 @@ + +import operator +from utils import x,y,euclidian_distance + +# Based on the excellent article by Tom Switzer +# http://tomswitzer.net/2010/12/2d-convex-hulls-chans-algorithm/ + +TURN_RIGHT, TURN_NONE, TURN_LEFT = (-1, 0, 1) + +def turn(p, q, r): + """Returns -1, 0, 1 if the sequence of points (p,q,r) forms a right, straight, or left turn.""" + qr = ( x(q) - x(p) ) * ( y(r) - y(p) ) + rq = ( x(r) - x(p) ) * ( y(q) - y(p) ) + # cmp(x,y) returns -1 if xy + return cmp( qr - rq, 0) + + +def keep_left(hull, point): + """ Remove the last elements of the given hull that does not form a left turn with the given point.""" + + while len(hull) > 1 and turn( hull[-2], hull[-1], point ) != TURN_LEFT: + hull.pop() + + # if the hull is empty or does not contains the point point + if len(hull) == 0 or hull[-1] != point: + # add at least the asked point + hull.append(point) + + return hull + + +def graham_scan(points): + """Returns points on convex hull of an array of points in counter clockwise order.""" + + # Sort from the furthest left point. + spots = sorted(points) + + # Browse the hull turning left from the furthest left point, + # this is thus the lower part of the convex hull. + lower_hull = reduce(keep_left, spots, []) + + # Going from the furthest right point, we get the upper hull. + upper_hull = reduce(keep_left, reversed(spots), []) + + # Merge the lower and the upper hull, + # omitting the extreme points that are in both hulls parts. + for p in upper_hull[1:-1]: + lower_hull.append( p ) + + return lower_hull + + +def right_tangent(hull, p): + """Return the index of the point in hull that the right tangent line from p to hull touches.""" + + def at( index, h = hull ): + """Secured index accessor to the hull: + if the given index is greater than the container length, then start from the beginning.""" + return index % len(h) + + # Consider the turn formed by a window sliding + # from the extremes points toward the center of the hull. + # + # >> sliding window << + # ___________________________ + # / \ + # hull = [ . . .(. . . . . . . . . . . . . . .). . . ] + # ^ ^ ^ ^ ^ + # | | | | | + # ileft | | icenter+1 iright + # | icenter + # icenter-1 + + # Remember that the points are ordered in the hull. + ileft, iright = 0, len(hull) + # With the last point. + l_prev = turn( p, hull[0], hull[-1] ) + # With the second point (if the hull contains only one point, then ileft stays at zero). + l_next = turn( p, hull[0], hull[ at(ileft+1) ] ) + + # While the right tangent is not found, + # and the sliding window is not null. + while ileft < iright: + # Index of the point in the middle of the window. + icenter = (ileft + iright) / 2 + # Consider the turn formed by the given point, the center point and... + # ... the point before the center, + c_prev = turn( p, hull[icenter], hull[ at(icenter-1) ] ) + # ... the point after the center, + c_next = turn( p, hull[icenter], hull[ at(icenter+1) ] ) + # ... the point on the left of the window. + c_side = turn( p, hull[ileft] , hull[icenter] ) + + # The right tangent is the middle of the window iff + # the turns formed with points around the center are to the LEFT (or straight) + # (i.e. are not to the right). + if c_prev != TURN_RIGHT and c_next != TURN_RIGHT: + return icenter + + # If the tangent touches the left point in the window. + elif ( c_side == TURN_LEFT and ( l_next == TURN_RIGHT or l_prev == l_next) ) \ + or ( c_side == TURN_RIGHT and c_prev == TURN_RIGHT ): + # Do not consider points at the RIGHT of the center. + iright = icenter + + # If the tangent touches the right point in the window. + else: + # Do not consider points at the LEFT of the center. + ileft = icenter + 1 + # Switch sides: if the turn toward the point before the center + # was to the right, search to the left and conversely. + l_prev = -1 * c_next + # Update the turn to the next left point. + l_next = turn( p, hull[ileft], hull[ at(ileft+1) ] ) + + return ileft + + +def min_hull_pt_pair(hulls): + """Returns the (hull, point) index pair that is minimal.""" + # Find an extreme point and the hull chunk to which it is related. + min_h_i, min_p_i = 0, 0 + for i,hull in enumerate(hulls): + # Find the index j of the minimal point in the ith hull + # itemgetter(1) will return the point, because enumerate produces (index,point) pairs + # and thus (index,point)[1] == point + j,pt = min(enumerate(hull), key=operator.itemgetter(1)) + # Minimize across the hulls + if pt < hull[min_p_i]: + min_h_i, min_p_i = i, j + # Return the index of the hull which holds the minimal point and the index of the point itself. + return (min_h_i, min_p_i) + + +def next_hull_pt_pair( hulls, (hi,pi) ): + """ Returns the (hull, point) index pair of the next point in the convex hull.""" + + # The current point itself + base_pt = hulls[hi][pi] + + # Indices of the next (hull,point) pair in hulls + next_hullpt = ( hi, (pi+1)%len(hulls[hi]) ) + + # Now search for a "next" point after the base point + # that forms a right turn along with its tangent point. + + # Loop over the indices of hulls, + # but do not consider the current index. + for nhi in (i for i in range(len(hulls)) if i != hi): + # Index of the right tangent point of the base point + rti = right_tangent( hulls[nhi], base_pt ) + + # The other points themselves + tangent_pt = hulls[nhi][rti] + next_pt = hulls[ next_hullpt[0] ][ next_hullpt[1] ] + + # How is the turn formed by the three points? + t = turn( base_pt, next_pt, tangent_pt ) + + # If it forms a right turn, it is on the convex hull. + # If it forms a left turn, it is not and we continue the loop. + # In the (rare) case the points are aligned, we consider the next point + # only if it is closer to the base point than the tangent point. + if t == TURN_RIGHT \ + or ( t == TURN_NONE and euclidian_distance(base_pt, next_pt) < euclidian_distance(base_pt, tangent_pt) ): + # save the indices of the hull and point + next_hullpt = (nhi, rti) + + return next_hullpt + + +def convex_hull(points): + """Returns the points on the convex hull of points in CCW order.""" + + # Increasing guesses for the hull size. + sizes = [] + for t in range(len(points)): + sizes.append( 2 ** (2 ** t) ) + + for guess in sizes: + hulls = [] + for i in range(0, len(points), guess): + # Split the points into chunks of (roughly) the guess. + chunk = points[i:i + guess] + # Find the corresponding convex hull of these chunks. + hulls.append( graham_scan(chunk) ) + + # Find the extreme point and initialize the list of (hull,point) with it. + hullpt_pairs = [min_hull_pt_pair(hulls)] + + # Ensure we stop after no more than "guess" iterations. + for __ in range(guess): + pair = next_hull_pt_pair(hulls, hullpt_pairs[-1]) + if pair == hullpt_pairs[0]: + # Return the points in sequence + return [hulls[h][i] for h,i in hullpt_pairs] + hullpt_pairs.append(pair) + + +if __name__ == "__main__": + import sys + import random + import utils + import uberplot + import matplotlib.pyplot as plot + + if len(sys.argv) > 1: + scale = 100 + nb = int(sys.argv[1]) + points = [ (scale*random.random(),scale*random.random()) for i in range(nb)] + else: + points = [ + (0,40), + (100,60), + (40,0), + (50,100), + (90,10), + (50,50), + ] + + fig = plot.figure() + + hull = convex_hull( points ) + edges = list(utils.tour(hull)) + + ax = fig.add_subplot(111) + ax.set_aspect('equal') + ax.scatter( [i[0] for i in points],[i[1] for i in points], facecolor="red") + uberplot.plot_segments( ax, edges, edgecolor = "blue" ) + plot.show() +