{-# 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: 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