Graham Scan Algorithm

graham_scan.hs

-- Read: http://en.wikipedia.org/wiki/Graham_scan
-- ================================================ --
-- ===========  Definitions/Imports  ============== --
-- ================================================ --

import Data.List

-- ch03/ex09.hs
data Direction = LeftTurn | RightTurn | Straight deriving (Show,Eq)
data Point2D = Point2D {x :: Double, y :: Double} deriving (Show,Ord,Eq)

distance :: Point2D -> Point2D -> Double
distance (Point2D x1 y1) (Point2D x2 y2) = sqrt((y2-y1)^2 + (x2-x1)^2)

dot :: Point2D -> Point2D -> Double
dot (Point2D x1 y1) (Point2D x2 y2) = x1*x2 + y1*y2

norm :: Point2D -> Double
norm (Point2D x y) = sqrt(x^2 + y^2)

minus :: Point2D -> Point2D -> Point2D
minus (Point2D x1 y1) (Point2D x2 y2) = Point2D (x1-x2) (y1-y2)

-- ch03/ex10.hs
direction :: Point2D -> Point2D -> Point2D -> Direction
direction (Point2D x1 y1) (Point2D x2 y2) (Point2D x3 y3)
              | cross_product_direction > 0 = LeftTurn
              | cross_product_direction < 0 = RightTurn
              | otherwise                   = Straight                    
              where cross_product_direction = (x2 - x1)*(y3 - y1) - (y2 - y1)*(x3 - x1)

-- ch03/ex11.hs
turns :: [Point2D] -> [Direction]
turns points 
      | length(points) >= 3 = turn:(turns (tail points)) 
      | otherwise           = []
      where [x,y,z] = take 3 points
            turn    = direction x y z

-- ch03/ex12.hs
compareHeight :: Point2D -> Point2D -> Ordering
compareHeight (Point2D x1 y1) (Point2D x2 y2)
              | y1 < y2   = LT
              | y1 > y2   = GT
              | otherwise = compare x1 x2

lowest :: [Point2D] -> Point2D
lowest points = head (sortBy compareHeight points)

projection :: Point2D -> Point2D -> Double
projection p1 p2 = let adjacent   = dot p1 p2
                       hypotenuse = (norm p1) * (norm p2)
                   in  adjacent / hypotenuse

cosine :: Point2D -> Double
cosine point = (x point) / (norm point)

angle :: Point2D -> Double
angle point = acos (cosine point)

--
-- sort the points by height (y coordinate)
-- choose the first sorted point as the origin point "P"
-- sort the current points by using the angle between the line from "P" to the given point and the x-axis
--
sortByAngle :: [Point2D] -> [Point2D]
sortByAngle points = let origin:rest        = sortBy compareHeight points
                         compareAngle p1 p2 = let origin_p1 = minus p1 origin
                                                  origin_p2 = minus p2 origin                         
                                              in  compare (angle origin_p1) (angle origin_p2)
                     in  origin:(sortBy compareAngle rest)

--
-- find the first right turn
-- remove the middle point of the 3 points which caused the right turn
-- return the list of points with the next right turn removed
--
removeNextRightTurn :: [Point2D] -> [Point2D]
removeNextRightTurn points 
      | length(points) >= 3 = if turn == LeftTurn 
                              then x:(removeNextRightTurn (tail points)) 
                              else [x,z]++(drop 3 points) 
      | otherwise           = points
      where [x,y,z] = take 3 points
            turn    = direction x y z

--
-- the convex hull will be the remaining sequence after all right turns are removed
--
grahamScan :: [Point2D] -> [Point2D]
grahamScan points = if elem RightTurn (turns points) 
                    then grahamScan $ removeNextRightTurn points 
                    else points

-- ====================================== --
-- ===========  Test Data  ============== --
-- ====================================== --

a = Point2D 9 2
b = Point2D 7 4
c = Point2D 8 7
d = Point2D 5 9
e = Point2D 3 8
f = Point2D 5 6
g = Point2D 3 4
h = Point2D 5 3
i = Point2D 4 2

points = [a,b,c,d,e,f,g,h,i]

-- sorted_points = [i,a,b,h,c,f,d,e,g]
sorted_points = sortByAngle points

-- convex_hull = [i,a,c,d,e,g]
convex_hull = grahamScan sorted_points