ConvexHull.frink

View or download ConvexHull.frink in plain text format


/** 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] :=
   {
      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]]


View or download ConvexHull.frink in plain text format


This is a program written in the programming language Frink.
For more information, view the Frink Documentation or see More Sample Frink Programs.

Alan Eliasen was born 17983 days, 22 hours, 30 minutes ago.