Difference between revisions of "Graham Scan Implementation"
Jump to navigation
Jump to 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) |
||
+ | import Data.Function (on) |
||
− | --Graham Scan exercise |
||
+ | -- define data `Direction` |
||
− | --Direction type |
||
− | data Direction = |
+ | data Direction = GoLeft | GoRight | GoStraight |
− | + | deriving (Show, Eq) |
|
− | | Straight |
||
− | deriving (Show, Eq) |
||
+ | -- determine direction change via point a->b->c |
||
− | --Point type |
||
+ | -- -- define 2D point |
||
− | data Point = Point (Double, Double) |
||
− | + | 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 |
||
− | --some points |
||
− | 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] |
||
+ | -- -- Helper functions |
||
− | -- Actually, I'd leave it as EQ, GT, LT. Then, actually, |
||
+ | -- -- Find the most button left point |
||
− | -- if you wanted to sort points rotationally around a single point, |
||
+ | buttonLeft :: [Pt] -> Pt |
||
− | -- sortBy (dir x) would actually work. --wasserman.louis@gmail.com |
||
+ | buttonLeft [] = Pt (1/0, 1/0) |
||
− | --Get direction of single set of line segments |
||
+ | buttonLeft [pt] = pt |
||
− | dir :: Point -> Point -> Point -> Direction |
||
+ | buttonLeft (pt:pts) = minY pt (buttonLeft pts) where |
||
− | dir (Point (ax, ay)) (Point (bx, by)) (Point (cx, cy)) = case sign of |
||
+ | minY (Pt (ax, ay)) (Pt (bx, by)) |
||
− | EQ -> Straight |
||
− | + | | ay > by = Pt (bx, by) |
|
− | + | | ay < by = Pt (ax, ay) |
|
− | + | | ax < bx = Pt (ax, ay) |
|
+ | | otherwise = Pt (bx, by) |
||
+ | -- -- Main |
||
− | --Get a list of Directions from a list of Points |
||
− | + | convex :: [Pt] -> [Pt] |
|
+ | convex [] = [] |
||
− | dirlist (x:y:z:xs) = dir x y z : dirlist (y:z:xs) |
||
− | + | 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> |
||
− | --Compare Y axes |
||
+ | spts = tail (sortBy (compare `on` compkey pt0) pts) where |
||
− | sortByY :: [Point] -> [Point] |
||
+ | compkey (Pt (ax, ay)) (Pt (bx, by)) = (atan2 (by - ay) (bx - ax), |
||
− | sortByY xs = sortBy lowestY xs |
||
+ | {-the secondary key make sure collinear points in order-} |
||
− | where lowestY (Point(x1,y1)) (Point (x2,y2)) = compare (y1,x1) (y2,x2) |
||
+ | abs (bx - ax)) |
||
− | --get COT of line defined by two points and the x-axis |
||
− | pointAngle :: Point -> Point -> Double |
||
− | pointAngle (Point (x1, y1)) (Point (x2, y2)) = (x2 - x1) / (y2 - y1) |
||
+ | -- Scan the points to find out convex |
||
− | --compare based on point angle |
||
+ | -- -- handle the case that all points are collinear |
||
− | pointOrdering :: Point -> Point -> Ordering |
||
+ | scan [p0] (p1:ps) |
||
− | pointOrdering a b = compare (pointAngle a b) 0.0 |
||
+ | | isTurned pz p0 p1 == GoStraight = [pz, p0] |
||
+ | where pz = last ps |
||
+ | scan (x:xs) (y:z:rsts) = case isTurned x y z of |
||
− | --Sort by angle |
||
+ | GoRight -> scan xs (x:z:rsts) |
||
− | sortByAngle :: [Point] -> [Point] |
||
+ | GoStraight -> scan (x:xs) (z:rsts) -- I choose to skip the collinear points |
||
− | sortByAngle ps = bottomLeft : sortBy (compareAngles bottomLeft) (tail (sortedPs)) |
||
− | + | GoLeft -> scan (y:x:xs) (z:rsts) |
|
+ | scan xs [z] = z : xs |
||
− | bottomLeft = head (sortedPs) |
||
+ | |||
+ | |||
+ | -- 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 |
||
− | --Compare angles |
||
+ | pts5 = [Pt (5,5), Pt (4,4), Pt (3,3), Pt (2,2), Pt (1,1)] |
||
− | compareAngles :: Point -> Point -> Point -> Ordering |
||
+ | pts6 = reverse pts5 |
||
− | compareAngles = comparing . pointAngle |
||
− | |||
− | --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