hgeometry 0.12.0.4 → 0.13

raw patch · 39 files changed

+2359/−427 lines, 39 filesdep ~hgeometry-combinatorialbinary-addedPVP ok

version bump matches the API change (PVP)

Dependency ranges changed: hgeometry-combinatorial

API changes (from Hackage documentation)

Files

• benchmark/Algorithms/Geometry/PolygonTriangulation/Bench.hs +1/−1
• changelog +18/−0
• changelog.org +18/−0
• docs/Data/Geometry/PlanarSubdivision/mySubdiv.jpg binary
• docs/Data/PlaneGraph/planegraph.png binary
• hgeometry.cabal +10/−8
• src/Algorithms/Geometry/DelaunayTriangulation/Types.hs +3/−3
• src/Algorithms/Geometry/PolyLineSimplification/ImaiIri.hs +13/−13
• src/Algorithms/Geometry/PolygonTriangulation.hs +14/−0
• src/Algorithms/Geometry/PolygonTriangulation/Triangulate.hs +3/−5
• src/Data/Geometry.hs +1/−1
• src/Data/Geometry/Ball.hs +38/−10
• src/Data/Geometry/BezierSpline.hs +508/−80
• src/Data/Geometry/Box/Corners.hs +3/−1
• src/Data/Geometry/Box/Internal.hs +6/−1
• src/Data/Geometry/HalfLine.hs +11/−1
• src/Data/Geometry/HalfSpace.hs +4/−4
• src/Data/Geometry/HyperPlane.hs +5/−1
• src/Data/Geometry/Interval.hs +1/−0
• src/Data/Geometry/Line.hs +11/−8
• src/Data/Geometry/Line/Internal.hs +17/−0
• src/Data/Geometry/LineSegment/Internal.hs +15/−4
• src/Data/Geometry/PlanarSubdivision.hs +24/−0
• src/Data/Geometry/PlanarSubdivision/Basic.hs +259/−55
• src/Data/Geometry/PlanarSubdivision/Dynamic.hs +530/−0
• src/Data/Geometry/PlanarSubdivision/Raw.hs +8/−0
• src/Data/Geometry/PlanarSubdivision/TreeRep.hs +110/−0
• src/Data/Geometry/Polygon.hs +4/−1
• src/Data/Geometry/Polygon/Bezier.hs +6/−3
• src/Data/Geometry/QuadTree/Cell.hs +3/−0
• src/Data/Geometry/Slab.hs +14/−0
• src/Data/Geometry/SubLine.hs +5/−0
• src/Data/Geometry/Transformation.hs +42/−175
• src/Data/Geometry/Transformation/Internal.hs +219/−0
• src/Data/Geometry/Triangle.hs +4/−0
• src/Data/PlaneGraph.hs +84/−5
• src/Data/PlaneGraph/AdjRep.hs +1/−1
• src/Data/PlaneGraph/Core.hs +268/−41
• src/Data/PlaneGraph/IO.hs +78/−5

@@ -9,7 +9,7 @@ import           Data.Ext import           Test.Tasty.Bench import qualified Data.Foldable as F-import           Data.Geometry.Ipe+import           Ipe import           Data.Geometry.LineSegment import           Data.Geometry.Polygon import           Data.Geometry.PlanarSubdivision

@@ -2,6 +2,24 @@  * Changelog +** 0.13++- Implementation of Logaritmic Method, wich allows us to transform a+  static data structure into an insertion only data structure+- Moved 'intersects' from the HasIntersectionWith class into a new+  class IsIntersectableWith. This allows separate (weaker) constraints+  for checking *if* geometries intersect rather than computing exact+  intersections.+- New BezierSpline features.+- "Zoom to fit" transformation.+- Many fixes related to PlaneGraph/PlanarSubdivison; i.e. bugs in+  which order the vertices/darts where reported when traversing a+  face. The polygon representing the outer boundary now is some area+  inside a bounding polygon.+- Fixed a bug in the DelaunayTriangulation.+- Preliminary implementations for updating planar subdivisions+  (e.g. subdividing edges).+ ** 0.12  - New website: https://hgeometry.org/

@@ -2,6 +2,24 @@  * Changelog +** 0.13++- Implementation of Logaritmic Method, wich allows us to transform a+  static data structure into an insertion only data structure+- Moved 'intersects' from the HasIntersectionWith class into a new+  class IsIntersectableWith. This allows separate (weaker) constraints+  for checking *if* geometries intersect rather than computing exact+  intersections.+- New BezierSpline features.+- "Zoom to fit" transformation.+- Many fixes related to PlaneGraph/PlanarSubdivison; i.e. bugs in+  which order the vertices/darts where reported when traversing a+  face. The polygon representing the outer boundary now is some area+  inside a bounding polygon.+- Fixed a bug in the DelaunayTriangulation.+- Preliminary implementations for updating planar subdivisions+  (e.g. subdividing edges).+ ** 0.12  - New website: https://hgeometry.org/

binary file changed (absent → 378223 bytes)

binary file changed (absent → 323390 bytes)

@@ -1,5 +1,6 @@+cabal-version:       2.4 name:                hgeometry-version:             0.12.0.4+version:             0.13 synopsis:            Geometric Algorithms, Data structures, and Data types. description:   HGeometry provides some basic geometry types, and geometric algorithms and@@ -8,7 +9,7 @@   asymptotic running time guarantees. Note that HGeometry is still highly experimental, don't be surprised to find bugs.  homepage:            https://fstaals.net/software/hgeometry-license:             BSD3+license:             BSD-3-Clause license-file:        LICENSE author:              Frank Staals maintainer:          frank@fstaals.net@@ -23,10 +24,8 @@                      changelog                      changelog.org -Extra-doc-files:     docs/Data/PlaneGraph/small.png-                     -- docs/**/*.png--cabal-version:       2.0+Extra-doc-files:     docs/**/*.png+                     docs/**/*.jpg source-repository head   type:     git   location: https://github.com/noinia/hgeometry@@ -99,6 +98,8 @@                     Data.Geometry.PlanarSubdivision.Raw                     Data.Geometry.PlanarSubdivision.Basic                     Data.Geometry.PlanarSubdivision.Merge+                    Data.Geometry.PlanarSubdivision.Dynamic+                    Data.Geometry.PlanarSubdivision.TreeRep                      Data.Geometry.Arrangement                     Data.Geometry.Arrangement.Internal@@ -156,7 +157,7 @@                     Algorithms.Geometry.Diameter.ConvexHull                      -- Algorithms.Geometry.Sweep-+                    Algorithms.Geometry.PolygonTriangulation                     Algorithms.Geometry.PolygonTriangulation.Types                     Algorithms.Geometry.PolygonTriangulation.Triangulate                     Algorithms.Geometry.PolygonTriangulation.MakeMonotone@@ -197,6 +198,7 @@    other-modules:                     Data.Geometry.Matrix.Internal+                    Data.Geometry.Transformation.Internal                      -- * Implementation Internals of Polygons                     Data.Geometry.Polygon.Core@@ -226,7 +228,7 @@   -- other-extensions:   build-depends:                 base                    >= 4.11      &&     < 5-              , hgeometry-combinatorial >= 0.12.0.3+              , hgeometry-combinatorial >= 0.13                , bifunctors              >= 4.1               , bytestring              >= 0.10

@@ -86,11 +86,11 @@ -- showDT :: (Show p, Show r)  => Triangulation p r -> IO () -- showDT = mapM_ print . edgesAsPoints -{- HLINT ignore edgesAsPoints -}+ -- | List add edges as point pairs. edgesAsPoints   :: Triangulation p r -> [(Point 2 r :+ p, Point 2 r :+ p)] edgesAsPoints t = let pts = _positions t-                       in map (\(u,v) -> (pts V.! u, pts V.! v)) . edgesAsVertices $ t+                   in map (bimap (pts V.!) (pts V.!)) . edgesAsVertices $ t  -- | List add edges as VertexID pairs. edgesAsVertices :: Triangulation p r -> [(VertexID,VertexID)]@@ -127,5 +127,5 @@ toPlaneGraph _ tr = PG.PlaneGraph $ g&PPG.vertexData .~ vtxData   where     g       = PPG.fromAdjacencyLists . V.toList . V.imap f $ tr^.neighbours-    f i adj = (VertexId i, VertexId <$> adj)+    f i adj = (VertexId i, C.leftElements $ VertexId <$> adj) -- report in CCW order     vtxData = (\(loc :+ p) -> VertexData loc p) <$> tr^.positions

@@ -25,8 +25,8 @@ import qualified Data.Vector as V import           Witherable -import Data.RealNumber.Rational-type R = RealNumber 5+-- import Data.RealNumber.Rational+-- type R = RealNumber 5  -------------------------------------------------------------------------------- @@ -122,17 +122,17 @@ --------------------------------------------------------------------------------  -tr :: Tree Int-tr = Node 0 [Node 1 [], Node 2 [Node 3 [], Node 2 [], Node 4 [Node 5 []]]]+-- tr :: Tree Int+-- tr = Node 0 [Node 1 [], Node 2 [Node 3 [], Node 2 [], Node 4 [Node 5 []]]] -poly :: PolyLine 2 Int R-poly = case fromPoints [origin :+ 0, Point2 1 1 :+ 1, Point2 2 2 :+ 2, Point2 3 3 :+ 3] of-         Just p -> p+-- poly :: PolyLine 2 Int R+-- poly = case fromPoints [origin :+ 0, Point2 1 1 :+ 1, Point2 2 2 :+ 2, Point2 3 3 :+ 3] of+--          Just p -> p -test = Seq.fromList [0..5]+-- test = Seq.fromList [0..5] -myTree :: Tree Int-myTree = Node {rootLabel = 0, subForest = [Node {rootLabel = 1, subForest = []}-                                       ,Node {rootLabel = 2, subForest = []}-                                       ,Node {rootLabel = 3, subForest = []}]-           }+-- myTree :: Tree Int+-- myTree = Node {rootLabel = 0, subForest = [Node {rootLabel = 1, subForest = []}+--                                        ,Node {rootLabel = 2, subForest = []}+--                                        ,Node {rootLabel = 3, subForest = []}]+--            }

