{-# LANGUAGE OverloadedStrings #-} -------------------------------------------------------------------------------- -- | -- Module : Data.Geometry.Polygon.Core -- Copyright : (C) Frank Staals -- License : see the LICENSE file -- Maintainer : Frank Staals -- -- A Polygon data type and some basic functions to interact with them. -- -------------------------------------------------------------------------------- module Data.Geometry.Polygon.Core ( PolygonType(..) , Polygon(..) , Vertices , _SimplePolygon, _MultiPolygon , SimplePolygon, MultiPolygon, SomePolygon -- * Construction , fromPoints , fromCircularVector , simpleFromPoints , simpleFromCircularVector , unsafeFromPoints , unsafeFromCircularVector , unsafeFromVector , toVector , toPoints , isSimple , size , polygonVertices, listEdges , outerBoundary, outerBoundaryVector , unsafeOuterBoundaryVector , outerBoundaryEdges , outerVertex, unsafeOuterVertex , outerBoundaryEdge , polygonHoles, polygonHoles' , holeList , area, signedArea , centroid , pickPoint , isTriangle , isCounterClockwise , toCounterClockWiseOrder, toCounterClockWiseOrder' , toClockwiseOrder, toClockwiseOrder' , reverseOuterBoundary , findDiagonal , withIncidentEdges, numberVertices -- * Testing for Reflex or Convex , isReflexVertex, isConvexVertex, isStrictlyConvexVertex , reflexVertices, convexVertices, strictlyConvexVertices -- * Specialized folds , maximumVertexBy , minimumVertexBy , findRotateTo , rotateLeft , rotateRight ) where import qualified Algorithms.Geometry.LineSegmentIntersection.BentleyOttmann as BO import Control.DeepSeq import Control.Lens (Getter, Lens', Prism', Traversal', lens, over, prism', to, toListOf, view, (%~), (&), (.~), (^.)) import Data.Aeson import Data.Bifoldable import Data.Bifunctor import Data.Bitraversable import Data.Ext import qualified Data.Foldable as F import Data.Geometry.Boundary import Data.Geometry.Box (IsBoxable (..), boundingBoxList') import Data.Geometry.Line import Data.Geometry.LineSegment import Data.Geometry.Point import Data.Geometry.Properties import Data.Geometry.Transformation import Data.Geometry.Triangle (Triangle (..), inTriangle) import Data.Geometry.Vector (Additive (zero, (^+^)), Affine ((.+^), (.-.)), (*^), (^*), (^/)) import qualified Data.List as List import qualified Data.List.NonEmpty as NonEmpty import Data.Maybe (catMaybes) import Data.Ord (comparing) import Data.Semigroup (sconcat) import Data.Semigroup.Foldable import Data.Util import Data.Vector (Vector) import qualified Data.Vector as V import Data.Vector.Circular (CircularVector) import qualified Data.Vector.Circular as CV import qualified Data.Vector.Circular.Util as CV -- import Data.RealNumber.Rational -------------------------------------------------------------------------------- {- $setup >>> import Data.RealNumber.Rational >>> import Data.Foldable >>> import Control.Lens.Extras >>> :{ -- import qualified Data.Vector.Circular as CV let simplePoly :: SimplePolygon () (RealNumber 10) simplePoly = fromPoints . map ext $ [ Point2 0 0 , Point2 10 0 , Point2 10 10 , Point2 5 15 , Point2 1 11 ] simpleTriangle :: SimplePolygon () (RealNumber 10) simpleTriangle = fromPoints . map ext $ [ Point2 0 0, Point2 2 0, Point2 1 1] multiPoly :: MultiPolygon () (RealNumber 10) multiPoly = MultiPolygon (fromPoints . map ext $ [Point2 (-1) (-1), Point2 3 (-1), Point2 2 2]) [simpleTriangle] :} -} -- | We distinguish between simple polygons (without holes) and polygons with holes. data PolygonType = Simple | Multi -- | Polygons are sequences of points and may or may not contain holes. -- -- Degenerate polygons (polygons with self-intersections or fewer than 3 points) -- are only possible if you use functions marked as unsafe. data Polygon (t :: PolygonType) p r where SimplePolygon :: Vertices (Point 2 r :+ p) -> SimplePolygon p r MultiPolygon :: SimplePolygon p r -> [SimplePolygon p r] -> MultiPolygon p r newtype Vertices a = Vertices (CircularVector a) deriving (Functor, Foldable, Foldable1, Traversable, NFData, Eq, Ord) -- | Prism to 'test' if we are a simple polygon -- -- >>> is _SimplePolygon simplePoly -- True _SimplePolygon :: Prism' (Polygon Simple p r) (Vertices (Point 2 r :+ p)) _SimplePolygon = prism' SimplePolygon (\(SimplePolygon vs) -> Just vs) -- | Prism to 'test' if we are a Multi polygon -- -- >>> is _MultiPolygon multiPoly -- True _MultiPolygon :: Prism' (Polygon Multi p r) (Polygon Simple p r, [Polygon Simple p r]) _MultiPolygon = prism' (uncurry MultiPolygon) (\(MultiPolygon vs hs) -> Just (vs,hs)) instance Functor (Polygon t p) where fmap = bimap id instance Bifunctor (Polygon t) where bimap = bimapDefault instance Bifoldable (Polygon t) where bifoldMap = bifoldMapDefault instance Bitraversable (Polygon t) where bitraverse f g p = case p of SimplePolygon vs -> SimplePolygon <$> bitraverseVertices f g vs MultiPolygon vs hs -> MultiPolygon <$> bitraverse f g vs <*> traverse (bitraverse f g) hs instance (NFData p, NFData r) => NFData (Polygon t p r) where rnf (SimplePolygon vs) = rnf vs rnf (MultiPolygon vs hs) = rnf (vs,hs) bitraverseVertices :: (Applicative f, Traversable t) => (p -> f q) -> (r -> f s) -> t (Point 2 r :+ p) -> f (t (Point 2 s :+ q)) bitraverseVertices f g = traverse (bitraverse (traverse g) f) -- | Polygon without holes. type SimplePolygon = Polygon Simple -- | Polygon with zero or more holes. type MultiPolygon = Polygon Multi -- | Either a simple or multipolygon type SomePolygon p r = Either (Polygon Simple p r) (Polygon Multi p r) type instance Dimension (SomePolygon p r) = 2 type instance NumType (SomePolygon p r) = r -- | Polygons are per definition 2 dimensional type instance Dimension (Polygon t p r) = 2 type instance NumType (Polygon t p r) = r instance (Show p, Show r) => Show (Polygon t p r) where show (SimplePolygon vs) = "SimplePolygon " <> show (F.toList vs) show (MultiPolygon vs hs) = "MultiPolygon (" <> show vs <> ") (" <> show hs <> ")" instance (Read p, Read r) => Read (SimplePolygon p r) where readsPrec d = readParen (d > app_prec) $ \r -> [ (unsafeFromPoints vs, t) | ("SimplePolygon", s) <- lex r, (vs, t) <- reads s ] where app_prec = 10 instance (Read p, Read r) => Read (MultiPolygon p r) where readsPrec d = readParen (d > app_prec) $ \r -> [ (MultiPolygon vs hs, t') | ("MultiPolygon", s) <- lex r , (vs, t) <- reads s , (hs, t') <- reads t ] where app_prec = 10 -- instance (Read p, Read r) => Show (Polygon t p r) where -- show (SimplePolygon vs) = "SimplePolygon (" <> show vs <> ")" -- show (MultiPolygon vs hs) = "MultiPolygon (" <> show vs <> ") (" <> show hs <> ")" instance (Eq p, Eq r) => Eq (Polygon t p r) where (SimplePolygon vs) == (SimplePolygon vs') = vs == vs' (MultiPolygon vs hs) == (MultiPolygon vs' hs') = vs == vs' && hs == hs' instance PointFunctor (Polygon t p) where pmap f (SimplePolygon vs) = SimplePolygon (fmap (first f) vs) pmap f (MultiPolygon vs hs) = MultiPolygon (pmap f vs) (map (pmap f) hs) instance Fractional r => IsTransformable (Polygon t p r) where transformBy = transformPointFunctor instance IsBoxable (Polygon t p r) where boundingBox = boundingBoxList' . toListOf (outerBoundaryVector.traverse.core) instance (ToJSON r, ToJSON p) => ToJSON (Polygon t p r) where toJSON = \case (SimplePolygon vs) -> object [ "tag" .= ("SimplePolygon" :: String) , "vertices" .= F.toList vs ] (MultiPolygon vs hs) -> object [ "tag" .= ("MultiPolygon" :: String) , "outerBoundary" .= getVertices vs , "holes" .= map getVertices hs ] where getVertices = view (outerBoundaryVector.to F.toList) instance (FromJSON r, Eq r, Num r, FromJSON p) => FromJSON (Polygon Simple p r) where parseJSON = withObject "Polygon" $ \o -> o .: "tag" >>= \case "SimplePolygon" -> pSimple o (_ :: String) -> fail "Not a SimplePolygon" where pSimple o = fromPoints <$> o .: "vertices" instance (FromJSON r, Eq r, Num r, FromJSON p) => FromJSON (Polygon Multi p r) where parseJSON = withObject "Polygon" $ \o -> o .: "tag" >>= \case "MultiPolygon" -> pMulti o (_ :: String) -> fail "Not a MultiPolygon" where pMulti o = (\vs hs -> MultiPolygon (fromPoints vs) (map fromPoints hs)) <$> o .: "outerBoundary" <*> o .: "holes" -- * Functions on Polygons -- | Getter access to the outer boundary vector of a polygon. -- -- >>> toList (simpleTriangle ^. outerBoundaryVector) -- [Point2 0 0 :+ (),Point2 2 0 :+ (),Point2 1 1 :+ ()] outerBoundaryVector :: forall t p r. Getter (Polygon t p r) (CircularVector (Point 2 r :+ p)) outerBoundaryVector = to g where g :: Polygon t p r -> CircularVector (Point 2 r :+ p) g (SimplePolygon (Vertices vs)) = vs g (MultiPolygon (SimplePolygon (Vertices vs)) _) = vs -- | Unsafe lens access to the outer boundary vector of a polygon. -- -- >>> toList (simpleTriangle ^. unsafeOuterBoundaryVector) -- [Point2 0 0 :+ (),Point2 2 0 :+ (),Point2 1 1 :+ ()] -- -- >>> simpleTriangle & unsafeOuterBoundaryVector .~ CV.singleton (Point2 0 0 :+ ()) -- SimplePolygon [Point2 0 0 :+ ()] unsafeOuterBoundaryVector :: forall t p r. Lens' (Polygon t p r) (CircularVector (Point 2 r :+ p)) unsafeOuterBoundaryVector = lens g s where g :: Polygon t p r -> CircularVector (Point 2 r :+ p) g (SimplePolygon (Vertices vs)) = vs g (MultiPolygon (SimplePolygon (Vertices vs)) _) = vs s :: Polygon t p r -> CircularVector (Point 2 r :+ p) -> Polygon t p r s SimplePolygon{} vs = SimplePolygon (Vertices vs) s (MultiPolygon _ hs) vs = MultiPolygon (SimplePolygon (Vertices vs)) hs -- | \( O(1) \) Lens access to the outer boundary of a polygon. outerBoundary :: forall t p r. Lens' (Polygon t p r) (SimplePolygon p r) outerBoundary = lens g s where g :: Polygon t p r -> SimplePolygon p r g poly@SimplePolygon{} = poly g (MultiPolygon simple _) = simple s :: Polygon t p r -> SimplePolygon p r -> Polygon t p r s SimplePolygon{} simple = simple s (MultiPolygon _ hs) simple = MultiPolygon simple hs -- | Lens access for polygon holes. -- -- >>> multiPoly ^. polygonHoles -- [SimplePolygon [Point2 0 0 :+ (),Point2 2 0 :+ (),Point2 1 1 :+ ()]] polygonHoles :: forall p r. Lens' (Polygon Multi p r) [Polygon Simple p r] polygonHoles = lens g s where g :: Polygon Multi p r -> [Polygon Simple p r] g (MultiPolygon _ hs) = hs s :: Polygon Multi p r -> [Polygon Simple p r] -> Polygon Multi p r s (MultiPolygon vs _) = MultiPolygon vs {- HLINT ignore polygonHoles' -} -- | \( O(1) \). Traversal lens for polygon holes. Does nothing for simple polygons. polygonHoles' :: Traversal' (Polygon t p r) [Polygon Simple p r] polygonHoles' = \f -> \case p@SimplePolygon{} -> pure p MultiPolygon vs hs -> MultiPolygon vs <$> f hs -- | /O(1)/ Access the i^th vertex on the outer boundary. Indices are modulo \(n\). -- -- >>> simplePoly ^. outerVertex 0 -- Point2 0 0 :+ () outerVertex :: Int -> Getter (Polygon t p r) (Point 2 r :+ p) outerVertex i = outerBoundaryVector . CV.item i -- | \( O(1) \) read and \( O(n) \) write. Access the i^th vertex on the outer boundary -- -- >>> simplePoly ^. unsafeOuterVertex 0 -- Point2 0 0 :+ () -- >>> simplePoly & unsafeOuterVertex 0 .~ (Point2 10 10 :+ ()) -- SimplePolygon [Point2 10 10 :+ (),Point2 10 0 :+ (),Point2 10 10 :+ (),Point2 5 15 :+ (),Point2 1 11 :+ ()] unsafeOuterVertex :: Int -> Lens' (Polygon t p r) (Point 2 r :+ p) unsafeOuterVertex i = unsafeOuterBoundaryVector . CV.item i -- | \( O(1) \) Get the n^th edge along the outer boundary of the polygon. The edge is half open. outerBoundaryEdge :: Int -> Polygon t p r -> LineSegment 2 p r outerBoundaryEdge i p = let u = p^.outerVertex i v = p^.outerVertex (i+1) in LineSegment (Closed u) (Open v) -- | Get all holes in a polygon holeList :: Polygon t p r -> [Polygon Simple p r] holeList SimplePolygon{} = [] holeList (MultiPolygon _ hs) = hs -- | \( O(1) \) Vertex count. Includes the vertices of holes. size :: Polygon t p r -> Int size (SimplePolygon (Vertices cv)) = F.length cv size (MultiPolygon b hs) = sum (map size (b:hs)) -- | \( O(n) \) The vertices in the polygon. No guarantees are given on the order in which -- they appear! polygonVertices :: Polygon t p r -> NonEmpty.NonEmpty (Point 2 r :+ p) polygonVertices p@SimplePolygon{} = toNonEmpty $ p^.outerBoundaryVector polygonVertices (MultiPolygon vs hs) = sconcat $ toNonEmpty (polygonVertices vs) NonEmpty.:| map polygonVertices hs -- FIXME: Get rid of 'Fractional r' constraint. -- | \( O(n \log n) \) Check if a polygon has any holes, duplicate points, or -- self-intersections. isSimple :: (Ord r, Fractional r) => Polygon p t r -> Bool isSimple p@SimplePolygon{} = null . BO.interiorIntersections . map ext $ listEdges p isSimple (MultiPolygon b []) = isSimple b isSimple MultiPolygon{} = False requireThree :: String -> [a] -> [a] requireThree _ lst@(_:_:_:_) = lst requireThree label _ = error $ "Data.Geometry.Polygon." ++ label ++ ": Polygons must have at least three points." -- | \( O(n) \) Creates a polygon from the given list of vertices. -- -- The points are placed in CCW order if they are not already. Overlapping -- edges and repeated vertices are allowed. -- fromPoints :: forall p r. (Eq r, Num r) => [Point 2 r :+ p] -> SimplePolygon p r fromPoints = fromCircularVector . CV.unsafeFromList . requireThree "fromPoints" -- | \( O(n) \) Creates a polygon from the given vector of vertices. -- -- The points are placed in CCW order if they are not already. Overlapping -- edges and repeated vertices are allowed. -- fromCircularVector :: forall p r. (Eq r, Num r) => CircularVector (Point 2 r :+ p) -> SimplePolygon p r fromCircularVector = toCounterClockWiseOrder . unsafeFromCircularVector -- | \( O(n \log n) \) Creates a simple polygon from the given list of vertices. -- -- The points are placed in CCW order if they are not already. Overlapping -- edges and repeated vertices are /not/ allowed and will trigger an exception. -- simpleFromPoints :: forall p r. (Ord r, Fractional r) => [Point 2 r :+ p] -> SimplePolygon p r simpleFromPoints = simpleFromCircularVector . CV.unsafeFromList . requireThree "simpleFromPoints" -- | \( O(n \log n) \) Creates a simple polygon from the given vector of vertices. -- -- The points are placed in CCW order if they are not already. Overlapping -- edges and repeated vertices are /not/ allowed and will trigger an exception. -- simpleFromCircularVector :: forall p r. (Ord r, Fractional r) => CircularVector (Point 2 r :+ p) -> SimplePolygon p r simpleFromCircularVector v = let p = fromCircularVector v hasInteriorIntersections = not . null . BO.interiorIntersections . map ext in if hasInteriorIntersections (listEdges p) then error "Data.Geometry.Polygon.simpleFromCircularVector: \ \Found self-intersections or repeated vertices." else p -- | \( O(n) \) Creates a simple polygon from the given list of vertices. -- -- pre: the input list constains no repeated vertices. unsafeFromPoints :: [Point 2 r :+ p] -> SimplePolygon p r unsafeFromPoints = unsafeFromCircularVector . CV.unsafeFromList -- | \( O(1) \) Creates a simple polygon from the given vector of vertices. -- -- pre: the input list constains no repeated vertices. unsafeFromCircularVector :: CircularVector (Point 2 r :+ p) -> SimplePolygon p r unsafeFromCircularVector = SimplePolygon . Vertices -- | \( O(1) \) Creates a simple polygon from the given vector of vertices. -- -- pre: the input list constains no repeated vertices. unsafeFromVector :: Vector (Point 2 r :+ p) -> SimplePolygon p r unsafeFromVector = unsafeFromCircularVector . CV.unsafeFromVector -- -- | Polygon points, from left to right. -- toList :: Polygon t p r -> [Point 2 r :+ p] -- toList (SimplePolygon c) = F.toList c -- toList (MultiPolygon s hs) = toList s ++ concatMap toList hs -- | \( O(n) \) -- Polygon points, from left to right. toVector :: Polygon t p r -> Vector (Point 2 r :+ p) toVector p@SimplePolygon{} = CV.toVector $ p^.outerBoundaryVector toVector (MultiPolygon s hs) = foldr (<>) (toVector s) (map toVector hs) -- | \( O(n) \) -- Polygon points, from left to right. toPoints :: Polygon t p r -> [Point 2 r :+ p] toPoints = V.toList . toVector -- | \( O(n) \) The edges along the outer boundary of the polygon. The edges are half open. outerBoundaryEdges :: Polygon t p r -> CircularVector (LineSegment 2 p r) outerBoundaryEdges = toEdges . (^.outerBoundaryVector) -- | \( O(n) \) Lists all edges. The edges on the outer boundary are given before the ones -- on the holes. However, no other guarantees are given on the order. listEdges :: Polygon t p r -> [LineSegment 2 p r] listEdges pg = let f = F.toList . outerBoundaryEdges in f pg <> concatMap f (holeList pg) -- | Pairs every vertex with its incident edges. The first one is its -- predecessor edge, the second one its successor edge (in terms of -- the ordering along the boundary). -- -- -- >>> mapM_ print . polygonVertices $ withIncidentEdges simplePoly -- Point2 0 0 :+ V2 (ClosedLineSegment (Point2 1 11 :+ ()) (Point2 0 0 :+ ())) (ClosedLineSegment (Point2 0 0 :+ ()) (Point2 10 0 :+ ())) -- Point2 10 0 :+ V2 (ClosedLineSegment (Point2 0 0 :+ ()) (Point2 10 0 :+ ())) (ClosedLineSegment (Point2 10 0 :+ ()) (Point2 10 10 :+ ())) -- Point2 10 10 :+ V2 (ClosedLineSegment (Point2 10 0 :+ ()) (Point2 10 10 :+ ())) (ClosedLineSegment (Point2 10 10 :+ ()) (Point2 5 15 :+ ())) -- Point2 5 15 :+ V2 (ClosedLineSegment (Point2 10 10 :+ ()) (Point2 5 15 :+ ())) (ClosedLineSegment (Point2 5 15 :+ ()) (Point2 1 11 :+ ())) -- Point2 1 11 :+ V2 (ClosedLineSegment (Point2 5 15 :+ ()) (Point2 1 11 :+ ())) (ClosedLineSegment (Point2 1 11 :+ ()) (Point2 0 0 :+ ())) withIncidentEdges :: Polygon t p r -> Polygon t (Two (LineSegment 2 p r)) r withIncidentEdges poly@SimplePolygon{} = unsafeFromCircularVector $ CV.zipWith3 f (CV.rotateLeft 1 vs) vs (CV.rotateRight 1 vs) where vs = poly ^. outerBoundaryVector f p c n = c&extra .~ Two (ClosedLineSegment p c) (ClosedLineSegment c n) withIncidentEdges (MultiPolygon vs hs) = MultiPolygon vs' hs' where vs' = withIncidentEdges vs hs' = map withIncidentEdges hs -- -- | Gets the i^th edge on the outer boundary of the polygon, that is the edge ---- with vertices i and i+1 with respect to the current focus. All indices -- -- modulo n. -- -- -- FIXME: Test that \poly -> fromEdges (toEdges poly) == poly -- | Given the vertices of the polygon. Produce a list of edges. The edges are -- half-open. toEdges :: CircularVector (Point 2 r :+ p) -> CircularVector (LineSegment 2 p r) toEdges vs = CV.zipWith (\p q -> LineSegment (Closed p) (Open q)) vs (CV.rotateRight 1 vs) -- | Compute the area of a polygon area :: Fractional r => Polygon t p r -> r area poly@SimplePolygon{} = abs $ signedArea poly area (MultiPolygon vs hs) = area vs - sum [area h | h <- hs] -- | Compute the signed area of a simple polygon. The the vertices are in -- clockwise order, the signed area will be negative, if the verices are given -- in counter clockwise order, the area will be positive. signedArea :: Fractional r => SimplePolygon p r -> r signedArea poly = signedArea2X poly / 2 -- | Compute the signed area times 2 of a simple polygon. The the vertices are in -- clockwise order, the signed area will be negative, if the verices are given -- in counter clockwise order, the area will be positive. signedArea2X :: Num r => SimplePolygon p r -> r signedArea2X poly = x where x = sum [ p^.core.xCoord * q^.core.yCoord - q^.core.xCoord * p^.core.yCoord | LineSegment' p q <- F.toList $ outerBoundaryEdges poly ] -- | Compute the centroid of a simple polygon. centroid :: Fractional r => SimplePolygon p r -> Point 2 r centroid poly = Point $ sum' xs ^/ (6 * signedArea poly) where xs = [ (toVec p ^+^ toVec q) ^* (p^.xCoord * q^.yCoord - q^.xCoord * p^.yCoord) | LineSegment' (p :+ _) (q :+ _) <- F.toList $ outerBoundaryEdges poly ] sum' = F.foldl' (^+^) zero -- | \( O(n) \) Pick a point that is inside the polygon. -- -- (note: if the polygon is degenerate; i.e. has <3 vertices, we report a -- vertex of the polygon instead.) -- -- pre: the polygon is given in CCW order pickPoint :: (Ord r, Fractional r) => Polygon p t r -> Point 2 r pickPoint pg | isTriangle pg = centroid $ pg^.outerBoundary | otherwise = let LineSegment' (p :+ _) (q :+ _) = findDiagonal pg in p .+^ (0.5 *^ (q .-. p)) -- | \( O(1) \) Test if the polygon is a triangle isTriangle :: Polygon p t r -> Bool isTriangle = \case p@SimplePolygon{} -> F.length (p^.outerBoundaryVector) == 3 MultiPolygon vs [] -> isTriangle vs MultiPolygon _ _ -> False -- | \( O(n) \) Find a diagonal of the polygon. -- -- pre: the polygon is given in CCW order findDiagonal :: (Ord r, Fractional r) => Polygon t p r -> LineSegment 2 p r findDiagonal pg = List.head . catMaybes . F.toList $ diags -- note that a diagonal is guaranteed to exist, so the usage of head is safe. where vs = pg^.outerBoundaryVector diags = CV.zipWith3 f (CV.rotateLeft 1 vs) vs (CV.rotateRight 1 vs) f u v w = case ccw (u^.core) (v^.core) (w^.core) of CCW -> Just $ findDiag u v w -- v is a convex vertex, so find a diagonal -- (either uw) or from v to a point inside the -- triangle CW -> Nothing -- v is a reflex vertex CoLinear -> Nothing -- colinear vertex!? -- we test if uw is a diagonal by figuring out if there is a vertex -- strictly inside the triangle t. If there is no such vertex then uw must -- be a diagonal (i.e. uw intersects the polygon boundary iff there is a -- vtx inside t). If there are vertices inside the triangle, we find the -- one z furthest from the line(segment) uw. It then follows that vz is a -- diagonal. Indeed this is pretty much the argument used to prove that any -- polygon can be triangulated. See BKOS Chapter 3 for details. findDiag u v w = let t = Triangle u v w uw = ClosedLineSegment u w in maybe uw (ClosedLineSegment v) . safeMaximumOn (distTo $ supportingLine uw) . filter (\(z :+ _) -> z `inTriangle` t == Inside) . F.toList . polygonVertices $ pg distTo l (z :+ _) = sqDistanceTo z l safeMaximumOn :: Ord b => (a -> b) -> [a] -> Maybe a safeMaximumOn f = \case [] -> Nothing xs -> Just $ List.maximumBy (comparing f) xs -- | \( O(n) \) Test if the outer boundary of the polygon is in clockwise or counter -- clockwise order. isCounterClockwise :: (Eq r, Num r) => Polygon t p r -> Bool isCounterClockwise = (\x -> x == abs x) . signedArea2X . view outerBoundary -- | \( O(n) \) Make sure that every edge has the polygon's interior on its -- right, by orienting the outer boundary into clockwise order, and -- the inner borders (i.e. any holes, if they exist) into -- counter-clockwise order. toClockwiseOrder :: (Eq r, Num r) => Polygon t p r -> Polygon t p r toClockwiseOrder p = toClockwiseOrder' p & polygonHoles'.traverse %~ toCounterClockWiseOrder' -- | \( O(n) \) Orient the outer boundary into clockwise order. Leaves any holes -- as they are. -- toClockwiseOrder' :: (Eq r, Num r) => Polygon t p r -> Polygon t p r toClockwiseOrder' pg | isCounterClockwise pg = reverseOuterBoundary pg | otherwise = pg -- | \( O(n) \) Make sure that every edge has the polygon's interior on its left, -- by orienting the outer boundary into counter-clockwise order, and -- the inner borders (i.e. any holes, if they exist) into clockwise order. toCounterClockWiseOrder :: (Eq r, Num r) => Polygon t p r -> Polygon t p r toCounterClockWiseOrder p = toCounterClockWiseOrder' p & polygonHoles'.traverse %~ toClockwiseOrder' -- | \( O(n) \) Orient the outer boundary into counter-clockwise order. Leaves -- any holes as they are. toCounterClockWiseOrder' :: (Eq r, Num r) => Polygon t p r -> Polygon t p r toCounterClockWiseOrder' p | not $ isCounterClockwise p = reverseOuterBoundary p | otherwise = p -- FIXME: Delete this function. -- | Reorient the outer boundary from clockwise order to counter-clockwise order or -- from counter-clockwise order to clockwise order. Leaves -- any holes as they are. -- reverseOuterBoundary :: Polygon t p r -> Polygon t p r reverseOuterBoundary p = p&unsafeOuterBoundaryVector %~ CV.reverse -- | assigns unique integer numbers to all vertices. Numbers start from 0, and -- are increasing along the outer boundary. The vertices of holes -- will be numbered last, in the same order. -- -- >>> numberVertices simplePoly -- SimplePolygon [Point2 0 0 :+ SP 0 (),Point2 10 0 :+ SP 1 (),Point2 10 10 :+ SP 2 (),Point2 5 15 :+ SP 3 (),Point2 1 11 :+ SP 4 ()] numberVertices :: Polygon t p r -> Polygon t (SP Int p) r numberVertices = snd . bimapAccumL (\a p -> (a+1,SP a p)) (,) 0 -- TODO: Make sure that this does not have the same issues as foldl vs foldl' -------------------------------------------------------------------------------- -- Specialized folds -- maximum and minimum probably aren't useful. Disabled for now. Lemmih, 2020-12-26. -- | \( O(n) \) Yield the maximum point of the polygon. Points are compared first by x-coordinate -- and then by y-coordinate. The maximum point will therefore be the right-most point in -- the polygon (and top-most if multiple points share the largest x-coordinate). -- -- Hole vertices are ignored since they cannot be the maximum. _maximum :: Ord r => Polygon t p r -> Point 2 r :+ p _maximum = F.maximumBy (comparing _core) . view outerBoundaryVector -- | \( O(n) \) Yield the maximum point of a polygon according to the given comparison function. maximumVertexBy :: (Point 2 r :+ p -> Point 2 r :+ p -> Ordering) -> Polygon t p r -> Point 2 r :+ p maximumVertexBy fn (SimplePolygon vs) = F.maximumBy fn vs maximumVertexBy fn (MultiPolygon b hs) = F.maximumBy fn $ map (maximumVertexBy fn) (b:hs) -- | \( O(n) \) Yield the maximum point of the polygon. Points are compared first by x-coordinate -- and then by y-coordinate. The minimum point will therefore be the left-most point in -- the polygon (and bottom-most if multiple points share the smallest x-coordinate). -- -- Hole vertices are ignored since they cannot be the minimum. _minimum :: Ord r => Polygon t p r -> Point 2 r :+ p _minimum = F.minimumBy (comparing _core) . view outerBoundaryVector -- | \( O(n) \) Yield the maximum point of a polygon according to the given comparison function. minimumVertexBy :: (Point 2 r :+ p -> Point 2 r :+ p -> Ordering) -> Polygon t p r -> Point 2 r :+ p minimumVertexBy fn (SimplePolygon vs) = F.minimumBy fn vs minimumVertexBy fn (MultiPolygon b hs) = F.minimumBy fn $ map (minimumVertexBy fn) (b:hs) -- | Rotate to the first point that matches the given condition. -- -- >>> toVector <$> findRotateTo (== (Point2 1 0 :+ ())) (unsafeFromPoints [Point2 0 0 :+ (), Point2 1 0 :+ (), Point2 1 1 :+ ()]) -- Just [Point2 1 0 :+ (),Point2 1 1 :+ (),Point2 0 0 :+ ()] -- >>> findRotateTo (== (Point2 7 0 :+ ())) $ unsafeFromPoints [Point2 0 0 :+ (), Point2 1 0 :+ (), Point2 1 1 :+ ()] -- Nothing findRotateTo :: (Point 2 r :+ p -> Bool) -> SimplePolygon p r -> Maybe (SimplePolygon p r) findRotateTo fn = fmap unsafeFromCircularVector . CV.findRotateTo fn . view outerBoundaryVector -------------------------------------------------------------------------------- -- Rotation -- | \( O(1) \) Rotate the polygon to the left by n number of points. rotateLeft :: Int -> SimplePolygon p r -> SimplePolygon p r rotateLeft n = over unsafeOuterBoundaryVector (CV.rotateLeft n) -- | \( O(1) \) Rotate the polygon to the right by n number of points. rotateRight :: Int -> SimplePolygon p r -> SimplePolygon p r rotateRight n = over unsafeOuterBoundaryVector (CV.rotateRight n) -------------------------------------------------------------------------------- -- Testing for reflex or convex -- | Test if a given vertex is a reflex vertex. -- -- \(O(1)\) isReflexVertex :: (Ord r, Num r) => Int -> Polygon Simple p r -> Bool isReflexVertex i pg = ccw' u v w == CW where u = pg^.outerVertex (i-1) v = pg^.outerVertex i w = pg^.outerVertex (i+1) -- | Test if a given vertex is a convex vertex (i.e. not a reflex vertex). -- -- \(O(1)\) isConvexVertex :: (Ord r, Num r) => Int -> Polygon Simple p r -> Bool isConvexVertex i = not . isReflexVertex i -- | Test if a given vertex is a strictly convex vertex. -- -- \(O(1)\) isStrictlyConvexVertex :: (Ord r, Num r) => Int -> Polygon t p r -> Bool isStrictlyConvexVertex i pg = ccw' u v w == CCW where u = pg^.outerVertex (i-1) v = pg^.outerVertex i w = pg^.outerVertex (i+1) -- | Computes all reflex vertices of the polygon. -- -- \(O(n)\) reflexVertices :: (Ord r, Num r) => Polygon t p r -> [Int :+ (Point 2 r :+ p)] reflexVertices p@(SimplePolygon _) = reflexVertices' p reflexVertices (numberVertices -> MultiPolygon vs hs) = map (\(_ :+ (p :+ SP i e)) -> i :+ (p :+ e)) $ reflexVertices' vs <> concatMap strictlyConvexVertices' hs -- | Computes all convex (i.e. non-reflex) vertices of the polygon. -- -- \(O(n)\) convexVertices :: (Ord r, Num r) => Polygon t p r -> [Int :+ (Point 2 r :+ p)] convexVertices = \case p@(SimplePolygon _) -> convexVertices' p (numberVertices -> MultiPolygon vs hs) -> map (\(_ :+ (p :+ SP i e)) -> i :+ (p :+ e)) $ convexVertices' vs <> concatMap reflexVertices' hs -- | Computes all strictly convex vertices of the polygon. -- -- \(O(n)\) strictlyConvexVertices :: (Ord r, Num r) => Polygon t p r -> [Int :+ (Point 2 r :+ p)] strictlyConvexVertices = \case p@(SimplePolygon _) -> convexVertices' p (numberVertices -> MultiPolygon vs hs) -> map (\(_ :+ (p :+ SP i e)) -> i :+ (p :+ e)) $ strictlyConvexVertices' vs <> concatMap reflexVertices' hs ---------------------------------------- -- | Return (the indices of) all reflex vertices, in increasing order -- along the boundary. -- -- \(O(n)\) reflexVertices' :: (Ord r, Num r) => SimplePolygon p r -> [Int :+ (Point 2 r :+ p)] reflexVertices' = filterReflexConvexWorker asReflex where asReflex u v w | ccw' (u^.extra) (v^.extra) (w^.extra) == CW = Just v | otherwise = Nothing -- | Return (the indices of) all strictly convex vertices, in -- increasing order along the boundary. -- -- \(O(n)\) strictlyConvexVertices' :: (Ord r, Num r) => SimplePolygon p r -> [Int :+ (Point 2 r :+ p)] strictlyConvexVertices' = filterReflexConvexWorker asStrictlyConvex where asStrictlyConvex u v w | ccw' (u^.extra) (v^.extra) (w^.extra) == CCW = Just v | otherwise = Nothing -- | Return (the indices of) all convex (= non-reflex) vertices, in increasing order -- along the boundary. -- -- \(O(n)\) convexVertices' :: (Ord r, Num r) => SimplePolygon p r -> [Int :+ (Point 2 r :+ p)] convexVertices' = filterReflexConvexWorker asConvex where asConvex u v w | ccw' (u^.extra) (v^.extra) (w^.extra) /= CW = Just v | otherwise = Nothing -- | Helper function to implement convexVertices, reflexVertices, and -- strictlyConvexVertices filterReflexConvexWorker :: (Ord r, Num r) => ( Int :+ (Point 2 r :+ p) -> Int :+ (Point 2 r :+ p) -> Int :+ (Point 2 r :+ p) -> Maybe (Int :+ (Point 2 r :+ p)) ) -> SimplePolygon p r -> [Int :+ (Point 2 r :+ p)] filterReflexConvexWorker g pg = catMaybes $ zip3RWith g (CV.rotateLeft 1 vs) vs (CV.rotateRight 1 vs) where vs = CV.withIndicesRight $ pg^.outerBoundaryVector zip3RWith f us' vs' ws' = zipWith3 f (F.toList us') (F.toList vs') (F.toList ws')