Wednesday, February 22, 2017

convex hulls

This code demonstrates the Graham scan, an O(n log n) method for finding the convex hull (smallest enclosing polygon) of a planar point set. It's a great example of exploiting order: it works by sorting the set by angle about a known extreme point, allowing the hull points to be found in linear time.

A variation on this algorithm, noted by A.M Andrew, sorts the set lexicographically and finds the upper and lower hull chains separately. This 'monotone chain' technique is often preferred, since it's easier to do robustly.

I've also written a couple JavaScript examples, the Graham scan and a slower (but constant-space and adaptable to higher dimensions) method called the Jarvis march.

{-# LANGUAGE NoMonomorphismRestriction #-}

import Control.Arrow ((***))
import Data.List (maximumBy, delete, sort, sortBy, unfoldr)
import Data.Ord (comparing)
import Graphics.UI.SDL as SDL
import Graphics.UI.SDL.Primitives (filledCircle, line)
import System.Random.Mersenne.Pure64 (newPureMT, randomDouble)

res = 250

main = withInit [InitVideo] $ do
  w  <- setVideoMode res res 32 []
  ps <- randPoints
  enableEvent SDLMouseMotion False
  setCaption "Graham Scan" "Graham Scan"
  fillRect w (Just $ Rect 0 0 res res) $ Pixel 0
  limn w $ map (round *** round) $ hull ps
  plot w ps
  pause

plot w ps = do
  mapM_ (f . (round *** round)) ps
  SDL.flip w
 where
  f (x,y) = filledCircle w x y 1 $ Pixel 0xFFFFFFFF

limn w ps = f $ ps ++ [head ps]
 where
  f [_] = return ()
  f ((a,b):(x,y):ps) = do
    line w a b x y $ Pixel 0xFF0000FF
    f $ (x,y) : ps

pause = do
  delay 128
  e <- pollEvent
  case e of
   KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
   _                              -> pause

hull qs = go (drop 2 ps) $ reverse $ take 2 ps
 where
  o  = bottomRightP qs
  ps = o : sortBy (ccw o) (delete o qs)

go [] qs = qs

go (p:ps) s@(a:b:qs)
  | ccw a b p /= GT = go (p:ps) $ b:qs
  | otherwise       = go ps $ p:s

ccw (ax, ay) (bx, by) (cx, cy)
  | d < 0 = LT
  | d > 0 = GT
  | True  = EQ
 where
  d = (bx - ax) * (cy - ay) - (by - ay) * (cx - ax)

bottomRightP = maximumBy (comparing snd) . sort

randPoints = fmap f newPureMT
 where
  f = uncurry zip . splitAt 20 . g
  g = map (* res) . unfoldr (Just . randomDouble)

No comments:

Post a Comment