Difference between revisions of "Graham Scan Implementation"
Jump to navigation
Jump to search
(Fixed bugs with the recursion) |
m (fixup Pt definition) |
||
(5 intermediate revisions by 2 users not shown) | |||
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 |
||
− | + | data Direction = LEFT | RIGHT | STRAIGHT |
|
+ | deriving (Show, Eq) |
||
− | data Direction = LeftTurn |
||
− | | RightTurn |
||
− | | Straight |
||
− | deriving (Show, Eq) |
||
+ | data Pt = Pt (Double, Double) |
||
− | --Point type |
||
+ | deriving (Show, Eq, Ord) |
||
− | data Point = Point (Double, Double) |
||
− | deriving (Show) |
||
+ | isTurned :: Pt -> Pt -> Pt -> Direction |
||
− | --some points |
||
+ | isTurned (Pt (ax, ay)) (Pt (bx, by)) (Pt (cx, cy)) = case sign of |
||
− | p0 = Point (2.1,2.0) |
||
+ | EQ -> STRAIGHT |
||
− | p1 = Point (4.2,2.0) |
||
+ | LT -> RIGHT |
||
− | p2 = Point (0.5,2.5) |
||
+ | GT -> LEFT |
||
− | p3 = Point (3.2,3.5) |
||
+ | where sign = compare ((bx - ax) * (cy - ay)) ((cx - ax) * (by - ay)) |
||
− | 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] |
||
+ | gScan :: [Pt] -> [Pt] |
||
− | -- Actually, I'd leave it as EQ, GT, LT. Then, actually, |
||
+ | gScan pts |
||
− | -- if you wanted to sort points rotationally around a single point, |
||
+ | | length pts >= 3 = scan [pt0] rests |
||
− | -- sortBy (dir x) would actually work. --wasserman.louis@gmail.com |
||
+ | | otherwise = pts |
||
− | --Get direction of single set of line segments |
||
+ | where |
||
− | dir :: Point -> Point -> Point -> Direction |
||
+ | -- Find the most bottom-left point pt0 |
||
− | dir (Point (ax, ay)) (Point (bx, by)) (Point (cx, cy)) = case sign of |
||
+ | pt0 = foldr bottomLeft (Pt (1/0, 1/0)) pts where |
||
− | EQ -> Straight |
||
− | + | bottomLeft pa pb = case ord of |
|
− | + | LT -> pa |
|
− | + | GT -> pb |
|
+ | EQ -> pa |
||
+ | where ord = (compare `on` (\ (Pt (x, y)) -> (y, x))) pa pb |
||
+ | -- Sort other points based on angle |
||
− | --Get a list of Directions from a list of Points |
||
+ | rests = tail (sortBy (compare `on` compkey pt0) pts) where |
||
− | dirlist :: [Point] -> [Direction] |
||
+ | compkey (Pt (x0, y0)) (Pt (x, y)) = (atan2 (y - y0) (x - x0), |
||
− | dirlist (x:y:z:xs) = dir x y z : dirlist (y:z:xs) |
||
+ | abs (x - x0)) |
||
− | dirlist _ = [] |
||
+ | -- Scan the points to find out convex |
||
− | --Compare Y axes |
||
+ | -- -- handle the case that all points are collinear |
||
− | sortByY :: [Point] -> [Point] |
||
+ | scan [p0] (p1:ps) |
||
− | sortByY xs = sortBy lowestY xs |
||
− | + | | isTurned pz p0 p1 == STRAIGHT = [pz, p0] |
|
+ | where pz = last ps |
||
− | --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 (x:xs) (y:z:rsts) = case isTurned x y z of |
||
− | --compare based on point angle |
||
+ | RIGHT -> scan xs (x:z:rsts) |
||
− | pointOrdering :: Point -> Point -> Ordering |
||
+ | STRAIGHT -> scan (x:xs) (z:rsts) -- skip collinear points |
||
− | pointOrdering a b = compare (pointAngle a b) 0.0 |
||
+ | LEFT -> scan (y:x:xs) (z:rsts) |
||
+ | scan xs [z] = z : xs |
||
− | --Sort by angle |
||
+ | </haskell> |
||
− | sortByAngle :: [Point] -> [Point] |
||
− | sortByAngle ps = bottomLeft : sortBy (compareAngles bottomLeft) (tail (sortedPs)) |
||
− | where sortedPs = sortByY ps |
||
− | bottomLeft = head (sortedPs) |
||
− | |||
+ | == Test == |
||
− | --Compare angles |
||
+ | <haskell> |
||
− | compareAngles :: Point -> Point -> Point -> Ordering |
||
+ | import Test.QuickCheck |
||
− | compareAngles = comparing . pointAngle |
||
+ | import Control.Monad |
||
+ | import Data.List |
||
+ | |||
− | --Graham Scan |
||
+ | convex' = map (\(Pt x) -> x) . convex . map Pt |
||
− | gscan :: [Point] -> [Point] |
||
+ | |||
− | gscan ps = scan (sortByAngle ps) |
||
+ | prop_convex n = forAll (replicateM n arbitrary) (\xs -> sort (convex' xs) == sort (convex' (convex' xs))) |
||
− | where scan (x:y:z:xs) = if dir x y z == LeftTurn |
||
+ | |||
− | then scan (x:z:xs) |
||
+ | </haskell> |
||
− | else x: scan (y:z:xs) |
||
+ | |||
− | scan [x,y] = [x,y] -- there's no shame in a pattern match |
||
+ | The results are: |
||
− | -- of this type! |
||
+ | |||
− | scan _ = [] |
||
+ | <haskell> |
||
+ | > quickCheck (prop_convex 0) |
||
+ | +++ OK, passed 100 tests. |
||
+ | > quickCheck (prop_convex 1) |
||
+ | +++ OK, passed 100 tests. |
||
+ | > quickCheck (prop_convex 2) |
||
+ | +++ OK, passed 100 tests. |
||
+ | > quickCheck (prop_convex 3) |
||
+ | +++ OK, passed 100 tests. |
||
+ | > quickCheck (prop_convex 4) |
||
+ | +++ OK, passed 100 tests. |
||
+ | > quickCheck (prop_convex 5) |
||
+ | +++ OK, passed 100 tests. |
||
+ | > quickCheck (prop_convex 10) |
||
+ | +++ OK, passed 100 tests. |
||
+ | > quickCheck (prop_convex 100) |
||
+ | +++ OK, passed 100 tests. |
||
+ | > quickCheck (prop_convex 1000) |
||
+ | +++ OK, passed 100 tests. |
||
+ | > quickCheck (prop_convex 10000) |
||
+ | +++ OK, passed 100 tests. |
||
</haskell> |
</haskell> |
Revision as of 14:36, 5 January 2015
Descriptions of this problem can be found in Real World Haskell, Chapter 3
import Data.List (sortBy)
import Data.Function (on)
data Direction = LEFT | RIGHT | STRAIGHT
deriving (Show, Eq)
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 -> STRAIGHT
LT -> RIGHT
GT -> LEFT
where sign = compare ((bx - ax) * (cy - ay)) ((cx - ax) * (by - ay))
gScan :: [Pt] -> [Pt]
gScan pts
| length pts >= 3 = scan [pt0] rests
| otherwise = pts
where
-- Find the most bottom-left point pt0
pt0 = foldr bottomLeft (Pt (1/0, 1/0)) pts where
bottomLeft pa pb = case ord of
LT -> pa
GT -> pb
EQ -> pa
where ord = (compare `on` (\ (Pt (x, y)) -> (y, x))) pa pb
-- Sort other points based on angle
rests = tail (sortBy (compare `on` compkey pt0) pts) where
compkey (Pt (x0, y0)) (Pt (x, y)) = (atan2 (y - y0) (x - x0),
abs (x - x0))
-- Scan the points to find out convex
-- -- handle the case that all points are collinear
scan [p0] (p1:ps)
| isTurned pz p0 p1 == STRAIGHT = [pz, p0]
where pz = last ps
scan (x:xs) (y:z:rsts) = case isTurned x y z of
RIGHT -> scan xs (x:z:rsts)
STRAIGHT -> scan (x:xs) (z:rsts) -- skip collinear points
LEFT -> scan (y:x:xs) (z:rsts)
scan xs [z] = z : xs
Test
import Test.QuickCheck
import Control.Monad
import Data.List
convex' = map (\(Pt x) -> x) . convex . map Pt
prop_convex n = forAll (replicateM n arbitrary) (\xs -> sort (convex' xs) == sort (convex' (convex' xs)))
The results are:
> quickCheck (prop_convex 0)
+++ OK, passed 100 tests.
> quickCheck (prop_convex 1)
+++ OK, passed 100 tests.
> quickCheck (prop_convex 2)
+++ OK, passed 100 tests.
> quickCheck (prop_convex 3)
+++ OK, passed 100 tests.
> quickCheck (prop_convex 4)
+++ OK, passed 100 tests.
> quickCheck (prop_convex 5)
+++ OK, passed 100 tests.
> quickCheck (prop_convex 10)
+++ OK, passed 100 tests.
> quickCheck (prop_convex 100)
+++ OK, passed 100 tests.
> quickCheck (prop_convex 1000)
+++ OK, passed 100 tests.
> quickCheck (prop_convex 10000)
+++ OK, passed 100 tests.