Essays/Convex Hull

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