-
Notifications
You must be signed in to change notification settings - Fork 3
/
hull.py
238 lines (187 loc) · 8.62 KB
/
hull.py
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
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
import operator
from utils import LOG,LOGN
from geometry import x,y,euclidian_distance
# Based on the excellent article by Tom Switzer <[email protected]>
# 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 x<y, 0 if x==y, +1 if x>y
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,
# but this is not the last possible tangent.
elif icenter+1 < iright:
# 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) ] )
# There is no more possible tangent, the hull contains a straight segment.
else:
return ileft
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.
for guess in ( 2**(2**t) for t in range(len(points)) ):
LOG( "Guess",guess)
hulls = []
for i in range(0, len(points), guess):
# LOG(".")
# 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):
LOG("*")
pair = next_hull_pt_pair(hulls, hullpt_pairs[-1])
if pair == hullpt_pairs[0]:
# Return the points in sequence
LOGN("o")
return [hulls[h][i] for h,i in hullpt_pairs]
hullpt_pairs.append(pair)
LOGN("x")
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()