--------------------------------------------------------------------------------
-- |
-- Module      :  Algorithms.Geometry.ConvexHull.DivideAndConquer
-- Copyright   :  (C) Frank Staals
-- License     :  see the LICENSE file
-- Maintainer  :  Frank Staals
--
-- \(O(n\log n)\) time divide and conquer algorithm to compute the convex hull
-- of a set of \(n\) points in \(\mathbb{R}^2\).
--
--------------------------------------------------------------------------------
module Algorithms.Geometry.ConvexHull.DivideAndConquer( convexHull
                                                      , upperHull
                                                      , lowerHull
                                                      ) where

import           Algorithms.DivideAndConquer
import           Control.Arrow ((&&&))
import           Data.Ext
import           Data.Geometry.Point
import           Data.Geometry.Polygon
import           Data.Geometry.Polygon.Convex
import           Data.List.NonEmpty (NonEmpty(..))
import qualified Data.List.NonEmpty as NonEmpty
import           Data.Util
--------------------------------------------------------------------------------

-- | \(O(n \log n)\) time ConvexHull using divide and conquer. The resulting polygon is
-- given in clockwise order.
convexHull           :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> ConvexPolygon p r
convexHull :: NonEmpty (Point 2 r :+ p) -> ConvexPolygon p r
convexHull (Point 2 r :+ p
p :| []) = SimplePolygon p r -> ConvexPolygon p r
forall p r. SimplePolygon p r -> ConvexPolygon p r
ConvexPolygon (SimplePolygon p r -> ConvexPolygon p r)
-> ([Point 2 r :+ p] -> SimplePolygon p r)
-> [Point 2 r :+ p]
-> ConvexPolygon p r
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Point 2 r :+ p] -> SimplePolygon p r
forall r p. [Point 2 r :+ p] -> SimplePolygon p r
unsafeFromPoints ([Point 2 r :+ p] -> ConvexPolygon p r)
-> [Point 2 r :+ p] -> ConvexPolygon p r
forall a b. (a -> b) -> a -> b
$ [Point 2 r :+ p
p]
convexHull NonEmpty (Point 2 r :+ p)
pts       = (NonEmpty (Point 2 r :+ p), NonEmpty (Point 2 r :+ p))
-> ConvexPolygon p r
forall r p.
(NonEmpty (Point 2 r :+ p), NonEmpty (Point 2 r :+ p))
-> ConvexPolygon p r
combine ((NonEmpty (Point 2 r :+ p), NonEmpty (Point 2 r :+ p))
 -> ConvexPolygon p r)
-> (NonEmpty (Point 2 r :+ p)
    -> (NonEmpty (Point 2 r :+ p), NonEmpty (Point 2 r :+ p)))
-> NonEmpty (Point 2 r :+ p)
-> ConvexPolygon p r
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
forall r p.
(Ord r, Num r) =>
NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
upperHull' (NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p))
-> (NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p))
-> NonEmpty (Point 2 r :+ p)
-> (NonEmpty (Point 2 r :+ p), NonEmpty (Point 2 r :+ p))
forall (a :: * -> * -> *) b c c'.
Arrow a =>
a b c -> a b c' -> a b (c, c')
&&& NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
forall r p.
(Ord r, Num r) =>
NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
lowerHull') (NonEmpty (Point 2 r :+ p)
 -> (NonEmpty (Point 2 r :+ p), NonEmpty (Point 2 r :+ p)))
-> (NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p))
-> NonEmpty (Point 2 r :+ p)
-> (NonEmpty (Point 2 r :+ p), NonEmpty (Point 2 r :+ p))
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Point 2 r :+ p) -> (Point 2 r :+ p) -> Ordering)
-> NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
forall a. (a -> a -> Ordering) -> NonEmpty a -> NonEmpty a
NonEmpty.sortBy (Point 2 r :+ p) -> (Point 2 r :+ p) -> Ordering
forall r p q.
Ord r =>
(Point 2 r :+ p) -> (Point 2 r :+ q) -> Ordering
incXdecY (NonEmpty (Point 2 r :+ p) -> ConvexPolygon p r)
-> NonEmpty (Point 2 r :+ p) -> ConvexPolygon p r
forall a b. (a -> b) -> a -> b
$ NonEmpty (Point 2 r :+ p)
pts
  where
    combine :: (NonEmpty (Point 2 r :+ p), NonEmpty (Point 2 r :+ p))
