Tuesday 8 June 2010

Graham Scan in Haskell

The last exercise in chapter 3 of Real World Haskell is quite more difficult than the previous ones, so it deserves a post with code. I have not used any function not yet seen at that point in the book, so the code is not very pretty, and I have followed quite closely the algorithm for Graham Scan as explained in Introduction to Algorithms, which can also be seen here.

When I learn a bit more about Haskell syntax I will try to pretify it, but for the moment here it goes in all its uglyness:


data Direction = LeftTurn
    | RightTurn
    | Straight
      deriving (Eq, Show)

turn (x1,y1) (x2,y2) (x3,y3)
     | cross_product == 0 = Straight
     | cross_product > 0 = LeftTurn
     | cross_product < 0 = RightTurn
     where
       cross_product = (x2 - x1)*(y3 - y1) - (y2 - y1)*(x3 - x1)

miny [x] = x
miny ((x1,y1):(x2,y2):xs)
     | y1 < y2 = miny ((x1,y1):xs)
     | y2 < y1 = miny ((x2,y2):xs)
     | x1 < x2 = miny ((x1,y1):xs)
     | x2 < x1 = miny ((x2,y2):xs)
     | otherwise = error "repeated points"

removepoint _ [] = []
removepoint (x,y) ((x1,y1):xs)
    | x == x1 && y == y1 = xs
    | otherwise = (x1,y1):removepoint (x,y) xs

sort_merge_polar origin [] = []
sort_merge_polar origin [x] = [x]
sort_merge_polar (x1,y1) lista =
    merge (sort_merge_polar (x1,y1) (take midpoint lista)) (sort_merge_polar (x1,y1) (drop midpoint lista))
        where
          midpoint = div (length lista) 2
          merge [] [] = []
          merge [] lisb = lisb
          merge lisa [] = lisa
          merge ((x2,y2):xs) ((x3,y3):ys)
              | turn (x1,y1) (x2,y2) (x3,y3) == LeftTurn = (x2,y2) : merge xs ((x3,y3):ys)
              | turn (x1,y1) (x2,y2) (x3,y3) == RightTurn = (x3,y3) : merge ((x2,y2):xs) ys
              | turn (x1,y1) (x2,y2) (x3,y3) == Straight = further_point : merge xs ys
              where further_point
                        | abs(x2 - x1) > abs(x3 - x1) = (x2,y2)
                        | abs(x2 - x1) < abs(x3 - x1) = (x3,y3)
                        | abs(y2 - y1) > abs(y3 - y1) = (x2,y2)
                        | abs(y2 - y1) < abs(y3 - y1) = (x3,y3)
                        | otherwise = error "repeated points.."

graham_scan (pa:pb:pc:ps) = result
    where
      p0 = miny (pa:pb:pc:ps)
      sorted_list = sort_merge_polar p0 (removepoint p0 (pa:pb:pc:ps))
      p1 = head sorted_list
      p2 = head (tail sorted_list)
      result = iterate_graham [p2,p1,p0] (tail (tail sorted_list))
      iterate_graham convex_hull [] = reverse(convex_hull)
      iterate_graham (top:ntop:xs) (pi:ps)
          | turn ntop top pi == RightTurn = iterate_graham (ntop:xs) (pi:ps)
          | turn ntop top pi == LeftTurn = iterate_graham (pi:top:ntop:xs) ps
          | turn ntop top pi == Straight = error "didn't I get rid of straight points relative to p0??"

graham_scan _ = error "we need at least three points"


As a basic test, if we compute the Convex Hull of the following points, we get the correct result. Other online solutions seem to choke on inputs like this, probably due to an incorrect implementation of the Graham Scan algorithm.

*Main> graham_scan [(10,0), (10,1),(-10,1),(-10,0),(-7,0),(-10,2),(-10,3),(-4,1),(-2,2),(-12,1)]
[(-10,0),(10,0),(10,1),(-10,3),(-12,1)]
*Main>

2 comments:

Jeff said...

Your code definitely looks like it works a whole lot better than mine!

As I was reading through the code, I had a few quick comments.

Could you replace "removePoint" with Data.List.deleteBy?

In "sort_merge_polar" you've got take and drop being used. I think you could use splitAt to do that in one go.

Is it possible that "minY" be expressed as a fold?

angelv said...

Hi Jeff,

I was not aware of deleteBy (I'm just reading the book sequentially), but I already came accross splitAt and folds, so when it is time to revisit my code I will use them.

Thanks