GrahamScan.hs

--------------------------------------------------------------------------------
-- |
-- Module      :  Algorithms.Geometry.ConvexHull.GrahamScan
-- Copyright   :  (C) Frank Staals
-- License     :  see the LICENSE file
-- Maintainer  :  Frank Staals
--------------------------------------------------------------------------------
module Algorithms.Geometry.ConvexHull.GrahamScan( convexHull
                                                , upperHull, upperHull'
                                                , lowerHull, lowerHull'

                                                , upperHullFromSorted, upperHullFromSorted'
                                                ) where

import           Control.Lens ((^.))
import           Data.Ext
import           Data.Geometry.Point
import           Data.Geometry.Polygon
import           Data.Geometry.Polygon.Convex (ConvexPolygon(..))
import qualified Data.List.NonEmpty as NonEmpty
import           Data.List.NonEmpty (NonEmpty(..))


-- | \(O(n \log n)\) time ConvexHull using Graham-Scan. The resulting polygon is
-- given in clockwise order.
convexHull            :: (Ord r, Num r)
                      => NonEmpty (Point 2 r :+ p) -> ConvexPolygon p r
convexHull (p :| []) = ConvexPolygon . unsafeFromPoints $ [p]
convexHull ps        = let ps' = NonEmpty.toList . NonEmpty.sortBy incXdecY $ ps
                           uh  = NonEmpty.tail . hull' $         ps'
                           lh  = NonEmpty.tail . hull' $ reverse ps'
                       in ConvexPolygon . unsafeFromPoints . reverse $ lh ++ uh

-- | Computes the upper hull. The upper hull is given from left to right.
--
-- Specifically. A pair of points defines an edge of the upper hull
-- iff all other points are strictly to the right of its supporting
-- line.
--
-- Note that this definition implies that the segment may be
-- vertical. Use 'upperHull'' if such an edge should not be reported.
--
-- running time: \(O(n\log n)\)
upperHull  :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
upperHull = NonEmpty.reverse . hull id

-- | Computes the upper hull, making sure that there are no vertical segments.
--
-- The upper hull is given from left to right
--
upperHull'  :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
upperHull' = NonEmpty.reverse . dropVertical . hull id

-- | Helper function to remove vertical segments from the hull.
--
-- Tests if the first two points are on a vertical line, if so removes
-- the first point.
dropVertical :: Eq r => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
dropVertical = \case
  h@(_ :| [])                                            -> h
  h@(p :| (q : rest)) | p^.core.xCoord == q^.core.xCoord -> q :| rest
                      | otherwise                        -> h


-- | Computes the upper hull. The upper hull is given from left to right.
--
-- Specifically. A pair of points defines an edge of the lower hull
-- iff all other points are strictly to the left of its supporting
-- line.
--
-- Note that this definition implies that the segment may be
-- vertical. Use 'lowerHull'' if such an edge should not be reported.
--
-- running time: \(O(n\log n)\)
lowerHull :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
lowerHull = hull reverse

-- | Computes the lower hull, making sure there are no vertical
-- segments. (Note that the only such segment could be the first
-- segment).
lowerHull' :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
lowerHull' = dropVertical . hull reverse

-- | Helper function so that that can compute both the upper or the lower hull, depending
-- on the function f
hull               :: (Ord r, Num r)
                   => ([Point 2 r :+ p] -> [Point 2 r :+ p])
                   -> NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
hull _ h@(_ :| []) = h
hull f pts         = hull' .  f
                   . NonEmpty.toList . NonEmpty.sortBy incXdecY $ pts

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


-- | Given a sequence of points that is sorted on increasing
-- x-coordinate and decreasing y-coordinate, computes the upper
-- hull, in *right to left order*.
--
-- Specifically. A pair of points defines an edge of the upper hull
-- iff all other points are strictly to the right of its supporting
-- line.
--
--
-- Note that In constrast to the 'upperHull' function, the result is
-- returned *from right to left* !!!
--
-- running time: \(O(n)\).
upperHullFromSorted :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
upperHullFromSorted = \case
  h@(_ :| [])  -> h
  pts          -> hull' $ NonEmpty.toList pts

-- | Computes the upper hull from a sorted input. Removes the last vertical segment.
--
--
-- running time: \(O(n)\).
upperHullFromSorted' :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)
upperHullFromSorted' = dropVertical . upperHullFromSorted


-- | Precondition: The list of input points is sorted
hull'          :: (Ord r, Num r) => [Point 2 r :+ p] -> NonEmpty (Point 2 r :+ p)
hull' (a:b:ps) = NonEmpty.fromList $ hull'' [b,a] ps
  where
    hull'' h []      = h
    hull'' h (p:ps') = hull'' (cleanMiddle (p:h)) ps'

    cleanMiddle h@[_,_]                         = h
    cleanMiddle h@(z:y:x:rest)
      | rightTurn (x^.core) (y^.core) (z^.core) = h
      | otherwise                               = cleanMiddle (z:x:rest)
    cleanMiddle _                               = error "cleanMiddle: too few points"
hull' _ = error
  "Algorithms.Geometry.ConvexHull.GrahamScan.hull' requires a list with at least \
  \two elements."

rightTurn       :: (Ord r, Num r) => Point 2 r -> Point 2 r -> Point 2 r -> Bool
rightTurn a b c = ccw a b c == CW