-> ConvexPolygon p r
combine (Point 2 r :+ p
l:|[Point 2 r :+ p]
uh,Point 2 r :+ p
_:|[Point 2 r :+ p]
lh) = SimplePolygon p r -> ConvexPolygon p r
forall p r. SimplePolygon p r -> ConvexPolygon p r
ConvexPolygon (SimplePolygon p r -> ConvexPolygon p r)
-> ([Point 2 r :+ p] -> SimplePolygon p r)
-> [Point 2 r :+ p]
-> ConvexPolygon p r
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Point 2 r :+ p] -> SimplePolygon p r
forall r p. [Point 2 r :+ p] -> SimplePolygon p r
unsafeFromPoints ([Point 2 r :+ p] -> ConvexPolygon p r)
-> [Point 2 r :+ p] -> ConvexPolygon p r
forall a b. (a -> b) -> a -> b
$ Point 2 r :+ p
l (Point 2 r :+ p) -> [Point 2 r :+ p] -> [Point 2 r :+ p]
forall a. a -> [a] -> [a]
: [Point 2 r :+ p]
uh [Point 2 r :+ p] -> [Point 2 r :+ p] -> [Point 2 r :+ p]
forall a. Semigroup a => a -> a -> a
<> [Point 2 r :+ p] -> [Point 2 r :+ p]
forall a. [a] -> [a]
reverse ([Point 2 r :+ p] -> [Point 2 r :+ p]
forall a. [a] -> [a]
init [Point 2 r :+ p]
lh)

----------------------------------------
-- * Computing a lower hull

-- | \(O(n \log n)\) time LowerHull using divide and conquer. The resulting Hull is
-- given from left to right, i.e. in counter clockwise order.
lowerHull :: (Ord r, Num r)
          => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
lowerHull :: NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
lowerHull = NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
forall r p.
(Ord r, Num r) =>
NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
lowerHull' (NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p))
-> (NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p))
-> NonEmpty (Point 2 r :+ p)
-> NonEmpty (Point 2 r :+ p)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Point 2 r :+ p) -> (Point 2 r :+ p) -> Ordering)
-> NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
forall a. (a -> a -> Ordering) -> NonEmpty a -> NonEmpty a
NonEmpty.sortBy (Point 2 r :+ p) -> (Point 2 r :+ p) -> Ordering
forall r p q.
Ord r =>
(Point 2 r :+ p) -> (Point 2 r :+ q) -> Ordering
incXdecY

lowerHull' :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
lowerHull' :: NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
lowerHull' = LH r p -> NonEmpty (Point 2 r :+ p)
forall r p. LH r p -> NonEmpty (Point 2 r :+ p)
unLH (LH r p -> NonEmpty (Point 2 r :+ p))
-> (NonEmpty (Point 2 r :+ p) -> LH r p)
-> NonEmpty (Point 2 r :+ p)
-> NonEmpty (Point 2 r :+ p)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Point 2 r :+ p) -> LH r p) -> NonEmpty (Point 2 r :+ p) -> LH r p
forall (f :: * -> *) s a.
(Foldable1 f, Semigroup s) =>
(a -> s) -> f a -> s
divideAndConquer1 (NonEmpty (Point 2 r :+ p) -> LH r p
forall r p. NonEmpty (Point 2 r :+ p) -> LH r p
LH (NonEmpty (Point 2 r :+ p) -> LH r p)
-> ((Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p))
-> (Point 2 r :+ p)
-> LH r p
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Point 2 r :+ p) -> [Point 2 r :+ p] -> NonEmpty (Point 2 r :+ p)
forall a. a -> [a] -> NonEmpty a
:|[]))

