Difference between revisions of "Graham Scan Implementation"

From HaskellWiki
Jump to: navigation, search
m (Testing: add graph of points)
Line 2: Line 2:
   
 
<haskell>
 
<haskell>
import Data.Ord (comparing)
 
 
import Data.List (sortBy)
 
import Data.List (sortBy)
--Graham Scan exercise
 
  +
import Data.Function (on)
   
--Direction type
 
  +
-- define data `Direction`
data Direction = LeftTurn
+
data Direction = GoLeft | GoRight | GoStraight
| RightTurn
+
deriving (Show, Eq)
| Straight
 
deriving (Show, Eq)
 
   
--Point type
 
  +
-- determine direction change via point a->b->c
data Point = Point (Double, Double)
 
  +
-- -- define 2D point
deriving (Show)
+
data Pt = Pt (Double, Double) deriving (Show, Eq, Ord)
  +
isTurned :: Pt -> Pt -> Pt -> Direction
  +
isTurned (Pt (ax, ay)) (Pt (bx, by)) (Pt (cx, cy)) = case sign of
  +
EQ -> GoStraight
  +
LT -> GoRight
  +
GT -> GoLeft
  +
where sign = compare ((bx - ax) * (cy - ay))
  +
((cx - ax) * (by - ay))
   
--some points
 
  +
-- implement Graham scan algorithm
p0 = Point (2.1,2.0)
 
p1 = Point (4.2,2.0)
 
p2 = Point (0.5,2.5)
 
p3 = Point (3.2,3.5)
 
p4 = Point (1.2,4.0)
 
p5 = Point (0.7,4.7)
 
p6 = Point (1.0,1.0)
 
p7 = Point (3.0,5.2)
 
p8 = Point (4.0,4.0)
 
p9 = Point (3.5,1.5)
 
pA = Point (0.5,1.0)
 
points = [p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,pA]
 
   
-- Actually, I'd leave it as EQ, GT, LT. Then, actually,
 
  +
-- -- Helper functions
-- if you wanted to sort points rotationally around a single point,
 
  +
-- -- Find the most button left point
-- sortBy (dir x) would actually work. --wasserman.louis@gmail.com
 
  +
buttonLeft :: [Pt] -> Pt
--Get direction of single set of line segments
 
  +
buttonLeft [] = Pt (1/0, 1/0)
dir :: Point -> Point -> Point -> Direction
 
  +
buttonLeft [pt] = pt
dir (Point (ax, ay)) (Point (bx, by)) (Point (cx, cy)) = case sign of
 
  +
buttonLeft (pt:pts) = minY pt (buttonLeft pts) where
EQ -> Straight
 
  +
minY (Pt (ax, ay)) (Pt (bx, by))
GT -> LeftTurn
 
  +
| ay > by = Pt (bx, by)
LT -> RightTurn
 
  +
| ay < by = Pt (ax, ay)
where sign = compare ((bx - ax) * (cy - ay)) ((by - ay) * (cx - ax))
 
  +
| ax < bx = Pt (ax, ay)
  +
| otherwise = Pt (bx, by)
   
--Get a list of Directions from a list of Points
 
  +
-- -- Main
dirlist :: [Point] -> [Direction]
+
convex :: [Pt] -> [Pt]
dirlist (x:y:z:xs) = dir x y z : dirlist (y:z:xs)
+
convex [] = []
dirlist _ = []
+
convex [pt] = [pt]
  +
convex [pt0, pt1] = [pt0, pt1]
  +
convex pts = scan [pt0] spts where
  +
-- Find the most buttonleft point pt0
  +
pt0 = buttonLeft pts
   
--Compare Y axes
 
  +
-- Sort other points `ptx` based on angle <pt0->ptx>
sortByY :: [Point] -> [Point]
 
  +
spts = tail (sortBy (compare `on` compkey pt0) pts) where
sortByY xs = sortBy lowestY xs
 
  +
compkey (Pt (ax, ay)) (Pt (bx, by)) = (atan2 (by - ay) (bx - ax),
where lowestY (Point(x1,y1)) (Point (x2,y2)) = compare (y1,x1) (y2,x2)
 
  +
{-the secondary key make sure collinear points in order-}
--get COT of line defined by two points and the x-axis
 
  +
abs (bx - ax))
pointAngle :: Point -> Point -> Double
 
pointAngle (Point (x1, y1)) (Point (x2, y2)) = (x2 - x1) / (y2 - y1)
 
   
--compare based on point angle
 
  +
-- Scan the points to find out convex
pointOrdering :: Point -> Point -> Ordering
 
  +
-- -- handle the case that all points are collinear
pointOrdering a b = compare (pointAngle a b) 0.0
 
  +
