hgeometry-0.5.0.0: src/Data/Geometry/Polygon/Convex.hs
{-# LANGUAGE ScopedTypeVariables #-}
{-|
Module : Data.Geometry.Polygon.Convex
Description: Convex Polygons
Copyright : (c) Frank Staals
License : See LICENCE file
-}
module Data.Geometry.Polygon.Convex( ConvexPolygon
, merge
, lowerTangent, upperTangent
, isLeftOf, isRightOf
, extremes
, maxInDirection
) where
import Control.Lens hiding ((:<), (:>))
import qualified Data.CircularSeq as C
import Data.CircularSeq (focus,CSeq)
import Data.Ext
import qualified Data.Foldable as F
import Data.Function (on, )
import Data.Geometry
import Data.Geometry.Polygon (fromPoints, cmpExtreme)
import Data.Maybe (fromJust)
import Data.Ord (comparing)
import qualified Data.Sequence as S
import Data.Sequence (viewl,viewr, ViewL(..), ViewR(..))
import Data.Geometry.Ipe
import Debug.Trace
--------------------------------------------------------------------------------
type ConvexPolygon = SimplePolygon
mainWith inFile outFile = do
ePage <- readSinglePageFile inFile
case ePage of
Left err -> error "" -- err
Right (page :: IpePage Rational) -> case page^..content.traverse._withAttrs _IpePath _asSimplePolygon.core of
[] -> error "No points found"
polies@(_:_) -> do
-- let out = [asIpe drawTriangulation dt, asIpe drawTree' emst]
-- print $ length $ edges' dt
-- print $ toPlaneGraph (Proxy :: Proxy DT) dt
-- writeIpeFile outFile . singlePageFromContent $ out
-- mapM_ (print . extremesNaive (v2 1 0)) polies
pure $ map (maxInDirection (v2 (-1) 0)) polies
-- | Finds the extreme points, minimum and maximum, in a given direction
--
-- pre: The input polygon is strictly convex.
--
-- running time: $O(\log n)$
--
--
extremes :: (Num r, Ord r) => Vector 2 r -> ConvexPolygon p r
-> (Point 2 r :+ p, Point 2 r :+ p)
extremes u p = (maxInDirection ((-1) *^ u) p, maxInDirection u p)
-- | Finds the extreme maximum point in the given direction. Based on
-- http://geomalgorithms.com/a14-_extreme_pts.html
--
--
-- pre: The input polygon is strictly convex.
--
-- running time: $O(\log^2 n)$
maxInDirection :: (Num r, Ord r) => Vector 2 r -> ConvexPolygon p r -> Point 2 r :+ p
maxInDirection u p = findMaxStart . C.rightElements $ p^.outerBoundary
where
findMaxStart s@(viewl -> (a:<r))
| isLocalMax r a r = a
| otherwise = findMax s
findMaxStart _ = error "absurd"
findMax s = let i = F.length s `div` 2
(ac,cb') = S.splitAt i s
(c :< cb) = viewl cb'
in findMax' ac c cb
findMax' ac c cb
| isLocalMax ac c cb = c
| otherwise = binSearch ac c cb
-- | Given the vertices [a..] c [..b] find the exteral vtx
binSearch ac@(viewl -> a:<r) c cb = case (isUpwards a r, isUpwards c cb, a >=. c) of
(True,False,_) -> findMax (ac |> c)
(True,True,True) -> findMax (ac |> c)
(True,True,False) -> findMax (c <| cb)
(False,True,_) -> findMax (c <| cb)
(False,False,False) -> findMax (ac |> c)
(False,False,True) -> findMax (c <| cb)
binSearch _ _ _ = error "maxInDirection, binSearch: empty chain"
isLocalMax (viewr -> _ :> l) c (viewl -> r :< _) = c >=. l && c >=. r
isLocalMax (viewr -> _ :> l) c _ = c >=. l
isLocalMax _ c (viewl -> r :< _) = c >=. r
isLocalMax _ _ _ = True
-- the Edge from a to b is upwards w.r.t b if a is not larger than b
isUpwards a (viewl -> b :< _) = cmpExtreme u a b /= GT
isUpwards _ _ = error "isUpwards: no edge endpoint"
p' >=. q = cmpExtreme u p' q /= LT
-- | Given a convex polygon poly, and a point outside the polygon, find the
-- right tangent of q and the polygon, i.e. the vertex v of the convex polygon
-- s.t. the polygon lies completely to the right of the line from q to v.
--
-- running time: $O(\log^2 n)$.
-- rightTangent :: ConvexPolygon p r -> Point 2 r -> Point 2 r :+ p
-- rightTangent poly q = undefined
-- TODO: same as maxInDirection, but with a slightly different ordering
-- * Merging Two convex Hulls
-- | Rotating Right <-> rotate clockwise
--
-- Merging two convex hulls, based on the paper:
--
-- Two Algorithms for Constructing a Delaunay Triangulation
-- Lee and Schachter
-- International Journal of Computer and Information Sciences, Vol 9, No. 3, 1980
--
-- : (combined hull, lower tangent that was added, upper tangent thtat was
-- added)
-- pre: - lp and rp are disjoint, and there is a vertical line separating
-- the two polygons.
-- - The vertices of the polygons are given in clockwise order
--
-- Running time: O(n+m), where n and m are the sizes of the two polygons respectively
merge :: (Num r, Ord r) => ConvexPolygon p r -> ConvexPolygon p r
-> (ConvexPolygon p r, LineSegment 2 p r, LineSegment 2 p r)
merge lp rp = (fromPoints $ r' ++ l', lt, ut)
where
lt@(ClosedLineSegment a b) = lowerTangent lp rp
ut@(ClosedLineSegment c d) = upperTangent lp rp
takeUntil p xs = let (xs',x:_) = break p xs in xs' ++ [x]
rightElems = F.toList . C.rightElements
r' = takeUntil (coreEq b) . rightElems . rotateTo' d $ rp^.outerBoundary
l' = takeUntil (coreEq c) . rightElems . rotateTo' a $ lp^.outerBoundary
rotateTo' :: Eq a => (a :+ b) -> CSeq (a :+ b) -> CSeq (a :+ b)
rotateTo' x = fromJust . C.findRotateTo (coreEq x)
coreEq :: Eq a => (a :+ b) -> (a :+ b) -> Bool
coreEq = (==) `on` (^.core)
-- | Compute the lower tangent of the two polgyons
--
-- pre: - polygons lp and rp have at least 1 vertex
-- - lp and rp are disjoint, and there is a vertical line separating
-- the two polygons.
-- - The vertices of the polygons are given in clockwise order
--
-- Running time: O(n+m), where n and m are the sizes of the two polygons respectively
lowerTangent :: (Num r, Ord r)
=> ConvexPolygon p r
-> ConvexPolygon p r
-> LineSegment 2 p r
lowerTangent (SimplePolygon l) (SimplePolygon r) = rotate xx yy zz zz''
where
xx = rightMost l
yy = leftMost r
zz = pred' yy
zz'' = succ' xx
rotate x y z z''
| focus z `isRightOf` (focus x, focus y) = rotate x z (pred' z) z''
-- rotate the right polygon CCW
| focus z'' `isRightOf` (focus x, focus y) = rotate z'' y z (succ' z'')
-- rotate the left polygon CW
| otherwise = ClosedLineSegment (focus x)
(focus y)
succ' :: CSeq a -> CSeq a
succ' = C.rotateR
pred' :: CSeq a -> CSeq a
pred' = C.rotateL
-- | Compute the upper tangent of the two polgyons
--
-- pre: - polygons lp and rp have at least 1 vertex
-- - lp and rp are disjoint, and there is a vertical line separating
-- the two polygons.
-- - The vertices of the polygons are given in clockwise order
--
-- Running time: O(n+m), where n and m are the sizes of the two polygons respectively
upperTangent :: (Num r, Ord r)
=> ConvexPolygon p r
-> ConvexPolygon p r
-> LineSegment 2 p r
upperTangent (SimplePolygon l) (SimplePolygon r) = rotate xx yy zz zz'
where
xx = rightMost l
yy = leftMost r
zz = succ' yy
zz' = pred' xx
rotate x y z z'
| focus z `isLeftOf` (focus x, focus y) = rotate x z (succ' z) z'
-- rotate the right polygon CW
| focus z' `isLeftOf` (focus x, focus y) = rotate z' y z (pred' z')
-- rotate the left polygon CCW
| otherwise = ClosedLineSegment (focus x)
(focus y)
isRightOf :: (Num r, Ord r)
=> Point 2 r :+ p -> (Point 2 r :+ p', Point 2 r :+ p'') -> Bool
a `isRightOf` (b,c) = ccw (b^.core) (c^.core) (a^.core) == CW
isLeftOf :: (Num r, Ord r)
=> Point 2 r :+ p -> (Point 2 r :+ p', Point 2 r :+ p'') -> Bool
a `isLeftOf` (b,c) = ccw (b^.core) (c^.core) (a^.core) == CCW
--------------------------------------------------------------------------------
-- | Rotate to the rightmost point
rightMost :: Ord r => CSeq (Point 2 r :+ p) -> CSeq (Point 2 r :+ p)
rightMost xs = let m = F.maximumBy (comparing (^.core.xCoord)) xs in rotateTo' m xs
-- | Rotate to the leftmost point
leftMost :: Ord r => CSeq (Point 2 r :+ p) -> CSeq (Point 2 r :+ p)
leftMost xs = let m = F.minimumBy (comparing (^.core.xCoord)) xs in rotateTo' m xs
-- -- | rotate right while p 'current' 'rightNeibhour' is true
-- rotateRWhile :: (a -> a -> Bool) -> C.CList a -> C.CList a
-- rotateRWhile p lst
-- | C.isEmpty lst = lst
-- | otherwise = go lst
-- where
-- go xs = let cur = focus xs
-- xs' = C.rotR xs
-- nxt = focus' xs'
-- in if p cur nxt then go xs' else xs