Graham Scan Implementation

From HaskellWiki
Revision as of 14:30, 5 January 2015 by Robturtle (talk | contribs) (make the implementation neat)
Jump to navigation Jump to search
The printable version is no longer supported and may have rendering errors. Please update your browser bookmarks and please use the default browser print function instead.

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.