Add convex hull algorithm via graham scan

Algorithms/ConvexHull.hs

@@ -0,0 +1,66 @@
+{-# OPTIONS_HADDOCK ignore-exports #-}
+
+module Algorithms.ConvexHull where
+
+import Data.List
+import Diagrams.TwoD.Types
+import Diagrams.TwoD.Vector
+import Util
+import LinearAlgebra.Vector
+
+
+-- |Find the point with the lowest Y coordinate.
+-- If the lowest y-coordinate exists in more than one point in the set,
+-- the point with the lowest x-coordinate out of the candidates is
+-- chosen.
+lowestYC :: [PT] -> PT
+lowestYC []          = error "lowestYC: empty list"
+lowestYC [a]         = a
+lowestYC (a:b:vs)
+  | ay >  by     = lowestYC (b:vs)
+  | ay == by &&
+    ax >  bx     = lowestYC (b:vs)
+  | otherwise    = lowestYC (a:vs)
+    where
+      (ax, ay) = unp2 a
+      (bx, by) = unp2 b
+
+
+-- |Sort the points in increasing order of their degree between
+-- P0 and the x-axis.
+grahamSort :: [PT] -- ^ the points to sort
+           -> [PT] -- ^ sorted points
+grahamSort xs = p0 : sortBy (\a b
+  -> noEqual a b                            .
+       compare
+       (getAngle (pt2Vec a - pt2Vec p0) xv) $
+       (getAngle (pt2Vec b - pt2Vec p0) xv))
+  (removeItem p0 xs)
+    where
+      xv = unitX
+      p0 = lowestYC xs
+      -- Have to account for corner cases when points are in
+      -- a straight line or have the same y coordinates. Eq is
+      -- not an option anyhow.
+      noEqual :: PT -> PT -> Ordering -> Ordering
+      noEqual a b EQ
+        | ay == by &&
+          ax < bx      = LT
+        | otherwise    = GT
+          where
+            (ax, ay) = unp2 a
+            (bx, by) = unp2 b
+      noEqual _ _ LT     = LT
+      noEqual _ _ GT     = GT
+
+
+-- |Get all points on a convex hull by using the graham scan
+-- algorithm.
+grahamGetCH :: [PT] -> [PT]
+grahamGetCH vs = f . grahamSort $ vs
+  where
+    f [] = []
+    f (x:y:z:xs)
+      | ccw x y z = x : f (y:z:xs)
+      | otherwise = f (x:z:xs)
+    f xs          = xs

Diagram.hs

@@ -9,6 +9,7 @@ module Diagram (t,
                 diagS,
                 whiteRect) where
 
+import Algorithms.ConvexHull
 import Class.Defaults
 import Diagrams.Prelude
 import Diagrams.Backend.Cairo
@@ -102,6 +103,19 @@ showCoordinates = Diag f
               dot = (circle $ t p :: Diagram Cairo R2) # fc black
 
 
+-- |Create a diagram which shows the points of the convex hull.
+showConvexHullPoints :: Diag
+showConvexHullPoints = Diag f
+  where
+    f p vt
+      = position (zip (filter (inRange (dX p) (dY p)) $ vtch)
+                      (repeat dot))   # moveTo (p2(xOffset p, yOffset p))
+            where
+              -- a dot itself is a diagram
+              dot = (circle $ t p :: Diagram Cairo R2) # fc red # lc red
+              vtch = grahamGetCH vt
+
+
 -- |Creates a Diagram that shows an XAxis which is bound
 -- by the dimensions given in xD from DiagProp.
 showXAxis :: Diag
@@ -132,6 +146,10 @@ diag p = case alg p of
   0 -> mkDiag
          (mconcat [showCoordinates, showXAxis, showYAxis, showWhiteRectB])
          p
+  1 -> mkDiag
+         (mconcat [showConvexHullPoints, showCoordinates,
+                   showXAxis, showYAxis, showWhiteRectB])
+         p
   _ -> mempty