Convex Hull

05_source
20161206
 ```1 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 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 ``` ``````#!/usr/bin/python import math import matplotlib.pyplot as plt import copy import random # only needed for random points import numpy as np # only needed for random points from normal distrib # convenience functions to get the x,y value from a tuple, # and avoid using the wrong tuple index # (sacrifice a bit of performance for legibility) def x(tp): return tp def y(tp): return tp # calculate the angle for a line between two points and the horizontal def calc_angle(pt1, pt2): dy=float( y(pt2) - y(pt1) ) dx=float( x(pt2) - x(pt1) ) angle=0.0 if dx!=0: angle=math.atan(dy/dx) else: if dy==0: angle=0.0 elif dy<0: angle=math.atan(float('-inf')) else: angle=math.atan(float('inf')) # if the calucated angle is negative, make it positive while (angle<0): angle=angle+2*math.pi; return angle # calculate the distance between two points def dist(pt1,pt2): dx=float( x(pt2) - x(pt1) ) dy=float( y(pt2) - y(pt1) ) return math.sqrt( dx*dx + dy*dy) # does moving from point pts0 -> pt1 -> pt2 describe a counterclockwise motion ? def clockwise(pt1, pt2, pt3): area = (x(pt2)-x(pt1))*(y(pt3)-y(pt1))-(y(pt2)-y(pt1))*(x(pt3)-x(pt1)) rv=0 # collinear if area<0: rv=-1 # clockwise elif area>0: rv=+1 # counter clockwise return rv ## MAIN ---------------------------------------------------------------------- # --- DATA ------------ points= [ (9,2), (5,6), (3,10), (11,13), (7,6), (15,18), (19,16), (13,4), (23,8), (17,12) ] # variation: random points #points=[] #for i in range(1,1000): points.append( (random.randint(1,100), random.randint(1,100)) ) # variation: sedgewick's points (approx) #points=[] # points= [(373, 234), (312, 235), (313, 128), (165, 164), (185, 91), (190, 34), # (452, 249), (404, 109), (237, 297), (262, 130), (124, 127), (334, 29), (459, 2)] # variation : random points, but normally distributed #points=list(zip( # [ int(r) for r in np.random.normal(100, 20, 1000)], # [ int(r) for r in np.random.normal(100, 20, 1000)] # ) ) # --------------------- # identify the point with the lowest y coordinate and put it first in the list of points # special case: in case of tie, decide by lowest x coordinate points=sorted( points, key=lambda a: a) # secondary sort on x points=sorted( points, key=lambda a: a) # primary sort on y (python sort is stable) # enrich the tuple list with the angle between the line through (point, first point) # and the horizontal points = [ (x(pt),y(pt),calc_angle(points,pt)) for pt in points ] # order the points by angle points=sorted(points, key=lambda a: a) # after sorting we don't need the angles anymore, let's drop them points = [ (x(pt),y(pt)) for pt in points ] # the core of the graham scan: # iterate over the points in ascending angle order # check if pt1-> pt2 -> pt3 makes a counterclockwise turn # if NOT, then discard pt2 hull=copy.deepcopy(points) # deep copy, because we are going to delete points i=0 while True: if clockwise(hull[i],hull[i+1],hull[i+2])<=0: # clockwise or collinear del hull[i+1] i-=1 # take a step back else: i+=1 # step forward if (i+2)>=len(hull): # are we done? break print "Points on the hull" ,hull # to close the polygon, (for drawing, and calculating distance) append the first point hull.append(hull) print "Perimeter distance: ", sum( [ dist( hull[i],hull[i+1] ) for i in range(len(hull)-1) ] ) # plotting the points and the hull plt.scatter(*zip(*points),s=20) # plot all the points plt.plot(*zip(*hull),linestyle='-') # plot the hull plt.savefig("convex_hull.png")``````

Notes by Data Munging Ninja. Generated on nini:sync/20151223_datamungingninja/convexhull at 2016-12-26 08:22