Graham’s Scan; from imperative to functional

This story actually happened in a few hours, but the ride was the must fun I have had in the last few months(professionally), so I thought in sharing it.

The task

I had to implement the Graham’s Scan in Scala. Given that the algorithm has been around since the 70’s, you’d think it would be a walk in the park. It was more like a game of rapid chess.

GRAHAM_SCAN

1.         Find p0 in Q with minimum y-coordinate (and minimum x-coordinate if there are ties).
2.         Sorted the remaining points of Q (that is, Q − {p0}) by polar angle in counterclockwise order with respect to p0.
3.         TOP [S] = 0                ▷ Lines 3-6 initialize the stack to contain, from bottom to top, first three points.
4.         PUSH (p0, S)
5.         PUSH (p1, S)
6.         PUSH (p2, S)
7.         for i = 3 to m             ▷ Perform test for each point p3, …, pm.
8.             do while the angle between NEXT_TO_TOP[S], TOP[S], and pi makes a nonleft turn  ▷ remove if not a vertex
9.                         do POP(S)
10.                 PUSH (S, pi)
11.       return S

As you can see, the “good stuff” happens in a stack.  All the versions describing the algorithm that I could found, were written this way. So instead of starting from the abstraction of the algorithm, I had to abstract an opinionated interpretation, to come up with a functional version.

First things first, I created a data set of points for plotting in Octave. Octave already has a convex hull function, so I could see the solution before hand.


X = [0.410036859415169 0.16990953347486892 0.7438936112344358 0.9879763068249436 0.09850362854534644 0.27745454018107196 0.2463477117635604 0.0369350463381225 0.7356804213480419 0.979020964716884 0.4104133424281955 0.9566561367393451 0.36405116056182674 0.6743515475246237 0.2014896354108069 0.5937635359142525 0.9990876238009762 0.4753053860797861 0.15127539785912514 0.3109492347362207 0.3710641858259751 0.4981026670484303 0.472172885898987 0.9804204746053471 0.20208956184140103 0.7386385686849415 0.3182368631449125 0.4529941453549946 0.3035123394853474 0.3785220082038615 0.6028255171683854 0.4754941518809386 0.02885621173133124 0.619495881160238 0.03310900644639769 0.2791084654454382 0.8674411544744653 0.393894922057599 0.271592427317065 0.31136184917695364];
Y = [0.0964872038572564 0.8409000674315462 0.5513623119954284 0.9591944778976642 0.3152016460947692 0.16362690913082834 0.31570455065654346 0.5538361485253371 0.7670638114116081 0.015668444677785942 0.3467907029082512 0.07399794456597841 0.12236562106445759 0.9981954727371612 0.9308695268083278 0.8560770018856956 0.6579161973082933 0.5959520081450651 0.13661844095029363 0.6514314696732396 0.19298645312277207 0.00994355691826554 0.7010697212377054 0.1634450128412981 0.3793185478977382 0.868604478585853 0.4590548622219143 0.18377356484707785 0.2866066023269005 0.8213505729786262 0.38257319264118383 0.4566803134679752 0.4748548161642837 0.26842322605193014 0.5464910173382869 0.6503654422352256 0.9950193150372821 0.730931140118583 0.5869777402729242 0.06419667711655552];

k = convhull(X,Y);

plot (X(k), Y(k), "r-", X, Y, "b+");

 

octave

There, a reliable sample from octave

Starting

First, I found a cool way of sorting the points by polar angle:

val coordList = origin :: coords.filterNot(_ == origin).
        sortBy(point => atan2(point._2 - origin._2, point._1 - origin._1))

Just to be sure, I plotted the sorted points without any further processing:

incomplete

Beautiful polar angles!

In another article, I found this amazing formula for deciding if a point makes a left turn in relation with a line created by the other two:

// isLeft(): tests if a point is Left|On|Right of an infinite line.
//    Input:  three points P0, P1, and P2
//    Return: > 0 for P2 left of the line through P0 and P1
//            = 0 for P2 on the line
//            < 0 for P2 right of the line
//    See: Algorithm 1 on Area of Triangles
inline float
isLeft( Point P0, Point P1, Point P2 )
{
    return (P1.x - P0.x)*(P2.y - P0.y) - (P2.x - P0.x)*(P1.y - P0.y);
}

So, P2 is the new point, P1 the last point added to the hull, and P0 the second to last.

Misinterpretation.

In my first version I was doing the following:

  1. Add the first three points to the hull.
  2. After every new point, check if it makes a left turn in relation with the two previous
  3. If it does, add it to the hull.
  4. If it doesn’t,  remove the last point from the solution and add the new one to the hull.

Sounds reasonable right?

I wrote the code and plotted the result in Octave.

 

foreigner

Oh noes! it doesn’t work!

Yes! No! It found all of the points of the hull, but it kept also the wrong ones. Is hard when a perfectly reasonable code doesn’t do you think it should do. Anyway, I knew that the algorithm was fine, I knew that the points were sorted in a proper way. So it had to be my interpretation of the algorithm. After a careful  reading of the pseudo code above, I realized that before adding a new point, I needed to scan the already found vertices in relation with the new point, and removing  P1 every time a non-left angle were found.

foldRight

The coolest moment for me was realizing that instead of the stack, I could just use foldRight. That gave away the solution.

case class Point(x: Double, y: Double) {
    //improves readability
    def goesLeft(p0: Point, p1: Point): Boolean =
      (p1.x - p0.x) * (this.y - p0.y) - (this.x - p0.x) * (p1.y - p0.y) > 0
  }

  def scan(coords: List[Point]): List[Point] = {
    //instead of using PUSH and POP, foldRight
    def addToHull(hull: List[Point], new_point: Point): List[Point] =
      //the new point is always added but filtering all the non-left turn vertices from the current hull
      new_point :: hull.foldRight(List.empty[Point]) {
      case (p1, currentHull@(p0 :: _)) =>
        if (new_point.goesLeft(p0, p1))
          p1 :: currentHull
        else
          currentHull
      case (p, currentHull) => p :: currentHull
    }
    val min = coords.minBy(_.y)
    //adding the min to close the polygon on octave
    min :: coords.sortBy(point => atan2(point.y - min.y, point.x - min.x))
      .foldLeft(List.empty[Point])(addToHull)
  }

good

It worked!

I created a project in GitHub, that also has a java 8 solution, not very functional.  What the test suites actually do is to generate random samples and a Octave file with the solution for plotting.

Advertisements

One thought on “Graham’s Scan; from imperative to functional

  1. Hey there, have you seen Lean Body Maximizer? (Go Google it) You will learn about the crimes we commit against ourselves. With “Lean Body Maremizix”, you will discover how to lose fat fast.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s