scan [p0] (p1:ps)
  +
| isTurned pz p0 p1 == GoStraight = [pz, p0]
  +
where pz = last ps
   
--Sort by angle
 
  +
scan (x:xs) (y:z:rsts) = case isTurned x y z of
sortByAngle :: [Point] -> [Point]
 
  +
GoRight -> scan xs (x:z:rsts)
sortByAngle ps = bottomLeft : sortBy (compareAngles bottomLeft) (tail (sortedPs))
 
  +
GoStraight -> scan (x:xs) (z:rsts) -- I choose to skip the collinear points
where sortedPs = sortByY ps
 
  +
GoLeft -> scan (y:x:xs) (z:rsts)
bottomLeft = head (sortedPs)
 
  +
scan xs [z] = z : xs
  +
  +
  +
-- Test data
  +
pts1 = [Pt (0,0), Pt (1,0), Pt (2,1), Pt (3,1), Pt (2.5,1), Pt (2.5,0), Pt (2,0)]
  +
-- (2,1)--(2.5,1)<>(3,1)
  +
-- -> |
  +
-- ---- v
  +
-- (0,0)------->(1,0)-- (2,0)<-(2.5,0)
  +
pts2 = [Pt (1,-2), Pt (-1, 2), Pt (-3, 6), Pt (-7, 3)]
  +
-- -- This case demonstrate a collinear points across the button-left point
  +
-- |-----(-3,6)
  +
-- |----- ^__
  +
-- <----- |---
  +
-- (-7,3) |
  +
-- |-- (-1,2)
  +
-- |---- <--
  +
-- |------------- |--|
  +
-- |------------>(1,-2)
  +
pts3 = [Pt (6,0), Pt (5.5,2), Pt (5,2), Pt (0,4), Pt (1,4), Pt (6,4), Pt (0,0)]
  +
-- -- This case demonstrates a consecutive inner points
  +
-- (0,4)<-(1,4)<--------------------(6,4)
  +
-- | |-----------^
  +
-- v (5,2)<(5.5,2)<-|
  +
-- (0,0)--------------------------->(6,0)
   
 
  +
pts4 = [Pt (0,0), Pt (0,4), Pt (1,4), Pt (2,3), Pt (6,5), Pt (6,0)]
  +
-- -- This case demonstrates a consecutive inner points
   
--Compare angles
 
  +
-- collinear points
compareAngles :: Point -> Point -> Point -> Ordering
 
  +
pts5 = [Pt (5,5), Pt (4,4), Pt (3,3), Pt (2,2), Pt (1,1)]
compareAngles = comparing . pointAngle
 
  +
pts6 = reverse pts5
 
--Graham Scan
 
gscan :: [Point] -> [Point]
 
gscan ps = scan (sortByAngle ps)
 
where scan (x:y:z:xs) = if dir x y z == LeftTurn
 
then scan (x:z:xs)
 
else x: scan (y:z:xs)
 
scan [x,y] = [x,y] -- there's no shame in a pattern match
 
-- of this type!
 
scan _ = []
 
 
</haskell>
 
</haskell>
 
== Testing ==
 
The above implementation has an error, which can be found using
 
 
<haskell>
 
import Test.QuickCheck
 
import Control.Monad
 
import Data.List
 
 
 