@@ -0,0 +1,14 @@+--------------------------------------------------------------------------------+-- |+-- Module      :  Algorithms.Geometry.PolygonTriangulation+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--------------------------------------------------------------------------------+module Algorithms.Geometry.PolygonTriangulation+  ( triangulate+  , triangulate'+  , computeDiagonals+  ) where++import Algorithms.Geometry.PolygonTriangulation.Triangulate

@@ -71,8 +71,6 @@                   . lefts . map (^._2.core)                   . filter (\mp -> mp^._2.extra == Inside) -- triangulate only the insides                   -- . filter (\f -> f^._1 /= outerFaceId')-                  . F.toList . rawFacePolygons $ monotoneP--    -- -- we alredy know we get the polgyons in *clockwise* order, so skip the-    -- -- check if it is counter clockwise-    -- toCounterClockWiseOrder'' = reverseOuterBoundary+                  . F.toList . internalFacePolygons $ monotoneP+    -- FIXME: we should already get all polygons in CCW order, so no+    -- need for the toClockwiseOrder' call

@@ -26,7 +26,7 @@ import Data.Geometry.LineSegment import Data.Geometry.Point import Data.Geometry.PolyLine hiding (fromPoints)-import Data.Geometry.Polygon hiding (fromPoints, maximumBy, minimumBy)+import Data.Geometry.Polygon hiding (fromPoints) import Data.Geometry.Properties import Data.Geometry.Transformation -- import Linear.Affine hiding (Point, Vector, origin)

@@ -31,6 +31,7 @@ import           Linear.Matrix import           Linear.V3 (V3(..)) + -------------------------------------------------------------------------------- -- * A d-dimensional ball @@ -192,16 +193,36 @@ newtype Touching p = Touching p deriving (Show,Eq,Ord,Functor,F.Foldable,T.Traversable)  -- | No intersection, one touching point, or two points-type instance IntersectionOf (Line 2 r) (Circle p r) = [ NoIntersection-                                                       , Touching (Point 2 r)-                                                       , (Point 2 r, Point 2 r)-                                                       ]+type instance IntersectionOf (Line d r) (Sphere d p r) = [ NoIntersection+                                                         , Touching (Point d r)+                                                         , (Point d r, Point d r)+                                                         ] +instance  {-# OVERLAPPABLE #-} (Ord r, Fractional r, Arity d)+          => Line d r `HasIntersectionWith` Sphere d q r where+  l `intersects` (Sphere (c :+ _) r) = let closest = pointClosestTo  c l+                                       in squaredEuclideanDist c closest <= r +instance {-# OVERLAPPING #-} (Ord r, Num r) => Line 2 r `HasIntersectionWith` Circle p r where+  (Line p' v) `intersects` (Circle (c :+ _) r) = discr >= 0+    where+      (Vector2 vx vy)   = v+      -- (px, py) is the vector/point after translating the circle s.t. it is centered at the+      -- origin+      (Vector2 px py) = p' .-. c++      -- let q lambda be the intersection point. We solve the following equation+      -- solving the equation (q_x)^2 + (q_y)^2 = r^2 then yields the equation+      -- L^2(vx^2 + vy^2) + L2(px*vx + py*vy) + px^2 + py^2 = 0+      -- where L = \lambda+      aa                   = vx^2 + vy^2+      bb                   = 2 * (px * vx + py * vy)+      cc                   = px^2 + py^2 - r^2+      discr                = bb^2 - 4*aa*cc+ instance (Ord r, Floating r) => Line 2 r `IsIntersectableWith` Circle p r where    nonEmptyIntersection = defaultNonEmptyIntersection-   (Line p' v) `intersect` (Circle (c :+ _) r) = case discr `compare` 0 of                                                 LT -> coRec NoIntersection                                                 EQ -> coRec . Touching $ q' (lambda (+))@@ -233,12 +254,19 @@  -- | A line segment may not intersect a circle, touch it, or intersect it -- properly in one or two points.-type instance IntersectionOf (LineSegment 2 p r) (Circle q r) = [ NoIntersection-                                                                , Touching (Point 2 r)-                                                                , Point 2 r-                                                                , (Point 2 r, Point 2 r)-                                                                ]+type instance IntersectionOf (LineSegment d p r) (Sphere d q r) = [ NoIntersection+                                                                  , Touching (Point d r)+                                                                  , Point d r+                                                                  , (Point d r, Point d r)+                                                                  ] +instance (Ord r, Fractional r, Arity d)+          => LineSegment d p r `HasIntersectionWith` Sphere d q r where+  seg `intersects` (Sphere (c :+ _) r) = let closest = pointClosestTo  c  (supportingLine seg)+                                         in case squaredEuclideanDist c closest `compare` r of+                                              LT -> True+                                              EQ -> closest `intersects` seg+                                              GT -> False  instance (Ord r, Floating r) => LineSegment 2 p r `IsIntersectableWith` Circle q r where 

@@ -8,40 +8,60 @@ -- Maintainer  :  Frank Staals -------------------------------------------------------------------------------- module Data.Geometry.BezierSpline(-    BezierSpline (BezierSpline)+    BezierSpline (BezierSpline, Bezier2, Bezier3)   , controlPoints   , fromPointSeq+  , endPoints+  , Data.Geometry.BezierSpline.reverse    , evaluate   , split+  , splitMany+  , splitMonotone+  , splitByPoints+  , extension+  , extend+  , growTo+  , merge   , subBezier   , tangent   , approximate   , parameterOf   , snap--  , pattern Bezier2, pattern Bezier3-+  , intersectB   , colinear-  , lineApproximate   , quadToCubic   ) where -import           Control.Lens                 hiding (Empty)-import qualified Data.Foldable                as F+import           Algorithms.Geometry.ConvexHull.GrahamScan+import           Algorithms.Geometry.SmallestEnclosingBall.RIC+import           Algorithms.Geometry.SmallestEnclosingBall.Types+import           Control.Lens hiding (Empty)+import           Data.Ext+import qualified Data.Foldable as F+import           Data.Geometry.Ball+import           Data.Geometry.Box.Internal import           Data.Geometry.Line+import           Data.Geometry.LineSegment hiding (endPoints) import           Data.Geometry.Point+import           Data.Geometry.PolyLine (PolyLine(..))+import           Data.Geometry.Polygon+import           Data.Geometry.Polygon.Convex hiding (merge) import           Data.Geometry.Properties import           Data.Geometry.Transformation-import           Data.Geometry.Vector-import           Data.LSeq                    (LSeq)-import qualified Data.LSeq                    as LSeq-import           Data.Sequence                (Seq (..))-import qualified Data.Sequence                as Seq-import           Data.Traversable             (fmapDefault, foldMapDefault)+import           Data.Geometry.Vector hiding (init)+import           Data.LSeq (LSeq)+import qualified Data.LSeq as LSeq+import           Data.List (sort)+import qualified Data.List.NonEmpty as NonEmpty+import           Data.Sequence (Seq(..))+import qualified Data.Sequence as Seq+import           Data.Traversable (fmapDefault,foldMapDefault) import           GHC.TypeNats-import qualified Test.QuickCheck              as QC+import qualified Test.QuickCheck as QC +-- import Debug.Trace+ --------------------------------------------------------------------------------  -- | Datatype representing a Bezier curve of degree \(n\) in \(d\)-dimensional space.@@ -49,7 +69,8 @@ -- makeLenses ''BezierSpline  -- | Bezier control points. With n degrees, there are n+1 control points.-controlPoints :: Iso (BezierSpline n1 d1 r1) (BezierSpline n2 d2 r2) (LSeq (1+n1) (Point d1 r1)) (LSeq (1+n2) (Point d2 r2))+controlPoints :: Iso (BezierSpline n1 d1 r1)     (BezierSpline n2 d2 r2)+                     (LSeq (1+n1) (Point d1 r1)) (LSeq (1+n2) (Point d2 r2)) controlPoints = iso _controlPoints BezierSpline  -- | Quadratic Bezier Spline@@ -66,19 +87,32 @@     Bezier3 p q r s = fromPointSeq . Seq.fromList $ [p,q,r,s] {-# COMPLETE Bezier3 #-} +-- | Constructs the Bezier Spline from a given sequence of points.+fromPointSeq :: Seq (Point d r) -> BezierSpline n d r+fromPointSeq = BezierSpline . LSeq.promise . LSeq.fromSeq++ deriving instance (Arity d, Eq r) => Eq (BezierSpline n d r)  type instance Dimension (BezierSpline n d r) = d type instance NumType   (BezierSpline n d r) = r - instance (Arity n, Arity d, QC.Arbitrary r) => QC.Arbitrary (BezierSpline n d r) where   arbitrary = fromPointSeq . Seq.fromList <$> QC.vector (fromIntegral . (1+) . natVal $ C @n) --- | Constructs the Bezier Spline from a given sequence of points.-fromPointSeq :: Seq (Point d r) -> BezierSpline n d r-fromPointSeq = BezierSpline . LSeq.promise . LSeq.fromSeq+{-+instance (Arity n, Arity d, QC.Arbitrary r, Ord r) => QC.Arbitrary (BezierSpline n d r) where+  arbitrary = fromPointSeq . Seq.fromList <$> allDifferent (fromIntegral . (1+) . natVal $ C @n) +-- | Generates a set of unique items.+allDifferent   :: (Ord a, QC.Arbitrary  a) => Int -> QC.Gen [a]+allDifferent n = take n . Set.toList . go maxattempts mempty <$> QC.infiniteList+  where+    maxattempts = 100+    go 0 s _                        = s -- too many attempts+    go t s (x:xs) | Set.size s == n = s+                  | otherwise       = go (t-1) (Set.insert x s) xs+-}  instance (Arity d, Show r) => Show (BezierSpline n d r) where   show (BezierSpline ps) =@@ -100,38 +134,79 @@ instance PointFunctor (BezierSpline n d) where   pmap f = over controlPoints (fmap f) +--------------------------------------------------------------------------------++-- | Convert a quadratic bezier to a cubic bezier.+quadToCubic :: Fractional r => BezierSpline 2 2 r -> BezierSpline 3 2 r+quadToCubic (Bezier2 a (Point b) c) =+  Bezier3 a (Point $ (1/3)*^ (toVec a ^+^ 2*^b)) (Point $ (1/3)*^ (2*^ b ^+^ toVec c)) c++--------------------------------------------------------------------------------++-- | Reverse a BezierSpline+reverse :: (Arity d, Ord r, Num r) => BezierSpline n d r -> BezierSpline n d r+reverse = controlPoints %~ LSeq.reverse++ -- | Evaluate a BezierSpline curve at time t in [0, 1] -- -- pre: \(t \in [0,1]\)-evaluate    :: (Arity d, Ord r, Num r) => BezierSpline n d r -> r -> Point d r+evaluate     :: (Arity d, Eq r, Num r) => BezierSpline n d r -> r -> Point d r+evaluate b 0 = fst $ endPoints b+evaluate b 1 = snd $ endPoints b evaluate b t = evaluate' (b^.controlPoints.to LSeq.toSeq)   where     evaluate' = \case       (p :<| Empty)  -> p       pts@(_ :<| tl) -> let (ini :|> _) = pts in evaluate' $ Seq.zipWith blend ini tl       _              -> error "evaluate: absurd"-     blend p q = p .+^ t *^ (q .-. p) --- | Tangent to the bezier spline at the starting point.+-- | Extract a tangent vector from the first to the second control point. tangent   :: (Arity d, Num r, 1 <= n) => BezierSpline n d r -> Vector d r-tangent b = b^?!controlPoints.ix 1  .-. b^?!controlPoints.ix 0+tangent b = b^?!controlPoints.ix 1 .-. b^?!controlPoints.ix 0 --- | Restrict a Bezier curve to th,e piece between parameters t < u in [0, 1].+-- | Return the endpoints of the Bezier spline.+endPoints   :: BezierSpline n d r -> (Point d r, Point d r)+endPoints b = let (p LSeq.:<| _) = b^.controlPoints+                  (_ LSeq.:|> q) = b^.controlPoints+              in (p,q)+++++-- | Restrict a Bezier curve to the piece between parameters t < u in [0, 1]. subBezier     :: (KnownNat n, Arity d, Ord r, Num r)               => r -> r -> BezierSpline n d r -> BezierSpline n d r subBezier t u = fst . split u . snd . split t ++-- | Compute the convex hull of the control polygon of a 2-dimensional Bezier curve.+--   Should also work in any dimension, but convex hull is not yet implemented.+convexHullB :: (Ord r, Fractional r) => BezierSpline n 2 r -> ConvexPolygon () r+convexHullB = convexHull . NonEmpty.fromList . fmap ext . F.toList . _controlPoints++--------------------------------------------------------------------------------+ -- | Split a Bezier curve at time t in [0, 1] into two pieces.-split :: forall n d r. (KnownNat n, Arity d, Ord r, Num r)-      => r -> BezierSpline n d r -> (BezierSpline n d r, BezierSpline n d r)-split t b | t < 0 || t > 1 = error "Split parameter out of bounds."-          | otherwise      = let n  = fromIntegral $ natVal (C @n)-                                 ps = collect t $ b^.controlPoints-                             in ( fromPointSeq . Seq.take (n + 1) $ ps-                                , fromPointSeq . Seq.drop n       $ ps-                                )+split                 :: forall n d r. (KnownNat n, Arity d, Ord r, Num r)+                      => r -> BezierSpline n d r -> (BezierSpline n d r, BezierSpline n d r)+split t b | t < 0     = error "split: t < 0" -- ++ show t ++ " < 0"+          | t > 1     = error "split: t > 1" -- ++ show t ++ " > 1"+          | otherwise = splitRaw t b ++-- | Split without parameter check. If t outside [0,1], will actually extend the curve+--   rather than split it.+splitRaw     :: forall n d r. (KnownNat n, Arity d, Ord r, Num r)+             => r -> BezierSpline n d r -> (BezierSpline n d r, BezierSpline n d r)+splitRaw t b = let n  = fromIntegral $ natVal (C @n)+                   ps = collect t $ b^.controlPoints+               in ( fromPointSeq . Seq.take (n + 1) $ ps+                  , fromPointSeq . Seq.drop (n + 0) $ ps+                  )++-- | implementation of splitRaw collect   :: (Arity d, Ord r, Num r) => r -> LSeq n (Point d r) -> Seq (Point d r) collect t = go . LSeq.toSeq   where@@ -142,48 +217,168 @@      blend p q = p .+^ t *^ (q .-. p) --- {-+-- | Split a Bezier curve into many pieces.+--   Todo: filter out duplicate parameter values!+splitMany :: forall n d r. (KnownNat n, Arity d, Ord r, Fractional r)+          => [r] -> BezierSpline n d r -> [BezierSpline n d r]+splitMany = splitManySorted . sort . map (restrict "splitMany" 0 1) --- -- | Merge to Bezier pieces. Assumes they can be merged into a single piece of the same degree--- --   (as would e.g. be the case for the result of a 'split' operation).--- --   Does not test whether this is the case!--- merge :: (Arity d, Ord r, Num r) => (Bezier d r, Bezier d r) -> Bezier d r+  where splitManySorted []       b' = [b']+        splitManySorted (t : ts) b' = let (a,c) = split t b'+                                      in a : splitManySorted (map (rescale t) ts) c+        rescale :: r -> r -> r+        rescale 1 _ = 1+        rescale t u = (u - t) / (1 - t) --- -} --- | Approximate Bezier curve by Polyline with given resolution.-approximate :: forall n d r. (KnownNat n, Arity d, Ord r, Fractional r)-            => r -> BezierSpline n d r -> [Point d r]-approximate eps b-    | squaredEuclideanDist p q < eps^2 = [p,q]-    | otherwise                        = let (b1, b2) = split 0.5 b-                                         in approximate eps b1 ++ tail (approximate eps b2)-  where-    p = b^.controlPoints.to LSeq.head-    q = b^.controlPoints.to LSeq.last+-- | Cut a Bezier curve into $x_i$-monotone pieces.+--   Can only be solved exactly for degree 4 or smaller.+--   Only gives rational result for degree 2 or smaller.+--   Currentlly implemented for degree 3.+splitMonotone :: (Arity d, Ord r, Enum r, Floating r) => Int -> BezierSpline 3 d r -> [BezierSpline 3 d r]+splitMonotone i b = splitMany (locallyExtremalParameters i b) b --- | Given a point on (or close to) a Bezier curve, return the corresponding parameter value.---   (For points far away from the curve, the function will return the parameter value of---   an approximate locally closest point to the input point.)-parameterOf      :: (Arity d, Ord r, Fractional r) => BezierSpline n d r -> Point d r -> r-parameterOf b p = binarySearch (qdA p . evaluate b) treshold (1 - treshold)-  where treshold = 0.0001+{-+type family RealTypeConstraint (n :: Nat) (r :: *) :: Constraint where+  RealTypeConstraint 1 r = (Fractional r)+  RealTypeConstraint 2 r = (Fractional r)+  RealTypeConstraint 3 r = (Floating r)+  RealTypeConstraint 4 r = (Floating r)+  RealTypeConstraint 5 r = (Floating r)+  RealTypeConstraint n r = TypeError ""+-} -binarySearch                                    :: (Ord r, Fractional r) => (r -> r) -> r -> r -> r-binarySearch f l r | abs (f l - f r) < treshold = m-                   | derivative f m  > 0        = binarySearch f l m-                   | otherwise                  = binarySearch f m r-  where m = (l + r) / 2-        treshold = 0.0001+-- | Report all parameter values at which the derivative of the $i$th coordinate is 0.+locallyExtremalParameters         :: (Arity d, Ord r, Enum r, Floating r)+                                  => Int -> BezierSpline 3 d r -> [r]+locallyExtremalParameters i curve =+  let [x1, x2, x3, x4] = map (view $ unsafeCoord i) $ F.toList $ _controlPoints curve+      a = 3 * x4 -  9 * x3 + 9 * x2 - 3 * x1+      b = 6 * x1 - 12 * x2 + 6 * x3+      c = 3 * x2 -  3 * x1+  in filter (\j -> 0 <= j && j <= 1) $ solveQuadraticEquation a b c -derivative     :: Fractional r => (r -> r) -> r -> r-derivative f x = (f (x + delta) - f x) / delta-  where delta = 0.00001 --- | Snap a point close to a Bezier curve to the curve.-snap   :: (Arity d, Ord r, Fractional r) => BezierSpline n d r -> Point d r -> Point d r-snap b = evaluate b . parameterOf b+-- | Subdivide a curve based on a sequence of points.+--   Assumes these points are all supposed to lie on the curve, and+--   snaps endpoints of pieces to these points.+--   (higher dimensions would work, but depends on convex hull)+splitByPoints :: (KnownNat n, Ord r, RealFrac r)+              => r -> [Point 2 r] -> BezierSpline n 2 r -> [BezierSpline n 2 r]+splitByPoints treshold points curve =+  let a      = fst $ endPoints curve+      b      = snd $ endPoints curve+      intern = filter (\p -> p /= a && p /= b) points+      times  = map (parameterOf treshold curve) intern+      tipos  = sort $ zip times intern+      pieces = splitMany (map fst tipos) curve+      stapts = a : map snd tipos+      endpts = map snd tipos ++ [b]+  in zipWith3 snapEndpoints stapts endpts pieces +--------------------------------------------------------------------------------++-- | Extend a Bezier curve to a parameter value t outside the interval [0,1].+--   For t < 0, returns a Bezier representation of the section of the underlying curve+--   from parameter value t until paramater value 0. For t > 1, the same from 1 to t.+--+-- pre: t outside [0,1]+extension :: forall n d r. (KnownNat n, Arity d, Ord r, Num r)+      => r -> BezierSpline n d r -> BezierSpline n d r+extension t b | t > 0 && t < 1        = error "extension: 0 < t < 1" -- ++ show t ++ " < 1"+              | t <= 0                = fst $ splitRaw t b+              | otherwise {- t >= 1-} = snd $ splitRaw t b++-- | Extend a Bezier curve to a parameter value t outside the interval [0,1].+--   For t < 0, returns a Bezier representation of the section of the underlying curve+--   from parameter value t until paramater value 1. For t > 1, the same from 0 to t.+--+-- pre: t outside [0,1]+extend :: forall n d r. (KnownNat n, Arity d, Ord r, Num r)+      => r -> BezierSpline n d r -> BezierSpline n d r+extend t b | t > 0 && t < 1         = error "extend: 0 < t < 1" -- ++ show t ++ " < 1"+           | t <= 0                 = snd $ splitRaw t b+           | otherwise {- t >= 1 -} = fst $ splitRaw t b+++-- | Extend a Bezier curve to a point not on the curve, but on / close+--   to the extended underlying curve.+growTo              :: (KnownNat n, Arity d, Ord r, Fractional r)+                    => r -> Point d r -> BezierSpline n d r -> BezierSpline n d r+growTo treshold p b =+  let t = extendedParameterOf treshold b p+      r | t < 0 = extend t b+        | t > 1 = extend t b+        | otherwise = b+  in r++{-++-- | Tries to fit a degree n Bezier curve through a list of points, with error parameter eps.+--   Either returns an appropriate curve, or fails.+fit :: r -> [Point 2 r] -> Maybe (Bezier n d r)+fit eps pts++-}+++--------------------------------------------------------------------------------++-- | Merge two Bezier pieces. Assumes they can be merged into a single piece of the same degree+--   (as would e.g. be the case for the result of a 'split' operation).+--   Does not test whether this is the case!+merge                :: (KnownNat n, Arity d, Ord r, Fractional r)+                     => r -> BezierSpline n d r -> BezierSpline n d r -> BezierSpline n d r+merge treshold b1 b2 = let (p1, q1) = endPoints b1+                           (p2, q2) = endPoints b2+                           result | q1 /= p2 = error "merge: something is wrong, maybe need to flip one of the curves?"+                                  | otherwise = snapEndpoints p1 q2 $ growTo treshold p1 b2+                       in result++-- need distance function between polyBeziers...+++--------------------------------------------------------------------------------+++-- | Approximate Bezier curve by Polyline with given resolution.  That+-- is, every point on the approximation will have distance at most res+-- to the Bezier curve.+approximate     :: (KnownNat n, Arity d, Ord r, Fractional r)+                => r -> BezierSpline n d r -> PolyLine d () r+approximate res = PolyLine . fmap ext . approximate' res++-- | implementation of approximate; returns the polyline as an LSeq+approximate'     :: (KnownNat n, Arity d, Ord r, Fractional r)+                 => r -> BezierSpline n d r -> LSeq 2 (Point d r)+approximate' res = LSeq.promise . LSeq.fromSeq . go+  where+    go b | flat res b = let (p,q) = endPoints b in Seq.fromList [p,q]+         | otherwise  = let (b1, b2) = split 0.5 b in go b1 <> Seq.drop 1 (go b2)++-- | Test whether a Bezier curve can be approximated by a single line segment,+--   given the resolution parameter.+flat :: (KnownNat n, Arity d, Ord r, Fractional r) => r -> BezierSpline n d r -> Bool+flat r b = let p = fst $ endPoints b+               q = snd $ endPoints b+               s = ClosedLineSegment (p :+ ()) (q :+ ())+               e t = sqDistanceToSeg (evaluate b t) s < r ^ 2+           in qdA p q < r ^ 2 || all e [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]++-- seems this is now covered by approximate+--+--+-- -- | Approximate curve as line segments where no point on the curve is further away+-- --   from the nearest line segment than the given tolerance.+-- lineApproximate :: (Ord r, Fractional r) => r -> BezierSpline 3 2 r -> [Point 2 r]+-- lineApproximate eps bezier+--   | colinear eps bezier =+--     [ bezier^.controlPoints.to LSeq.head+--     , bezier^.controlPoints.to LSeq.last ]+--   | otherwise =+--     let (b1, b2) = split 0.5 bezier+--     in lineApproximate eps b1 ++ tail (lineApproximate eps b2)+ -- If both control points are on the same side of the straight line from the start and end -- points then the curve is guaranteed to be within 3/4 of the distance from the straight line -- to the furthest control point.@@ -206,18 +401,251 @@           | sameSide  = 9/16  * maxDist           | otherwise = 16/81 * maxDist --- | Approximate curve as line segments where no point on the curve is further away---   from the nearest line segment than the given tolerance.-lineApproximate :: (Ord r, Fractional r) => r -> BezierSpline 3 2 r -> [Point 2 r]-lineApproximate eps bezier-  | colinear eps bezier =-    [ bezier^.controlPoints.to LSeq.head-    , bezier^.controlPoints.to LSeq.last ]-  | otherwise =-    let (b1, b2) = split 0.5 bezier-    in lineApproximate eps b1 ++ tail (lineApproximate eps b2)+-------------------------------------------------------------------------------- --- | Convert a quadratic bezier to a cubic bezier.-quadToCubic :: (Fractional r) => BezierSpline 2 2 r -> BezierSpline 3 2 r-quadToCubic (Bezier2 a b c) =-  Bezier3 a ((1/3)*^(Point (toVec a ^+^ 2*^toVec b))) ((1/3)*^(Point (2*^ toVec b ^+^ toVec c))) c+-- general d depends on convex hull+-- parameterOf :: (Arity d, Ord r, Fractional r) => BezierSpline n d r -> Point d r -> r+--+-- | Given a point on (or within distance treshold to) a Bezier curve, return the parameter value+--   of some point on the curve within distance treshold from p.+--   For points farther than treshold from the curve, the function will attempt to return the+--   parameter value of an approximate locally closest point to the input point, but no guarantees.+parameterOf :: (KnownNat n, Ord r, RealFrac r) => r -> BezierSpline n 2 r -> Point 2 r -> r+parameterOf treshold b p | closeEnough treshold p $ fst $ endPoints b = 0+                         | closeEnough treshold p $ snd $ endPoints b = 1+                         | otherwise = parameterInterior treshold b p++-- parameterInterior is slow, look into algebraic solution?++-- general d depends on convex hull+parameterInterior :: (KnownNat n, Ord r, RealFrac r) => r -> BezierSpline n 2 r -> Point 2 r -> r+parameterInterior treshold b p | sqrad (F.toList $ view controlPoints b) < (0.5 * treshold)^2 = 0.5+                               | otherwise =+  let (b1, b2) = split 0.5 b+      recurse1 =       0.5 * parameterInterior treshold b1 p+      recurse2 = 0.5 + 0.5 * parameterInterior treshold b2 p+      chb1     = _simplePolygon $ convexHullB b1+      chb2     = _simplePolygon $ convexHullB b2+      in1      = sqDistanceToPolygon p chb1 < treshold^2+      in2      = sqDistanceToPolygon p chb2 < treshold^2+      result |     in1 &&     in2 = betterFit b p recurse1 recurse2+             |     in2 && not in2 = recurse1+             | not in2 &&     in2 = recurse2+             | sqDistanceToPolygon p chb1 < sqDistanceToPolygon p chb2 = recurse1+             | otherwise                                               = recurse2+  in result++-- | Given a point on (or close to) the extension of a Bezier curve, return the corresponding+--   parameter value, which might also be smaller than 0 or larger than 1.+--   (For points far away from the curve, the function will return the parameter value of+--   an approximate locally closest point to the input point.)+--+--   This implementation is not robust: might return a locally closest point on the curve+--   even though the point lies on another part of the curve. For points on the actual+--   curve, use parameterOf instead.+extendedParameterOf      :: (Arity d, KnownNat n, Ord r, Fractional r)+                         => r -> BezierSpline n d r -> Point d r -> r+extendedParameterOf treshold b p | p == fst (endPoints b) = 0+                                 | p == snd (endPoints b) = 1+                                 | otherwise = binarySearch treshold (qdA p . evaluate b) (-100) 100++----------------------------------------+-- * Stuff to implement parameterOf and extendedParameterOf++betterFit         :: (KnownNat n, Arity d, Ord r, Fractional r)+                  => BezierSpline n d r -> Point d r -> r -> r -> r+betterFit b p t u =+  let q = evaluate b t+      r = evaluate b u+  in if qdA q p < qdA r p then t else u++sqDistanceToPolygon :: (Ord r, Fractional r) => Point 2 r -> SimplePolygon p r -> r+sqDistanceToPolygon point poly | insidePolygon point poly = 0+                               | otherwise = minimum $ map (sqDistanceToSeg point) $ listEdges poly+++--------------------------------------------------------------------------------+++++--------------------------------------------------------------------------------++-- | Given two Bezier curves, list all intersection points.+--   Not exact, since for degree >= 3 there is no closed form.+--   (In principle, this algorithm works in any dimension+--   but this requires convexHull, area/volume, and intersect.)+intersectB :: (KnownNat n, Ord r, RealFrac r) => r -> BezierSpline n 2 r -> BezierSpline n 2 r -> [Point 2 r]+intersectB treshold a b+  | a == b    = [fst $ endPoints b, snd $ endPoints b] -- should really return the whole curve+  | otherwise = let [a1, _a2, _a3, a4] = F.toList $ _controlPoints a+                    [b1, _b2, _b3, b4] = F.toList $ _controlPoints b+                in    intersectPointsPoints     treshold [a1, a4] [b1, b4]+                   ++ intersectPointsInterior   treshold [a1, a4] b+                   ++ intersectPointsInterior   treshold [b1, b4] a+                   ++ intersectInteriorInterior treshold [a1, a4, b1, b4] a b+++closeEnough :: (Arity d, Ord r, Fractional r) => r -> Point d r -> Point d r -> Bool+closeEnough treshold p q = qdA p q < treshold ^ 2++intersectPointsPoints :: (Ord r, Fractional r) => r -> [Point 2 r] -> [Point 2 r] -> [Point 2 r]+intersectPointsPoints treshold ps = filter (\q -> any (closeEnough treshold q) ps)++intersectPointsInterior :: (KnownNat n, Ord r, RealFrac r) => r -> [Point 2 r] -> BezierSpline n 2 r -> [Point 2 r]+intersectPointsInterior treshold ps b =+  let [b1, _b2, _b3, b4] = F.toList $ _controlPoints b+      nearc p = closeEnough treshold (snap treshold b p) p+      near1 = closeEnough treshold b1+      near4 = closeEnough treshold b4+  in filter (\p -> nearc p && not (near1 p) && not (near4 p)) ps+++intersectInteriorInterior :: (KnownNat n, Ord r, RealFrac r) => r -> [Point 2 r] -> BezierSpline n 2 r -> BezierSpline n 2 r -> [Point 2 r]+intersectInteriorInterior treshold forbidden a b =+  let cha      = _simplePolygon $ convexHullB a+      chb      = _simplePolygon $ convexHullB b+      (a1, a2) = split 0.5 a+      (b1, b2) = split 0.5 b+      points   = F.toList (view controlPoints a)+              ++ F.toList (view controlPoints b)+      approx   = average points+      done | not (cha `intersectsP` chb) = True+           | sqrad points < treshold^2   = True+           | otherwise                   = False+      result | not (cha `intersectsP` chb)        = []+             | any (closeEnough treshold approx) forbidden = []+             | otherwise                          = [approx]+      recurse = intersectInteriorInterior treshold forbidden a1 b1+             ++ intersectInteriorInterior treshold forbidden a1 b2+             ++ intersectInteriorInterior treshold forbidden a2 b1+             ++ intersectInteriorInterior treshold forbidden a2 b2+  in if done then result else recurse++sqrad :: (Ord r, RealFrac r) => [Point 2 r] -> r+sqrad points | length points < 2 = error "sqrad: not enough points"+sqrad points | otherwise =+  let rationalPoints :: [Point 2 Rational] -- smallestEnclosingDisk fails on Floats+      rationalPoints = map (traverse %~ realToFrac) points+      (a : b : cs) = map (:+ ()) rationalPoints+      diskResult   = smallestEnclosingDisk' a b cs+  in realToFrac $ view squaredRadius $ view enclosingDisk $ diskResult++average :: (Functor t, Foldable t, Arity d, Fractional r) => t (Point d r) -> Point d r+average ps = origin .+^ foldr1 (^+^) (fmap toVec ps) ^/ realToFrac (length ps)++{-+type instance IntersectionOf (BezierSpline n 2 r) (BezierSpline n 2 r) = [ NoIntersection+                                                                                   , [Point 2 r]+                                                                                   , BezierSpline n 2 r+                                                                                   ]+++instance (KnownNat n, Ord r, Fractional r) => (BezierSpline n 2 r) `IsIntersectableWith` (BezierSpline n 2 r) where+  nonEmptyIntersection = defaultNonEmptyIntersection+  a `intersect` b = a `intersectB` b+-}+++-- function to test whether two convex polygons intersect+-- for speed, first test bounding boxes+-- maybe would be faster to directly compare bounding boxes of points, rather than+-- call convex hull first?+intersectsP :: (Ord r, Fractional r) => SimplePolygon p r -> SimplePolygon p r -> Bool+intersectsP p q | not $ boundingBox p `intersects` boundingBox q = False+                | otherwise = or [a `intersects` b | a <- p & listEdges, b <- q & listEdges]+                           || (any (flip insidePolygon p) $ map _core $ F.toList $ polygonVertices q)+                           || (any (flip insidePolygon q) $ map _core $ F.toList $ polygonVertices p)+  -- first test bounding box?+++{-++instance (Arity d, Floating r) => IsBoxable (BezierSpline 3 d r) where+  boundingBox b = foldr1 (<>) $ map (\i -> boundingBox (extremal True i b) <> boundingBox (extremal False i b)) [1 .. d]++-- | Find extremal points on curve in the $i$th dimension.+extremal :: Floating r => Bool -> Int -> BezierSpline 3 d r -> Point d r+extremal pos i b =+  let [p1, _, _, p4] = F.toList $ view controlPoints b+      ps = map evaluate $ locallyExtremalParameters i b+      candidates = [p1, p4] ++ ps+      result | pos     = maximumBy (unsafeCoord i . snd) candidates+             | not pos = minimumBy (unsafeCoord i . snd) candidates+  in result++-}+++--------------------------------------------------------------------------------++snapEndpoints           :: (KnownNat n, Arity d, Ord r, Fractional r)+                        => Point d r -> Point d r -> BezierSpline n d r -> BezierSpline n d r+snapEndpoints p q curve =+  let points = F.toList $ _controlPoints curve+      middle = tail . init $ points+      new    = [p] ++ middle ++ [q]+  in  fromPointSeq $ Seq.fromList new+++-- | Snap a point close to a Bezier curve to the curve.+snap   :: (KnownNat n, Ord r, RealFrac r) => r -> BezierSpline n 2 r -> Point 2 r -> Point 2 r+snap treshold b = evaluate b . parameterOf treshold b++--------------------------------------------------------------------------------+-- * Helper functions++-- | Solve equation of the form ax^2 + bx + c = 0.+--   If there are multiple solutions, report in ascending order.+--   Attempt at a somewhat robust implementation.+solveQuadraticEquation :: (Ord r, Enum r, Floating r) => r -> r -> r -> [r]+solveQuadraticEquation 0 0 0 = [0..] -- error "infinite solutions"+solveQuadraticEquation _ 0 0 = [0]+solveQuadraticEquation 0 _ 0 = [0]+solveQuadraticEquation 0 0 _ = []+solveQuadraticEquation a b 0 = sort [0, -b / a]+solveQuadraticEquation a 0 c | (-c / a) <  0 = []+                             | (-c / a) == 0 = [0]+                             | (-c / a) >  0 = [sqrt (-c / a)]+solveQuadraticEquation 0 b c = [-c / b]+solveQuadraticEquation a b c | almostzero a || almostzero (a / b) || almostzero (a / c) = solveQuadraticEquation 0 b c+solveQuadraticEquation a b c =+  let d = b^2 - 4 * a * c+      result | d == 0 = [-b / (2 * a)]+             | d >  0 = [(-b - sqrt d) / (2 * a), (-b + sqrt d) / (2 * a)]+             | otherwise = []+  in result+  -- trace ("soving equation " ++ show a ++ "x^2 + " ++ show b ++ "x + " ++ show c ++ " = 0") $ result++-- | Test whether a floating point number is close enough to zero, taking rounding errors into account.+almostzero :: (Floating r, Ord r) => r -> Bool+almostzero x = abs x < epsilon++-- | Treshold for rounding errors in almostzero test.+--   TODO: Should be different depending on the type.+epsilon :: Floating r => r+epsilon = 0.0001++++-- | This function tests whether a value lies within bounds of a given interval.+--   If not, graciously continues with value snapped to interval.+--   This should never happen, but apparently it sometimes does?+restrict :: (Ord r) => String -> r -> r -> r -> r+restrict f l r x | l > r = error $ f <> ": restrict [l,r] is not an interval" --error $ f ++ ": restrict: [" ++ show l ++ ", " ++ show r ++ "] is not an interval"+                 --   | x < l = trace (f ++ ": restricting " ++ show x ++ " to [" ++ show l ++ ", " ++ show r ++ "]") l+                 --   | x > r = trace (f ++ ": restricting " ++ show x ++ " to [" ++ show l ++ ", " ++ show r ++ "]") r+                 | otherwise = x+++binarySearch                                    :: (Ord r, Fractional r)+                                                => r -> (r -> r) -> r -> r -> r+binarySearch treshold f l r+    | abs (f l - f r) < treshold = restrict "binarySearch" l r   m+    | derivative f m  > 0        = restrict "binarySearch" l r $ binarySearch treshold f l m+    | otherwise                  = restrict "binarySearch" l r $ binarySearch treshold f m r+  where m = (l + r) / 2++derivative     :: Fractional r => (r -> r) -> r -> r+derivative f x = (f (x + delta) - f x) / delta+  where delta = 0.0000001

@@ -23,7 +23,9 @@  -------------------------------------------------------------------------------- --- | A Quadrant data type+-- | A data type rperesenting the corners of a box.  the order of the+-- Corners is 'northWest, northEast, southEast, southWest', i.e. in+-- clockwise order starting from the topleft. data Corners a = Corners { _northWest  :: !a                          , _northEast  :: !a                          , _southEast  :: !a

@@ -22,7 +22,7 @@ import qualified Data.Foldable as F import           Data.Geometry.Point import           Data.Geometry.Properties-import           Data.Geometry.Transformation+import           Data.Geometry.Transformation.Internal import           Data.Geometry.Vector import qualified Data.Geometry.Vector as V import qualified Data.List.NonEmpty as NE@@ -106,6 +106,8 @@  type instance IntersectionOf (Box d p r) (Box d q r) = '[ NoIntersection, Box d () r] +instance (Ord r, Arity d) => Box d p r `HasIntersectionWith` Box d q r+ instance (Ord r, Arity d) => Box d p r `IsIntersectableWith` Box d q r where   nonEmptyIntersection = defaultNonEmptyIntersection @@ -144,6 +146,9 @@   type instance IntersectionOf (Point d r) (Box d p r) = '[ NoIntersection, Point d r]++instance (Arity d, Ord r) => Point d r `HasIntersectionWith` Box d p r where+  intersects = inBox  instance (Arity d, Ord r) => Point d r `IsIntersectableWith` Box d p r where   nonEmptyIntersection = defaultNonEmptyIntersection

@@ -120,6 +120,7 @@                                                           , Point d r                                                           ] +instance (Ord r, Fractional r) => HalfLine 2 r `HasIntersectionWith` Line 2 r  instance (Ord r, Fractional r) => HalfLine 2 r `IsIntersectableWith` Line 2 r where   nonEmptyIntersection = defaultNonEmptyIntersection@@ -130,6 +131,7 @@     :& RNil  +instance (Ord r, Fractional r) => HalfLine 2 r `HasIntersectionWith` HalfLine 2 r  instance (Ord r, Fractional r) => HalfLine 2 r `IsIntersectableWith` HalfLine 2 r where   nonEmptyIntersection = defaultNonEmptyIntersection@@ -155,6 +157,8 @@            )       :& RNil +instance (Ord r, Fractional r) => LineSegment 2 () r `HasIntersectionWith` HalfLine 2 r+ instance (Ord r, Fractional r) => LineSegment 2 () r `IsIntersectableWith` HalfLine 2 r where   nonEmptyIntersection = defaultNonEmptyIntersection @@ -173,9 +177,11 @@       :& RNil  +instance (Ord r, Fractional r, Arity d) => Point d r `HasIntersectionWith` HalfLine d r where+  intersects = onHalfLine+ instance (Ord r, Fractional r, Arity d) => Point d r `IsIntersectableWith` HalfLine d r where   nonEmptyIntersection = defaultNonEmptyIntersection-  intersects = onHalfLine   p `intersect` hl | p `intersects` hl = coRec p                    | otherwise         = coRec NoIntersection @@ -188,6 +194,8 @@                                                               , Point 2 r                                                               , LineSegment 2 () r                                                               ]+instance (Ord r, Fractional r)+         => HalfLine 2 r `HasIntersectionWith` Boundary (Rectangle p r)  instance (Ord r, Fractional r)          => HalfLine 2 r `IsIntersectableWith` Boundary (Rectangle p r) where@@ -207,6 +215,8 @@                       (True,False)  -> coRec $ ClosedLineSegment (ext o) p                       (True,True)   -> coRec s)     :& RNil+instance (Ord r, Fractional r)+         => HalfLine 2 r `HasIntersectionWith` Rectangle p r  instance (Ord r, Fractional r)          => HalfLine 2 r `IsIntersectableWith` Rectangle p r where

@@ -98,19 +98,19 @@  type instance IntersectionOf (Point d r) (HalfSpace d r) = [NoIntersection, Point d r] -instance (Num r, Ord r, Arity d) => Point d r `IsIntersectableWith` HalfSpace d r where-  nonEmptyIntersection = defaultNonEmptyIntersection-+instance (Num r, Ord r, Arity d) => Point d r `HasIntersectionWith` HalfSpace d r where   q `intersects` h = q `inHalfSpace` h /= Outside +instance (Num r, Ord r, Arity d) => Point d r `IsIntersectableWith` HalfSpace d r where+  nonEmptyIntersection = defaultNonEmptyIntersection   q `intersect` h | q `intersects` h = coRec q                   | otherwise        = coRec NoIntersection  - type instance IntersectionOf (Line d r) (HalfSpace d r) =     [NoIntersection, HalfLine d r, Line d r] +instance (Fractional r, Ord r) => Line 2 r `HasIntersectionWith` HalfSpace 2 r  instance (Fractional r, Ord r) => Line 2 r `IsIntersectableWith` HalfSpace 2 r where   nonEmptyIntersection = defaultNonEmptyIntersection

@@ -48,9 +48,11 @@  type instance IntersectionOf (Point d r) (HyperPlane d r) = [NoIntersection, Point d r] +instance (Num r, Eq r, Arity d) => Point d r `HasIntersectionWith` HyperPlane d r where+  q `intersects` (HyperPlane p n) = n `dot` (q .-. p) == 0+ instance (Num r, Eq r, Arity d) => Point d r `IsIntersectableWith` HyperPlane d r where   nonEmptyIntersection = defaultNonEmptyIntersection-  q `intersects` (HyperPlane p n) = n `dot` (q .-. p) == 0    q `intersect` h | q `intersects` h = coRec q                   | otherwise        = coRec NoIntersection@@ -103,6 +105,8 @@                                    GT -> Above  type instance IntersectionOf (Line 3 r) (Plane r) = [NoIntersection, Point 3 r, Line 3 r]++instance (Eq r, Fractional r) => Line 3 r `HasIntersectionWith` Plane r  instance (Eq r, Fractional r) => Line 3 r `IsIntersectableWith` Plane r where   nonEmptyIntersection = defaultNonEmptyIntersection

@@ -129,6 +129,7 @@  type instance IntersectionOf (Interval a r) (Interval a r) = [NoIntersection, Interval a r] +instance Ord r => Interval a r `HasIntersectionWith` Interval a r instance Ord r => Interval a r `IsIntersectableWith` Interval a r where    nonEmptyIntersection = defaultNonEmptyIntersection

@@ -48,25 +48,27 @@ type instance IntersectionOf (Point d r) (Line d r) = [NoIntersection, Point d r]  -instance (Eq r, Fractional r, Arity d) => Point d r `IsIntersectableWith` Line d r where-  nonEmptyIntersection = defaultNonEmptyIntersection+instance (Eq r, Fractional r, Arity d)      => Point d r `HasIntersectionWith` Line d r where   intersects = onLine+instance {-# OVERLAPPING #-} (Ord r, Num r) => Point 2 r `HasIntersectionWith` Line 2 r where+  intersects = onLine2++instance (Eq r, Fractional r, Arity d)      => Point d r `IsIntersectableWith` Line d r where+  nonEmptyIntersection = defaultNonEmptyIntersection   p `intersect` l | p `intersects` l = coRec p                   | otherwise        = coRec NoIntersection -instance {-# OVERLAPPING #-} (Ord r, Num r)-        => Point 2 r `IsIntersectableWith` Line 2 r where+instance {-# OVERLAPPING #-} (Ord r, Num r) => Point 2 r `IsIntersectableWith` Line 2 r where   nonEmptyIntersection = defaultNonEmptyIntersection-  intersects = onLine2   p `intersect` l | p `intersects` l = coRec p                   | otherwise        = coRec NoIntersection - type instance IntersectionOf (Line 2 r) (Boundary (Rectangle p r)) =   [ NoIntersection, Point 2 r, (Point 2 r, Point 2 r) , LineSegment 2 () r] - instance (Ord r, Fractional r)+         => Line 2 r `HasIntersectionWith` Boundary (Rectangle p r)+instance (Ord r, Fractional r)          => Line 2 r `IsIntersectableWith` Boundary (Rectangle p r) where   nonEmptyIntersection = defaultNonEmptyIntersection @@ -102,7 +104,8 @@ type instance IntersectionOf (Line 2 r) (Rectangle p r) =   [ NoIntersection, Point 2 r, LineSegment 2 () r] -+instance (Ord r, Fractional r)+         => Line 2 r `HasIntersectionWith` Rectangle p r instance (Ord r, Fractional r)          => Line 2 r `IsIntersectableWith` Rectangle p r where   nonEmptyIntersection = defaultNonEmptyIntersection

@@ -171,6 +171,8 @@                                                      , Line 2 r                                                      ] +instance (Eq r, Fractional r) => Line 2 r `HasIntersectionWith` Line 2 r+ instance (Eq r, Fractional r) => Line 2 r `IsIntersectableWith` Line 2 r where  @@ -234,6 +236,21 @@     :& (H $ \(Point2 _ b)   -> Just (vy / vx,b))     :& (H $ \_              -> Nothing)    -- l is a vertical line (through x=0)     :& RNil++-- -- | get values a,b,c s.t. the input line is described by ax + by + c = 0+-- toLinearFunction' :: Line 2 r -> (r,r,r)+-- toLinearFunction' ()++-- | Given a point p and a line l, computes the point q on l closest to p.+pointClosestTo              :: (Fractional r, Arity d) => Point d r -> Line d r -> Point d r+pointClosestTo p (Line a m) = a .+^ (t0 *^ m)+  where+    -- see https://monkeyproofsolutions.nl/wordpress/how-to-calculate-the-shortest-distance-between-a-point-and-a-line/+    t0 = numerator / divisor+    numerator = (p .-. a) `dot` m+    divisor  = m `dot` m++   -- | Result of a side test

@@ -43,7 +43,7 @@ import           Data.Geometry.Point import           Data.Geometry.Properties import           Data.Geometry.SubLine-import           Data.Geometry.Transformation+import           Data.Geometry.Transformation.Internal import           Data.Geometry.Vector import           Data.Ord (comparing) import           Data.Vinyl@@ -234,16 +234,23 @@   instance {-# OVERLAPPING #-} (Ord r, Num r)+         => Point 2 r `HasIntersectionWith` LineSegment 2 p r where+  intersects = onSegment2++instance {-# OVERLAPPING #-} (Ord r, Num r)          => Point 2 r `IsIntersectableWith` LineSegment 2 p r where   nonEmptyIntersection = defaultNonEmptyIntersection-  intersects = onSegment2   p `intersect` seg | p `intersects` seg = coRec p                     | otherwise          = coRec NoIntersection + instance {-# OVERLAPPABLE #-} (Ord r, Fractional r, Arity d)+         => Point d r `HasIntersectionWith` LineSegment d p r where+  intersects = onSegment++instance {-# OVERLAPPABLE #-} (Ord r, Fractional r, Arity d)          => Point d r `IsIntersectableWith` LineSegment d p r where   nonEmptyIntersection = defaultNonEmptyIntersection-  intersects = onSegment   p `intersect` seg | p `intersects` seg = coRec p                     | otherwise          = coRec NoIntersection @@ -265,7 +272,8 @@   -- work in higher dimensions that might allow us to drop the   -- Fractional constraint -+instance (Ord r, Fractional r) =>+         LineSegment 2 p r `HasIntersectionWith` LineSegment 2 p r  instance (Ord r, Fractional r) =>          LineSegment 2 p r `IsIntersectableWith` LineSegment 2 p r where@@ -276,6 +284,9 @@       :& H coRec       :& H (coRec . subLineToSegment)       :& RNil++instance (Ord r, Fractional r) =>+         LineSegment 2 p r `HasIntersectionWith` Line 2 r where  instance (Ord r, Fractional r) =>          LineSegment 2 p r `IsIntersectableWith` Line 2 r where

@@ -22,6 +22,7 @@ import qualified Data.List.NonEmpty as NonEmpty import           Data.Geometry.PlanarSubdivision.Basic import           Data.Geometry.PlanarSubdivision.Merge+import           Data.Geometry.PlanarSubdivision.TreeRep import           Data.Geometry.Polygon import           Data.Proxy @@ -156,3 +157,26 @@ -- mySubDiv = fromSimplePolygons (Id Test) --                               0 --                               (NonEmpty.fromList [simplePg' :+ 1, trianglePG :+ 2])+++++-- type R = Int+-- data MyWorld++-- mySubDiv :: PlanarSubdivision MyWorld Int (Int,Int) String R+-- mySubDiv = undefined++--     faceData xs = FaceData (Seq.fromList xs)++++-- fromTreeRep                                 :: TreeRep v e f r -> PlanarSubdivision s v e f r+-- fromTreeRep (PlanarSD of' (InnerSD ajs fs)) = undefined+++-- fromInnerRep                     :: forall s v e f r. (Ord r, Fractional r)+--                                  => InnerRep v e f r -> PlanarSubdivision s v e () r+-- fromInnerRep f (InnerSD ajs fs) = fromConnectedSegments (Proxy @s) segs+--   where+--     segs = adjs

@@ -1,6 +1,4 @@ {-# LANGUAGE TemplateHaskell #-}-{-# LANGUAGE OverloadedStrings #-}-{-# LANGUAGE PartialTypeSignatures #-} {-# LANGUAGE ScopedTypeVariables #-} -------------------------------------------------------------------------------- -- |@@ -35,14 +33,16 @@                                             , components, component                                             , vertices', vertices                                             , edges', edges-                                            , faces', faces, internalFaces+                                            , faces', internalFaces', faces, internalFaces                                             , darts'-                                            -- , traverseVertices, traverseDarts, traverseFaces+                                            , traverseVertices, traverseDarts, traverseFaces+                                            , mapVertices, mapDarts, mapFaces                                              , headOf, tailOf, twin, endPoints                                              , incidentEdges, incomingEdges, outgoingEdges-                                            , nextIncidentEdge+                                            , nextIncidentEdge, prevIncidentEdge+                                            , nextIncidentEdgeFrom, prevIncidentEdgeFrom                                             , neighboursOf                                              , leftFace, rightFace@@ -58,8 +58,10 @@                                             , faceDataOf                                              , edgeSegment, edgeSegments-                                            , rawFacePolygon, rawFaceBoundary-                                            , rawFacePolygons+                                            , faceBoundary+                                            , internalFacePolygon, internalFacePolygons+                                            , outerFacePolygon, outerFacePolygon'+                                            , facePolygons                                              , VertexId(..), FaceId(..), Dart, World(..) @@ -68,9 +70,14 @@                                             , dataVal                                              , dartMapping, Raw(..)++                                            , asLocalD, asLocalV, asLocalF+                                            , Incident (incidences)+                                            , common, commonVertices, commonDarts, commonFaces                                             ) where  import           Control.Lens hiding (holes, holesOf, (.=))+import           Data.Bifunctor (first, second) import           Data.Coerce import           Data.Ext import qualified Data.Foldable as F@@ -92,6 +99,7 @@                                 , HasDataOf(..)                                 ) import qualified Data.Sequence as Seq+import qualified Data.Set as Set import qualified Data.Vector as V import qualified Data.Vector.Mutable as MV import           GHC.Generics (Generic)@@ -136,6 +144,7 @@   boundingBox = boundingBoxList' . V.toList . _components  +-- | Lens to access a particular component of the planar subdivision. component    :: ComponentId s -> Lens' (PlanarSubdivision s v e f r)                                        (Component s r) component ci = components.singular (ix $ unCI ci)@@ -333,6 +342,10 @@ faces' ps = let n = numFaces ps             in V.fromList $ map (FaceId . VertexId) [0..n-1] +-- | \( O(n) \). Vector of all primal faces.+internalFaces'    :: PlanarSubdivision s v e f r -> V.Vector (FaceId' s)+internalFaces' = V.tail . faces'+ -- | \( O(n) \). Vector of all primal faces with associated data. faces    :: PlanarSubdivision s v e f r -> V.Vector (FaceId' s, FaceData (Dart s) f) faces ps = (\fi -> (fi,ps^.faceDataOf fi)) <$> faces' ps@@ -416,8 +429,8 @@                      in (\d -> g^.dataOf d) <$> ds  --- | Given a dart d that points into some vertex v, report the next--- dart e in the cyclic order around v.+-- | Given a dart d that points into some vertex v, report the next dart in the+-- cyclic (counterclockwise) order around v. -- -- running time: \(O(1)\) nextIncidentEdge      :: Dart s -> PlanarSubdivision s v e f r -> Dart s@@ -425,6 +438,40 @@                             d''      = PG.nextIncidentEdge d' g                         in g^.dataOf d'' +-- | Given a dart d that points into some vertex v, report the+-- previous dart in the cyclic (counterclockwise) order around v.+--+-- running time: \(O(1)\)+--+-- >>> prevIncidentEdge (dart 1 "+1") smallG+-- Dart (Arc 3) +1+prevIncidentEdge      :: Dart s -> PlanarSubdivision s v e f r -> Dart s+prevIncidentEdge d ps = let (_,d',g) = asLocalD d ps+                            d''      = PG.prevIncidentEdge d' g+                        in g^.dataOf d''++-- | Given a dart d that points away from some vertex v, report the+-- next dart in the cyclic (counterclockwise) order around v.+--+--+-- running time: \(O(1)\)+--+nextIncidentEdgeFrom      :: Dart s -> PlanarSubdivision s v e f r -> Dart s+nextIncidentEdgeFrom d ps = let (_,d',g) = asLocalD d ps+                                d''      = PG.nextIncidentEdgeFrom d' g+                            in g^.dataOf d''++-- | Given a dart d that points into away from vertex v, report the previous dart in the+-- cyclic (counterclockwise) order around v.+--+-- running time: \(O(1)\)+--+prevIncidentEdgeFrom      :: Dart s -> PlanarSubdivision s v e f r -> Dart s+prevIncidentEdgeFrom d ps = let (_,d',g) = asLocalD d ps+                                d''      = PG.prevIncidentEdgeFrom d' g+                            in g^.dataOf d''++ -- | All edges incident to vertex v in incoming direction -- (i.e. pointing into v) in counterclockwise order around v. --@@ -467,10 +514,11 @@                      fi       = PG.rightFace d' g                 in g^.dataOf fi --- | The darts on the outer boundary of the face, for internal faces--- the darts are in clockwise order. For the outer face the darts are--- in counterclockwise order, and the darts from various components are in no particular order.---+-- | The darts on the outer boundary of this face. The darts are+-- reported in order along the face. This means that for internal+-- faces the darts are reported in *clockwise* order along the+-- boundary, whereas for the outer face the darts are reported in+-- counter clockwise order. -- -- running time: \(O(k)\), where \(k\) is the output size. outerBoundaryDarts      :: FaceId' s -> PlanarSubdivision s v e f r  -> V.Vector (Dart s)@@ -478,6 +526,7 @@   where     single (_,f',g) = (\d -> g^.dataOf d) <$> PG.boundary f' g + -- | Get the local face and component from a given face. asLocalF                          :: FaceId' s -> PlanarSubdivision s v e f r                                   -> NonEmpty (ComponentId s, FaceId' (Wrap s), Component s r)@@ -487,8 +536,9 @@   where     toLocalF d = let (ci,d',c) = asLocalD d ps in (ci,PG.leftFace d' c,c) --- | The vertices of the outer boundary of the face, for internal faces in--- clockwise order, for the outer face in counter clockwise order.+-- | The vertices of the outer boundary of the face, for internal+-- faces in clockwise order, for the outer face in counter clockwise+-- order. -- -- -- running time: \(O(k)\), where \(k\) is the output size.@@ -526,7 +576,10 @@ asLocalV (VertexId v) ps = let (Raw ci v' _) = ps^?!rawVertexData.ix v                            in (ci,v',ps^.component ci) --- | Note that using the setting part of this lens may be very expensive!!+-- | Lens to access the vertex data+--+-- Note that using the setting part of this lens may be very+-- expensive!!  (O(n)) vertexDataOf               :: VertexId' s                            -> Lens' (PlanarSubdivision s v e f r ) (VertexData r v) vertexDataOf (VertexId vi) = lens get' set''@@ -538,10 +591,17 @@                  in ps&rawVertexData.ix vi.dataVal                .~ (x^.vData)                       &component ci.PG.vertexDataOf wvdi.location .~ (x^.location) ++-- | Get the location of a vertex in the planar subdivision.+--+-- Note that the setting part of this lens may be very expensive!+-- Moreover, use with care (as this may destroy planarity etc.) locationOf   :: VertexId' s -> Lens' (PlanarSubdivision s v e f r ) (Point 2 r) locationOf v = vertexDataOf v.location  +-- | Lens to get the face data of a particular face. Note that the+-- setting part of this lens may be very expensive! (O(n)) faceDataOf    :: FaceId' s -> Lens' (PlanarSubdivision s v e f r)                                     (FaceData (Dart s) f) faceDataOf fi = lens getF setF@@ -562,33 +622,54 @@   type DataOf (PlanarSubdivision s v e f r) (FaceId' s) = f   dataOf f = faceDataOf f.fData +-- | Traverse the vertices of the planar subdivision+traverseVertices :: Applicative g+                 => (VertexId' s -> v -> g v')+                 -> PlanarSubdivision s v e f r -> g (PlanarSubdivision s v' e f r)+traverseVertices h = traverseOf rawVertexData (traverseWith VertexId h) --- -- | Traverse the vertices--- ----- traverseVertices   :: Applicative m---                    => (VertexId' s -> v -> m v')---                    -> PlanarSubdivision s v e f r---                    -> m (PlanarSubdivision s v' e f r)--- traverseVertices f = itraverseOf (vertexData.itraversed) (\i -> f (VertexId i))+-- | Traverse the darts of the Planar subdivision+traverseDarts   :: Applicative g+                => (Dart s -> e -> g e')+                -> PlanarSubdivision s v e f r -> g (PlanarSubdivision s v e' f r)+traverseDarts h = traverseOf rawDartData (traverseWith toEnum h) --- -- | Traverses the darts--- ----- traverseDarts   :: Applicative m---                 => (Dart s -> e -> m e')---                 -> PlanarSubdivision s v e f r---                 -> m (PlaneGraph s v e' f r)--- traverseDarts f = traverseOf (dart) (PG.traverseDarts f) +-- | Traverse the faces of the planar subdivision.+traverseFaces   :: Applicative g+                => (FaceId' s -> f -> g f')+                -> PlanarSubdivision s v e f r -> g (PlanarSubdivision s v e f' r)+traverseFaces h = traverseOf rawFaceData (traverseFaces' h)+  where+    traverseFaces' h' = itraverse (\i -> traverse (h' (FaceId . VertexId $ i))) --- -- | Traverses the faces--- ----- traverseFaces   :: Applicative m---                 => (FaceId' s  -> f -> m f')---                 -> PlaneGraph s v e f r---                 -> m (PlaneGraph s v e f' r)--- traverseFaces f = traverseOf graph (PG.traverseFaces f)+-- | Helper function to implement traver(vertertices|darts|faces)+traverseWith         :: Applicative g+                     => (Int -> w s)+                     -> (w s -> v -> g v')+                     -> V.Vector (Raw ci i v)+                     -> g (V.Vector (Raw ci i v'))+traverseWith mkIdx h = itraverse (\i -> traverse (h $ mkIdx i)) +-------------------------------------------------------------------------------- +-- | Map with index over all faces+mapFaces   :: (FaceId' s -> t -> f')+           -> PlanarSubdivision s v e t r -> PlanarSubdivision s v e f' r+mapFaces h = runIdentity . traverseFaces (\i x -> Identity $ h i x)++-- | Map with index over all vertices+mapVertices   :: (VertexId' s -> t -> v')+              -> PlanarSubdivision s t e f r -> PlanarSubdivision s v' e f r+mapVertices h = runIdentity . traverseVertices (\i x -> Identity $ h i x)++-- | Map with index over all darts+mapDarts   :: (Dart s -> t -> e')+           -> PlanarSubdivision s v t f r -> PlanarSubdivision s v e' f r+mapDarts h = runIdentity . traverseDarts (\i x -> Identity $ h i x)++--------------------------------------------------------------------------------+ -- | Getter for the data at the endpoints of a dart -- -- running time: \(O(1)\)@@ -621,7 +702,9 @@ edgeSegments ps = (\d -> (d,edgeSegment d ps)) <$> edges' ps  --- | Given a dart and the subdivision constructs the line segment representing it+-- | Given a dart and the subdivision constructs the line segment+-- representing it. The segment \(\overline{uv})\) is has \(u\) as its+-- tail and \(v\) as its head. -- -- \(O(1)\) edgeSegment      :: Dart s -> PlanarSubdivision s v e f r -> LineSegment 2 v r :+ e@@ -629,45 +712,106 @@                    in ClosedLineSegment p q :+ ps^.dataOf d  --- | Generates the darts incident to a face, starting with the given dart.+-- | Given a dart d, generates the darts on (the current component of)+-- the boundary of the the face that is to the right of the given+-- dart. The darts are reported in order along the face. This means+-- that for --+-- - (the outer boundary of an) internal faces the darts are reported+--   in *clockwise* order along the boundary,+-- - the "inner" boundary of a face, i.e. the boundary of ahole, the+--   darts are reported in *counter clockwise* order. --+-- Note that this latter case means that in the darts of a a component+-- of the outer face are reported in counter clockwise order.+-- -- \(O(k)\), where \(k\) is the number of darts reported boundary'     :: Dart s -> PlanarSubdivision s v e f r -> V.Vector (Dart s) boundary' d ps = let (_,d',g) = asLocalD d ps                  in (\d'' -> g^.dataOf d'') <$> PG.boundary' d' g ---- | Constructs the outer boundary of the face+-- | The outerboundary of the face as a simple polygon. For internal+-- faces the polygon that is reported has its vertices stored in CCW+-- order (as expected). ----- \(O(k)\), where \(k\) is the complexity of the outer boundary of the face-rawFaceBoundary      :: FaceId' s -> PlanarSubdivision s v e f r -> SimplePolygon v r :+ f-rawFaceBoundary i ps = unsafeFromPoints pts :+ (ps^.dataOf i)+-- pre: FaceId refers to an internal face.+--+-- \(O(k)\), where \(k\) is the complexity of the outer boundary of+-- the face+faceBoundary      :: FaceId' s -> PlanarSubdivision s v e f r -> SimplePolygon v r :+ f+faceBoundary i ps = unsafeFromPoints (reverse pts) :+ (ps^.dataOf i)   where     d   = V.head $ outerBoundaryDarts i ps     pts = (\d' -> PG.vtxDataToExt $ ps^.vertexDataOf (headOf d' ps))        <$> V.toList (boundary' d ps)-+    -- for internal faces boundary' produces the boundary darts in+    -- clockwise order. Hence, we reverse the sequence of points we+    -- obtain to get the points/vertices in CCW order, so that we can+    -- construct a simplepolygon out of them. --- | Constructs the boundary of the given face+-- | Constructs the boundary of the given face. -- -- \(O(k)\), where \(k\) is the complexity of the face-rawFacePolygon      :: FaceId' s -> PlanarSubdivision s v e f r-                    -> SomePolygon v r :+ f-rawFacePolygon i ps = case F.toList $ holesOf i ps of+internalFacePolygon      :: FaceId' s -> PlanarSubdivision s v e f r+                         -> SomePolygon v r :+ f+internalFacePolygon i ps = case F.toList $ holesOf i ps of                         [] -> Left  res                               :+ x                         hs -> Right (MultiPolygon res $ map toHole hs) :+ x   where-    res :+ x = rawFaceBoundary i ps-    toHole d = rawFaceBoundary (leftFace d ps) ps ^. core+    res :+ x = faceBoundary i ps+    toHole d = faceBoundary (leftFace d ps) ps ^. core+-- TODO: Verify that holes are in the right orientation. --- | Lists all *internal* faces of the planar subdivision.-rawFacePolygons    :: PlanarSubdivision s v e f r-                   -> V.Vector (FaceId' s, SomePolygon v r :+ f)-rawFacePolygons ps = fmap (\(i,_) -> (i,rawFacePolygon i ps)) . internalFaces $ ps +-- | Returns a sufficiently large, rectangular, polygon that contains+-- the entire planar subdivision. Each component corresponds to a hole+-- in this polygon.+outerFacePolygon    :: (Num r, Ord r)+                    => PlanarSubdivision s v e f r -> MultiPolygon (Maybe v) r :+ f+outerFacePolygon ps = outerFacePolygon' outer ps & core %~ first (either (const Nothing) Just)+  where+    outer = rectToPolygon . grow 1 . boundingBox $ ps+    rectToPolygon = unsafeFromPoints . reverse . F.toList . corners +-- | Given a sufficiently large outer boundary, draw the outerface as+-- a polygon with a hole.+outerFacePolygon'          :: SimplePolygon v' r+                           -> PlanarSubdivision s v e f r -> MultiPolygon (Either v' v) r :+ f+outerFacePolygon' outer ps = MultiPolygon (first Left outer) holePgs :+ ps^.dataOf i+  where+    i       = outerFaceId ps+    holePgs = map getBoundary . F.toList $ holesOf i ps+    -- get the bondary of a hole. Note that for holes, the function+    -- 'boundary' promisses to report the darts, and therefore the+    -- vertices in CCW order. Hence, we can directly construct a SimplePolygon out of it.+    getBoundary d = unsafeFromPoints . fmap (second Right) $ faceBoundary' (twin d)+    faceBoundary' d = (\d' -> PG.vtxDataToExt $ ps^.vertexDataOf (headOf d' ps))+                      <$> V.toList (boundary' d ps) +-- | Procuces a polygon for each *internal* face of the planar+-- subdivision.+internalFacePolygons    :: PlanarSubdivision s v e f r+                        -> V.Vector (FaceId' s, SomePolygon v r :+ f)+internalFacePolygons ps = fmap (\(i,_) -> (i,internalFacePolygon i ps)) . internalFaces $ ps+++-- | Procuces a polygon for each face of the planar subdivision.+facePolygons    :: (Num r, Ord r)+                => PlanarSubdivision s v e f r+                -> V.Vector (FaceId' s, SomePolygon (Maybe v) r :+ f)+facePolygons ps = V.cons (outerFaceId ps, first Right $ outerFacePolygon ps) ifs+  where+    ifs = wrapJust <$> internalFacePolygons ps+    g :: Bifunctor g => g a b -> g (Maybe a) b+    g = first Just++    wrapJust                 :: (FaceId' s, SomePolygon v r :+ f)+                             -> (FaceId' s, SomePolygon (Maybe v) r :+ f)+    wrapJust (i,(spg :+ f)) = (i,bimap g g spg :+ f)++++-- | Mapping between the internal and extenral darts dartMapping    :: PlanarSubdivision s v e f r -> V.Vector (Dart (Wrap s), Dart s) dartMapping ps = ps^.component (ComponentId 0).PG.dartData @@ -683,3 +827,63 @@ --          $ trianglePG  -- trianglePG = fromPoints . map ext $ [origin, Point2 10 0, Point2 10 10]+++++++++++++++--------------------------------------------------------------------------------+++-- | A class for describing which features (vertex, edge, face) of a planar subdivision+--   can be incident to each other.+class Incident s a b where+  incidences :: PlanarSubdivision s v e f r -> a -> [b]++instance Incident s (VertexId' s) (Dart s) where+  incidences psd i = V.toList (incidentEdges i psd) ++ map twin (V.toList $ incidentEdges i psd)++instance Incident s (VertexId' s) (FaceId' s) where+  incidences psd i = map ((flip leftFace) psd) $ V.toList $ incidentEdges i psd++instance Incident s (Dart s) (VertexId' s) where+  incidences psd i = [headOf i psd, tailOf i psd]++instance Incident s (Dart s) (FaceId' s) where+  incidences psd i = [leftFace i psd, rightFace i psd]++instance Incident s (FaceId' s) (VertexId' s) where+  incidences psd i = V.toList $ boundaryVertices i psd++instance Incident s (FaceId' s) (Dart s) where+  incidences psd i = V.toList (outerBoundaryDarts i psd) ++ map twin (V.toList $ outerBoundaryDarts i psd)++-- | Given two features (vertex, edge, or face) of a subdivision,+--   report all features of a given type that are incident to both.+common :: (Incident s a c, Incident s b c, Ord c) => PlanarSubdivision s v e f r -> a -> b -> [c]+common psd a b = Set.toList $ Set.intersection (Set.fromList $ incidences psd a) (Set.fromList $ incidences psd b)++-- | Given two features (edge or face) of a subdivision, report all+-- vertices that are incident to both.+commonVertices :: (Incident s a (VertexId' s), Incident s b (VertexId' s)) => PlanarSubdivision s v e f r -> a -> b -> [VertexId' s]+commonVertices = common++-- | Given two features (vertex or face) of a subdivision, report all+--   edges that are incident to both.  Returns both darts of each+--   qualifying edge.+commonDarts :: (Incident s a (Dart s), Incident s b (Dart s)) => PlanarSubdivision s v e f r -> a -> b -> [Dart s]+commonDarts = common++-- | Given two features (vertex or edge) of a subdivision, report all+-- faces that are incident to both.+commonFaces :: (Incident s a (FaceId' s), Incident s b (FaceId' s)) => PlanarSubdivision s v e f r -> a -> b -> [FaceId' s]+commonFaces = common

@@ -0,0 +1,530 @@+module Data.Geometry.PlanarSubdivision.Dynamic +  ( splitEdge, unSplitEdge+  , sproutIntoFace+  , splitFace+  ) where++import Control.Lens++import Data.Vector (Vector, toList, (//), empty)+import qualified Data.Vector as V+import Data.List (sort, sortOn, findIndex)++import Data.Functor.Identity+import Data.Ext+import Data.Geometry hiding (Vector, head, imap)+import Data.Geometry.PlanarSubdivision+import Data.Geometry.PlanarSubdivision.Raw++import Data.PlanarGraph (Dart (Dart), Arc (Arc), VertexId (VertexId), FaceId (FaceId), Direction (Positive, Negative))+import Data.PlaneGraph (PlaneGraph)+import qualified Data.PlaneGraph as PG+import Data.PlaneGraph.AdjRep hiding (id, vData, faces)+import qualified Data.PlaneGraph.AdjRep as AR (id, vData, fData, faces, Face (..))++import           Data.List.NonEmpty (NonEmpty(..))+import qualified Data.List.NonEmpty as NonEmpty++import Debug.Trace++import Data.Geometry.PlanarSubdivision.Basic++tracingOn = False++tr :: Show a => String -> a -> a+tr s a | tracingOn = trace ("\9608 " ++ s ++ ": " ++ show a) a+       | otherwise = a+++-- TO DO:+-- ADD EDGE JOINING TWO COMPONENTS+-- CREATE NEW COMPONENT (SINGLE VERTEX)+-- DELETIONS+++-- | Splits a given edge of a planar subdivision by inserting a new vertex on the edges.+--   Increases #vertices and #edges by 1.+splitEdge +  :: (Show v, Show e, Show f, Show r)+  => VertexId' s +  -> VertexId' s +  -> Point 2 r +  -> v +  -> (e -> (e, e)) +  -> PlanarSubdivision s v e f r +  -> PlanarSubdivision s v e f r++splitEdge a b p v f d = +  let (_, la, _) = asLocalV a d+      (_, lb, _) = asLocalV b d+      v' = (freeVertexId d, v)+      fd = freeDart d+      f' (Dart i Positive, e) = ((Dart i Positive, fst $ f e), (fd, snd $ f e))+      f' (Dart i Negative, e) = ((twin fd, fst $ f e), (Dart i Negative, snd $ f e))+  in  tr "splitEdge" $ d & components' %~ fmap (splitEdgeInPlaneGraph la lb p v' f')++-- | Sprouts a new edge from a given vertex into the interior of a given (incident) face.+--   Increases #vertices and #edges by 1.+sproutIntoFace+  :: (Show v, Show e, Show f, Show r)+  => VertexId' s +  -> FaceId' s +  -> Point 2 r +  -> v                       +  -> (e, e)+  -> PlanarSubdivision s v e f r +  -> PlanarSubdivision s v e f r++sproutIntoFace a f p v (e1, e2) d =+  let [ea] = tr "[ea]" $ filter (\e -> headOf e d == a && leftFace e d == f) $ commonDarts d a f+      (_, la, _) = asLocalV a d+      (_, lc, _) = asLocalV (tailOf ea d) d+      v' = (freeVertexId d, v)+      fd = freeDart d+      e1' = (fd, e1)+      e2' = (twin fd, e2)+  in  tr "sproutIntoFace" $ d & components' %~ fmap (sproutIntoFaceInPlaneGraph la lc p v' (e1', e2'))++-- | Inserts a new edge between two given vertices, adjacent to a common face.+--   Increases #edges and #faces by 1.+splitFace+  :: (Show v, Show e, Show f, Show r)+  => VertexId' s +  -> VertexId' s +  -> (e, e)                       +  -> (f -> (f, f)) +  -> PlanarSubdivision s v e f r +  -> PlanarSubdivision s v e f r++splitFace a b e g d =+  let (ca, _, _) = asLocalV a d+      (cb, _, _) = asLocalV b d+  in if ca == cb then splitFaceSameComponent a b e g d+                 else splitFaceDifferentComponents a b e g d++splitFaceSameComponent a b e g d =+  let fs   = commonFaces d a b+      f | length fs == 1 = tr "f(a)" $ headTrace "splitFaceSameComponent f" fs+        | otherwise = tr "f(b)" $ headTrace "splitFaceSameComponent f" $ filter (not . isOuterFace) fs+      [ea] = tr "[ea]" $ filter (\e -> headOf e d == a && leftFace e d == f) $ commonDarts d a f+      [eb] = tr "[eb]" $ filter (\e -> headOf e d == b && leftFace e d == f) $ commonDarts d b f+      (_, la, _) = asLocalV a d+      (_, lb, _) = asLocalV b d+      (_, lc, _) = asLocalV (tailOf ea d) d+      (_, ld, _) = asLocalV (tailOf eb d) d+      (_, lf, _) :| [] = asLocalF f d+      fd = freeDart d+      e' = ((fd, fst e), (twin fd, snd e))+      tf = freeFaceId d+      g' (ef, x) = ((ef, fst $ g x), (tf, snd $ g x))+  in tr "splitFaceSameComponent" $ d & components' %~ fmap (splitFaceInPlaneGraph (tr "la" la) (tr "lb" lb) (tr "lc" lc) (tr "ld" ld) (tr "lf" lf) e' g')++splitFaceDifferentComponents = undefined+++-- | Splits a given edge of a planar subdivision by inserting a new vertex on the edges.+--   Increases #vertices and #edges by 1.+unSplitEdge +  :: (Show v, Show e, Show f, Show r)+  => VertexId' s +  -> ((e, e) -> e)+  -> PlanarSubdivision s v e f r +  -> PlanarSubdivision s v e f r++unSplitEdge b f d = +  let [a, c] = tr "[a, c]" $ toList $ neighboursOf b d+      (_, la, _) = asLocalV a d+      (_, lb, _) = asLocalV b d+      (_, lc, _) = asLocalV c d+      [dab] = filter (\e -> tailOf e d == a) $ commonDarts d a b+      [dcb] = filter (\e -> tailOf e d == c) $ commonDarts d b c+      f' ((di, ei), (dj, ej)) | di == dab = (     dab, f (ei, ej))+                              | di == dcb = (twin dab, f (ei, ej))+                              | otherwise = error "you shouldn't call f' on any other dart"+      -- no longer used: vertex id b and dart id dcb+  in  tr "unSplitEdge" $ d & components' %~ fmap (unSplitEdgeInPlaneGraph la lb lc f')+-- globally, need to restore VertexId and DartIds ???++++++-- nodig:++freeVertexId :: PlanarSubdivision s v e f r -> VertexId' s+freeDart :: PlanarSubdivision s v e f r -> Dart s+freeFaceId :: PlanarSubdivision s v e f r -> FaceId' s++freeVertexId = VertexId . numVertices+freeDart     = flip Dart Positive . Arc . numEdges+freeFaceId   = FaceId . VertexId . numFaces++components' :: (Show v, Show e, Show f, Show r) => Lens' (PlanarSubdivision s v e f r) (Vector (Component' s v e f r))+type Component' s v e f r = PlaneGraph (Wrap s) (VertexId' s, v) (Dart s, e) (FaceId' s, f) r+components' = lens getComponents' setComponents'++getComponents' :: PlanarSubdivision s v e f r -> Vector (Component' s v e f r)+getComponents' p = fmap (addExtraData p) $ p ^. components++addExtraData :: PlanarSubdivision s v e f r -> Component s r -> Component' s v e f r+addExtraData p c = c & PG.vertexData  . traverse %~ (\i -> (i, p ^. dataOf i))+                     & PG.rawDartData . traverse %~ (\i -> (i, p ^. dataOf i))+                     & PG.faceData    . traverse %~ (\i -> (i, p ^. dataOf i))++setComponents' :: (Show v, Show e, Show f, Show r) => PlanarSubdivision s v e f r -> Vector (Component' s v e f r) -> PlanarSubdivision s v e f r+setComponents' p cs = p & components .~ fmap remExtraData cs+                        & rawVertexData .~ (tr "rawVertexData" . vectorise $ getRawVertexData cs)+                        & rawDartData   .~ (tr "rawDartData"   . vectorise $ getRawEdgeData cs)+                        & rawFaceData   .~ (tr "rawFaceData"   . vectorise $ getRawFaceData cs)++getRawVertexData :: Vector (Component' s v e f r) +                 -> [(VertexId' s, Raw s (VertexId' (Wrap s)) v)]+getRawVertexData = concat . imap (\ci g -> map (\(li, VertexData _ (gi, v)) -> (gi, Raw (toEnum ci) li v)) $ toList $ PG.vertices g) . toList++--getEdgeData :: Vector (Component' s v e f r) -> [(Dart s, (Dart s, e))]+--getEdgeData = map (\(a, b) -> (a, (a, b))) . concatMap (toList . (^. PG.rawDartData)) . toList++getRawEdgeData :: Vector (Component' s v e f r)+               -> [(Dart s, Raw s (Dart (Wrap s)) e)]+getRawEdgeData = concat . imap (\ci g -> map (\(li, (gi, e)) -> (gi, Raw (toEnum ci) li e)) $ toList $ PG.darts g) . toList+++--getFaceData :: Vector (Component' s v e f r) -> [(FaceId' s, f)]+--getFaceData = concatMap (toList . (^. PG.faceData)) . toList+++-- data RawFace	s f+-- _faceIdx :: !(Maybe (ComponentId s, FaceId' (Wrap s)))	 +-- _faceDataVal :: !(FaceData (Dart s) f)++-- | Something in this implementation is not right. It makes asLocalF produce an error.+getRawFaceData :: Vector (Component' s v e f r)+               -> [(FaceId' s, RawFace s f)]+getRawFaceData = concat . imap (\ci g -> map (bla ci) $ toList $ PG.faces g) . toList+  where+    bla ci (li, (gi, f)) | isOuterFace gi = (gi, RawFace Nothing (FaceData Empty f))+                         | otherwise      = (gi, RawFace (Just (toEnum ci, li)) (FaceData Empty f))+-- holes are always empty! (where to get them from?)++isOuterFace :: FaceId' s -> Bool+isOuterFace i = fromEnum i == 0++remExtraData :: Component' s v e f r -> Component s r+remExtraData c = c & PG.vertexData  . traverse %~ fst+                   & PG.rawDartData . traverse %~ fst+                   & PG.faceData    . traverse %~ fst+++vectorise :: (Enum i, Show i) => [(i, a)] -> Vector a+vectorise vs = V.replicate (length vs) undefined // map (\(i, a) -> (fromEnum i, a)) vs+++++------------------+-- PLANE GRAPHS --+------------------+++-- INSERTIONS --+++splitEdgeInPlaneGraph +  :: (Show v, Show e, Show f, Show r) +  => VertexId' s +  -> VertexId' s +  -> Point 2 r +  -> v +  -> (e -> (e, e)) +  -> PlaneGraph s v e f r +  -> PlaneGraph s v e f r+-- LET OP! TEST OF a EN b WEL VOORKOMEN!+splitEdgeInPlaneGraph a b p v f +  = tr "splitEdgeInPlaneGraph" +  . PG.fromAdjRep undefined +  . splitEdgeInAdjRep (fromEnum a) (fromEnum b) p v f +  . PG.toAdjRep++sproutIntoFaceInPlaneGraph+  :: (Show v, Show e, Show f, Show r) +  => VertexId' s +  -> VertexId' s +  -> Point 2 r +  -> v +  -> (e, e)+  -> PlaneGraph s v e f r +  -> PlaneGraph s v e f r+sproutIntoFaceInPlaneGraph a c p v e g =+  let ai = fromEnum a+      ci = fromEnum c+  in tr "splitEdgeInPlaneGraph" +   $ PG.fromAdjRep undefined +   $ sproutInAdjRep ai ci p v e+   $ PG.toAdjRep g+++-- PG.toAdjRep :: PlaneGraph s v e f r -> Gr (Vtx v e r) (Face f)+-- PG.fromAdjRep :: proxy s -> Gr (Vtx v e r) (Face f) -> PlaneGraph s v e f r+++splitFaceInPlaneGraph+  :: (Show v, Show e, Show f, Show r)+  => VertexId' s             -- index van vertex a+  -> VertexId' s             -- index van vertex b+  -> VertexId' s             -- index van vertex c+  -> VertexId' s             -- index van vertex d+  -> FaceId' s               -- index van te splitsen face+  -> (e, e)                  -- extra data voor nieuwe edge ab+  -> (f -> (f, f))           -- functie om face data in twee stukken te knippen+  -> PlaneGraph s v e f r -- input graaf+  -> PlaneGraph s v e f r -- output graaf++splitFaceInPlaneGraph a b c d f e h g = +  let ai = fromEnum a+      bi = fromEnum b+      ci = fromEnum c+      di = fromEnum d+      fi = fromEnum $ tr "fi" $ traceShow (g ^. dataOf f) $ PG.tailOf (PG.boundaryDart f g) g+      fj = fromEnum $ tr "fj" $ PG.headOf (PG.boundaryDart f g) g+      -- ^ boundaryDart seems not working either+  in tr "splitFaceInPlaneGraph" +   $ PG.fromAdjRep undefined +   $ splitFaceInAdjRep ai bi ci di fi fj e h +   $ PG.toAdjRep g+++-- DELETIONS --+++unSplitEdgeInPlaneGraph +  :: (Show v, Show e, Show f, Show r) +  => VertexId' s +  -> VertexId' s +  -> VertexId' s +  -> ((e, e) -> e) +  -> PlaneGraph s v e f r +  -> PlaneGraph s v e f r++unSplitEdgeInPlaneGraph a b c f +  = tr "unSplitEdgeInPlaneGraph" +  . PG.fromAdjRep undefined +  . unSplitEdgeInAdjRep (fromEnum a) (fromEnum b) (fromEnum c) f +  . PG.toAdjRep+++-------------+-- ADJREPS --+-------------++-- Gr   +-- adjacencies :: [v]  +-- faces :: [f]   ++-- Vtx  +-- id :: Int  +-- loc :: Point 2 r   +-- adj :: [(Int, e)] +-- vData :: v   ++-- Face   +-- incidentEdge :: (Int, Int)   +-- fData :: f++--deriving instance (Show v, Show f) => Show (Gr v f)+--deriving instance (Show v, Show e, Show r) => Show (Vtx v e r)+--deriving instance Show f => Show (Face f)+++-- instance {-# OVERLAPS #-} Show (VertexId s Primal) where show i = 'v' : show (fromEnum i)+-- instance {-# OVERLAPS #-} Show (FaceId   s Primal) where show i = 'f' : show (fromEnum i)+-- instance {-# OVERLAPS #-} Show (Dart s, v) where +--   show (Dart (Arc s) Positive, _) = 'd' : show (fromEnum s) ++ "+"+--   show (Dart (Arc s) Negative, _) = 'd' : show (fromEnum s) ++ "-"++-- instance Show f => Show (Face f) where show f = (show $ AR.fData f) ++ "~>" ++ (show $ incidentEdge f)+-- instance (Show e, Show r) => Show (Vtx v e r) where show v = (show $ AR.id v) ++ "~>" ++ (show $ adj v)+-- instance (Show v, Show f) => Show (Gr v f) where show g = "Gr " ++ (show $ adjacencies g) ++ " " ++ (show $ AR.faces g)++-- ik heb:+splitEdgeInAdjRep +  :: (Show v, Show e, Show f, Show r)+  => Int                     -- index van vertex a+  -> Int                     -- index van vertex b+  -> Point 2 r               -- locatie voor nieuwe vertex c+  -> v                       -- extra data voor vertex c+  -> (e -> (e, e))           -- functie om edge data in twee stukken te knippen+  -> Gr (Vtx v e r) (Face f) -- input graaf+  -> Gr (Vtx v e r) (Face f) -- output graaf++splitEdgeInAdjRep a b p v f g = +  let n  = length $ adjacencies g+      -- first find vertices a and b+      oa = headTrace "splitEdgeInAdjRep oa" $ filter ((== a) . AR.id) $ adjacencies g+      ob = headTrace "splitEdgeInAdjRep ob" $ filter ((== b) . AR.id) $ adjacencies g+      os = filter ((lift (&&) (/= a) (/= b)) . AR.id) $ adjacencies g+      -- find edge data+      e1 = snd $ headTrace "splitEdgeInAdjRep e1" $ filter ((== b) . fst) $ adj oa+      e2 = snd $ headTrace "splitEdgeInAdjRep e2" $ filter ((== a) . fst) $ adj ob+      -- create new adjacencies to c in a and b+      na = oa {adj = replace ((== b) . fst) (const (n, fst $ f e1)) $ adj oa}+      nb = ob {adj = replace ((== a) . fst) (const (n, fst $ f e2)) $ adj ob}+      -- create new vertex c+      nc = Vtx {AR.id = n, loc = p, adj = [(a, snd $ f e2), (b, snd $ f e1)], AR.vData = v}+      -- update faces (only if incidentEdge happens to point to ab)+      nf = replace ((== (a, b)) . incidentEdge) (\f -> f {incidentEdge = (a, n)}) +         $ replace ((== (b, a)) . incidentEdge) (\f -> f {incidentEdge = (b, n)}) +         $ AR.faces g+  in tr "splitEdgeInAdjRep" $ (tr "original" g) {adjacencies = sortOn AR.id $ na : nb : nc : os, AR.faces = nf}+  ++sproutInAdjRep+  :: (Show v, Show e, Show f, Show r)+  => Int                     -- index van vertex a+  -> Int                     -- index van vertex c (andere kant van edge a)+  -> Point 2 r               -- locatie voor nieuwe vertex c+  -> v                       -- extra data voor vertex c+  -> (e, e)                  -- extra data voor nieuwe edge+  -> Gr (Vtx v e r) (Face f) -- input graaf+  -> Gr (Vtx v e r) (Face f) -- output graaf++sproutInAdjRep a c p v e g =+  let n  = length $ adjacencies g+      -- first find vertex a+      oa = tr "oa" $ headTrace "sproutInAdjRep oa" $ filter ((== a) . AR.id) $ adjacencies g+      os = tr "os" $ filter ((/= a) . AR.id) $ adjacencies g+      -- need to find index of c+      fj (Just x) = x+      fj Nothing  = error "splitFaceInAdjRep got Nothing"      +      ci = tr "ci" $ fj $ findIndex ((== c) . fst) $ adj oa+      -- create new adjacency to new vertex z in a+      na = tr "na" $ oa {adj = take ci (adj oa) ++ (n, fst e) : drop ci (adj oa)}+      -- create new vertex z+      nz = Vtx {AR.id = n, loc = p, adj = [(a, snd e)], AR.vData = v}+  in tr "splitFaceInAdjRep" $ (tr "original" g) {adjacencies = sortOn AR.id $ na : nz : os}++splitFaceInAdjRep +  :: (Show v, Show e, Show f, Show r)+  => Int                     -- index van vertex a+  -> Int                     -- index van vertex b+  -> Int                     -- index van vertex c (andere kant van edge a)+  -> Int                     -- index van vertex d (andere kant van edge b)+  -> Int                     -- index van face edge start+  -> Int                     -- index van face edge eind+  -> (e, e)                  -- extra data voor nieuwe edge ab+  -> (f -> (f, f))           -- functie om face data in twee stukken te knippen+  -> Gr (Vtx v e r) (Face f) -- input graaf+  -> Gr (Vtx v e r) (Face f) -- output graaf++-- is it easier to split a vertex than a face?++splitFaceInAdjRep a b c d u v e f g =+  let +      -- first find vertices a and b+      oa = tr "oa" $ headTrace "splitFaceInAdjRep oa" $ filter ((== a) . AR.id) $ adjacencies g+      ob = tr "ob" $ headTrace "splitFaceInAdjRep ob" $ filter ((== b) . AR.id) $ adjacencies g+      os = tr "os" $ filter ((lift (&&) (/= a) (/= b)) . AR.id) $ adjacencies g+      -- insert new adjacency between a and b+      fj (Just x) = x+      fj Nothing  = error "splitFaceInAdjRep got Nothing"      +      -- need to find indices c and d!+      ci = tr "ci" $ fj $ findIndex ((== c) . fst) $ adj oa+      di = tr "di" $ fj $ findIndex ((== d) . fst) $ adj ob+      -- insert new adjacencies to each other in a and b+      na = tr "na" $ oa {adj = take ci (adj oa) ++ (b, fst e) : drop ci (adj oa)}+      nb = tr "nb" $ ob {adj = take di (adj ob) ++ (a, snd e) : drop di (adj ob)}+      -- find the face that is incident to both a and b+      i  = tr "i"  $ fj $ findIndex ((== (u, v)) . incidentEdge) $ AR.faces g+      fd = tr "fd" $ AR.fData $ AR.faces g !! i+      ef = tr "ef" $ take i (AR.faces g) ++ drop (i + 1) (AR.faces g)+      f1 = tr "f1" $ AR.Face {incidentEdge = (a, b), AR.fData = fst $ f fd}+      f2 = tr "f2" $ AR.Face {incidentEdge = (b, a), AR.fData = snd $ f fd}+  in tr "splitFaceInAdjRep" $ (tr "original" g) {adjacencies = sortOn AR.id $ na : nb : os, AR.faces = ef ++ [f1, f2]}+  +++++unSplitEdgeInAdjRep +  :: (Show v, Show e, Show f, Show r)+  => Int                     -- index van vertex a+  -> Int                     -- index van vertex b (te verwijderen)+  -> Int                     -- index van vertex c+  -> ((e, e) -> e)           -- functie om edge data te mergen+  -> Gr (Vtx v e r) (Face f) -- input graaf+  -> Gr (Vtx v e r) (Face f) -- output graaf++unSplitEdgeInAdjRep a b c f g = +  let n  = length $ adjacencies g+      -- first find vertices a, b and c+      oa = head $ filter ((== a) . AR.id) $ adjacencies g+      ob = head $ filter ((== b) . AR.id) $ adjacencies g+      oc = head $ filter ((== c) . AR.id) $ adjacencies g+      os = filter ((\i -> i /= a && i /= b && i /= c) . AR.id) $ adjacencies g+      -- find edge data+      eab = snd $ head $ filter ((== b) . fst) $ adj oa+      eba = snd $ head $ filter ((== a) . fst) $ adj ob+      ebc = snd $ head $ filter ((== c) . fst) $ adj ob+      ecb = snd $ head $ filter ((== b) . fst) $ adj oc+      -- create new adjacencies between a and c+      na = oa {adj = replace ((== b) . fst) (const (c, f (eab, ebc))) $ adj oa}+      nc = oc {adj = replace ((== b) . fst) (const (a, f (ecb, eba))) $ adj oc}+      nv = sortOn AR.id $ na : nc : os+      -- update faces (only if incidentEdge happens to point to ab or bc)+      nf = replace ((== (a, b)) . incidentEdge) (\f -> f {incidentEdge = (a, c)}) +         $ replace ((== (b, a)) . incidentEdge) (\f -> f {incidentEdge = (c, a)}) +         $ replace ((== (b, c)) . incidentEdge) (\f -> f {incidentEdge = (a, c)}) +         $ replace ((== (c, b)) . incidentEdge) (\f -> f {incidentEdge = (c, a)}) +         $ AR.faces g+      -- restore consecutive numbering: replace vertex n-1 by b+      ng = replaceIndex (n - 1) b $ (tr "original" g) {adjacencies = nv, AR.faces = nf}+  in tr "unSplitEdgeInAdjRep" $ ng++-- Gr   +-- adjacencies :: [v]  +-- faces :: [f]   ++-- Vtx  +-- id :: Int  +-- loc :: Point 2 r   +-- adj :: [(Int, e)] +-- vData :: v   ++-- Face   +-- incidentEdge :: (Int, Int)   +-- fData :: f++replaceIndex :: Int -> Int -> Gr (Vtx v e r) (Face f) -> Gr (Vtx v e r) (Face f)+replaceIndex i j g = g { adjacencies = map (replaceIndexAdjacency i j) $ adjacencies g+                       , AR.faces    = map (replaceIndexFace      i j) $ AR.faces    g+                       }++replaceIndexAdjacency :: Int -> Int -> Vtx v e r -> Vtx v e r+replaceIndexAdjacency i j v = v { AR.id = if AR.id v == i then j else AR.id v+                                , adj   = replace ((== i) . fst) (set _1 j) $ adj v+                                }++replaceIndexFace :: Int -> Int -> Face f -> Face f+replaceIndexFace i j f | fst (incidentEdge f) == i = f {incidentEdge = incidentEdge f & set _1 j}+                       | snd (incidentEdge f) == i = f {incidentEdge = incidentEdge f & set _2 j}+                       | otherwise = f+++-------------+-- HELPERS --+-------------++replace :: (a -> Bool) -> (a -> a) -> [a] -> [a]+replace f g = map $ replace' f g++replace' :: (a -> Bool) -> (a -> a) -> a -> a+replace' f g x | f x = g x+               | otherwise = x++lift :: (a -> b -> c) -> (d -> a) -> (d -> b) -> d -> c+lift f g h x = f (g x) (h x)++++headTrace :: String -> [a] -> a+headTrace s xs | null xs   = error $ s ++ ": head of empty list"+               | otherwise = head xs

@@ -39,6 +39,14 @@ instance (ToJSON ia, ToJSON a) => ToJSON (Raw s ia a) where   toEncoding = genericToEncoding defaultOptions +instance FunctorWithIndex i (Raw ci i) where+  imap f (Raw ci i x) = Raw ci i (f i x)+instance FoldableWithIndex i (Raw ci i) where+  ifoldMap f (Raw _ i x) = f i x+instance TraversableWithIndex i (Raw ci i) where+  itraverse f (Raw ci i x) = Raw ci i <$> f i x++ -- | get the dataVal of a Raw dataVal :: Lens (Raw s ia a) (Raw s ia b) a b dataVal = lens (\(Raw _ _ x) -> x) (\(Raw c i _) y -> Raw c i y)

@@ -0,0 +1,110 @@+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.PlanarSubdivision.TreeRep+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Data types that help encode/decode a planegraph as a JSON/YAML file.+--+--------------------------------------------------------------------------------+module Data.Geometry.PlanarSubdivision.TreeRep( PlanarSD(..)+                                              , Vtx(..)+                                              , myTreeRep+                                              ) where++-- FIXME; uncomment myTreeRep++import Data.Aeson+import Data.PlaneGraph.AdjRep (Vtx(..))+import GHC.Generics (Generic)++import Data.Geometry.Point+import Data.RealNumber.Rational++--------------------------------------------------------------------------------++++-- | Specify the planar subdivison as a tree of components+data PlanarSD v e f r = PlanarSD+  { outerFace :: f           -- ^ outer face+  , inner     :: InnerSD v e f r+  } deriving (Show,Eq,Functor,Generic)++instance (ToJSON r,   ToJSON v, ToJSON e, ToJSON f)     => ToJSON   (PlanarSD v e f r) where+  toEncoding = genericToEncoding defaultOptions+instance (FromJSON r, FromJSON v, FromJSON e, FromJSON f) => FromJSON (PlanarSD v e f r)+++data InnerSD v e f r = InnerSD+  { adjs     :: [Vtx v e r] -- ^ list of vertices and edges in the+                                -- components incident to the outer+                                -- face+  , faces    :: [(f, [InnerSD v e f r])] -- ^ for each internal+                 -- face in the component described by adjs its data,+                 -- and possible holes+  } deriving (Show,Eq,Functor,Generic)++instance (ToJSON r,   ToJSON v, ToJSON e, ToJSON f)     => ToJSON   (InnerSD v e r f) where+  toEncoding = genericToEncoding defaultOptions+instance (FromJSON r, FromJSON v, FromJSON e, FromJSON f) => FromJSON (InnerSD v e r f)++++--------------------------------------------------------------------------------++-- | This represents the following Planar subdivision. Note that the+-- graph is undirected, the arrows are just to indicate what the+-- Positive direction of the darts is.+--+-- ![mySubDiv](docs/Data/Geometry/PlanarSubdivision/mySubDiv.jpg)+myTreeRep :: PlanarSD Int () String (RealNumber 3)+myTreeRep = PlanarSD "f_infty" (InnerSD ads fs)+  where+    fs = [ ("f_1", [])+         , ("f_2", [f5, f6])+         , ("f_3", [])+         , ("f_4", [f7])+         ]++    f5 = InnerSD [ vtx 16 (Point2 3    8) [e 17, e 18]+                 , vtx 17 (Point2 0    7) [e 16, e 18]+                 , vtx 18 (Point2 (-1) 4) [e 16, e 17]+                 ] [("f_5",[])]++    f6 = InnerSD [ vtx 15 (Point2 3   3) [e 14, e 13]+                 , vtx 13 (Point2 6   4) [e 14, e 15]+                 , vtx 14 (Point2 3   6) [e 13, e 15]+                 ] [("f_6",[])]++    f7 = InnerSD [ vtx 19 (Point2 0   9) [e 20, e 23]+                 , vtx 20 (Point2 0   4) [e 19, e 21]+                 , vtx 21 (Point2 15  2) [e 20, e 22]+                 , vtx 22 (Point2 17  5) [e 21, e 23]+                 , vtx 23 (Point2 15  8) [e 19, e 22]+                 ] [("f_7",[f8])]++    f8 = InnerSD [ vtx 24 (Point2 14  6) [e 25, e 26]+                 , vtx 25 (Point2 13  8) [e 24, e 26]+                 , vtx 26 (Point2 12  5) [e 24, e 25]+                 ] [("f_8",[])]++    ads = [ vtx 0 (Point2 0    0)    [e 1, e 4]+          , vtx 1 (Point2 10   2)    [e 0, e 5]+          , vtx 2 (Point2 9    9)    [e 1, e 7, e 3]+          , vtx 3 (Point2 0    10)   [e 2, e 4]+          , vtx 4 (Point2 (-4) 5)    [e 0, e 3]+          , vtx 5 (Point2 15   3)    [e 1, e 6]+          , vtx 6 (Point2 20   6)    [e 5, e 7]+          , vtx 7 (Point2 10   14)   [e 2, e 6, e 8]+          , vtx 8 (Point2 4    13)   [e 7, e 3]+          , vtx 9 (Point2 4    (-4)) [e 10, e 11]+          , vtx 10 (Point2 8   (-4)) [e 11, e 9]+          , vtx 11 (Point2 11  (-2)) [e 10, e 12]+          , vtx 12 (Point2 7   (-1)) [e 9, e 11]+          ]++    e i = (i,())++    vtx i p as = Vtx i p as i

@@ -1,3 +1,4 @@+{-# OPTIONS_GHC -fno-warn-orphans #-} -------------------------------------------------------------------------------- -- | -- Module      :  Data.Geometry.Polygon@@ -121,9 +122,11 @@  type instance IntersectionOf (Point 2 r) (Polygon t p r) = [NoIntersection, Point 2 r] +instance (Fractional r, Ord r) => Point 2 r `HasIntersectionWith` Polygon t p r where+  q `intersects` pg = q `inPolygon` pg /= Outside+ instance (Fractional r, Ord r) => Point 2 r `IsIntersectableWith` Polygon t p r where   nonEmptyIntersection = defaultNonEmptyIntersection-  q `intersects` pg = q `inPolygon` pg /= Outside   q `intersect` pg | q `intersects` pg = coRec q                    | otherwise         = coRec NoIntersection 

@@ -9,10 +9,13 @@  import           Control.Lens import           Data.Ext-import           Data.Geometry.BezierSpline (BezierSpline, lineApproximate, pattern Bezier3)+import           Data.Geometry.BezierSpline (BezierSpline, pattern Bezier3)+import qualified Data.Geometry.BezierSpline as Bezier import           Data.Geometry.Point+import           Data.Geometry.PolyLine(points) import           Data.Geometry.Polygon import qualified Data.Vector.Circular       as CV+import qualified Data.Foldable as F  data PathJoin r   = JoinLine@@ -46,8 +49,8 @@   where     f :: (Point 2 r :+ PathJoin r, Point 2 r :+ PathJoin r) -> CV.CircularVector (Point 2 r :+ ())     f (a :+ JoinLine, _) = CV.singleton (ext a)-    f (a :+ JoinCurve b c, d :+ _) =-      CV.unsafeFromList $ map ext $ init (lineApproximate eps (Bezier3 a b c d))+    f (a :+ JoinCurve b c, d :+ _) = let poly = Bezier.approximate eps (Bezier3 a b c d)+                                     in CV.unsafeFromList . init . F.toList $ poly^.points  approximateSome :: (Ord r, Fractional r) => r -> SomePolygon (PathJoin r) r -> SomePolygon () r approximateSome eps (Left p)  = Left $ approximate eps p

@@ -48,6 +48,9 @@  type instance IntersectionOf (Point 2 r) (Cell r) = '[ NoIntersection, Point 2 r] +instance (Ord r, Fractional r) => Point 2 r `HasIntersectionWith` Cell r where+  p `intersects` c = p `intersects` toBox c+ instance (Ord r, Fractional r) => Point 2 r `IsIntersectableWith` Cell r where   nonEmptyIntersection = defaultNonEmptyIntersection   p `intersect` c = p `intersect` toBox c

@@ -63,6 +63,8 @@   '[Rectangle (a,a) r]  +instance Ord r => Slab o a r `HasIntersectionWith` Slab o a r+ instance Ord r => Slab o a r `IsIntersectableWith` Slab o a r where   nonEmptyIntersection = defaultNonEmptyIntersection @@ -71,6 +73,9 @@      :& H (\i''            -> coRec (Slab i'' :: Slab o a r))      :& RNil +instance Slab Horizontal a r `HasIntersectionWith` Slab Vertical a r where+  _ `intersects` _ = True+ instance Slab Horizontal a r `IsIntersectableWith` Slab Vertical a r where   nonEmptyIntersection _ _ _ = True @@ -105,6 +110,9 @@   [NoIntersection, Line 2 r, LineSegment 2 a r]  instance (Fractional r, Ord r, HasBoundingLines o) =>+         Line 2 r `HasIntersectionWith` Slab o a r++instance (Fractional r, Ord r, HasBoundingLines o) =>          Line 2 r `IsIntersectableWith` Slab o a r where   nonEmptyIntersection = defaultNonEmptyIntersection @@ -132,6 +140,9 @@   [NoIntersection, SubLine 2 () s r]  instance (Fractional r, Ord r, HasBoundingLines o) =>+         SubLine 2 a r r `HasIntersectionWith` Slab o a r++instance (Fractional r, Ord r, HasBoundingLines o) =>          SubLine 2 a r r `IsIntersectableWith` Slab o a r where    nonEmptyIntersection = defaultNonEmptyIntersection@@ -151,6 +162,9 @@  type instance IntersectionOf (LineSegment 2 p r) (Slab o a r) =   [NoIntersection, LineSegment 2 () r]++instance (Fractional r, Ord r, HasBoundingLines o) =>+         LineSegment 2 a r `HasIntersectionWith` Slab o a r  instance (Fractional r, Ord r, HasBoundingLines o) =>          LineSegment 2 a r `IsIntersectableWith` Slab o a r where

@@ -141,6 +141,8 @@                                                                    , Point 2 r                                                                    , SubLine 2 p s r                                                                    ]+instance (Ord r, Fractional r) =>+         SubLine 2 p r r `HasIntersectionWith` SubLine 2 p r r  {- HLINT ignore "Redundant bracket" -} instance (Ord r, Fractional r) =>@@ -164,6 +166,9 @@       s'' = asProperInterval . first (^.extra)           $ s'&start.core .~ toOffset' (s'^.start.extra.core) l               &end.core   .~ toOffset' (s'^.end.extra.core)   l++instance (Ord r, Fractional r) =>+         SubLine 2 p (UnBounded r) r `HasIntersectionWith` SubLine 2 p (UnBounded r) r  instance (Ord r, Fractional r) =>          SubLine 2 p (UnBounded r) r `IsIntersectableWith` SubLine 2 p (UnBounded r) r where

@@ -1,5 +1,3 @@-{-# LANGUAGE Unsafe #-}-{-# LANGUAGE UndecidableInstances #-} -------------------------------------------------------------------------------- -- | -- Module      :  Data.Geometry.Transformation@@ -7,186 +5,55 @@ -- License     :  see the LICENSE file -- Maintainer  :  Frank Staals ---------------------------------------------------------------------------------module Data.Geometry.Transformation where--import           Control.Lens (iso,set,Iso,imap)-import           Data.Geometry.Matrix-import           Data.Geometry.Matrix.Internal (mkRow)-import           Data.Geometry.Point-import           Data.Geometry.Properties-import           Data.Geometry.Vector-import qualified Data.Geometry.Vector as V-import           Data.Proxy-import           GHC.TypeLits--{- $setup->>> import Data.Geometry.LineSegment->>> import Data.Ext--}------------------------------------------------------------------------------------- * Transformations---- | A type representing a Transformation for d dimensional objects-newtype Transformation d r = Transformation { _transformationMatrix :: Matrix (d + 1) (d + 1) r }---- | Transformations and Matrices are isomorphic.-transformationMatrix :: Iso (Transformation d r)       (Transformation d       s)-                            (Matrix (d + 1) (d + 1) r) (Matrix (d + 1) (d + 1) s)-transformationMatrix = iso _transformationMatrix Transformation+module Data.Geometry.Transformation+  ( Transformation(Transformation)+  , transformationMatrix+  , (|.|), identity, inverseOf -deriving instance (Show r, Arity (d + 1)) => Show (Transformation d r)-deriving instance (Eq r, Arity (d + 1))   => Eq (Transformation d r)-deriving instance (Ord r, Arity (d + 1))  => Ord (Transformation d r)-deriving instance Arity (d + 1)           => Functor (Transformation d)-deriving instance Arity (d + 1)           => Foldable (Transformation d)-deriving instance Arity (d + 1)           => Traversable (Transformation d)+  , IsTransformable(..)+  , transformAllBy+  , transformPointFunctor -type instance NumType (Transformation d r) = r+  , translation, scaling, uniformScaling --- | Compose transformations (right to left)-(|.|) :: (Num r, Arity (d + 1)) => Transformation d r -> Transformation d r -> Transformation d r-(Transformation f) |.| (Transformation g) = Transformation $ f `multM` g+  , translateBy, scaleBy, scaleUniformlyBy +  , rotateTo --- if it exists?+  , skewX, rotation, reflection, reflectionV, reflectionH --- | Compute the inverse transformation------ >>> inverseOf $ translation (Vector2 (10.0) (5.0))--- Transformation {_transformationMatrix = Matrix (Vector3 (Vector3 1.0 0.0 (-10.0)) (Vector3 0.0 1.0 (-5.0)) (Vector3 0.0 0.0 1.0))}-inverseOf :: (Fractional r, Invertible (d + 1) r)-          => Transformation d r -> Transformation d r-inverseOf = Transformation . inverse' . _transformationMatrix+  , fitToBox+  , fitToBoxTransform+  ) where +import           Control.Lens+import           Data.Ext+import           Data.Geometry.Box (Rectangle, IsBoxable)+import qualified Data.Geometry.Box as Box+import           Data.Geometry.Properties+import           Data.Geometry.Point+import           Data.Geometry.Transformation.Internal+import           Data.Geometry.Vector ----------------------------------------------------------------------------------- * Transformable geometry objects --- | A class representing types that can be transformed using a transformation-class IsTransformable g where-  transformBy :: Transformation (Dimension g) (NumType g) -> g -> g---- | Apply a transformation to a collection of objects.------ >>> transformAllBy (uniformScaling 2) [Point1 1, Point1 2, Point1 3]--- [Point1 2.0,Point1 4.0,Point1 6.0]-transformAllBy :: (Functor c, IsTransformable g)-               => Transformation (Dimension g) (NumType g) -> c g -> c g-transformAllBy t = fmap (transformBy t)---- | Apply transformation to a PointFunctor, ie something that contains---   points. Polygons, triangles, line segments, etc, are all PointFunctors.------ >>> transformPointFunctor (uniformScaling 2) $ OpenLineSegment (Point1 1 :+ ()) (Point1 2 :+ ())--- OpenLineSegment (Point1 2.0 :+ ()) (Point1 4.0 :+ ())-transformPointFunctor   :: ( PointFunctor g, Fractional r, d ~ Dimension (g r)-                           , Arity d, Arity (d + 1)-                           ) => Transformation d r -> g r -> g r-transformPointFunctor t = pmap (transformBy t)--instance (Fractional r, Arity d, Arity (d + 1))-         => IsTransformable (Point d r) where-  transformBy t = Point . transformBy t . toVec--instance (Fractional r, Arity d, Arity (d + 1))-         => IsTransformable (Vector d r) where-  transformBy (Transformation m) v = f $ m `mult` snoc v 1-    where-      f u   = (/ V.last u) <$> V.init u-------------------------------------------------------------------------------------- * Common transformations---- | Create translation transformation from a vector.------ >>> transformBy (translation $ Vector2 1 2) $ Point2 2 3--- Point2 3.0 5.0-translation   :: (Num r, Arity d, Arity (d + 1))-              => Vector d r -> Transformation d r-translation v = Transformation . Matrix $ imap transRow (snoc v 1)---- | Create scaling transformation from a vector.------ >>> transformBy (scaling $ Vector2 2 (-1)) $ Point2 2 3--- Point2 4.0 (-3.0)-scaling   :: (Num r, Arity d, Arity (d + 1))-          => Vector d r -> Transformation d r-scaling v = Transformation . Matrix $ imap mkRow (snoc v 1)---- | Create scaling transformation from a scalar that is applied---   to all dimensions.------ >>> transformBy (uniformScaling 5) $ Point2 2 3--- Point2 10.0 15.0--- >>> uniformScaling 5 == scaling (Vector2 5 5)--- True--- >>> uniformScaling 5 == scaling (Vector3 5 5 5)--- True-uniformScaling :: (Num r, Arity d, Arity (d + 1)) => r -> Transformation d r-uniformScaling = scaling . pure-------------------------------------------------------------------------------------- * Functions that execute transformations---- | Translate a given point.------ >>> translateBy (Vector2 1 2) $ Point2 2 3--- Point2 3.0 5.0-translateBy :: ( IsTransformable g, Num (NumType g)-               , Arity (Dimension g), Arity (Dimension g + 1)-               ) => Vector (Dimension g) (NumType g) -> g -> g-translateBy = transformBy . translation---- | Scale a given point.------ >>> scaleBy (Vector2 2 (-1)) $ Point2 2 3--- Point2 4.0 (-3.0)-scaleBy :: ( IsTransformable g, Num (NumType g)-           , Arity (Dimension g), Arity (Dimension g + 1)-           ) => Vector (Dimension g) (NumType g) -> g -> g-scaleBy = transformBy . scaling----- | Scale a given point uniformly in all dimensions.------ >>> scaleUniformlyBy 5 $ Point2 2 3--- Point2 10.0 15.0-scaleUniformlyBy :: ( IsTransformable g, Num (NumType g)-                    , Arity (Dimension g), Arity (Dimension g + 1)-                    ) => NumType g -> g -> g-scaleUniformlyBy = transformBy  . uniformScaling----- | Row in a translation matrix--- transRow     :: forall n r. ( Arity n, Arity (n- 1), ((n - 1) + 1) ~ n---                             , Num r) => Int -> r -> Vector n r--- transRow i x = set (V.element (Proxy :: Proxy (n-1))) x $ mkRow i 1--transRow     :: forall n r. (Arity n, Arity (n + 1), Num r)-             => Int -> r -> Vector (n + 1) r-transRow i x = set (V.element (Proxy :: Proxy n)) x $ mkRow i 1------------------------------------------------------------------------------------- * 3D Rotations---- | Given three new unit-length basis vectors (u,v,w) that map to (x,y,z),--- construct the appropriate rotation that does this.-------rotateTo                 :: Num r => Vector 3 (Vector 3 r) -> Transformation 3 r-rotateTo (Vector3 u v w) = Transformation . Matrix $ Vector4 (snoc u        0)-                                                             (snoc v        0)-                                                             (snoc w        0)-                                                             (Vector4 0 0 0 1)------------------------------------------------------------------------------------- * 2D Transformations+-- | Given a rectangle r and a geometry g with its boundingbox,+-- transform the g to fit r.+fitToBox     :: forall g r q.+                ( IsTransformable g, IsBoxable g, NumType g ~ r, Dimension g ~ 2+                , Ord r, Fractional r+                ) => Rectangle q r -> g -> g+fitToBox r g = transformBy (fitToBoxTransform r g) g --- | Skew transformation that keeps the y-coordinates fixed and shifts--- the x coordinates.-skewX        :: Num r => r -> Transformation 2 r-skewX lambda = Transformation . Matrix $ Vector3 (Vector3 1 lambda 0)-                                                 (Vector3 0 1      0)-                                                 (Vector3 0 0      1)+-- | Given a rectangle r and a geometry g with its boundingbox,+-- compute a transformation can fit g to r.+fitToBoxTransform     :: forall g r q. ( IsTransformable g, IsBoxable g+                                       , NumType g ~ r, Dimension g ~ 2+                                       , Ord r, Fractional r+                      ) => Rectangle q r -> g -> Transformation 2 r+fitToBoxTransform r g = translation v2 |.| uniformScaling lam |.| translation v1+  where+    b = Box.boundingBox g+    v1  :: Vector 2 r+    v1  = negate <$> b^.to Box.minPoint.core.vector+    v2  = r^.to Box.minPoint.core.vector+    lam = minimum $ (/) <$> Box.size r <*> Box.size b

@@ -0,0 +1,219 @@+{-# LANGUAGE UndecidableInstances #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.Transformation.Internal+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--------------------------------------------------------------------------------+module Data.Geometry.Transformation.Internal where++import           Control.Lens (iso,set,Iso,imap)+import           Data.Geometry.Matrix+import           Data.Geometry.Matrix.Internal (mkRow)+import           Data.Geometry.Point+import           Data.Geometry.Properties+import           Data.Geometry.Vector+import qualified Data.Geometry.Vector as V+import           Data.Proxy+import           GHC.TypeLits++{- $setup+>>> import Data.Geometry.LineSegment+>>> import Data.Ext+-}++--------------------------------------------------------------------------------+-- * Transformations++-- | A type representing a Transformation for d dimensional objects+newtype Transformation d r = Transformation { _transformationMatrix :: Matrix (d + 1) (d + 1) r }++-- | Transformations and Matrices are isomorphic.+transformationMatrix :: Iso (Transformation d r)       (Transformation d       s)+                            (Matrix (d + 1) (d + 1) r) (Matrix (d + 1) (d + 1) s)+transformationMatrix = iso _transformationMatrix Transformation++deriving instance (Show r, Arity (d + 1)) => Show (Transformation d r)+deriving instance (Eq r, Arity (d + 1))   => Eq (Transformation d r)+deriving instance (Ord r, Arity (d + 1))  => Ord (Transformation d r)+deriving instance Arity (d + 1)           => Functor (Transformation d)+deriving instance Arity (d + 1)           => Foldable (Transformation d)+deriving instance Arity (d + 1)           => Traversable (Transformation d)++type instance NumType (Transformation d r) = r++-- | Compose transformations (right to left)+(|.|) :: (Num r, Arity (d + 1)) => Transformation d r -> Transformation d r -> Transformation d r+(Transformation f) |.| (Transformation g) = Transformation $ f `multM` g++-- | Identity transformation; i.e. the transformation which does not change anything.+identity :: (Num r, Arity (d + 1)) => Transformation d r+identity = Transformation identityMatrix++-- if it exists?++-- | Compute the inverse transformation+--+-- >>> inverseOf $ translation (Vector2 (10.0) (5.0))+-- Transformation {_transformationMatrix = Matrix (Vector3 (Vector3 1.0 0.0 (-10.0)) (Vector3 0.0 1.0 (-5.0)) (Vector3 0.0 0.0 1.0))}+inverseOf :: (Fractional r, Invertible (d + 1) r)+          => Transformation d r -> Transformation d r+inverseOf = Transformation . inverse' . _transformationMatrix++--------------------------------------------------------------------------------+-- * Transformable geometry objects++-- | A class representing types that can be transformed using a transformation+class IsTransformable g where+  transformBy :: Transformation (Dimension g) (NumType g) -> g -> g++-- | Apply a transformation to a collection of objects.+--+-- >>> transformAllBy (uniformScaling 2) [Point1 1, Point1 2, Point1 3]+-- [Point1 2.0,Point1 4.0,Point1 6.0]+transformAllBy :: (Functor c, IsTransformable g)+               => Transformation (Dimension g) (NumType g) -> c g -> c g+transformAllBy t = fmap (transformBy t)++-- | Apply transformation to a PointFunctor, ie something that contains+--   points. Polygons, triangles, line segments, etc, are all PointFunctors.+--+-- >>> transformPointFunctor (uniformScaling 2) $ OpenLineSegment (Point1 1 :+ ()) (Point1 2 :+ ())+-- OpenLineSegment (Point1 2.0 :+ ()) (Point1 4.0 :+ ())+transformPointFunctor   :: ( PointFunctor g, Fractional r, d ~ Dimension (g r)+                           , Arity d, Arity (d + 1)+                           ) => Transformation d r -> g r -> g r+transformPointFunctor t = pmap (transformBy t)++instance (Fractional r, Arity d, Arity (d + 1))+         => IsTransformable (Point d r) where+  transformBy t = Point . transformBy t . toVec++instance (Fractional r, Arity d, Arity (d + 1))+         => IsTransformable (Vector d r) where+  transformBy (Transformation m) v = f $ m `mult` snoc v 1+    where+      f u   = (/ V.last u) <$> V.init u+++--------------------------------------------------------------------------------+-- * Common transformations++-- | Create translation transformation from a vector.+--+-- >>> transformBy (translation $ Vector2 1 2) $ Point2 2 3+-- Point2 3.0 5.0+translation   :: (Num r, Arity d, Arity (d + 1))+              => Vector d r -> Transformation d r+translation v = Transformation . Matrix $ imap transRow (snoc v 1)++-- | Create scaling transformation from a vector.+--+-- >>> transformBy (scaling $ Vector2 2 (-1)) $ Point2 2 3+-- Point2 4.0 (-3.0)+scaling   :: (Num r, Arity d, Arity (d + 1))+          => Vector d r -> Transformation d r+scaling v = Transformation . Matrix $ imap mkRow (snoc v 1)++-- | Create scaling transformation from a scalar that is applied+--   to all dimensions.+--+-- >>> transformBy (uniformScaling 5) $ Point2 2 3+-- Point2 10.0 15.0+-- >>> uniformScaling 5 == scaling (Vector2 5 5)+-- True+-- >>> uniformScaling 5 == scaling (Vector3 5 5 5)+-- True+uniformScaling :: (Num r, Arity d, Arity (d + 1)) => r -> Transformation d r+uniformScaling = scaling . pure+++--------------------------------------------------------------------------------+-- * Functions that execute transformations++-- | Translate a given point.+--+-- >>> translateBy (Vector2 1 2) $ Point2 2 3+-- Point2 3.0 5.0+translateBy :: ( IsTransformable g, Num (NumType g)+               , Arity (Dimension g), Arity (Dimension g + 1)+               ) => Vector (Dimension g) (NumType g) -> g -> g+translateBy = transformBy . translation++-- | Scale a given point.+--+-- >>> scaleBy (Vector2 2 (-1)) $ Point2 2 3+-- Point2 4.0 (-3.0)+scaleBy :: ( IsTransformable g, Num (NumType g)+           , Arity (Dimension g), Arity (Dimension g + 1)+           ) => Vector (Dimension g) (NumType g) -> g -> g+scaleBy = transformBy . scaling+++-- | Scale a given point uniformly in all dimensions.+--+-- >>> scaleUniformlyBy 5 $ Point2 2 3+-- Point2 10.0 15.0+scaleUniformlyBy :: ( IsTransformable g, Num (NumType g)+                    , Arity (Dimension g), Arity (Dimension g + 1)+                    ) => NumType g -> g -> g+scaleUniformlyBy = transformBy  . uniformScaling+++-- | Row in a translation matrix+-- transRow     :: forall n r. ( Arity n, Arity (n- 1), ((n - 1) + 1) ~ n+--                             , Num r) => Int -> r -> Vector n r+-- transRow i x = set (V.element (Proxy :: Proxy (n-1))) x $ mkRow i 1++transRow     :: forall n r. (Arity n, Arity (n + 1), Num r)+             => Int -> r -> Vector (n + 1) r+transRow i x = set (V.element (Proxy :: Proxy n)) x $ mkRow i 1++--------------------------------------------------------------------------------+-- * 3D Rotations++-- | Given three new unit-length basis vectors (u,v,w) that map to (x,y,z),+-- construct the appropriate rotation that does this.+--+--+rotateTo                 :: Num r => Vector 3 (Vector 3 r) -> Transformation 3 r+rotateTo (Vector3 u v w) = Transformation . Matrix $ Vector4 (snoc u        0)+                                                             (snoc v        0)+                                                             (snoc w        0)+                                                             (Vector4 0 0 0 1)++--------------------------------------------------------------------------------+-- * 2D Transformations++-- | Skew transformation that keeps the y-coordinates fixed and shifts+-- the x coordinates.+skewX        :: Num r => r -> Transformation 2 r+skewX lambda = Transformation . Matrix $ Vector3 (Vector3 1 lambda 0)+                                                 (Vector3 0 1      0)+                                                 (Vector3 0 0      1)++-- | Create a matrix that corresponds to a rotation by 'a' radians counter-clockwise+--   around the origin.+rotation :: Floating r => r -> Transformation 2 r+rotation a = Transformation . Matrix $ Vector3 (Vector3 (cos a) (- sin a) 0)+                                               (Vector3 (sin a) (  cos a) 0)+                                               (Vector3 0       0         1)++-- | Create a matrix that corresponds to a reflection in a line through the origin+--   which makes an angle of 'a' radians with the positive 'x'-asis, in counter-clockwise+--   orientation.+reflection :: Floating r => r -> Transformation 2 r+reflection a = rotation a |.| reflectionV |.| rotation (-a)++-- | Vertical reflection+reflectionV :: Num r => Transformation 2 r+reflectionV = Transformation . Matrix $ Vector3 (Vector3 1   0  0)+                                                (Vector3 0 (-1) 0)+                                                (Vector3 0   0  1)++-- | Horizontal reflection+reflectionH :: Num r => Transformation 2 r+reflectionH = Transformation . Matrix $ Vector3 (Vector3 (-1) 0  0)+                                                (Vector3   0  1  0)+                                                (Vector3   0  0  1)

@@ -182,6 +182,8 @@ type instance IntersectionOf (Line 2 r) (Triangle 2 p r) =   [ NoIntersection, Point 2 r, LineSegment 2 () r ] +instance (Fractional r, Ord r) => Line 2 r `HasIntersectionWith` Triangle 2 p r+ instance (Fractional r, Ord r) => Line 2 r `IsIntersectableWith` Triangle 2 p r where    nonEmptyIntersection = defaultNonEmptyIntersection @@ -206,6 +208,8 @@  type instance IntersectionOf (Line 3 r) (Triangle 3 p r) =   [ NoIntersection, Point 3 r, LineSegment 3 () r ]++instance (Fractional r, Ord r) => Line 3 r `HasIntersectionWith` Triangle 3 p r  {- HLINT ignore "Use const" -} instance (Fractional r, Ord r) => Line 3 r `IsIntersectableWith` Triangle 3 p r where

@@ -11,7 +11,8 @@ -- embedding. -- ---------------------------------------------------------------------------------module Data.PlaneGraph( PlaneGraph(PlaneGraph), graph+module Data.PlaneGraph( -- $setup+                        PlaneGraph(PlaneGraph), graph                       , PlanarGraph                       , VertexData(VertexData), vData, location, vtxDataToExt                       , fromSimplePolygon, fromConnectedSegments@@ -31,7 +32,7 @@                       , incidentEdges, incomingEdges, outgoingEdges                       , neighboursOf                       , nextIncidentEdge, prevIncidentEdge-+                      , nextIncidentEdgeFrom, prevIncidentEdgeFrom                        , leftFace, rightFace                       , nextEdge, prevEdge@@ -44,15 +45,93 @@                       , vertexData, faceData, dartData, rawDartData                        , edgeSegment, edgeSegments-                      , rawFacePolygon, rawFaceBoundary-                      , rawFacePolygons+                      , faceBoundary, internalFacePolygon+                      , outerFacePolygon, outerFacePolygon'+                      , facePolygons, facePolygons'                        , VertexId(..), FaceId(..), Dart, World(..), VertexId', FaceId' -                       , withEdgeDistances                       , writePlaneGraph, readPlaneGraph                       ) where  import           Data.PlaneGraph.IO import           Data.PlaneGraph.Core+++--------------------------------------------------------------------------------++-- $setup+-- >>> import Data.Proxy+-- >>> import Data.PlaneGraph.AdjRep(Gr(Gr),Face(Face),Vtx(Vtx))+-- >>> import Data.PlaneGraph.IO(fromAdjRep)+-- >>> import Data.PlanarGraph.Dart(Dart(..),Arc(..))+-- >>> :{+-- let dart i s = Dart (Arc i) (read s)+--     small :: Gr (Vtx Int String Int) (Face String)+--     small = Gr [ Vtx 0 (Point2 0 0) [ (2,"0->2")+--                                     , (1,"0->1")+--                                     , (3,"0->3")+--                                     ] 0+--                , Vtx 1 (Point2 2 2) [ (0,"1->0")+--                                     , (2,"1->2")+--                                     , (3,"1->3")+--                                     ] 1+--                , Vtx 2 (Point2 2 0) [ (0,"2->0")+--                                     , (1,"2->1")+--                                     ] 2+--                , Vtx 3 (Point2 (-1) 4) [ (0,"3->0")+--                                        , (1,"3->1")+--                                        ] 3+--                ]+--                [ Face (2,1) "OuterFace"+--                , Face (0,1) "A"+--                , Face (1,0) "B"+--                ]+--     smallG = fromAdjRep (Proxy :: Proxy ()) small+-- :}+--+--+-- This represents the following graph. Note that the graph is undirected, the+-- arrows are just to indicate what the Positive direction of the darts is.+--+-- ![myGraph](docs/Data/PlaneGraph/small.png)+--+--+-- Here is also a slightly larger example graph:+-- ![myGraph](docs/Data/PlaneGraph/planegraph.png)+--+-- >>> import Data.RealNumber.Rational+-- >>> data MyWorld+-- >>> :{+-- let myPlaneGraph :: PlaneGraph MyWorld Int () String (RealNumber 5)+--     myPlaneGraph = fromAdjRep (Proxy @MyWorld) myPlaneGraphAdjrep+--     myPlaneGraphAdjrep :: Gr (Vtx Int () (RealNumber 5)) (Face String)+--     myPlaneGraphAdjrep = Gr [ vtx 0 (Point2 0   0   ) [e 9, e 5, e 1, e 2]+--                             , vtx 1 (Point2 4   4   ) [e 0, e 5, e 12]+--                             , vtx 2 (Point2 3   7   ) [e 0, e 3]+--                             , vtx 3 (Point2 0   5   ) [e 4, e 2]+--                             , vtx 4 (Point2 3   8   ) [e 3, e 13]+--                             , vtx 5 (Point2 8   1   ) [e 0, e 6, e 8, e 1]+--                             , vtx 6 (Point2 6   (-1)) [e 5, e 9]+--                             , vtx 7 (Point2 9   (-1)) [e 8, e 11]+--                             , vtx 8 (Point2 12  1   ) [e 7, e 12, e 5]+--                             , vtx 9 (Point2 8   (-5)) [e 0, e 10, e 6]+--                             , vtx 10 (Point2 12 (-3)) [e 9, e 11]+--                             , vtx 11 (Point2 14 (-1)) [e 10, e 7]+--                             , vtx 12 (Point2 10 4   ) [e 1, e 8, e 13, e 14]+--                             , vtx 13 (Point2 9  6   ) [e 4, e 14, e 12]+--                             , vtx 14 (Point2 8  5   ) [e 13, e 12]+--                             ]+--                             [ Face (0,9) "OuterFace"+--                             , Face (0,5) "A"+--                             , Face (0,1) "B"+--                             , Face (0,2) "C"+--                             , Face (14,13) "D"+--                             , Face (1,12) "E"+--                             , Face (5,8) "F"+--                             ]+--       where+--         e i = (i,())+--         vtx i p es = Vtx i p es i+-- :}

@@ -27,7 +27,7 @@                                           -- edge. Adjacencies are given in                                           -- arbitrary order                      , vData :: v-                     } deriving (Generic, Functor)+                     } deriving (Generic, Show, Eq, Functor)  instance (ToJSON r,   ToJSON v, ToJSON e)     => ToJSON   (Vtx v e r) where   toEncoding = genericToEncoding defaultOptions

@@ -25,7 +25,7 @@                             , vertices', vertices                            , edges', edges-                           , faces', faces, internalFaces, faces''+                           , faces', internalFaces', faces, internalFaces, faces''                            , darts', darts                            , traverseVertices, traverseDarts, traverseFaces @@ -34,6 +34,7 @@                            , incidentEdges, incomingEdges, outgoingEdges                            , neighboursOf                            , nextIncidentEdge, prevIncidentEdge+                           , nextIncidentEdgeFrom, prevIncidentEdgeFrom                              , leftFace, rightFace@@ -47,41 +48,43 @@                            , vertexData, faceData, dartData, rawDartData                             , edgeSegment, edgeSegments-                           , rawFacePolygon, rawFaceBoundary-                           , rawFacePolygons+                           , faceBoundary, internalFacePolygon+                           , outerFacePolygon, outerFacePolygon'+                           , facePolygons, facePolygons', internalFacePolygons                             , VertexId(..), FaceId(..), Dart, World(..), VertexId', FaceId' -                            , withEdgeDistances                            -- , writePlaneGraph, readPlaneGraph                            ) where  -import           Control.Lens              hiding (holes, holesOf, (.=))+import           Control.Lens hiding (holes, holesOf, (.=)) import           Data.Aeson+import           Data.Bifunctor (first) import           Data.Ext-import qualified Data.Foldable             as F-import           Data.Function             (on)+import qualified Data.Foldable as F+import           Data.Function (on) import           Data.Geometry.Box import           Data.Geometry.Interval-import           Data.Geometry.Line        (cmpSlope, supportingLine)+import           Data.Geometry.Line (cmpSlope, supportingLine) import           Data.Geometry.LineSegment hiding (endPoints) import           Data.Geometry.Point import           Data.Geometry.Polygon import           Data.Geometry.Properties-import qualified Data.List.NonEmpty        as NonEmpty-import qualified Data.Map                  as M-import           Data.Ord                  (comparing)+import qualified Data.List.NonEmpty as NonEmpty+import qualified Data.Map as M+import           Data.Ord (comparing) import           Data.PlanarGraph          (Arc (..), Dart (..), Direction (..), FaceId (..),                                             FaceId', HasDataOf (..), PlanarGraph, VertexId (..),                                             VertexId', World (..), dual, planarGraph, twin)-import qualified Data.PlanarGraph          as PG+import qualified Data.PlanarGraph as PG import           Data.Util-import qualified Data.Vector               as V-import           Data.Vector.Circular      (CircularVector)-import           GHC.Generics              (Generic)+import qualified Data.Vector as V+import           Data.Vector.Circular (CircularVector)+import           GHC.Generics (Generic) + --------------------------------------------------------------------------------  -- $setup@@ -119,7 +122,45 @@ -- arrows are just to indicate what the Positive direction of the darts is. -- -- ![myGraph](docs/Data/PlaneGraph/small.png)-+--+--+-- Here is also a slightly larger example graph:+-- ![myGraph](docs/Data/PlaneGraph/planegraph.png)+--+-- >>> import Data.RealNumber.Rational+-- >>> data MyWorld+-- >>> :{+-- let myPlaneGraph :: PlaneGraph MyWorld Int () String (RealNumber 5)+--     myPlaneGraph = fromAdjRep (Proxy @MyWorld) myPlaneGraphAdjrep+--     myPlaneGraphAdjrep :: Gr (Vtx Int () (RealNumber 5)) (Face String)+--     myPlaneGraphAdjrep = Gr [ vtx 0 (Point2 0   0   ) [e 9, e 5, e 1, e 2]+--                             , vtx 1 (Point2 4   4   ) [e 0, e 5, e 12]+--                             , vtx 2 (Point2 3   7   ) [e 0, e 3]+--                             , vtx 3 (Point2 0   5   ) [e 4, e 2]+--                             , vtx 4 (Point2 3   8   ) [e 3, e 13]+--                             , vtx 5 (Point2 8   1   ) [e 0, e 6, e 8, e 1]+--                             , vtx 6 (Point2 6   (-1)) [e 5, e 9]+--                             , vtx 7 (Point2 9   (-1)) [e 8, e 11]+--                             , vtx 8 (Point2 12  1   ) [e 7, e 12, e 5]+--                             , vtx 9 (Point2 8   (-5)) [e 0, e 10, e 6]+--                             , vtx 10 (Point2 12 (-3)) [e 9, e 11]+--                             , vtx 11 (Point2 14 (-1)) [e 10, e 7]+--                             , vtx 12 (Point2 10 4   ) [e 1, e 8, e 13, e 14]+--                             , vtx 13 (Point2 9  6   ) [e 4, e 14, e 12]+--                             , vtx 14 (Point2 8  5   ) [e 13, e 12]+--                             ]+--                             [ Face (0,9) "OuterFace"+--                             , Face (0,5) "A"+--                             , Face (0,1) "B"+--                             , Face (0,2) "C"+--                             , Face (14,13) "D"+--                             , Face (1,12) "E"+--                             , Face (5,8) "F"+--                             ]+--       where+--         e i = (i,())+--         vtx i p es = Vtx i p es i+-- :}  -------------------------------------------------------------------------------- -- * Vertex Data@@ -131,6 +172,7 @@                                             ,Functor,Foldable,Traversable) makeLenses ''VertexData +-- | Convert to an Ext vtxDataToExt                  :: VertexData r v -> Point 2 r :+ v vtxDataToExt (VertexData p v) = p :+ v @@ -235,6 +277,8 @@ -- -- >>> numVertices smallG -- 4+-- >>> numVertices myPlaneGraph+-- 15 numVertices :: PlaneGraph s v e f r  -> Int numVertices = PG.numVertices . _graph @@ -242,6 +286,7 @@ -- -- >>> numDarts smallG -- 10+-- numDarts :: PlaneGraph s v e f r  -> Int numDarts = PG.numDarts . _graph @@ -256,6 +301,8 @@ -- -- >>> numFaces smallG -- 3+-- >>> numFaces myPlaneGraph+-- 7 numFaces :: PlaneGraph s v e f r  -> Int numFaces = PG.numFaces . _graph @@ -281,9 +328,12 @@ darts' = PG.darts' . _graph  -- | Get all darts together with their data+--+-- darts :: PlaneGraph s v e f r  -> V.Vector (Dart s, e) darts = PG.darts . _graph + -- | Enumerate all edges. We report only the Positive darts edges' :: PlaneGraph s v e f r  -> V.Vector (Dart s) edges' = PG.edges' . _graph@@ -327,16 +377,30 @@ faces' :: PlaneGraph s v e f r  -> V.Vector (FaceId' s) faces' = PG.faces' . _graph ++-- | face Ids of all internal faces in the plane graph+--+-- running time: \(O(n)\)+internalFaces'   :: (Ord r, Fractional r) => PlaneGraph s v e f r  -> V.Vector (FaceId' s)+internalFaces' g = let i = outerFaceId g in V.filter (/= i) $ faces' g+ -- | All faces with their face data. -- -- >>> mapM_ print $ faces smallG -- (FaceId 0,"OuterFace") -- (FaceId 1,"A") -- (FaceId 2,"B")+-- >>> mapM_ print $ faces myPlaneGraph+-- (FaceId 0,"OuterFace")+-- (FaceId 1,"A")+-- (FaceId 2,"B")+-- (FaceId 3,"C")+-- (FaceId 4,"E")+-- (FaceId 5,"F")+-- (FaceId 6,"D") faces :: PlaneGraph s v e f r  -> V.Vector (FaceId' s, f) faces = PG.faces . _graph - -- | Reports the outerface and all internal faces separately. -- running time: \(O(n)\) faces''   :: (Ord r, Fractional r)@@ -384,6 +448,11 @@ -- -- >>> incidentEdges (VertexId 1) smallG -- [Dart (Arc 1) -1,Dart (Arc 4) +1,Dart (Arc 3) +1]+-- >>> mapM_ print $ incidentEdges (VertexId 5) myPlaneGraph+-- Dart (Arc 1) -1+-- Dart (Arc 7) +1+-- Dart (Arc 10) +1+-- Dart (Arc 4) -1 incidentEdges   :: VertexId' s -> PlaneGraph s v e f r -> V.Vector (Dart s) incidentEdges v = PG.incidentEdges v . _graph @@ -417,31 +486,69 @@ -- -- >>> neighboursOf (VertexId 1) smallG -- [VertexId 0,VertexId 2,VertexId 3]+-- >>> neighboursOf (VertexId 5) myPlaneGraph+-- [VertexId 0,VertexId 6,VertexId 8,VertexId 1] neighboursOf   :: VertexId' s -> PlaneGraph s v e f r                -> V.Vector (VertexId' s) neighboursOf v = PG.neighboursOf v . _graph  -- | Given a dart d that points into some vertex v, report the next dart in the--- cyclic order around v in clockwise direction.+-- cyclic (counterclockwise) order around v. -- -- running time: \(O(1)\) -- -- >>> nextIncidentEdge (dart 1 "+1") smallG--- Dart (Arc 2) +1+-- Dart (Arc 4) +1+-- >>> nextIncidentEdge (dart 1 "+1") myPlaneGraph+-- Dart (Arc 7) +1+-- >>> nextIncidentEdge (dart 17 "-1") myPlaneGraph+-- Dart (Arc 15) -1 nextIncidentEdge   :: Dart s -> PlaneGraph s v e f r -> Dart s nextIncidentEdge d = PG.nextIncidentEdge d . _graph --- | Given a dart d that points into some vertex v, report the next dart in the--- cyclic order around v (in clockwise order)+-- | Given a dart d that points into some vertex v, report the previous dart in the+-- cyclic (counterclockwise) order around v. -- -- running time: \(O(1)\) -- -- >>> prevIncidentEdge (dart 1 "+1") smallG--- Dart (Arc 0) +1+-- Dart (Arc 3) +1+-- >>> prevIncidentEdge (dart 1 "+1") myPlaneGraph+-- Dart (Arc 4) -1+-- >>> prevIncidentEdge (dart 7 "-1") myPlaneGraph+-- Dart (Arc 1) -1 prevIncidentEdge   :: Dart s -> PlaneGraph s v e f r -> Dart s prevIncidentEdge d = PG.prevIncidentEdge d . _graph  +-- | Given a dart d that points away from some vertex v, report the+-- next dart in the cyclic (counterclockwise) order around v.+--+--+-- running time: \(O(1)\)+--+-- >>> nextIncidentEdgeFrom (dart 1 "+1") smallG+-- Dart (Arc 2) +1+-- >>> nextIncidentEdgeFrom (dart 1 "+1") myPlaneGraph+-- Dart (Arc 2) +1+-- >>> nextIncidentEdgeFrom (dart 4 "+1") myPlaneGraph+-- Dart (Arc 15) +1+nextIncidentEdgeFrom   :: Dart s -> PlaneGraph s v e f r -> Dart s+nextIncidentEdgeFrom d = PG.nextIncidentEdgeFrom d . _graph++-- | Given a dart d that points into away from vertex v, report the previous dart in the+-- cyclic (counterclockwise) order around v.+--+-- running time: \(O(1)\)+--+-- >>> prevIncidentEdgeFrom (dart 1 "+1") smallG+-- Dart (Arc 0) +1+-- >>> prevIncidentEdgeFrom (dart 4 "+1") myPlaneGraph+-- Dart (Arc 2) -1+prevIncidentEdgeFrom   :: Dart s -> PlaneGraph s v e f r -> Dart s+prevIncidentEdgeFrom d = PG.prevIncidentEdgeFrom d . _graph++ -- | The face to the left of the dart -- -- running time: \(O(1)\).@@ -494,20 +601,34 @@ prevEdge d = PG.prevEdge d . _graph  --- | The darts bounding this face, for internal faces in clockwise order, for--- the outer face in counter clockwise order.---+-- | The darts bounding this face. The darts are reported in order+-- along the face. This means that for internal faces the darts are+-- reported in *clockwise* order along the boundary, whereas for the+-- outer face the darts are reported in counter clockwise order. -- -- running time: \(O(k)\), where \(k\) is the output size. --+-- >>> boundary (FaceId $ VertexId 2) smallG -- around face B+-- [Dart (Arc 2) +1,Dart (Arc 3) -1,Dart (Arc 1) -1]+-- >>> boundary (FaceId $ VertexId 0) smallG -- around outer face+-- [Dart (Arc 0) +1,Dart (Arc 4) -1,Dart (Arc 3) +1,Dart (Arc 2) -1] -- boundary   :: FaceId' s -> PlaneGraph s v e f r  -> V.Vector (Dart s) boundary f = PG.boundary f . _graph --- | Generates the darts incident to a face, starting with the given dart.+-- | Given a dart d, generates the darts bounding the face that is to+-- the right of the given dart. The darts are reported in order along+-- the face. This means that for internal faces the darts are reported+-- in *clockwise* order along the boundary, whereas for the outer face+-- the darts are reported in counter clockwise order. --+-- running time: \(O(k)\), where \(k\) is the number of darts reported ----- \(O(k)\), where \(k\) is the number of darts reported+-- >>> boundary' (dart 2 "+1") smallG -- around face B+-- [Dart (Arc 2) +1,Dart (Arc 3) -1,Dart (Arc 1) -1]+-- >>> boundary' (dart 0 "+1") smallG -- around outer face+-- [Dart (Arc 0) +1,Dart (Arc 4) -1,Dart (Arc 3) +1,Dart (Arc 2) -1]+-- boundary'   :: Dart s -> PlaneGraph s v e f r -> V.Vector (Dart s) boundary' d = PG.boundary' d . _graph @@ -519,8 +640,24 @@ -- | The vertices bounding this face, for internal faces in clockwise order, for -- the outer face in counter clockwise order. ----- -- running time: \(O(k)\), where \(k\) is the output size.+--+-- >>> boundaryVertices (FaceId $ VertexId 2) smallG -- around B+-- [VertexId 0,VertexId 3,VertexId 1]+-- >>> boundaryVertices (FaceId $ VertexId 0) smallG -- around outerface+-- [VertexId 0,VertexId 2,VertexId 1,VertexId 3]+-- >>> mapM_ print $ boundaryVertices (FaceId $ VertexId 0) myPlaneGraph+-- VertexId 0+-- VertexId 9+-- VertexId 10+-- VertexId 11+-- VertexId 7+-- VertexId 8+-- VertexId 12+-- VertexId 13+-- VertexId 4+-- VertexId 3+-- VertexId 2 boundaryVertices   :: FaceId' s -> PlaneGraph s v e f r                    -> V.Vector (VertexId' s) boundaryVertices f = PG.boundaryVertices f . _graph@@ -529,9 +666,18 @@ -------------------------------------------------------------------------------- -- * Access data ++-- | Lens to access the vertex data+--+-- Note that using the setting part of this lens may be very+-- expensive!!  (O(n)) vertexDataOf   :: VertexId' s -> Lens' (PlaneGraph s v e f r ) (VertexData r v) vertexDataOf v = graph.PG.dataOf v +-- | Get the location of a vertex in the plane graph+--+-- Note that the setting part of this lens may be very expensive!+-- Moreover, use with care (as this may destroy planarity etc.) locationOf   :: VertexId' s -> Lens' (PlaneGraph s v e f r ) (Point 2 r) locationOf v = vertexDataOf v.location @@ -652,6 +798,27 @@ -- (Dart (Arc 2) +1,ClosedLineSegment (Point2 0 0 :+ 0) (Point2 (-1) 4 :+ 3) :+ "0->3") -- (Dart (Arc 4) +1,ClosedLineSegment (Point2 2 2 :+ 1) (Point2 2 0 :+ 2) :+ "1->2") -- (Dart (Arc 3) +1,ClosedLineSegment (Point2 2 2 :+ 1) (Point2 (-1) 4 :+ 3) :+ "1->3")+-- >>> mapM_ print $ edgeSegments myPlaneGraph+-- (Dart (Arc 0) +1,ClosedLineSegment (Point2 0 0 :+ 0) (Point2 8 (-5) :+ 9) :+ ())+-- (Dart (Arc 1) +1,ClosedLineSegment (Point2 0 0 :+ 0) (Point2 8 1 :+ 5) :+ ())+-- (Dart (Arc 2) +1,ClosedLineSegment (Point2 0 0 :+ 0) (Point2 4 4 :+ 1) :+ ())+-- (Dart (Arc 3) +1,ClosedLineSegment (Point2 0 0 :+ 0) (Point2 3 7 :+ 2) :+ ())+-- (Dart (Arc 4) +1,ClosedLineSegment (Point2 4 4 :+ 1) (Point2 8 1 :+ 5) :+ ())+-- (Dart (Arc 15) +1,ClosedLineSegment (Point2 4 4 :+ 1) (Point2 10 4 :+ 12) :+ ())+-- (Dart (Arc 5) +1,ClosedLineSegment (Point2 3 7 :+ 2) (Point2 0 5 :+ 3) :+ ())+-- (Dart (Arc 6) +1,ClosedLineSegment (Point2 0 5 :+ 3) (Point2 3 8 :+ 4) :+ ())+-- (Dart (Arc 18) +1,ClosedLineSegment (Point2 3 8 :+ 4) (Point2 9 6 :+ 13) :+ ())+-- (Dart (Arc 7) +1,ClosedLineSegment (Point2 8 1 :+ 5) (Point2 6 (-1) :+ 6) :+ ())+-- (Dart (Arc 10) +1,ClosedLineSegment (Point2 8 1 :+ 5) (Point2 12 1 :+ 8) :+ ())+-- (Dart (Arc 12) +1,ClosedLineSegment (Point2 6 (-1) :+ 6) (Point2 8 (-5) :+ 9) :+ ())+-- (Dart (Arc 8) +1,ClosedLineSegment (Point2 9 (-1) :+ 7) (Point2 12 1 :+ 8) :+ ())+-- (Dart (Arc 14) +1,ClosedLineSegment (Point2 9 (-1) :+ 7) (Point2 14 (-1) :+ 11) :+ ())+-- (Dart (Arc 9) +1,ClosedLineSegment (Point2 12 1 :+ 8) (Point2 10 4 :+ 12) :+ ())+-- (Dart (Arc 11) +1,ClosedLineSegment (Point2 8 (-5) :+ 9) (Point2 12 (-3) :+ 10) :+ ())+-- (Dart (Arc 13) +1,ClosedLineSegment (Point2 12 (-3) :+ 10) (Point2 14 (-1) :+ 11) :+ ())+-- (Dart (Arc 16) +1,ClosedLineSegment (Point2 10 4 :+ 12) (Point2 9 6 :+ 13) :+ ())+-- (Dart (Arc 17) +1,ClosedLineSegment (Point2 10 4 :+ 12) (Point2 8 5 :+ 14) :+ ())+-- (Dart (Arc 19) +1,ClosedLineSegment (Point2 9 6 :+ 13) (Point2 8 5 :+ 14) :+ ()) edgeSegments    :: PlaneGraph s v e f r -> V.Vector (Dart s, LineSegment 2 v r :+ e) edgeSegments ps = fmap withSegment . edges $ ps   where@@ -671,31 +838,87 @@     seg = let (p,q) = bimap vtxDataToExt vtxDataToExt $ ps^.endPointsOf d           in ClosedLineSegment p q --- | The polygon describing the face++-- | The boundary of the face as a simple polygon. For internal faces+-- the polygon that is reported has its vertices stored in CCW order+-- (as expected). ----- runningtime: \(O(k)\), where \(k\) is the size of the face.+-- pre: FaceId refers to an internal face. --+-- For the other face this prodcuces a polygon in CW order (this may+-- lead to unexpected results.) ---rawFaceBoundary      :: FaceId' s -> PlaneGraph s v e f r-                    -> SimplePolygon v r :+ f-rawFaceBoundary i ps = pg :+ (ps^.dataOf i)+-- runningtime: \(O(k)\), where \(k\) is the size of the face.+faceBoundary      :: FaceId' s -> PlaneGraph s v e f r -> SimplePolygon v r :+ f+faceBoundary i ps = pg :+ (ps^.dataOf i)   where-    pg = unsafeFromPoints . F.toList . fmap (\j -> ps^.graph.dataOf j.to vtxDataToExt)+    pg = unsafeFromVector . V.reverse . fmap (\j -> ps^.graph.dataOf j.to vtxDataToExt)        . boundaryVertices i $ ps+    -- polygons are stored in CCW order, the boundaryVertices of+    -- internal faces are reported in CW order we reverse them. --- | Alias for rawFace Boundary+--------------------------------------------------------------------------------++-- | The boundary of the face as a simple polygon. For internal faces+-- the polygon that is reported has its vertices stored in CCW order+-- (as expected). --+-- pre: FaceId refers to an internal face.+--+-- For the other face this prodcuces a polygon in CW order (this may+-- lead to unexpected results.)+-- -- runningtime: \(O(k)\), where \(k\) is the size of the face.-rawFacePolygon :: FaceId' s -> PlaneGraph s v e f r -> SimplePolygon v r :+ f-rawFacePolygon = rawFaceBoundary+internalFacePolygon :: FaceId' s -> PlaneGraph s v e f r -> SimplePolygon v r :+ f+internalFacePolygon = faceBoundary --- | Lists all faces of the plane graph.-rawFacePolygons    :: PlaneGraph s v e f r-                   -> V.Vector (FaceId' s, SimplePolygon v r :+ f)-rawFacePolygons ps = fmap (\i -> (i,rawFacePolygon i ps)) . faces' $ ps+-- | Given the outerFaceId and the graph, construct a sufficiently+-- large rectangular multipolygon ith a hole containing the boundary+-- of the outer face.+outerFacePolygon      :: (Num r, Ord r)+                       => FaceId' s -> PlaneGraph s v e f r -> MultiPolygon (Maybe v) r :+ f+outerFacePolygon i pg =+    outerFacePolygon' i outer pg & core %~ first (either (const Nothing) Just)+  where+    outer = rectToPolygon . grow 1 . boundingBox $ pg+    rectToPolygon = unsafeFromPoints . reverse . F.toList . corners +-- | Given the outerface id, and a sufficiently large outer boundary,+-- draw the outerface as a polygon with a hole.+outerFacePolygon'            :: FaceId' s -> SimplePolygon v' r+                             -> PlaneGraph s v e f r -> MultiPolygon (Either v' v) r :+ f+outerFacePolygon' i outer pg = MultiPolygon (first Left outer) [hole] :+ pg^.dataOf i+  where+    hole = reverseOuterBoundary . first Right . view core $ faceBoundary i pg+    -- if we call faceBoundary on the outerface we get a polygon in+    -- the wrong orientation. So reverse it.+ -------------------------------------------------------------------------------- +-- | Given the outerFace Id, construct polygons for all faces. We+-- construct a polygon with a hole for the outer face.+--+facePolygons      :: (Num r, Ord r) => FaceId' s -> PlaneGraph s v e f r+                  -> ( (FaceId' s, MultiPolygon (Maybe v) r :+ f)+                     , V.Vector (FaceId' s, SimplePolygon v r :+ f)+                     )+facePolygons i ps = ((i, outerFacePolygon i ps), facePolygons' i ps)++-- | Given the outerFace Id, lists all internal faces of the plane+-- graph with their boundaries.+facePolygons'      :: FaceId' s -> PlaneGraph s v e f r+                   ->  V.Vector (FaceId' s, SimplePolygon v r :+ f)+facePolygons' i ps = fmap (\j -> (j,internalFacePolygon j ps)) . V.filter (/= i) . faces' $ ps+++-- | lists all internal faces of the plane graph with their+-- boundaries.+internalFacePolygons    :: (Ord r, Fractional r)+                        => PlaneGraph s v e f r ->  V.Vector (FaceId' s, SimplePolygon v r :+ f)+internalFacePolygons pg = facePolygons' (outerFaceId pg) pg++--------------------------------------------------------------------------------+ -- | Labels the edges of a plane graph with their distances, as specified by -- the distance function. withEdgeDistances     :: (Point 2 r ->  Point 2 r -> a)@@ -703,3 +926,7 @@ withEdgeDistances f g = g&graph.PG.dartData %~ fmap (\(d,x) -> (d,len d :+ x))   where     len d = uncurry f . over both (^.location) $ endPointData d g++++--------------------------------------------------------------------------------

@@ -21,13 +21,18 @@ import qualified Data.PlanarGraph.AdjRep as PGA import qualified Data.PlanarGraph.IO as PGIO import           Data.PlaneGraph.Core-import           Data.PlaneGraph.AdjRep (Face,Vtx(Vtx),Gr(Gr))+import           Data.PlaneGraph.AdjRep import           Data.Proxy import qualified Data.Vector as V import qualified Data.Vector.Mutable as MV import           Data.Yaml (ParseException) import           Data.Yaml.Util ++import Data.RealNumber.Rational+-- import Data.PlanarGraph.Dart+-- import Data.PlaneGraph.AdjRep+ --------------------------------------------------------------------------------  -- $setup@@ -68,10 +73,10 @@ -- * Reading and Writing the Plane Graph  -- | Reads a plane graph from a bytestring-readPlaneGraph   :: (FromJSON v, FromJSON e, FromJSON f, FromJSON r)-                 => proxy s -> B.ByteString-                 -> Either ParseException (PlaneGraph s v e f r)-readPlaneGraph _ = decodeYaml+readPlaneGraph :: forall s v e f r. (FromJSON v, FromJSON e, FromJSON f, FromJSON r)+               => B.ByteString+               -> Either ParseException (PlaneGraph s v e f r)+readPlaneGraph = decodeYaml  -- | Writes a plane graph to a bytestring writePlaneGraph :: (ToJSON v, ToJSON e, ToJSON f, ToJSON r)@@ -132,3 +137,71 @@                            -- hence, no need to pick a secondary comparison  --------------------------------------------------------------------------------++-- smallG = fromAdjRep (Proxy :: Proxy ()) small+--   where+--     small :: Gr (Vtx Int String Int) (Face String)+--     small = Gr [ Vtx 0 (Point2 0 0) [ (2,"0->2")+--                                     , (1,"0->1")+--                                     , (3,"0->3")+--                                     ] 0+--                , Vtx 1 (Point2 2 2) [ (0,"1->0")+--                                     , (2,"1->2")+--                                     , (3,"1->3")+--                                     ] 1+--                , Vtx 2 (Point2 2 0) [ (0,"2->0")+--                                     , (1,"2->1")+--                                     ] 2+--                , Vtx 3 (Point2 (-1) 4) [ (0,"3->0")+--                                        , (1,"3->1")+--                                        ] 3+--                ]+--                [ Face (2,1) "OuterFace"+--                , Face (0,1) "A"+--                , Face (1,0) "B"+--                ]++-- dart i s = Dart (Arc i) (read s)++data MyWorld++-- ![myGraph](docs/Data/PlaneGraph/planegraph.png)+myPlaneGraph :: PlaneGraph MyWorld Int () String (RealNumber 5)+myPlaneGraph = fromAdjRep (Proxy @MyWorld) myPlaneGraphAdjrep++myPlaneGraphAdjrep :: Gr (Vtx Int () (RealNumber 5)) (Face String)+myPlaneGraphAdjrep = Gr [ vtx 0 (Point2 0   0   ) [e 9, e 5, e 1, e 2]+                        , vtx 1 (Point2 4   4   ) [e 0, e 5, e 12]+                        , vtx 2 (Point2 3   7   ) [e 0, e 3]+                        , vtx 3 (Point2 0   5   ) [e 4, e 2]+                        , vtx 4 (Point2 3   8   ) [e 3, e 13]+                        , vtx 5 (Point2 8   1   ) [e 0, e 6, e 8, e 1]+                        , vtx 6 (Point2 6   (-1)) [e 5, e 9]+                        , vtx 7 (Point2 9   (-1)) [e 8, e 11]+                        , vtx 8 (Point2 12  1   ) [e 7, e 12, e 5]+                        , vtx 9 (Point2 8   (-5)) [e 0, e 10, e 6]+                        , vtx 10 (Point2 12 (-3)) [e 9, e 11]+                        , vtx 11 (Point2 14 (-1)) [e 10, e 7]+                        , vtx 12 (Point2 10 4   ) [e 1, e 8, e 13, e 14]+                        , vtx 13 (Point2 9  6   ) [e 4, e 14, e 12]+                        , vtx 14 (Point2 8  5   ) [e 13, e 12]+                        ]+                        [ Face (0,9) "OuterFace"+                        , Face (0,5) "A"+                        , Face (0,1) "B"+                        , Face (0,2) "C"+                        , Face (14,13) "D"+                        , Face (1,12) "E"+                        , Face (5,8) "F"+                        ]+  where+    e i = (i,())+    vtx i p es = Vtx i p es i+++++-- myPlaneGraph' :: IO (PlaneGraph MyWorld () () () (RealNumber 5))+-- myPlaneGraph' = let err x  = error $ show x+--                 in either err id . readPlaneGraph+--                 <$> B.readFile "docs/Data/PlaneGraph/myPlaneGraph.yaml"