Essays/Convex Hull
< Essays
Jump to navigation
Jump to search
Convex Hull
The convex hull of a set of 2D points is a subset of the points that define a convex polygon that contains all the points in the set.
The following verbs generate the convex hull.
The first verb is akin to a bubble sort, repeatedly processing edges and discarding ones found not to be on the hull.
The second verb is like a quicksort. It partitions the points repeatedly, discarding points as they are found not to be on the hull. This is the QuickHull algorithm.
For random data, the second verb is faster, but both versions are pretty fast, since the Eddy-Floyd step (removing points inside a quadrilateral) discards most of the points before the main processing starts.
Iterative version:
convexhull =: 3 : 0 pts =. y NB. Find points with max/in x and y; make a ccw quadrilateral out of them: 4 edges, closed ccwminmax =. (, {.) pts {~ 0 2 1 3 { , (i.!.0 <./ , >./)"1 |: pts NB. Calculate outside-the-edge masks for each point, and discard points inside all edges (Eddy-Floyd) NB. The determinant of (y,x) in ccw order in left-handed system is positive, so we take polygon-point to get - determinant for cw outmask =. ccwminmax 0&(>!.0) @: (2&((-/ . *)\)"2) @: (-"1"_ 1) pts assert. *./ 1 >: +/"1 outmask NB. each point can be outside only one edge NB. Each surviving point will be outside exactly one edge; associate the point with that edge outpts =. ((=/~ i. 4) -:"1/ outmask) <@# pts NB. Find the convex hull corresponding to each edge. The two endpoints PQ of the edge are known to NB. be on the hull. Sort the points outside the edge into the order they will be encountered in a CCW sweep from P. NB. For each edge AB, we calculate the winding of the triangle ABC (C is the next point). If ABC is cw, delete point B. NB. Repeat the procedure until no points are deleted. NB. Taking, say, the first edge, which goes from ymin to xmin, the points must have a smaller x than the edge point but they could have equal y. NB. So calculate the slope as dy/dx, sort decending. For the second edge, the points must have larger y, so NB. calculate -dx/dy, sort decending. For third, dy/dx descending; for fourth, -dx/dy descending sortpts =. outpts \:&.> 1 _1 1 _1 *&.> %/"1&.> 0 _1 0 _1 |."1&.> outpts -"1&.> <"1 }: ccwminmax hullpts =. (<"1 }. ccwminmax) ({.@] , }.@] (#~ 0&(<!.0)) 3&((-/ . *)@:(}. -"1 {.)\)@:,~)^:(1<#@])^:_&.> sortpts NB. Append each quadrilateral point with the following hull points and run together to form the result NB. If a point is two vertices of the Floyd rectangle it might appear twice, so delete duplicates ~.!.0 ; (<"1 }: ccwminmax) ,&.> hullpts )
Partitioning version:
NB. convex hull, recursive and faster for random points NB. y is a table of y,x value (or x, y if you want to think of it that way) NB. result is y,x points of the convex hull, in CCW order in left-handed coordinates convexhull =: 3 : 0 NB. Find points with min/max y. This is the starting edge with endpoints on the hull NB. minmax =. (({~ (i.!.0 <./ , >./)) {."1) y NB. For a better partition, this is to pick the min/max distance from a point not inside the hull. NB. More work, so it's about the same speed. minmax =. (({~ (i.!.0 <./ , >./)) +/"1@:*:) (-"1 <./) y NB. Classify each point on the correct side of the edge or its reverse NB. The determinant of (y,x) in ccw order in left-handed system is positive NB. So, - determinant for ABy means point goes with AB, + with BA outmask =. y 0&(>!.0) @: (-/ . *) @: (-"1/) minmax NB. Create AB,&<BA, and associate the points with each NB. Find the convex hull corresponding to each edge, and run together to produce the hull NB. The processing might produce sequences of collinear points, so remove them (#~ 0 <!.0 (3) ((-/ . *) @: (}. -"1 {.))\ {:,],{.) ; ((,&< |.) minmax) hullsection&.> outmask (# ; (#~ -.)~) y ) NB. x is a table of two points on the hull AB, in ccw order NB. y is a table, a set of points that are all known to be outside the chord described by x NB. result is 0{x, followed by any other points on the hull, in ccw order. NB. 1{x is not part of the result hullsection =: 4 : 0 NB. If y has 0 or 1 point, we're done if. 2 > #y do. ({.x) , y else. NB. Find the point of y P with the largest (negative) perpendicular distance from the chord AB. Create AP,:PB appb =. x ([ (0 1 |."0 2 ,:"1) ] {~ [: (i.!.0 <./) [: -/ . * -"1"2 1) y NB. For each point of y Q, look at APQ and PBQ to produce the mask of which edges the point is outside outmask =. 0 >!.0 appb ([: -/ . * -"1"3 1) y NB. Partition the points according to which edge they are outside. Points not outside an edge are discarded NB. Recur to produce the hull points for each partition NB. Result is the hull points for each partition ; (<"2 appb) hullsection&.> ((=/~ i. 2) -:"1/ outmask) <@# y end. )
Contributed by Henry Rich