gscan' = map (\(Point x) -> x) . gscan . map Point
 
 
prop_gscan n = forAll (replicateM n arbitrary) (\xs -> sort (gscan' xs) == sort (gscan' (gscan' xs)))
 
 
</haskell>
 
 
This leads to a counterexample for lists of 5 elements
 
 
<haskell>
 
> quickCheck (prop_gscan 3)
 
+++ OK, passed 100 tests.
 
> quickCheck (prop_gscan 4)
 
+++ OK, passed 100 tests.
 
> quickCheck (prop_gscan 5)
 
*** Failed! Falsifiable (after 10 tests):
 
[(-2.4444173639942894,13.729627457461254),(65.72912810263666,5.955962930412828),(-34.288098030422404,80.6230134460068),(-6.446932942713564,-11.632835144720378),(2.861905401095031,1.1159493836896193)]
 
 
> let f = [(-2.4444173639942894,13.729627457461254),(65.72912810263666,5.955962930412828),(-34.288098030422404,80.6230134460068),(-6.446932942713564,-11.632835144720378),(2.861905401095031,1.1159493836896193)]
 
> gscan' f
 
[(-6.446932942713564,-11.632835144720378),(-34.288098030422404,80.6230134460068),(-2.4444173639942894,13.729627457461254),(65.72912810263666,5.955962930412828)]
 
> gscan' (gscan' f)
 
[(-6.446932942713564,-11.632835144720378),(-34.288098030422404,80.6230134460068),(65.72912810263666,5.955962930412828)]
 
> gscan' (gscan' (gscan' f))
 
[(-6.446932942713564,-11.632835144720378),(-34.288098030422404,80.6230134460068),(65.72912810263666,5.955962930412828)]
 
</haskell>
 
 
The versions of <tt>f</tt> with differing numbers of gscan applied are:
 
[[File:Graham_scan_implementation_failure.png]]
 

Revision as of 16:55, 11 April 2014

Descriptions of this problem can be found in Real World Haskell, Chapter 3

import Data.List (sortBy)
import Data.Function (on)

-- define data `Direction`
data Direction = GoLeft | GoRight | GoStraight
                 deriving (Show, Eq)

-- determine direction change via point a->b->c
-- -- define 2D point
data Pt = Pt (Double, Double) deriving (Show, Eq, Ord)
isTurned :: Pt -> Pt -> Pt -> Direction
isTurned (Pt (ax, ay)) (Pt (bx, by)) (Pt (cx, cy)) = case sign of
                    EQ -> GoStraight
                    LT -> GoRight
                    GT -> GoLeft
                    where sign = compare ((bx - ax) * (cy - ay)) 
                                         ((cx - ax) * (by - ay))

-- implement Graham scan algorithm

-- -- Helper functions
-- -- Find the most button left point
buttonLeft :: [Pt] -> Pt
buttonLeft [] = Pt (1/0, 1/0)
buttonLeft [pt] = pt
buttonLeft (pt:pts) = minY pt (buttonLeft pts) where
    minY (Pt (ax, ay)) (Pt (bx, by))
        | ay > by = Pt (bx, by)
        | ay < by = Pt (ax, ay)
        | ax < bx = Pt (ax, ay)
        | otherwise = Pt (bx, by)

-- -- Main
convex :: [Pt] -> [Pt]
convex []   = []
convex [pt] = [pt]
convex [pt0, pt1] = [pt0, pt1]
convex pts = scan [pt0] spts where
    -- Find the most buttonleft point pt0
    pt0 = buttonLeft pts

    -- Sort other points `ptx` based on angle <pt0->ptx>
    spts = tail (sortBy (compare `on` compkey pt0) pts) where
        compkey (Pt (ax, ay)) (Pt (bx, by)) = (atan2 (by - ay) (bx - ax),
                                               {-the secondary key make sure collinear points in order-}
                                               abs (bx - ax))

    -- Scan the points to find out convex
    -- -- handle the case that all points are collinear
    scan [p0] (p1:ps)
        | isTurned pz p0 p1 == GoStraight = [pz, p0]
        where pz = last ps

    scan (x:xs) (y:z:rsts) = case isTurned x y z of
                             GoRight    -> scan xs (x:z:rsts)
                             GoStraight -> scan (x:xs) (z:rsts) -- I choose to skip the collinear points
                             GoLeft     -> scan (y:x:xs) (z:rsts)
    scan xs [z] =  z : xs
   
                         
-- Test data
pts1 = [Pt (0,0), Pt (1,0), Pt (2,1), Pt (3,1), Pt (2.5,1), Pt (2.5,0), Pt (2,0)]
--                           (2,1)--(2.5,1)<>(3,1)
--                         ->          |
--                     ----            v
-- (0,0)------->(1,0)--      (2,0)<-(2.5,0)
pts2 = [Pt (1,-2), Pt (-1, 2), Pt (-3, 6), Pt (-7, 3)]
-- -- This case demonstrate a collinear points across the button-left point
--                |-----(-3,6)
--           |-----         ^__
--      <-----                |---
--  (-7,3)                       |
--     |--                      (-1,2)
--       |----                       <--
--           |-------------            |--|
--                        |------------>(1,-2)
pts3 = [Pt (6,0), Pt (5.5,2), Pt (5,2), Pt (0,4), Pt (1,4), Pt (6,4), Pt (0,0)]
-- -- This case demonstrates a consecutive inner points
-- (0,4)<-(1,4)<--------------------(6,4)
--   |                    |-----------^
--   v                  (5,2)<(5.5,2)<-|
-- (0,0)--------------------------->(6,0)

pts4 = [Pt (0,0), Pt (0,4), Pt (1,4), Pt (2,3), Pt (6,5), Pt (6,0)]
-- -- This case demonstrates a consecutive inner points

-- collinear points
pts5 = [Pt (5,5), Pt (4,4), Pt (3,3), Pt (2,2), Pt (1,1)]
pts6 = reverse pts5