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