/** This program calculates the convex hull of a list of points. It uses Graham's scan algorithm which is approximately O(n log n) with respect to the number of points. You wil generally call this with ConvexHull.hull[points] where points is an array of [x,y] points. To graph your solution, call: ConvexHull.plotHull[points].show[] where points is an array of [x,y] points. This somewhat follows the discussion at: http://www.geeksforgeeks.org/convex-hull-set-2-graham-scan/ It should be noted that the sample implementation at that site is broken in several ways. */ class ConvexHull { // This functions calculates the convex hull of an array of points specified // as [x,y] coordinates. class hull[points] := { points = toArray[points] // Points can come in as a set or enumeration size = length[points] // Find the smallest y value, smallest right if there's a tie ymin = points@0@1 minIndex = 0 for i = 1 to size-1 { [x, y] = points@i if ((y < ymin) or (ymin == y and x < points@minIndex@0)) { ymin = y minIndex = i } } // Place the smallest y value at index 0 by swapping temp = points@0 points@0 = points@minIndex points@minIndex = temp // Sort items 1 through size-1 with respect to point 0. // A point p1 comes betwen p2 in sorted output if p2 has a larger polar // angle (in counterclockwise direction) than p1. // TODO: Implement sort with sub-region begin and end specified result = sort[slice[points,1,undef], getFunction["compare", 3], points@0] result.pushFirst[points@0] points = result // If two or more points make the same angle with p0, remove all but // the farthest. The sort above made the farthest last in the list. m = 1 for i=1 to size-1 { // Keep removing element i while i and i+1 are colinear with p0 while i < size-1 and orientation[points@0, points@i, points@(i+1)] == 0 i = i + 1 points@m = points@i m = m + 1 // Update size of modified array } // If modified array has fewer than 3 points, no convex hull is possible if m < 3 return undef // Create an empty stack and push first three points onto it. s = sliceLength[points, 0, 3] // Process remaining n-3 points for i = 2 to m-1 { // Keep removing top while the angle formed by points next-to-top, // top, and points@i makes a non-left turn. The first clause is // apparently necessary to fix a bug in the original while length[s] >= 2 and orientation[s@(length[s]-2), s@(length[s]-1), points@i] != 2 s.pop[] s.push[points@i] } return s } // This is a comparison function that compares points with respect to the // first point p0. (Which is passed as an extra "data" argument to the sort // function) class compare[p1, p2, p0] := { o = orientation[p0, p1, p2] // If points are colinear, put nearest point first if o == 0 return -(distanceSquared[p0, p2] <=> distanceSquared[p0, p1]) return o==2 ? -1 : 1 } /** Returns the squared distance between two points. */ class distanceSquared[p1, p2] := { [x1,y1] = p1 [x2,y2] = p2 // TODO: This doesn't make sense with disparate units return (x1-x2)^2 + (y1-y2)^2 } /** Finds the orientation of the ordered triplet of points [p,q,r]. Returns: 0 if p, q, and r are colinear 1 if clockwise 2 if counterclockwise */ class orientation[p, q, r] := { [px, py] = p [qx, qy] = q [rx, ry] = r v = (qy - py) * (rx - qx) - (qx - px) * (ry - qy) units = px py if v == 0 units // This makes units come out right return 0 // Colinear else return v > 0 units ? 1 : 2 // Clockwise or counterclockwise } /** This is a function to calculate and plot a hull of points. It returns a graphics object which you can display or print with its show[] or print[] methods. */ class plotHull[points] := { g = new graphics for [x,y] = points g.fillEllipseCenter[x,y,.1, .1] hull = hull[points] ph = new polygon[hull] g.color[0,0,1] g.add[ph] g.color[1,0,0] // Draw the actually used hull points in red for [x,y] = hull g.fillEllipseCenter[x,y,.2, .2] return g } } /* // Testing stuff. p = [ [0,3], [1,1], [2,2], [4,4], [0,0], [1,2], [3,1], [3,3] ] ConvexHull.plotHull[p].show[] println[ConvexHull.hull[p]] // Sample from website above (in reverse order) p = [ [0,3], [1,1], [2,2], [4,4], [0,0], [1,2], [3,1], [3,3], [2,4] ] ConvexHull.plotHull[p].show[] println[ConvexHull.hull[p]] println["Colinear3"] p = [ [0,0], [1,1], [2,2]] ConvexHull.plotHull[p].show[] println[ConvexHull.hull[p]] println["Colinear4"] p = [ [0,0], [1,1], [2,2], [3,3]] ConvexHull.plotHull[p].show[] println[ConvexHull.hull[p]] mean = 5 sd = 1 d = new array for p = 1 to 2000 d.push[ [randomGaussian[mean, sd], randomGaussian[mean, sd]] ] g = ConvexHull.plotHull[d] g.show[] */ // Testing with units of measure. Not sure if this can be made to make sense //p2 = [ [0,3s], [1,1s], [2,2s], [4,4s], [0,0s], [1,2s], [3,1s], [3,3s], [2,4s] ] //println[ConvexHull.hull[p2]]