Difference between revisions of "Graham Scan Implementation"

From HaskellWiki
Jump to: navigation, search
(Add test code)
(make the implementation neat)
Line 5: Line 5:
 
import Data.Function (on)
 
import Data.Function (on)
   
-- define data `Direction`
+
data Direction = LEFT | RIGHT | STRAIGHT
data Direction = GoLeft | GoRight | GoStraight
+
deriving (Show, Eq)
deriving (Show, Eq)
+
  +
data Pt = Pt (Double, Double)
   
-- 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 -> Pt -> Pt -> Direction
 
isTurned (Pt (ax, ay)) (Pt (bx, by)) (Pt (cx, cy)) = case sign of
 
isTurned (Pt (ax, ay)) (Pt (bx, by)) (Pt (cx, cy)) = case sign of
EQ -> GoStraight
+
EQ -> STRAIGHT
LT -> GoRight
+
LT -> RIGHT
GT -> GoLeft
+
GT -> LEFT
where sign = compare ((bx - ax) * (cy - ay))
+
where sign = compare ((bx - ax) * (cy - ay)) ((cx - ax) * (by - 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>
 
  +
gScan :: [Pt] -> [Pt]
spts = tail (sortBy (compare `on` compkey pt0) pts) where
 
  +
gScan pts
compkey (Pt (ax, ay)) (Pt (bx, by)) = (atan2 (by - ay) (bx - ax),
 
  +
| length pts >= 3 = scan [pt0] rests
{-the secondary key make sure collinear points in order-}
 
  +
| otherwise = pts
abs (bx - ax))
 
  +
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
   
-- Scan the points to find out convex
+
-- Sort other points based on angle
-- -- handle the case that all points are collinear
+
rests = tail (sortBy (compare `on` compkey pt0) pts) where
scan [p0] (p1:ps)
+
compkey (Pt (x0, y0)) (Pt (x, y)) = (atan2 (y - y0) (x - x0),
| isTurned pz p0 p1 == GoStraight = [pz, p0]
+
abs (x - x0))
where pz = last ps
 
   
scan (x:xs) (y:z:rsts) = case isTurned x y z of
 
  +
-- Scan the points to find out convex
GoRight -> scan xs (x:z:rsts)
 
  +
-- -- handle the case that all points are collinear
GoStraight -> scan (x:xs) (z:rsts) -- I choose to skip the collinear points
 
  +
scan [p0] (p1:ps)
GoLeft -> scan (y:x:xs) (z:rsts)
 
  +
| isTurned pz p0 p1 == STRAIGHT = [pz, p0]
scan xs [z] = z : xs
+
where pz = last ps
 
 
-- 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)]
 
  +
scan (x:xs) (y:z:rsts) = case isTurned x y z of
-- -- This case demonstrates a consecutive inner points
 
  +
RIGHT -> scan xs (x:z:rsts)
  +
STRAIGHT -> scan (x:xs) (z:rsts) -- skip collinear points
  +
LEFT -> scan (y:x:xs) (z:rsts)
   
-- collinear points
 
  +
scan xs [z] = z : xs
pts5 = [Pt (5,5), Pt (4,4), Pt (3,3), Pt (2,2), Pt (1,1)]
 
pts6 = reverse pts5
 
 
</haskell>
 
</haskell>
   

Revision as of 14:30, 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)

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.