newtype LH r p = LH { LH r p -> NonEmpty (Point 2 r :+ p)
unLH :: NonEmpty (Point 2 r :+ p) } deriving (LH r p -> LH r p -> Bool
(LH r p -> LH r p -> Bool)
-> (LH r p -> LH r p -> Bool) -> Eq (LH r p)
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
forall r p. (Eq r, Eq p) => LH r p -> LH r p -> Bool
/= :: LH r p -> LH r p -> Bool
$c/= :: forall r p. (Eq r, Eq p) => LH r p -> LH r p -> Bool
== :: LH r p -> LH r p -> Bool
$c== :: forall r p. (Eq r, Eq p) => LH r p -> LH r p -> Bool
Eq,Int -> LH r p -> ShowS
[LH r p] -> ShowS
LH r p -> String
(Int -> LH r p -> ShowS)
-> (LH r p -> String) -> ([LH r p] -> ShowS) -> Show (LH r p)
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
forall r p. (Show r, Show p) => Int -> LH r p -> ShowS
forall r p. (Show r, Show p) => [LH r p] -> ShowS
forall r p. (Show r, Show p) => LH r p -> String
showList :: [LH r p] -> ShowS
$cshowList :: forall r p. (Show r, Show p) => [LH r p] -> ShowS
show :: LH r p -> String
$cshow :: forall r p. (Show r, Show p) => LH r p -> String
showsPrec :: Int -> LH r p -> ShowS
$cshowsPrec :: forall r p. (Show r, Show p) => Int -> LH r p -> ShowS
Show)

instance (Num r, Ord r) => Semigroup (LH r p) where
  (LH NonEmpty (Point 2 r :+ p)
lh) <> :: LH r p -> LH r p -> LH r p
<> (LH NonEmpty (Point 2 r :+ p)
rh) = NonEmpty (Point 2 r :+ p) -> LH r p
forall r p. NonEmpty (Point 2 r :+ p) -> LH r p
LH (NonEmpty (Point 2 r :+ p) -> LH r p)
-> NonEmpty (Point 2 r :+ p) -> LH r p
forall a b. (a -> b) -> a -> b
$ (NonEmpty (Point 2 r :+ p)
 -> NonEmpty (Point 2 r :+ p)
 -> Two ((Point 2 r :+ p) :+ [Point 2 r :+ p]))
-> NonEmpty (Point 2 r :+ p)
-> NonEmpty (Point 2 r :+ p)
-> NonEmpty (Point 2 r :+ p)
forall p.
(NonEmpty p -> NonEmpty p -> Two (p :+ [p]))
-> NonEmpty p -> NonEmpty p -> NonEmpty p
hull NonEmpty (Point 2 r :+ p)
-> NonEmpty (Point 2 r :+ p)
-> Two ((Point 2 r :+ p) :+ [Point 2 r :+ p])
forall r (f :: * -> *) p.
(Ord r, Num r, Foldable1 f) =>
f (Point 2 r :+ p)
-> f (Point 2 r :+ p) -> Two ((Point 2 r :+ p) :+ [Point 2 r :+ p])
lowerTangent' NonEmpty (Point 2 r :+ p)
lh NonEmpty (Point 2 r :+ p)
rh

----------------------------------------
-- * Computing an upper hull

-- | \(O(n \log n)\) time UpperHull using divide and conquer. The resulting Hull is
-- given from left to right, i.e. in clockwise order.
upperHull :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
upperHull :: NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
upperHull = NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
forall r p.
(Ord r, Num r) =>
NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
upperHull' (NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p))
-> (NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p))
-> NonEmpty (Point 2 r :+ p)
-> NonEmpty (Point 2 r :+ p)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Point 2 r :+ p) -> (Point 2 r :+ p) -> Ordering)
-> NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
forall a. (a -> a -> Ordering) -> NonEmpty a -> NonEmpty a
NonEmpty.sortBy (Point 2 r :+ p) -> (Point 2 r :+ p) -> Ordering
forall r p q.
Ord r =>
(Point 2 r :+ p) -> (Point 2 r :+ q) -> Ordering
incXdecY

upperHull' :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
upperHull' :: NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
upperHull' = UH r p -> NonEmpty (Point 2 r :+ p)
forall r p. UH r p -> NonEmpty (Point 2 r :+ p)
unUH (UH r p -> NonEmpty (Point 2 r :+ p))
-> (NonEmpty (Point 2 r :+ p) -> UH r p)
-> NonEmpty (Point 2 r :+ p)
-> NonEmpty (Point 2 r :+ p)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Point 2 r :+ p) -> UH r p) -> NonEmpty (Point 2 r :+ p) -> UH r p
forall (f :: * -> *) s a.
(Foldable1 f, Semigroup s) =>
(a -> s) -> f a -> s
divideAndConquer1 (NonEmpty (Point 2 r :+ p) -> UH r p
forall r p. NonEmpty (Point 2 r :+ p) -> UH r p
UH (NonEmpty (Point 2 r :+ p) -> UH r p)
-> ((Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p))
-> (Point 2 r :+ p)
-> UH r p
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Point 2 r :+ p) -> [Point 2 r :+ p] -> NonEmpty (Point 2 r :+ p)
forall a. a -> [a] -> NonEmpty a
:|[]))

