1module Algo.Quickhull (quickhull) where
2
3import Data.Vector.Unboxed as V
4
5quickhull :: (Vector Double, Vector Double) -> (Vector Double, Vector Double)
6{-# NOINLINE quickhull #-}
7quickhull (xs, ys) = xs' `seq` ys' `seq` (xs',ys')
8    where
9      (xs',ys') = V.unzip
10                $ hsplit points pmin pmax V.++ hsplit points pmax pmin
11
12      imin = V.minIndex xs
13      imax = V.maxIndex xs
14
15      points = V.zip xs ys
16      pmin   = points V.! imin
17      pmax   = points V.! imax
18
19
20      hsplit points p1 p2
21        | V.length packed < 2 = p1 `V.cons` packed
22        | otherwise = hsplit packed p1 pm V.++ hsplit packed pm p2
23        where
24          cs     = V.map (\p -> cross p p1 p2) points
25          packed = V.map fst
26                 $ V.filter (\t -> snd t > 0)
27                 $ V.zip points cs
28
29          pm     = points V.! V.maxIndex cs
30
31      cross (x,y) (x1,y1) (x2,y2) = (x1-x)*(y2-y) - (y1-y)*(x2-x)
32
33