Graham’s Scan in Scala

Sometimes my job throws an interesting problem my way. ¬†This week I was presented with a very odd geometry problem ūüôā

I needed to generate KML files from geographic data and one of my requirements was to represent certain geographic areas as polygons, the vertices of which would be supplied (along with the rest of the data) by another OSGi service running elsewhere. ¬†Seems fairly straightforward — In KML, the vertices of a polygon are usually specified as follows:

<Placemark>
  <name>LinearRing.kml</name>
  <Polygon>
    <outerBoundaryIs>
      <LinearRing>
        <coordinates>
          -122.365662,37.826988,0
          -122.365202,37.826302,0
          -122.364581,37.82655,0
          -122.365038,37.827237,0
          -122.365662,37.826988,0
        </coordinates>
      </LinearRing>
    </outerBoundaryIs>
  </Polygon>
</Placemark>

The coordinates tag requires at least four longitude/latitude/altitude triples to be specified, with the last coordinate being the same as the first. ¬† Here is where the problem comes in — The order in which these coordinates are specified matters (they must be specified in counter-clockwise order). ¬†To mix up the order of the coordinates would have unpredictable results (e.g. crazy geometric shapes) when the data is later displayed via Google Earth (or some other application that supports KML files). ¬†However, the area vertices are indeed fed to my KML generator in no particular order (and the services providing the data cannot be changed to guarantee a particular ordering).

So… how do I put the points in order? ¬†“Surely this is a solved problem.” ¬†I thought, turning to the all-knowing internet. ¬†A bit of searching turned up an algorithm called Graham’s Scan. ¬†Basically, this algorithm takes a bag of random coordinates and generates a convex hull with vertices defined in counter-clockwise order (Note: This may not be suitable if you’re trying to faithfully recreate complex geometries, fortunately I’m mostly concerned with rectangular areas). ¬†Roughly, the algorithm works as follows:

  1. Find the coordinate with the lowest y-value.
  2. Sort the remaining points by the polar angle between the line defined by the current point and the point found in step 1, and the x-axis.
  3. Go through the list, evaluating 3 points at a time (you’re concerned with the angle of the turn made by the two resulting line segments). ¬†Eliminate any non counter-clockwise turns (i.e., non-convex corners).

I found several example implementations for this algorithm in various languages: C++, Java, etc. ¬†Since I’m coding this up in Scala, I wasn’t too happy with any of those, and I couldn’t find an example in Scala to rip off draw inspiration from. ¬†However, I did manage to find an implementation in Haskell¬†which I used as a rough guide. ¬†Anyway, here’s my attempt at Graham’s Scan in Scala:

type Coordinate = (Double, Double) //Longitude, Latitude

  protected def processArea(coords: List[Coordinate]): List[Coordinate] = {
      //returns > 0 if points form a counter clockwise turn, 
      // < 0 if clockwise, and 0 if collinear
      def ccw(p1: Coordinate, p2: Coordinate, p3: Coordinate) =
        (p2._1 - p1._1)*(p3._2 - p1._2) - (p2._2 - p1._2)*(p3._1 - p1._1)

      //Scan the List of coordinates and find our vertices
      def scan(theCoords: List[Coordinate]): List[Coordinate] = theCoords match {
        case xs if xs.isEmpty => List()
        case xs if xs.size == 2 => theCoords
        case x::y::z::xs if ccw(x,y,z) > 0 => x::scan(y::z::xs)
        case x::y::z::xs => scan(x::z::xs)
      }
      //find the coordinate with the lowest latitude
      val origin = coords.minBy(_._2)
      //sort the rest of the points according to their polar angle (the angle between the line
      //defined by the origin and the current point, and the x-axis)   
      val coordList = origin :: coords.filterNot(_ == origin).
        sortBy(point => atan2(point._2 - origin._2, point._1 - origin._1))
      //do the graham scan
      scan(coordList)
    }

I think you’ll find that’s a bit shorter than some of the imperative language implementations out there ūüôā