newtype UH r p = UH { UH r p -> NonEmpty (Point 2 r :+ p)
unUH :: NonEmpty (Point 2 r :+ p) }

instance (Num r, Ord r) => Semigroup (UH r p) where
  (UH NonEmpty (Point 2 r :+ p)
lh) <> :: UH r p -> UH r p -> UH r p
<> (UH NonEmpty (Point 2 r :+ p)
rh) = NonEmpty (Point 2 r :+ p) -> UH r p
forall r p. NonEmpty (Point 2 r :+ p) -> UH r p
UH (NonEmpty (Point 2 r :+ p) -> UH r p)
-> NonEmpty (Point 2 r :+ p) -> UH r p
forall a b. (a -> b) -> a -> b
$ (NonEmpty (Point 2 r :+ p)
 -> NonEmpty (Point 2 r :+ p)
 -> Two ((Point 2 r :+ p) :+ [Point 2 r :+ p]))
-> NonEmpty (Point 2 r :+ p)
-> NonEmpty (Point 2 r :+ p)
-> NonEmpty (Point 2 r :+ p)
forall p.
(NonEmpty p -> NonEmpty p -> Two (p :+ [p]))
-> NonEmpty p -> NonEmpty p -> NonEmpty p
hull NonEmpty (Point 2 r :+ p)
-> NonEmpty (Point 2 r :+ p)
-> Two ((Point 2 r :+ p) :+ [Point 2 r :+ p])
forall r (f :: * -> *) p.
(Ord r, Num r, Foldable1 f) =>
f (Point 2 r :+ p)
-> f (Point 2 r :+ p) -> Two ((Point 2 r :+ p) :+ [Point 2 r :+ p])
upperTangent' NonEmpty (Point 2 r :+ p)
lh NonEmpty (Point 2 r :+ p)
rh

----------------------------------------

-- | The function that does the actual merging part
hull               :: (NonEmpty p -> NonEmpty p -> Two (p :+ [p]))
                   -> NonEmpty p -> NonEmpty p -> NonEmpty p
hull :: (NonEmpty p -> NonEmpty p -> Two (p :+ [p]))
-> NonEmpty p -> NonEmpty p -> NonEmpty p
hull NonEmpty p -> NonEmpty p -> Two (p :+ [p])
tangent NonEmpty p
lh NonEmpty p
rh = let Two (p
l :+ [p]
lh') (p
r :+ [p]
rh') = NonEmpty p -> NonEmpty p -> Two (p :+ [p])
tangent (NonEmpty p -> NonEmpty p
forall a. NonEmpty a -> NonEmpty a
NonEmpty.reverse NonEmpty p
lh) NonEmpty p
rh
                     in [p] -> NonEmpty p
forall a. [a] -> NonEmpty a
NonEmpty.fromList ([p] -> NonEmpty p) -> [p] -> NonEmpty p
forall a b. (a -> b) -> a -> b
$ [p] -> [p]
forall a. [a] -> [a]
reverse [p]
lh' [p] -> [p] -> [p]
forall a. Semigroup a => a -> a -> a
<> [p
l,p
r] [p] -> [p] -> [p]
forall a. Semigroup a => a -> a -> a
<> [p]
rh'

--------------------------------------------------------------------------------

incXdecY  :: Ord r => Point 2 r :+ p -> Point 2 r :+ q -> Ordering
incXdecY :: (Point 2 r :+ p) -> (Point 2 r :+ q) -> Ordering
incXdecY (Point2 r
px r
py :+ p
_) (Point2 r
qx r
qy :+ q
_) =
  r -> r -> Ordering
forall a. Ord a => a -> a -> Ordering
compare r
px r
qx Ordering -> Ordering -> Ordering
forall a. Semigroup a => a -> a -> a
<> r -> r -> Ordering
forall a. Ord a => a -> a -> Ordering
compare r
qy r
py