/** 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.

    These points can be obtained or passed to Frink's graphics polygon class:
    https://frinklang.org/#PolygonMethods

    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.  It returns the result as an array of [x,y] points.
   // These can be obtained or passed to Frink's graphics polygon class:
   // https://frinklang.org/#PolygonMethods
   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[rest[points], 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.
      // TODO:  Should first[array] return an array?  It currently returns an
      // EnumeratingExpression which led to confusion
      s = toArray[first[points, 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
         // TODO:  Factor out multiple calls to length[s]
         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]]