hgeometry 0.6.0.0 → 0.7.0.0

raw patch · 116 files changed

+8425/−1073 lines, 116 filesdep +aesondep +approximate-equalitydep +cairodep −Framesdep ~QuickCheckdep ~basedep ~bytestringnew-component:exe:hgeometry-viewerPVP ok

version bump matches the API change (PVP)

Dependencies added: aeson, approximate-equality, cairo, cairo-canvas, colour, dlist, fingertree, gi-cairo, gi-gdk, gi-gtk, haskell-gi-base, quickcheck-instances, reactive-banana, reactive-banana-gi-gtk, reflection, transformers, yaml

Dependencies removed: Frames

Dependency ranges changed: QuickCheck, base, bytestring, containers, criterion, fixed-vector, hexpat, lens, linear, optparse-applicative, semigroupoids, semigroups, singletons, text, vector, vinyl

API changes (from Hackage documentation)

Files

• README.md +1/−1
• benchmark/Algorithms/Geometry/ConvexHull/Bench.hs +72/−0
• benchmark/Algorithms/Geometry/ConvexHull/GrahamFam.hs +103/−0
• benchmark/Algorithms/Geometry/ConvexHull/GrahamFam6.hs +103/−0
• benchmark/Algorithms/Geometry/ConvexHull/GrahamFamPeano.hs +107/−0
• benchmark/Algorithms/Geometry/ConvexHull/GrahamFixed.hs +104/−0
• benchmark/Algorithms/Geometry/ConvexHull/GrahamV2.hs +95/−0
• benchmark/Benchmark/Util.hs +7/−0
• benchmark/Benchmarks.hs +3/−66
• benchmark/Data/Geometry/IntervalTreeBench.hs +79/−0
• benchmark/Data/Geometry/Vector/VectorFamily6.hs +258/−0
• benchmark/WSPDBench.hs +26/−0
• changelog.md +60/−0
• examples/BAPC2014/Armybase.lhs +2/−2
• examples/Demo/Delaunay.hs +0/−2
• examples/Demo/ExpectedPairwiseDistance.hs +4/−4
• examples/Demo/GPXParser.hs +1/−1
• examples/Demo/TriangulateWorld.hs +106/−0
• examples/Main.hs +4/−0
• hgeometry.cabal +174/−30
• interactive/RenderCanvas.hs +177/−0
• interactive/RenderUtil.hs +81/−0
• interactive/Viewer.hs +172/−0
• src/Algorithms/Geometry/ConvexHull/DivideAndConqueror.hs +2/−1
• src/Algorithms/Geometry/ConvexHull/GrahamScan.hs +2/−2
• src/Algorithms/Geometry/DelaunayTriangulation/Types.hs +12/−39
• src/Algorithms/Geometry/EuclideanMST/EuclideanMST.hs +2/−2
• src/Algorithms/Geometry/LineSegmentIntersection.hs +16/−0
• src/Algorithms/Geometry/LineSegmentIntersection/BentleyOttmann.hs +60/−46
• src/Algorithms/Geometry/LineSegmentIntersection/Naive.hs +1/−1
• src/Algorithms/Geometry/LineSegmentIntersection/Types.hs +8/−0
• src/Algorithms/Geometry/LowerEnvelope/DualCH.hs +49/−0
• src/Algorithms/Geometry/PolygonTriangulation/MakeMonotone.hs +310/−0
• src/Algorithms/Geometry/PolygonTriangulation/Triangulate.hs +70/−0
• src/Algorithms/Geometry/PolygonTriangulation/TriangulateMonotone.hs +178/−0
• src/Algorithms/Geometry/PolygonTriangulation/Types.hs +138/−0
• src/Algorithms/Geometry/Sweep.hs +151/−0
• src/Algorithms/Geometry/WellSeparatedPairDecomposition/Types.hs +1/−0
• src/Algorithms/Geometry/WellSeparatedPairDecomposition/WSPD.hs +10/−13
• src/Control/Monad/State/Persistent.hs +1/−1
• src/Data/BalBST.hs +27/−4
• src/Data/CircularSeq.hs +6/−1
• src/Data/Ext.hs +23/−2
• src/Data/Geometry.hs +6/−7
• src/Data/Geometry/Ball.hs +11/−2
• src/Data/Geometry/Boundary.hs +2/−2
• src/Data/Geometry/Box.hs +7/−0
• src/Data/Geometry/Box/Internal.hs +21/−6
• src/Data/Geometry/HalfLine.hs +3/−2
• src/Data/Geometry/Interval.hs +2/−2
• src/Data/Geometry/Ipe.hs +6/−0
• src/Data/Geometry/Ipe/Attributes.hs +27/−38
• src/Data/Geometry/Ipe/FromIpe.hs +109/−12
• src/Data/Geometry/Ipe/IpeOut.hs +46/−18
• src/Data/Geometry/Ipe/ParserPrimitives.hs +7/−1
• src/Data/Geometry/Ipe/PathParser.hs +31/−24
• src/Data/Geometry/Ipe/Reader.hs +42/−76
• src/Data/Geometry/Ipe/Types.hs +25/−14
• src/Data/Geometry/Ipe/Writer.hs +39/−23
• src/Data/Geometry/KDTree.hs +7/−6
• src/Data/Geometry/Line.hs +66/−17
• src/Data/Geometry/Line/Internal.hs +32/−5
• src/Data/Geometry/LineSegment.hs +27/−14
• src/Data/Geometry/PlanarSubdivision.hs +83/−108
• src/Data/Geometry/PlanarSubdivision/Basic.hs +598/−0
• src/Data/Geometry/PlanarSubdivision/Draw.hs +17/−0
• src/Data/Geometry/Point.hs +57/−29
• src/Data/Geometry/PolyLine.hs +6/−1
• src/Data/Geometry/Polygon.hs +111/−11
• src/Data/Geometry/Polygon/Convex.hs +20/−13
• src/Data/Geometry/Properties.hs +16/−12
• src/Data/Geometry/SegmentTree/Generic.hs +0/−2
• src/Data/Geometry/Slab.hs +17/−2
• src/Data/Geometry/SubLine.hs +37/−10
• src/Data/Geometry/Transformation.hs +32/−27
• src/Data/Geometry/Triangle.hs +17/−12
• src/Data/Geometry/Vector.hs +37/−25
• src/Data/Geometry/Vector/VectorFamily.hs +156/−0
• src/Data/Geometry/Vector/VectorFamilyPeano.hs +291/−0
• src/Data/Geometry/Vector/VectorFixed.hs +63/−43
• src/Data/OrdSeq.hs +190/−0
• src/Data/Permutation.hs +4/−0
• src/Data/PlanarGraph.hs +371/−114
• src/Data/PlaneGraph.hs +532/−21
• src/Data/PlaneGraph/Draw.hs +29/−0
• src/Data/Range.hs +11/−3
• src/Data/Seq.hs +9/−8
• src/Data/Sequence/Util.hs +21/−2
• src/Data/SlowSeq.hs +206/−0
• src/Data/Util.hs +0/−1
• src/Data/Yaml/Util.hs +22/−0
• src/Test/QuickCheck/HGeometryInstances.hs +94/−0
• test/Algorithms/Geometry/ConvexHull/ConvexHullSpec.hs +41/−0
• test/Algorithms/Geometry/LineSegmentIntersection/BentleyOttmannSpec.hs +52/−1
• test/Algorithms/Geometry/LineSegmentIntersection/selfIntersections.ipe +313/−0
• test/Algorithms/Geometry/LowerEnvelope/LowerEnvSpec.hs +104/−0
• test/Algorithms/Geometry/LowerEnvelope/manual.ipe +299/−0
• test/Algorithms/Geometry/PolygonTriangulation/MakeMonotoneSpec.hs +58/−0
• test/Algorithms/Geometry/PolygonTriangulation/TriangulateMonotoneSpec.hs +73/−0
• test/Algorithms/Geometry/PolygonTriangulation/monotone.ipe +364/−0
• test/Algorithms/Geometry/PolygonTriangulation/simplepolygon6.ipe +297/−0
• test/Algorithms/Geometry/WellSeparatedPairDecomposition/WSPDSpec.hs +8/−7
• test/Data/EdgeOracleSpec.hs +2/−2
• test/Data/Geometry/IntervalSpec.hs +17/−12
• test/Data/Geometry/Ipe/ReaderSpec.hs +1/−1
• test/Data/Geometry/KDTreeSpec.hs +2/−2
• test/Data/Geometry/LineSegmentSpec.hs +16/−0
• test/Data/Geometry/LineSpec.hs +48/−0
• test/Data/Geometry/PlanarSubdivisionSpec.hs +221/−0
• test/Data/Geometry/Polygon/Convex/ConvexSpec.hs +3/−2
• test/Data/Geometry/SubLineSpec.hs +50/−0
• test/Data/OrdSeqSpec.hs +58/−0
• test/Data/PlanarGraphSpec.hs +6/−6
• test/Data/RangeSpec.hs +14/−17
• test/QuickCheck/Instances.hs +0/−91
• test/Util.hs +9/−1

@@ -47,7 +47,7 @@ * A KD-Tree. The base tree is static.  HGeometry also includes a datastructure/data type for planar graphs. In-particular, it has a `EdgeOracle' data structure, that can be built in $O(n)$+particular, it has a `EdgeOracle` data structure, that can be built in \(O(n)\) time that can test if the graph contains an edge in constant time.  Numeric Types

@@ -0,0 +1,72 @@+module Algorithms.Geometry.ConvexHull.Bench where++import qualified Algorithms.Geometry.ConvexHull.DivideAndConqueror as DivideAndConqueror+import qualified Algorithms.Geometry.ConvexHull.GrahamScan as GrahamScan++-- | copies of the convex hull algo with different point types+import qualified Algorithms.Geometry.ConvexHull.GrahamV2   as GV+import qualified Algorithms.Geometry.ConvexHull.GrahamFam  as GFam+import qualified Algorithms.Geometry.ConvexHull.GrahamFamPeano  as GPeano+import qualified Algorithms.Geometry.ConvexHull.GrahamFam6  as GFam6+import qualified Algorithms.Geometry.ConvexHull.GrahamFixed as GFix+++import           Benchmark.Util+import           Control.DeepSeq+import           Criterion.Main+import           Criterion.Types+import           Data.Ext+import           Data.Geometry.Point+import           Data.List.NonEmpty (NonEmpty(..))+import qualified Data.List.NonEmpty as NonEmpty+import           Data.Proxy+import           Test.QuickCheck+import           Test.QuickCheck.HGeometryInstances ()++--------------------------------------------------------------------------------++main :: IO ()+main = defaultMainWith cfg [ benchmark ]+  where+    cfg = defaultConfig { reportFile = Just "bench.html" }++benchmark :: Benchmark+benchmark = bgroup "convexHullBench"+    [ env (genPts (Proxy :: Proxy Int) 10000) benchBuild+    ]++--------------------------------------------------------------------------------++genPts     :: (Ord r, Arbitrary r) => proxy r -> Int -> IO (NonEmpty (Point 2 r :+ ()))+genPts _ n = generate (NonEmpty.fromList <$> vectorOf n arbitrary)++-- | Benchmark building the convexHull+benchBuild    :: (Ord r, Num r, NFData r) => NonEmpty (Point 2 r :+ ()) -> Benchmark+benchBuild ps = bgroup "build" [ bgroup (show n) (build $ take' n ps)+                               | n <- sizes' ps+                               ]+  where+    take' n = NonEmpty.fromList . NonEmpty.take n+    sizes' _ = [2000]++    build pts = [ bench "sort"                 $ nf NonEmpty.sort pts+                , bench "sort_Linear.V2"       $ nf NonEmpty.sort ptsV2+                , bench "sort_FamPeano"        $ nf NonEmpty.sort ptsFamPeano+                , bench "sort_Family"          $ nf NonEmpty.sort ptsFam+                , bench "sort_Family6"         $ nf NonEmpty.sort ptsFam6+                , bench "sort_Fixed"           $ nf NonEmpty.sort ptsFix++                , bench "grahamScan"           $ nf GrahamScan.convexHull pts+                , bench "grahamScan_Linear.V2" $ nf GV.convexHull         ptsV2+                , bench "grahamScan_FamPeano"  $ nf GPeano.convexHull     ptsFamPeano+                , bench "grahamScan_Family"    $ nf GFam.convexHull       ptsFam+                , bench "grahamScan_Fixed"     $ nf GFix.convexHull       ptsFix++                , bench "Div&Conq"             $ nf DivideAndConqueror.convexHull pts+                ]+      where+        ptsV2       = fmap (GV.fromP) pts+        ptsFamPeano = fmap (GPeano.fromP) pts+        ptsFam      = fmap (GFam.fromP) pts+        ptsFam6     = fmap (GFam6.fromP) pts+        ptsFix      = fmap (GFix.fromP) pts

@@ -0,0 +1,103 @@+{-# LANGUAGE UndecidableInstances #-}+module Algorithms.Geometry.ConvexHull.GrahamFam( convexHull+                                               , upperHull+                                               , lowerHull, fromP+                                               ) where++import           Control.DeepSeq+import           Control.Lens ((^.))+import           Data.Ext+import           Data.Geometry.Point+import qualified Data.Geometry.Vector.VectorFamily as VF+import           Data.List.NonEmpty (NonEmpty(..))+import qualified Data.List.NonEmpty as NonEmpty+import           Data.Monoid+import           GHC.TypeLits+++newtype MyPoint d r = MyPoint (VF.Vector d r)++deriving instance (VF.Arity d, Eq r)  => Eq (MyPoint d r)+deriving instance (VF.Arity d, Ord r) => Ord (MyPoint d r)+deriving instance (VF.Arity d, Show r) => Show (MyPoint d r)+deriving instance (NFData (VF.Vector d r)) => NFData (MyPoint d r)++pattern MyPoint2 x y = MyPoint (VF.Vector2 x y)+++-- instance (NFData r, Arity d) => NFData (MyPoint d r)  where+--   rnf (MyPoint x y) = rnf (x,y)+--   rnf (MyP p)       = rnf p++toP                    :: MyPoint 2 r :+ e -> Point 2 r :+ e+toP (MyPoint2 x y :+ e) = Point2 x y :+ e++fromP                   :: Point 2 r :+ e -> MyPoint 2 r :+ e+fromP (Point2 x y :+ e) = MyPoint2 x y :+ e+++subt :: Num r => MyPoint 2 r -> MyPoint 2 r -> MyPoint 2 r+(MyPoint2 x y) `subt` (MyPoint2 a b) = MyPoint2 (x-a) (y-b)++newtype ConvexPolygon p r = ConvexPolygon [Point 2 r :+ p] deriving (Show,Eq,NFData)++-- | \(O(n \log n)\) time ConvexHull using Graham-Scan. The resulting polygon is+-- given in clockwise order.+convexHull            :: (Ord r, Num r)+                      => NonEmpty (MyPoint 2 r :+ p) -> ConvexPolygon p r+convexHull (p :| []) = ConvexPolygon $ [toP p]+convexHull ps        = let ps' = NonEmpty.toList . NonEmpty.sortBy incXdecY $ ps+                           uh  = NonEmpty.tail . hull' $         ps'+                           lh  = NonEmpty.tail . hull' $ reverse ps'+                       in ConvexPolygon . map toP . reverse $ lh ++ uh++upperHull  :: (Ord r, Num r) => NonEmpty (MyPoint 2 r :+ p) -> NonEmpty (MyPoint 2 r :+ p)+upperHull = hull id+++lowerHull :: (Ord r, Num r) => NonEmpty (MyPoint 2 r :+ p) -> NonEmpty (MyPoint 2 r :+ p)+lowerHull = hull reverse+++-- | Helper function so that that can compute both the upper or the lower hull, depending+-- on the function f+hull               :: (Ord r, Num r)+                   => ([MyPoint 2 r :+ p] -> [MyPoint 2 r :+ p])+                   -> NonEmpty (MyPoint 2 r :+ p) -> NonEmpty (MyPoint 2 r :+ p)+hull _ h@(_ :| []) = h+hull f pts         = hull' .  f+                   . NonEmpty.toList . NonEmpty.sortBy incXdecY $ pts++incXdecY  :: Ord r => (MyPoint 2 r) :+ p -> (MyPoint 2 r) :+ q -> Ordering+incXdecY (MyPoint2 px py :+ _) (MyPoint2 qx qy :+ _) =+  compare px qx <> compare qy py+++-- | Precondition: The list of input points is sorted+hull'          :: (Ord r, Num r) => [MyPoint 2 r :+ p] -> NonEmpty (MyPoint 2 r :+ p)+hull' (a:b:ps) = NonEmpty.fromList $ hull'' [b,a] ps+  where+    hull'' h []      = h+    hull'' h (p:ps') = hull'' (cleanMiddle (p:h)) ps'++    cleanMiddle h@[_,_]                         = h+    cleanMiddle h@(z:y:x:rest)+      | rightTurn (x^.core) (y^.core) (z^.core) = h+      | otherwise                               = cleanMiddle (z:x:rest)+    cleanMiddle _                               = error "cleanMiddle: too few points"++rightTurn       :: (Ord r, Num r) => MyPoint 2 r -> MyPoint 2 r -> MyPoint 2 r -> Bool+rightTurn a b c = ccwP a b c == CW++++ccwP :: (Ord r, Num r) => MyPoint 2 r -> MyPoint 2 r -> MyPoint 2 r -> CCW+ccwP p q r = case z `compare` 0 of+              LT -> CW+              GT -> CCW+              EQ -> CoLinear+     where++       MyPoint2 ux uy = q `subt` p+       MyPoint2 vx vy = r `subt` p+       z              = ux * vy - uy * vx

@@ -0,0 +1,103 @@+{-# LANGUAGE UndecidableInstances #-}+module Algorithms.Geometry.ConvexHull.GrahamFam6( convexHull+                                                , upperHull+                                                , lowerHull, fromP+                                                ) where++import           Control.DeepSeq+import           Control.Lens ((^.))+import           Data.Ext+import           Data.Geometry.Point+import qualified Data.Geometry.Vector.VectorFamily6 as VF+import           Data.List.NonEmpty (NonEmpty(..))+import qualified Data.List.NonEmpty as NonEmpty+import           Data.Monoid+import           GHC.TypeLits+++newtype MyPoint d r = MyPoint (VF.Vector d r)++deriving instance (VF.Arity d, Eq r)  => Eq (MyPoint d r)+deriving instance (VF.Arity d, Ord r) => Ord (MyPoint d r)+deriving instance (VF.Arity d, Show r) => Show (MyPoint d r)+deriving instance (NFData (VF.Vector d r)) => NFData (MyPoint d r)++pattern MyPoint2 x y = MyPoint (VF.Vector2 x y)+++-- instance (NFData r, Arity d) => NFData (MyPoint d r)  where+--   rnf (MyPoint x y) = rnf (x,y)+--   rnf (MyP p)       = rnf p++toP                    :: MyPoint 2 r :+ e -> Point 2 r :+ e+toP (MyPoint2 x y :+ e) = Point2 x y :+ e++fromP                   :: Point 2 r :+ e -> MyPoint 2 r :+ e+fromP (Point2 x y :+ e) = MyPoint2 x y :+ e+++subt :: Num r => MyPoint 2 r -> MyPoint 2 r -> MyPoint 2 r+(MyPoint2 x y) `subt` (MyPoint2 a b) = MyPoint2 (x-a) (y-b)++newtype ConvexPolygon p r = ConvexPolygon [Point 2 r :+ p] deriving (Show,Eq,NFData)++-- | \(O(n \log n)\) time ConvexHull using Graham-Scan. The resulting polygon is+-- given in clockwise order.+convexHull            :: (Ord r, Num r)+                      => NonEmpty (MyPoint 2 r :+ p) -> ConvexPolygon p r+convexHull (p :| []) = ConvexPolygon $ [toP p]+convexHull ps        = let ps' = NonEmpty.toList . NonEmpty.sortBy incXdecY $ ps+                           uh  = NonEmpty.tail . hull' $         ps'+                           lh  = NonEmpty.tail . hull' $ reverse ps'+                       in ConvexPolygon . map toP . reverse $ lh ++ uh++upperHull  :: (Ord r, Num r) => NonEmpty (MyPoint 2 r :+ p) -> NonEmpty (MyPoint 2 r :+ p)+upperHull = hull id+++lowerHull :: (Ord r, Num r) => NonEmpty (MyPoint 2 r :+ p) -> NonEmpty (MyPoint 2 r :+ p)+lowerHull = hull reverse+++-- | Helper function so that that can compute both the upper or the lower hull, depending+-- on the function f+hull               :: (Ord r, Num r)+                   => ([MyPoint 2 r :+ p] -> [MyPoint 2 r :+ p])+                   -> NonEmpty (MyPoint 2 r :+ p) -> NonEmpty (MyPoint 2 r :+ p)+hull _ h@(_ :| []) = h+hull f pts         = hull' .  f+                   . NonEmpty.toList . NonEmpty.sortBy incXdecY $ pts++incXdecY  :: Ord r => (MyPoint 2 r) :+ p -> (MyPoint 2 r) :+ q -> Ordering+incXdecY (MyPoint2 px py :+ _) (MyPoint2 qx qy :+ _) =+  compare px qx <> compare qy py+++-- | Precondition: The list of input points is sorted+hull'          :: (Ord r, Num r) => [MyPoint 2 r :+ p] -> NonEmpty (MyPoint 2 r :+ p)+hull' (a:b:ps) = NonEmpty.fromList $ hull'' [b,a] ps+  where+    hull'' h []      = h+    hull'' h (p:ps') = hull'' (cleanMiddle (p:h)) ps'++    cleanMiddle h@[_,_]                         = h+    cleanMiddle h@(z:y:x:rest)+      | rightTurn (x^.core) (y^.core) (z^.core) = h+      | otherwise                               = cleanMiddle (z:x:rest)+    cleanMiddle _                               = error "cleanMiddle: too few points"++rightTurn       :: (Ord r, Num r) => MyPoint 2 r -> MyPoint 2 r -> MyPoint 2 r -> Bool+rightTurn a b c = ccwP a b c == CW++++ccwP :: (Ord r, Num r) => MyPoint 2 r -> MyPoint 2 r -> MyPoint 2 r -> CCW+ccwP p q r = case z `compare` 0 of+              LT -> CW+              GT -> CCW+              EQ -> CoLinear+     where++       MyPoint2 ux uy = q `subt` p+       MyPoint2 vx vy = r `subt` p+       z              = ux * vy - uy * vx

@@ -0,0 +1,107 @@+{-# LANGUAGE UndecidableInstances #-}+module Algorithms.Geometry.ConvexHull.GrahamFamPeano( convexHull+                                                    , upperHull+                                                    , lowerHull, fromP+                                                    ) where++import           Control.DeepSeq+import           Control.Lens ((^.))+import           Data.Ext+import           Data.Geometry.Point+import qualified Data.Vector.Fixed.Cont as V+import qualified Data.Geometry.Vector.VectorFamilyPeano as VF+import           Data.List.NonEmpty (NonEmpty(..))+import qualified Data.List.NonEmpty as NonEmpty+import           Data.Monoid+import           GHC.TypeLits+import qualified Linear.V2 as V2+++newtype MyPoint d r = MyPoint (VF.VectorFamily d r)++deriving instance (VF.ImplicitArity d, Eq r)  => Eq (MyPoint d r)+deriving instance (VF.ImplicitArity d, Ord r) => Ord (MyPoint d r)+deriving instance (VF.ImplicitArity d, Show r) => Show (MyPoint d r)+deriving instance (NFData (VF.VectorFamily d r)) => NFData (MyPoint d r)++pattern Vector2 x y = VF.VectorFamily (V2.V2 x y)++pattern MyPoint2 x y = MyPoint (Vector2 x y)+++-- instance (NFData r, Arity d) => NFData (MyPoint d r)  where+--   rnf (MyPoint x y) = rnf (x,y)+--   rnf (MyP p)       = rnf p++toP                    :: MyPoint VF.Two r :+ e -> Point 2 r :+ e+toP (MyPoint2 x y :+ e) = Point2 x y :+ e++fromP                   :: Point 2 r :+ e -> MyPoint VF.Two r :+ e+fromP (Point2 x y :+ e) = MyPoint2 x y :+ e+++subt :: Num r => MyPoint VF.Two r -> MyPoint VF.Two r -> MyPoint VF.Two r+(MyPoint2 x y) `subt` (MyPoint2 a b) = MyPoint2 (x-a) (y-b)++newtype ConvexPolygon p r = ConvexPolygon [Point 2 r :+ p] deriving (Show,Eq,NFData)++-- | \(O(n \log n)\) time ConvexHull using Graham-Scan. The resulting polygon is+-- given in clockwise order.+convexHull            :: (Ord r, Num r)+                      => NonEmpty (MyPoint VF.Two r :+ p) -> ConvexPolygon p r+convexHull (p :| []) = ConvexPolygon $ [toP p]+convexHull ps        = let ps' = NonEmpty.toList . NonEmpty.sortBy incXdecY $ ps+                           uh  = NonEmpty.tail . hull' $         ps'+                           lh  = NonEmpty.tail . hull' $ reverse ps'+                       in ConvexPolygon . map toP . reverse $ lh ++ uh++upperHull  :: (Ord r, Num r) => NonEmpty (MyPoint VF.Two r :+ p) -> NonEmpty (MyPoint VF.Two r :+ p)+upperHull = hull id+++lowerHull :: (Ord r, Num r) => NonEmpty (MyPoint VF.Two r :+ p) -> NonEmpty (MyPoint VF.Two r :+ p)+lowerHull = hull reverse+++-- | Helper function so that that can compute both the upper or the lower hull, depending+-- on the function f+hull               :: (Ord r, Num r)+                   => ([MyPoint VF.Two r :+ p] -> [MyPoint VF.Two r :+ p])+                   -> NonEmpty (MyPoint VF.Two r :+ p) -> NonEmpty (MyPoint VF.Two r :+ p)+hull _ h@(_ :| []) = h+hull f pts         = hull' .  f+                   . NonEmpty.toList . NonEmpty.sortBy incXdecY $ pts++incXdecY  :: Ord r => (MyPoint VF.Two r) :+ p -> (MyPoint VF.Two r) :+ q -> Ordering+incXdecY (MyPoint2 px py :+ _) (MyPoint2 qx qy :+ _) =+  compare px qx <> compare qy py+++-- | Precondition: The list of input points is sorted+hull'          :: (Ord r, Num r) => [MyPoint VF.Two r :+ p] -> NonEmpty (MyPoint VF.Two r :+ p)+hull' (a:b:ps) = NonEmpty.fromList $ hull'' [b,a] ps+  where+    hull'' h []      = h+    hull'' h (p:ps') = hull'' (cleanMiddle (p:h)) ps'++    cleanMiddle h@[_,_]                         = h+    cleanMiddle h@(z:y:x:rest)+      | rightTurn (x^.core) (y^.core) (z^.core) = h+      | otherwise                               = cleanMiddle (z:x:rest)+    cleanMiddle _                               = error "cleanMiddle: too few points"++rightTurn       :: (Ord r, Num r) => MyPoint VF.Two r -> MyPoint VF.Two r -> MyPoint VF.Two r -> Bool+rightTurn a b c = ccwP a b c == CW++++ccwP :: (Ord r, Num r) => MyPoint VF.Two r -> MyPoint VF.Two r -> MyPoint VF.Two r -> CCW+ccwP p q r = case z `compare` 0 of+              LT -> CW+              GT -> CCW+              EQ -> CoLinear+     where++       MyPoint2 ux uy = q `subt` p+       MyPoint2 vx vy = r `subt` p+       z              = ux * vy - uy * vx

@@ -0,0 +1,104 @@+{-# LANGUAGE UndecidableInstances #-}+module Algorithms.Geometry.ConvexHull.GrahamFixed( convexHull+                                                 , upperHull+                                                 , lowerHull, fromP+                                                 ) where++import           Control.DeepSeq+import           Control.Lens ((^.))+import           Data.Ext+import           Data.Geometry.Point+import           Data.Vector.Fixed (Arity)+import qualified Data.Geometry.Vector.VectorFixed as VF+import           Data.List.NonEmpty (NonEmpty(..))+import qualified Data.List.NonEmpty as NonEmpty+import           Data.Monoid+import           GHC.TypeLits+++newtype MyPoint d r = MyPoint (VF.Vector d r)++deriving instance (Arity d, Eq r)  => Eq (MyPoint d r)+deriving instance (Arity d, Ord r) => Ord (MyPoint d r)+deriving instance (Arity d, Show r) => Show (MyPoint d r)+deriving instance (NFData (VF.Vector d r)) => NFData (MyPoint d r)++pattern MyPoint2 x y = MyPoint (VF.Vector2 x y)+++-- instance (NFData r, Arity d) => NFData (MyPoint d r)  where+--   rnf (MyPoint x y) = rnf (x,y)+--   rnf (MyP p)       = rnf p++toP                    :: MyPoint 2 r :+ e -> Point 2 r :+ e+toP (MyPoint2 x y :+ e) = Point2 x y :+ e++fromP                   :: Point 2 r :+ e -> MyPoint 2 r :+ e+fromP (Point2 x y :+ e) = MyPoint2 x y :+ e+++subt :: Num r => MyPoint 2 r -> MyPoint 2 r -> MyPoint 2 r+(MyPoint2 x y) `subt` (MyPoint2 a b) = MyPoint2 (x-a) (y-b)++newtype ConvexPolygon p r = ConvexPolygon [Point 2 r :+ p] deriving (Show,Eq,NFData)++-- | \(O(n \log n)\) time ConvexHull using Graham-Scan. The resulting polygon is+-- given in clockwise order.+convexHull            :: (Ord r, Num r)+                      => NonEmpty (MyPoint 2 r :+ p) -> ConvexPolygon p r+convexHull (p :| []) = ConvexPolygon $ [toP p]+convexHull ps        = let ps' = NonEmpty.toList . NonEmpty.sortBy incXdecY $ ps+                           uh  = NonEmpty.tail . hull' $         ps'+                           lh  = NonEmpty.tail . hull' $ reverse ps'+                       in ConvexPolygon . map toP . reverse $ lh ++ uh++upperHull  :: (Ord r, Num r) => NonEmpty (MyPoint 2 r :+ p) -> NonEmpty (MyPoint 2 r :+ p)+upperHull = hull id+++lowerHull :: (Ord r, Num r) => NonEmpty (MyPoint 2 r :+ p) -> NonEmpty (MyPoint 2 r :+ p)+lowerHull = hull reverse+++-- | Helper function so that that can compute both the upper or the lower hull, depending+-- on the function f+hull               :: (Ord r, Num r)+                   => ([MyPoint 2 r :+ p] -> [MyPoint 2 r :+ p])+                   -> NonEmpty (MyPoint 2 r :+ p) -> NonEmpty (MyPoint 2 r :+ p)+hull _ h@(_ :| []) = h+hull f pts         = hull' .  f+                   . NonEmpty.toList . NonEmpty.sortBy incXdecY $ pts++incXdecY  :: Ord r => (MyPoint 2 r) :+ p -> (MyPoint 2 r) :+ q -> Ordering+incXdecY (MyPoint2 px py :+ _) (MyPoint2 qx qy :+ _) =+  compare px qx <> compare qy py+++-- | Precondition: The list of input points is sorted+hull'          :: (Ord r, Num r) => [MyPoint 2 r :+ p] -> NonEmpty (MyPoint 2 r :+ p)+hull' (a:b:ps) = NonEmpty.fromList $ hull'' [b,a] ps+  where+    hull'' h []      = h+    hull'' h (p:ps') = hull'' (cleanMiddle (p:h)) ps'++    cleanMiddle h@[_,_]                         = h+    cleanMiddle h@(z:y:x:rest)+      | rightTurn (x^.core) (y^.core) (z^.core) = h+      | otherwise                               = cleanMiddle (z:x:rest)+    cleanMiddle _                               = error "cleanMiddle: too few points"++rightTurn       :: (Ord r, Num r) => MyPoint 2 r -> MyPoint 2 r -> MyPoint 2 r -> Bool+rightTurn a b c = ccwP a b c == CW++++ccwP :: (Ord r, Num r) => MyPoint 2 r -> MyPoint 2 r -> MyPoint 2 r -> CCW+ccwP p q r = case z `compare` 0 of+              LT -> CW+              GT -> CCW+              EQ -> CoLinear+     where++       MyPoint2 ux uy = q `subt` p+       MyPoint2 vx vy = r `subt` p+       z              = ux * vy - uy * vx

@@ -0,0 +1,95 @@+{-# Language DeriveGeneric #-}+module Algorithms.Geometry.ConvexHull.GrahamV2( convexHull+                                              , upperHull+                                              , lowerHull, fromP+                                              ) where+++import           Control.DeepSeq+import           Control.Lens ((^.))+import           Data.Ext+import           Data.Geometry.Point+import           Data.List.NonEmpty (NonEmpty(..))+import qualified Data.List.NonEmpty as NonEmpty+import           Data.Monoid+import           GHC.Generics+import qualified Linear.V2 as V2++++newtype MyPoint r = MKPoint (V2.V2 r) deriving (Show,Eq,Ord,Generic)+-- data MyPoint r = MyPoint !r !r deriving (Show,Eq,Ord,Generic)++pattern MyPoint x y = MKPoint (V2.V2 x y)++instance NFData r => NFData (MyPoint r)+++toP (MyPoint x y :+ e) = Point2 x y :+ e+fromP (Point2 x y :+ e) = MyPoint x y :+ e++(MyPoint x y) `subt` (MyPoint a b) = MyPoint (x-a) (y-b)+++newtype ConvexPolygon p r = ConvexPolygon [Point 2 r :+ p] deriving (Show,Eq,NFData)++-- | \(O(n \log n)\) time ConvexHull using Graham-Scan. The resulting polygon is+-- given in clockwise order.+convexHull            :: (Ord r, Num r)+                      => NonEmpty (MyPoint r :+ p) -> ConvexPolygon p r+convexHull (p :| []) = ConvexPolygon $ [toP p]+convexHull ps        = let ps' = NonEmpty.toList . NonEmpty.sortBy incXdecY $ ps+                           uh  = NonEmpty.tail . hull' $         ps'+                           lh  = NonEmpty.tail . hull' $ reverse ps'+                       in ConvexPolygon . map toP . reverse $ lh ++ uh++upperHull  :: (Ord r, Num r) => NonEmpty (MyPoint r :+ p) -> NonEmpty (MyPoint r :+ p)+upperHull = hull id+++lowerHull :: (Ord r, Num r) => NonEmpty (MyPoint r :+ p) -> NonEmpty (MyPoint r :+ p)+lowerHull = hull reverse+++-- | Helper function so that that can compute both the upper or the lower hull, depending+-- on the function f+hull               :: (Ord r, Num r)+                   => ([MyPoint r :+ p] -> [MyPoint r :+ p])+                   -> NonEmpty (MyPoint r :+ p) -> NonEmpty (MyPoint r :+ p)+hull _ h@(_ :| []) = h+hull f pts         = hull' .  f+                   . NonEmpty.toList . NonEmpty.sortBy incXdecY $ pts++incXdecY  :: Ord r => (MyPoint r) :+ p -> (MyPoint r) :+ q -> Ordering+incXdecY (MyPoint px py :+ _) (MyPoint qx qy :+ _) =+  compare px qx <> compare qy py+++-- | Precondition: The list of input points is sorted+hull'          :: (Ord r, Num r) => [MyPoint r :+ p] -> NonEmpty (MyPoint r :+ p)+hull' (a:b:ps) = NonEmpty.fromList $ hull'' [b,a] ps+  where+    hull'' h []      = h+    hull'' h (p:ps') = hull'' (cleanMiddle (p:h)) ps'++    cleanMiddle h@[_,_]                         = h+    cleanMiddle h@(z:y:x:rest)+      | rightTurn (x^.core) (y^.core) (z^.core) = h+      | otherwise                               = cleanMiddle (z:x:rest)+    cleanMiddle _                               = error "cleanMiddle: too few points"++rightTurn       :: (Ord r, Num r) => MyPoint r -> MyPoint r -> MyPoint r -> Bool+rightTurn a b c = ccwP a b c == CW++++ccwP :: (Ord r, Num r) => MyPoint r -> MyPoint r -> MyPoint r -> CCW+ccwP p q r = case z `compare` 0 of+              LT -> CW+              GT -> CCW+              EQ -> CoLinear+     where++       MyPoint ux uy = q `subt` p+       MyPoint vx vy = r `subt` p+       z             = ux * vy - uy * vx

@@ -0,0 +1,7 @@+module Benchmark.Util where++++-- | Generates different size benchmarks+sizes    :: Foldable f => f a -> [Int]+sizes xs = let n = length xs in (\i -> n*i `div` 100) <$> [5,10..100]

@@ -1,69 +1,6 @@ module Main where -import           Control.DeepSeq-import           Control.Lens-import           Criterion.Main-import           Data.Ext-import           Data.Geometry.Interval-import qualified Data.Geometry.IntervalTree as IT-import           Data.Geometry.SegmentTree (I(..))-import qualified Data.Geometry.SegmentTree as SegTree-import qualified Data.List.NonEmpty as NonEmpty-import           QuickCheck.Instances-import           Test.QuickCheck----- | generates n random intervals-genIntervals                  :: (Ord r, Arbitrary r)-                              => proxy r -> Int -> IO [Interval () r]-genIntervals _ n | n <= 0     = error "genIntervals: need n > 0"-                 | otherwise  = generate (vectorOf n arbitrary)--genQueries                      :: (Ord r, Arbitrary r)-                                => proxy r -> Int -> IO [r]-genQueries _ n | n <= 0     = error "genQueries: need n > 0"-               | otherwise  = generate (vectorOf n arbitrary)----- genQuerySetup     :: (Ord r, Arbitrary r)---                   => proxy r -> Int -> IO (Int,IT.IntervalTree (I (Interval () r)) r, [r])--- genQuerySetup p n = (\is qs -> (n, IT.fromIntervals . fmap I $ is, qs))---                  <$> genIntervals p n---                  <*> genQueries   p n---main = defaultMain [-  bgroup "IntervalTree" [ env (genIntervals (I (5 :: Int)) (100000 :: Int)) benchBuild-                         -- env (genIntervals (I (5 :: Int)) (100000 :: Int)) benchQueryIT-                        ]-  ]--benchBuild    :: (Ord r, NFData r) => [Interval () r] -> Benchmark-benchBuild is = bgroup "build" [ bench (show n) $ nf IT.fromIntervals (take n is')-                               | n <- sizes is-                               ]-  where-    is' = I <$> is--benchQueryIT    :: (Ord r, Arbitrary r, NFData r) => [Interval () r] -> Benchmark-benchQueryIT is = bgroup "queries"-    [ env (setup n) (\(t,qs) ->-                        bench ("queries on size" ++ show n) $ whnf (queryAll t) qs)-    | n <- sizes is-    ]-  where-    is'        = I <$> is-    r          = is^.to head.start.core-    setup n    = (IT.fromIntervals (take n is'),) <$> genQueries (I r) 100000-    queryAll t = map (flip IT.search t)----- benchQueryIT          :: Ord r---                       => (Int, IT.IntervalTree (I (Interval () r)) r, [r]) -> Benchmark--- benchQueryIT (n,t,qs) = bgroup "queries" [ bench "query" $ whnf (flip IT.search t) q---                                          | q <- qs---                                          ]-+import qualified Algorithms.Geometry.ConvexHull.Bench as M -sizes    :: [a] -> [Int]-sizes xs = let n = length xs in (\i -> n*i `div` 100) <$> [5,10..100]+main :: IO ()+main = M.main

@@ -0,0 +1,79 @@+module Data.Geometry.IntervalTreeBench where++import           Benchmark.Util+import           Control.DeepSeq+import           Control.Lens+import           Criterion.Main+import           Criterion.Types+import           Data.Ext+import           Data.Geometry.Interval+import qualified Data.Geometry.IntervalTree as IT+import           Data.Geometry.SegmentTree (I(..))+import qualified Data.Geometry.SegmentTree as SegTree+import qualified Data.List.NonEmpty as NonEmpty+import           Debug.Trace+import           Test.QuickCheck+import           Test.QuickCheck.HGeometryInstances ()++--------------------------------------------------------------------------------++main :: IO ()+main = defaultMainWith cfg [ intervalBench ]+  where+    cfg = defaultConfig { reportFile = Just "bench.html" }++intervalBench :: Benchmark+intervalBench = bgroup "IntervalTree"+    [ -- env (genIntervals (I (5 :: Int)) 1000) benchBuild+      env (genIntervals (I (5 :: Int)) 100) benchQueryIT+    ]++--------------------------------------------------------------------------------++-- | generates n random intervals+genIntervals                  :: (Ord r, Arbitrary r)+                              => proxy r -> Int -> IO [Interval () r]+genIntervals _ n | n <= 0     = error "genIntervals: need n > 0"+                 | otherwise  = generate (vectorOf n arbitrary)++genQueries                      :: (Ord r, Arbitrary r)+                                => proxy r -> Int -> IO [r]+genQueries _ n | n <= 0     = error "genQueries: need n > 0"+               | otherwise  = generate (vectorOf n arbitrary)+++-- genQuerySetup     :: (Ord r, Arbitrary r)+--                   => proxy r -> Int -> IO (Int,IT.IntervalTree (I (Interval () r)) r, [r])+-- genQuerySetup p n = (\is qs -> (n, IT.fromIntervals . fmap I $ is, qs))+--                  <$> genIntervals p n+--                  <*> genQueries   p n+++-- | Benchmark building the interval tree+benchBuild    :: (Ord r, NFData r) => [Interval () r] -> Benchmark+benchBuild is = bgroup "build" [ bench (show n) $ nf IT.fromIntervals (take n is')+                               | n <- sizes is+                               ]+  where+    is' = I <$> is++benchQueryIT    :: (Ord r, Arbitrary r, NFData r) => [Interval () r] -> Benchmark+benchQueryIT is = bgroup "queries"+    [ env (setup' n) (\(t,qs) ->+                        bench ("queries on size" ++ show n) $ whnf (queryAll t) qs)+    | n <- sizes is+    ]+  where+    is'        = I <$> is+    r          = is^.to head.start.core+    setup' n  = traceShow "setup" $ setup n++    setup n    = (IT.fromIntervals (take n is'),) <$> genQueries (I r) 100000+    queryAll t = map (flip IT.search t)+++-- benchQueryIT          :: Ord r+--                       => (Int, IT.IntervalTree (I (Interval () r)) r, [r]) -> Benchmark+-- benchQueryIT (n,t,qs) = bgroup "queries" [ bench "query" $ whnf (flip IT.search t) q+--                                          | q <- qs+--                                          ]

@@ -0,0 +1,258 @@+{-# LANGUAGE ScopedTypeVariables  #-}+{-# LANGUAGE UndecidableInstances #-}+module Data.Geometry.Vector.VectorFamily6 where++import           Control.Applicative (liftA2)+import           Control.DeepSeq+import           Control.Lens hiding (element)+-- import           Data.Aeson (ToJSON(..),FromJSON(..))+import qualified Data.Foldable as F+import qualified Data.Geometry.Vector.VectorFixed as FV+import           Data.Maybe (fromMaybe)+import           Data.Proxy+import           Data.Semigroup+import           Data.Traversable (foldMapDefault,fmapDefault)+import qualified Data.Vector.Fixed as V+import           Data.Vector.Fixed.Cont (Peano(..), PeanoNum(..), Fun(..))+import           GHC.TypeLits+import           Linear.Affine (Affine(..))+import           Linear.Metric+import qualified Linear.V2 as L2+import qualified Linear.V3 as L3+import qualified Linear.V4 as L4+import           Linear.Vector++--------------------------------------------------------------------------------+-- * d dimensional Vectors+++type One = S Z+type Two = S One+type Three = S Two+type Four = S Three+type Many d = S (S (S (S (S d))))+++type family FromPeano (d :: PeanoNum) :: Nat where+  FromPeano Z     = 0+  FromPeano (S d) = 1 + FromPeano d+++data SingPeano (d :: PeanoNum) where+  SZ :: SingPeano Z+  SS :: !(SingPeano d) -> SingPeano (S d)++class ImplicitPeano (d :: PeanoNum) where+  implicitPeano :: SingPeano d+instance ImplicitPeano Z where+  implicitPeano = SZ+instance ImplicitPeano d => ImplicitPeano (S d) where+  implicitPeano = SS implicitPeano++-- | Mapping between the implementation type, and the actual implementation.+type family VectorFamilyF (d :: PeanoNum) :: * -> * where+  VectorFamilyF Z        = Const ()+  VectorFamilyF One      = Identity+  VectorFamilyF Two      = L2.V2+  VectorFamilyF Three    = L3.V3+  VectorFamilyF Four     = L4.V4+  VectorFamilyF (Many d) = FV.Vector (FromPeano (Many d))+++-- | Datatype representing d dimensional vectors. The default implementation is+-- based n VectorFixed. However, for small vectors we automatically select a+-- more efficient representation.+newtype VectorFamily (d :: PeanoNum) (r :: *) =+  VectorFamily { _unVF :: VectorFamilyF d r }++type ImplicitArity d = (ImplicitPeano d, V.Arity (FromPeano d))++++instance (Eq r, ImplicitArity d) => Eq (VectorFamily d r) where+  (VectorFamily u) == (VectorFamily v) = case (implicitPeano :: SingPeano d) of+        SZ                         -> u == v+        (SS SZ)                    -> u == v+        (SS (SS SZ))               -> u == v+        (SS (SS (SS SZ)))          -> u == v+        (SS (SS (SS (SS SZ))))     -> u == v+        (SS (SS (SS (SS (SS _))))) -> u == v+  {-# INLINE (==) #-}++instance (Ord r, ImplicitArity d) => Ord (VectorFamily d r) where+  (VectorFamily u) `compare` (VectorFamily v) = case (implicitPeano :: SingPeano d) of+        SZ                         -> u `compare` v+        (SS SZ)                    -> u `compare` v+        (SS (SS SZ))               -> u `compare` v+        (SS (SS (SS SZ)))          -> u `compare` v+        (SS (SS (SS (SS SZ))))     -> u `compare` v+        (SS (SS (SS (SS (SS _))))) -> u `compare` v+  {-# INLINE compare #-}+++instance ImplicitArity d => Functor (VectorFamily d) where+  fmap f = VectorFamily . g f . _unVF+    where g = case (implicitPeano :: SingPeano d) of+                SZ                         -> fmap+                (SS SZ)                    -> fmap+                (SS (SS SZ))               -> fmap+                (SS (SS (SS SZ)))          -> fmap+                (SS (SS (SS (SS SZ))))     -> fmap+                (SS (SS (SS (SS (SS _))))) -> fmap+  {-# INLINE fmap #-}+++instance ImplicitArity d => Foldable (VectorFamily d) where+  foldMap f = g f . _unVF+    where g = case (implicitPeano :: SingPeano d) of+                SZ                         -> foldMap+                (SS SZ)                    -> foldMap+                (SS (SS SZ))               -> foldMap+                (SS (SS (SS SZ)))          -> foldMap+                (SS (SS (SS (SS SZ))))     -> foldMap+                (SS (SS (SS (SS (SS _))))) -> foldMap+  {-# INLINE foldMap #-}++instance ImplicitArity d => Traversable (VectorFamily d) where+  traverse f = fmap VectorFamily . g f . _unVF+    where g = case (implicitPeano :: SingPeano d) of+                SZ                         -> traverse+                (SS SZ)                    -> traverse+                (SS (SS SZ))               -> traverse+                (SS (SS (SS SZ)))          -> traverse+                (SS (SS (SS (SS SZ))))     -> traverse+                (SS (SS (SS (SS (SS _))))) -> traverse+  {-# INLINE traverse #-}++instance ImplicitArity d => Applicative (VectorFamily d) where+  pure = VectorFamily . case (implicitPeano :: SingPeano d) of+                SZ                         -> pure+                (SS SZ)                    -> pure+                (SS (SS SZ))               -> pure+                (SS (SS (SS SZ)))          -> pure+                (SS (SS (SS (SS SZ))))     -> pure+                (SS (SS (SS (SS (SS _))))) -> pure+  {-# INLINE pure #-}+  liftA2 f (VectorFamily u) (VectorFamily v) = VectorFamily $+      case (implicitPeano :: SingPeano d) of+                SZ                         -> liftA2 f u v+                (SS SZ)                    -> liftA2 f u v+                (SS (SS SZ))               -> liftA2 f u v+                (SS (SS (SS SZ)))          -> liftA2 f u v+                (SS (SS (SS (SS SZ))))     -> liftA2 f u v+                (SS (SS (SS (SS (SS _))))) -> liftA2 f u v+  {-# INLINE liftA2 #-}+++++type instance V.Dim (VectorFamily d)  = FromPeano d+++++instance ImplicitArity d => V.Vector (VectorFamily d) r where+  construct = fmap VectorFamily $ case (implicitPeano :: SingPeano d) of+                SZ                         -> Fun $ Const ()+                (SS SZ)                    -> V.construct+                (SS (SS SZ))               -> Fun L2.V2+                (SS (SS (SS SZ)))          -> Fun L3.V3+                (SS (SS (SS (SS SZ))))     -> Fun L4.V4+                (SS (SS (SS (SS (SS _))))) -> V.construct+  {-# INLINE construct #-}+  inspect (VectorFamily v) ff@(Fun f) = case (implicitPeano :: SingPeano d) of+                SZ                         -> f+                (SS SZ)                    -> V.inspect v ff+                (SS (SS SZ))               -> let (L2.V2 x y) = v     in f x y+                (SS (SS (SS SZ)))          -> let (L3.V3 x y z) = v   in f x y z+                (SS (SS (SS (SS SZ))))     -> let (L4.V4 x y z w) = v in f x y z w+                (SS (SS (SS (SS (SS _))))) -> V.inspect v ff+  {-# INLINE inspect #-}+  -- basicIndex (VectorFamily v) i = case (implicitPeano :: SingPeano d) of+  --               SZ                         -> err+  --               (SS SZ)                    -> if i == 0 then runIdentity v else err+  --               (SS (SS SZ))               -> let (L2.V2 x y) = v     in f x y+  --               (SS (SS (SS SZ)))          -> let (L3.V3 x y z) = v   in f x y z+  --               (SS (SS (SS (SS SZ))))     -> let (L4.V4 x y z w) = v in f x y z w+  --               (SS (SS (SS (SS (SS _))))) -> V.basicIndex v i+  --   where+  --     err = error "VectorFamily: basicIndex out of range"+  -- {-# INLINE basicIndex #-}+++instance (ImplicitArity d, Show r) => Show (VectorFamily d r) where+  show v = mconcat [ "Vector", show $ F.length v , " "+                   , show $ F.toList v ]++deriving instance (NFData (VectorFamilyF d r)) => NFData (VectorFamily d r)+++type instance Index   (VectorFamily d r) = Int+type instance IxValue (VectorFamily d r) = r++--------------------------------------------------------------------------------+++newtype Vector (d :: Nat) (r :: *) = MKVector { _unV :: VectorFamily (Peano d) r }++type instance V.Dim (Vector d)  = d+++type instance Index   (Vector d r) = Int+type instance IxValue (Vector d r) = r++type Arity d = ImplicitArity (Peano d)++deriving instance (Eq r,  Arity d) => Eq  (Vector d r)+deriving instance (Ord r, Arity d) => Ord (Vector d r)++deriving instance Arity d => Functor     (Vector d)+deriving instance Arity d => Foldable    (Vector d)+deriving instance Arity d => Traversable (Vector d)++instance (Arity d, Show r) => Show (Vector d r) where+  show v = mconcat [ "Vector", show $ F.length v , " "+                   , show $ F.toList v ]+++deriving instance (NFData (VectorFamily (Peano d) r)) => NFData (Vector d r)+++++--------------------------------------------------------------------------------+-- * Convenience "constructors"++pattern Vector   :: VectorFamilyF (Peano d) r -> Vector d r+pattern Vector v = MKVector (VectorFamily v)++pattern Vector1   :: r -> Vector 1 r+pattern Vector1 x = (Vector (Identity x))++pattern Vector2     :: r -> r -> Vector 2 r+pattern Vector2 x y = (Vector (L2.V2 x y))++pattern Vector3        :: r -> r -> r -> Vector 3 r+pattern Vector3 x y z  = (Vector (L3.V3 x y z))++pattern Vector4         :: r -> r -> r -> r -> Vector 4 r+pattern Vector4 x y z w = (Vector (L4.V4 x y z w))++--------------------------------------------------------------------------------++-- -- destruct            :: (Vec d r, Vec (d + 1) r, 1 <= (d + 1))+-- --                     => Vector (d + 1) r -> (r, Vector d r)+-- -- destruct (Vector v) = (V.head v, Vector $ V.tail v)+++-- -- -- vectorFromList :: Arity d => [a] -> Maybe (Vector d a)+-- -- vectorFromList = fmap Vector . V.fromListM++-- -- vectorFromListUnsafe :: V.Arity d => [a] -> Vector d a+-- -- vectorFromListUnsafe = Vector . V.fromList++ --------------------------------------------------------------------------------++-- | Cross product of two three-dimensional vectors+cross       :: Num r => Vector 3 r -> Vector 3 r -> Vector 3 r+(Vector u) `cross` (Vector v) = Vector $ u `L3.cross` v

@@ -0,0 +1,26 @@+{-# LANGUAGE PartialTypeSignatures #-}+module WSPDBench where++import Control.DeepSeq+import Control.Lens+import Criterion.Main+import Data.Ext+import Data.Geometry.Point+import Demo.ExpectedPairwiseDistance+++readInput'     :: FilePath -> Int -> IO [Point 2 Double :+ _]+readInput' fp k = take k <$> readInput fp+++benchWSPD :: Benchmark+benchWSPD = bgroup "Well-Separated Pair Decomposition"+    [ env (readInput' "pco.9420_convert.txt" 200) $ \pts -> bgroup "pco"+        [ bench "exact"     $ nf pairwiseDist                  pts+        , bench "wspd 0.05" $ nf (approxPairwiseDistance 0.05) pts+        , bench "wspd 0.10" $ nf (approxPairwiseDistance 0.10) pts+        , bench "wspd 0.20" $ nf (approxPairwiseDistance 0.20) pts+        ]+    ]++-- main = defaultMain [ benchWSPD ]

@@ -0,0 +1,60 @@+#### 0.7 ###++- Compatibility with GHC 8.0-8.4+- Implemented more Algorithms and Data Structures. This includes+  * Polygon triangulation+- A new implementation of PlanarSubdivision that now also supports disconnected+  subdivsions.+- Performance improvements by changing to a different Vector+  implementation. For low dimensional vectors (of dimension at most four) we+  now essentially use the types from+  [linear](https://hackage.haskell.org/package/linear), this gives significant+  speedups on several small benchmarks.+- bugfixes.++#### 0.6 ###++- Implemented more Algorithms and Data Structures. This includes+  * Bentley-Ottmannn line-segment intersection,+  * Well-Separated Pair decompositions,+  * extremal point/tangents for Convex hulls,+  * Minkowski sum for convex polygons,+  * one dimensional segment trees,+  * one dimensional interval trees, and a+  * KD-tree.+- Several bug fixes, including a very stupid bug in Box+- Separate ConvexPolygon type.+- More thorough testing for some of the algorithms.+- Started work on a proper representation for planar subdivsions. This includes+  a representation of planar graphs that support querying if two vertices are+  connected by an edge in $O(1)$ time.+- Dropped support for GHC 7.8++#### 0.5 ###++- Implemented several algorithms, including Delaunay Triangulation, EMST, and+Douglas Peucker.+- Revamped the data types for Intersections++#### 0.4 ###++- Major rewrite from scratch, providing much stronger type-level+  guarantees. Incompatible with older versions.+- Convex Hull and Smallest enclosing disk algorithms.+- HGeometry now includes some very experimental and preliminary support for+  reading and writing Ipe7 files.++#### 0.2 & 0.3 ###++- Internal releases.++### 0.1.1 ###++- Fixed a bug in point on n the line segment test+- Generalized the types of inCircle, inDisc, onCircle, onDisc etc. We now need+  only that the type representing precision model implements the typeclass+  `Num` instead of `Floating'.++### 0.1 ###++- Initial release.

@@ -60,7 +60,7 @@   -> triangArea :: Triangle p Int -> Half+> triangArea :: Triangle 2 p Int -> Half > triangArea = Half . doubleArea  @@ -193,7 +193,7 @@   > findLargestTriang        :: Point 2 Int -> Point 2 Int->                          -> Unimodal (Array Int) (Point 2 Int) -> Triangle () Int+>                          -> Unimodal (Array Int) (Point 2 Int) -> Triangle 2 () Int > findLargestTriang p q us = triang . ternarySearchArray area' $ us >   where >     triang v = Triangle (ext p) (ext q) (ext v)

@@ -4,7 +4,6 @@ import           Algorithms.Geometry.DelaunayTriangulation.DivideAndConqueror import           Algorithms.Geometry.DelaunayTriangulation.Types import           Algorithms.Geometry.EuclideanMST.EuclideanMST-import           Control.Applicative import           Control.Lens import           Data.Data import           Data.Ext@@ -12,7 +11,6 @@ import           Data.Geometry.Ipe import qualified Data.List.NonEmpty as NonEmpty import           Data.Semigroup-import           Data.Traversable import           Options.Applicative  

@@ -20,7 +20,7 @@ import           Data.Data import           Data.Semigroup import qualified Data.Set as Set-import           GHC.TypeLits (natVal,KnownNat)+import           GHC.TypeLits import           Options.Applicative hiding ((<>))  @@ -56,7 +56,7 @@ -- -- running time: $O(n(1/eps)^d + n\log n)$, where $n$ is the number of points approxExpectedPairwiseDistance          :: (Floating r, Ord r-                                           , AlwaysTrueWSPD d, Index' 0 d+                                           , Arity d, Arity (d+1), 1 <= d                                          , Show r, Show p)                                          => r -> Int -> [Point d r :+ p] -> r approxExpectedPairwiseDistance eps k pts =@@ -73,7 +73,7 @@ -- | $(1+\eps)$-approximation of the sum of the pairwise distances. -- -- running time: $O(n(1/eps)^d + n\log n)$, where $n$ is the number of points-approxPairwiseDistance         :: (Floating r, Ord r, AlwaysTrueWSPD d, Index' 0 d+approxPairwiseDistance         :: (Floating r, Ord r, Arity d, Arity (d+1), 1 <= d                                   , Show r, Show p)                                => r -> [Point d r :+ p] -> r approxPairwiseDistance _   []  = 0@@ -179,7 +179,7 @@   -- | Computes all pairs of points that are uncovered by the WSPD with separation s-uncovered         :: (Floating r, Ord r, AlwaysTrueWSPD d, Ord p)+uncovered         :: (Floating r, Ord r, Arity d, Arity (d+1), Ord p)                   => [Point d r :+ p] -> r -> SplitTree d p r a -> [(Point d r :+ p, Point d r :+ p)] uncovered pts s t = Set.toList $ allPairs `Set.difference` covered   where

@@ -85,7 +85,7 @@ extract = (\(Text s) -> s) . head . eChildren  readTime' :: String -> UTCTime-readTime' = readTime defaultTimeLocale "%0C%y-%m-%dT%TZ"+readTime' = parseTimeOrError True defaultTimeLocale "%0C%y-%m-%dT%TZ"  -- instance ReadGPX Position where --   parseGPX x@(Element "Position" _ _) = (\l l' -> Position $ point2 l l') <$> lat <*> lon

@@ -0,0 +1,106 @@+{-# LANGUAGE ScopedTypeVariables #-}+module Demo.TriangulateWorld where++import Algorithms.Geometry.LineSegmentIntersection (hasSelfIntersections)+import Algorithms.Geometry.PolygonTriangulation.Triangulate (triangulate)+import Algorithms.Geometry.PolygonTriangulation.MakeMonotone (makeMonotone)+import Data.Maybe(mapMaybe)+import Control.Lens+import Data.Data+import Data.Ext+import Data.Geometry.Ipe+import Data.Geometry.Polygon+import Data.Geometry.PlanarSubdivision+import Data.Semigroup+import Options.Applicative+import qualified Data.Foldable as F+++--------------------------------------------------------------------------------++data Options = Options { _inPath    :: FilePath+                       , _outFile   :: FilePath+                       }+               deriving Data++options :: ParserInfo Options+options = info (helper <*> parser)+               (  progDesc "Triangulate all polygons in the input file."+               <> header   "trianguldateWorld"+               )+  where+    parser = Options+          <$> strOption (help "Input file (in ipe7 xml format)"+                         <> short 'i'+                        )+          <*> strOption (help "Output File (in ipe7 xml format)"+                         <> short 'o'+                        )++++-- runExcept'   :: (Show e) => ExceptT e IO () -> IO ()+-- runExcept' m = runExceptT m >>= \case+--                 Left e   -> print e+--                 Right () -> pure ()++-- mainWith                          :: Options -> IO ()+-- mainWith (Options inFile outFile) = runExcept' $ do+--     (page :: IpePage Rational) <- readSinglePageFile inFile+--     let polies = page^..content.traverse._withAttrs _IpePath _asSimplePolygon+--     let out = undefinedL+--     lift $  writeIpeFile outFile . singlePageFromContent $ out++data PX = PX++mainWith                          :: Options -> IO ()+mainWith (Options inFile outFile) = do+    ePage <- readSinglePageFile inFile+    case ePage of+      Left err                         -> print err+      Right (page :: IpePage Rational) -> runPage page+  where+    runPage page = do+      let polies  = page^..content.to flattenGroups.traverse._withAttrs _IpePath _asSimplePolygon+          polies' = filter (not . hasSelfIntersections . (^.core)) polies+          subdivs = map (\(pg :+ _) -> triangulate (Identity PX) pg) polies'+          yMonotones = tail . mapMaybe (^?_2.core._Left)+                     . concatMap (F.toList.rawFacePolygons) $ subdivs+          ofs = map (\s -> rawFaceBoundary (outerFaceId s) s) subdivs+          segs    = map (^._2.core) . concatMap (F.toList . edgeSegments) $ subdivs+          out     = [ asIpeObject pg a+                    | pg :+ a <- polies+                    ] <>+                    [ asIpeObject s mempty+                    | s <- segs+                    ] <>+                    [ asIpeObject pg mempty+                    | pg <- yMonotones ]+      mapM_ print . map (\pg -> pg^.core.to polygonVertices.to length) $ polies'+      writeIpeFile outFile . singlePageFromContent $ out+++-- mainWith                          :: Options -> IO ()+-- mainWith (Options inFile outFile) = do+--     ePage <- readSinglePageFile inFile+--     case ePage of+--       Left err                         -> print err+--       Right (page :: IpePage Rational) -> runPage page+--   where+--     runPage page = do+--       let orig = page^.content+--           all' = page^.content.to flattenGroups+--       writeIpeFile outFile . singlePageFromContent $ orig <> all'+++-- type ValT = EitherT IO++-- flattenGroups :: [IpeObject r] -> [IpeObject r]+-- flattenGroups = concatMap flattenGroups'++-- flattenGroups'                              :: IpeObject r -> [IpeObject r]+-- flattenGroups' (IpeGroup (Group gs :+ ats)) =+--       map (applyAts ats) . concatMap flattenGroups' $ gs+--     where+--       applyAts ats = id+-- flattenGroups' o                            = [o]

@@ -12,6 +12,7 @@ import qualified Demo.MinDisk as MinDisk import qualified Demo.Delaunay as Delaunay import qualified Demo.ExpectedPairwiseDistance as ExpPWD+import qualified Demo.TriangulateWorld as TriangulateWorld   @@ -24,6 +25,7 @@              | MinDisk                   MinDisk.Options              | Delaunay                  Delaunay.Options              | ExpectedPairwiseDistance  ExpPWD.Options+             | TriangulateWorld          TriangulateWorld.Options              deriving Data  parser :: Parser Options@@ -33,6 +35,7 @@     <> command' MinDisk                        MinDisk.options     <> command' Delaunay                       Delaunay.options     <> command' ExpectedPairwiseDistance       ExpPWD.options+    <> command' TriangulateWorld               TriangulateWorld.options     )  @@ -44,6 +47,7 @@   MinDisk opts                        -> MinDisk.mainWith opts   Delaunay opts                       -> Delaunay.mainWith opts   ExpectedPairwiseDistance opts       -> ExpPWD.mainWith opts+  TriangulateWorld opts               -> TriangulateWorld.mainWith opts  -------------------------------------------------------------------------------- 

@@ -2,7 +2,7 @@ -- documentation, see http://haskell.org/cabal/users-guide/  name:                hgeometry-version:             0.6.0.0+version:             0.7.0.0 synopsis:            Geometric Algorithms, Data structures, and Data types. description:   HGeometry provides some basic geometry types, and geometric algorithms and@@ -17,17 +17,20 @@ maintainer:          frank@fstaals.net -- copyright: -tested-with:         GHC >= 7.10.2+tested-with:         GHC >= 8.2  category:            Geometry build-type:          Simple -extra-source-files:  README.md-                     resources/basic.isy+data-files:          resources/basic.isy+                     test/Algorithms/Geometry/LineSegmentIntersection/manual.ipe+                     test/Algorithms/Geometry/LineSegmentIntersection/selfIntersections.ipe+                     test/Algorithms/Geometry/LowerEnvelope/manual.ipe+                     test/Algorithms/Geometry/PolygonTriangulation/monotone.ipe+                     test/Algorithms/Geometry/PolygonTriangulation/simplepolygon6.ipe+                     test/Algorithms/Geometry/SmallestEnclosingDisk/manual.ipe                      test/Data/Geometry/pointInPolygon.ipe                      test/Data/Geometry/Polygon/Convex/convexTests.ipe-                     test/Algorithms/Geometry/SmallestEnclosingDisk/manual.ipe-                     test/Algorithms/Geometry/LineSegmentIntersection/manual.ipe                      examples/BAPC2014/sample.in                      examples/BAPC2014/sample.out                      examples/BAPC2014/testdata.in@@ -38,7 +41,11 @@                      examples/BAPC2012/sampleG.out  -cabal-version:       >=1.10++extra-source-files:  README.md+                     changelog.md++cabal-version:       2.0 source-repository head   type:     git   location: https://github.com/noinia/hgeometry@@ -49,8 +56,18 @@   default:     False   manual:      True +flag with-quickcheck+  description: Include QuickCheck instances+  default:     True+  manual:      False++flag interactive+  description: Build  interactive parts+  default:     False+  manual:      True+ library-  ghc-options: -Wall -fno-warn-unticked-promoted-constructors -fno-warn-type-defaults+  ghc-options: -O2 -Wall -fno-warn-unticked-promoted-constructors -fno-warn-type-defaults    exposed-modules:                     -- * Generic Geometry@@ -63,6 +80,9 @@                     -- * Basic Geometry Types                     Data.Geometry.Vector                     Data.Geometry.Vector.VectorFixed+                    Data.Geometry.Vector.VectorFamily+                    Data.Geometry.Vector.VectorFamilyPeano+                     -- Data.Geometry.Vector.Vinyl                     Data.Geometry.Interval                     Data.Geometry.Interval.Util@@ -89,6 +109,8 @@                      Data.Geometry.KDTree                     Data.Geometry.PlanarSubdivision+                    Data.Geometry.PlanarSubdivision.Basic+                    Data.Geometry.PlanarSubdivision.Draw                      -- * Algorithms                     Algorithms.Util@@ -97,6 +119,8 @@                     Algorithms.Geometry.ConvexHull.GrahamScan                     Algorithms.Geometry.ConvexHull.DivideAndConqueror +                    Algorithms.Geometry.LowerEnvelope.DualCH+                     Algorithms.Geometry.SmallestEnclosingBall.Types                     Algorithms.Geometry.SmallestEnclosingBall.RandomizedIncrementalConstruction                     Algorithms.Geometry.SmallestEnclosingBall.Naive@@ -114,6 +138,14 @@                      Algorithms.Geometry.Diameter +                    Algorithms.Geometry.Sweep++                    Algorithms.Geometry.PolygonTriangulation.Types+                    Algorithms.Geometry.PolygonTriangulation.Triangulate+                    Algorithms.Geometry.PolygonTriangulation.MakeMonotone+                    Algorithms.Geometry.PolygonTriangulation.TriangulateMonotone++                    Algorithms.Geometry.LineSegmentIntersection                     Algorithms.Geometry.LineSegmentIntersection.Naive                     Algorithms.Geometry.LineSegmentIntersection.BentleyOttmann                     Algorithms.Geometry.LineSegmentIntersection.Types@@ -146,16 +178,20 @@                      Data.CircularList.Util                     Data.BalBST+                    Data.OrdSeq+                    Data.SlowSeq                     Data.Util                      -- * Planar Graphs                     Data.Permutation                     Data.PlanarGraph                     Data.PlaneGraph+                    Data.PlaneGraph.Draw                      -- * Other                     System.Random.Shuffle                     Control.Monad.State.Persistent+                    Data.Yaml.Util   @@ -167,12 +203,11 @@    -- other-extensions:   build-depends:--                  Frames           >= 0.1.3.0-                , base             >= 4.8       &&     < 5+                  base             >= 4.9       &&     < 5                 , bifunctors       >= 4.1                 , bytestring       >= 0.10                 , containers       >= 0.5.5+                , dlist            >= 0.7                 , contravariant    >= 1.4                 , lens             >= 4.2                 , linear           >= 1.10@@ -180,8 +215,11 @@                 , semigroups       >= 0.18                 , singletons       >= 2.0                 , text             >= 1.1.1.0-                , vinyl            >= 0.5       &&     < 0.6+                , vinyl            >= 0.6                 , deepseq          >= 1.1+                , fingertree       >= 0.1+                , colour           >= 2.3.3+                , reflection       >= 2.1                 -- , validation       >= 0.4 @@ -189,18 +227,22 @@                -- , tranformers      > 0.3                 , vector           >= 0.11-               , fixed-vector     >= 0.6.4.0+               , fixed-vector     >= 1.0                , data-clist       >= 0.0.7.2                 , hexpat           >= 0.20.9+               , aeson            >= 1.0+               , yaml             >= 0.8++                , mtl                , random                , template-haskell --               , time-               , directory-               , optparse-applicative+  if flag(with-quickcheck)+    build-depends: QuickCheck              >= 2.5+                 , quickcheck-instances    >= 0.3+    exposed-modules: Test.QuickCheck.HGeometryInstances     hs-source-dirs: src@@ -235,6 +277,80 @@                     , FlexibleContexts                     , MultiParamTypeClasses +executable hgeometry-viewer+  if !flag(interactive)+    buildable: False+  main-is: Viewer.hs++  if flag(interactive)+    build-depends: base+                 , hgeometry+                 , lens+                 , containers+                 , vinyl+                 , semigroups+                 , optparse-applicative   >= 0.13.0.0+                 , text+                 , hexpat+                 , bytestring+                 , directory+                 , time+                 , random+                 , vector+                 , colour+                 , cairo-canvas           >= 0.1.0.0+                 -- , sdl2                   >= 2.2.0+                 , gi-gtk                 >= 3.0.15+                 , reactive-banana-gi-gtk >= 0.2.0.0+                 , cairo                  >= 0.13.3.1+                 -- , gi-glib+                 , gi-cairo+                 , gi-gdk+                 -- , gi-gdkpixbuf+                 , gi-gtk+                 , transformers+                 , linear+                 , haskell-gi-base+                 , reactive-banana        >= 1.1.0.1++++  hs-source-dirs: interactive++  other-modules: RenderCanvas+                 RenderUtil+                 -- ConvexHull++  default-language: Haskell2010++  default-extensions: TypeFamilies+                    , GADTs+                    , KindSignatures+                    , DataKinds+                    , TypeOperators+                    , ConstraintKinds+                    , PolyKinds+                    , RankNTypes++                    , PatternSynonyms+                    , ViewPatterns++                    , StandaloneDeriving+                    , GeneralizedNewtypeDeriving+                    , DeriveFunctor+                    , DeriveFoldable+                    , DeriveTraversable++                    , DeriveDataTypeable+                    , AutoDeriveTypeable++                    , FlexibleInstances+                    , FlexibleContexts+                    , MultiParamTypeClasses++++ executable hgeometry-examples   if !flag(examples)     buildable: False@@ -246,7 +362,6 @@                  , lens                  , containers                  , vinyl-                 , Frames                  , semigroups                  , optparse-applicative   >= 0.13.0.0                  , text@@ -255,7 +370,7 @@                  , directory                  , time                  , random-+                 , QuickCheck    hs-source-dirs: examples @@ -264,7 +379,7 @@                  Demo.MinDisk                  Demo.Delaunay                  Demo.ExpectedPairwiseDistance-+                 Demo.TriangulateWorld                  Demo.GPXParser    default-language: Haskell2010@@ -300,7 +415,8 @@   type:          exitcode-stdio-1.0   ghc-options:   -threaded   main-is:       doctests.hs-  build-depends: base, doctest >= 0.8+  build-depends: base+               , doctest >= 0.8    default-language:    Haskell2010 @@ -311,32 +427,43 @@   main-is:              Spec.hs   ghc-options:   -O2 -fno-warn-unticked-promoted-constructors +  build-tool-depends: hspec-discover:hspec-discover++   other-modules: Data.RangeSpec                  Data.EdgeOracleSpec                  Data.PlanarGraphSpec+                 Data.OrdSeqSpec                  Data.Geometry.Ipe.ReaderSpec                  Data.Geometry.PolygonSpec+                 Data.Geometry.LineSegmentSpec                  Data.Geometry.PointSpec                  Data.Geometry.Polygon.Convex.ConvexSpec                  Data.Geometry.KDTreeSpec                  Data.Geometry.IntervalSpec                  Data.Geometry.BoxSpec-+                 Data.Geometry.LineSpec+                 Data.Geometry.SubLineSpec+                 Data.Geometry.PlanarSubdivisionSpec                   Algorithms.Geometry.SmallestEnclosingDisk.RISpec                  Algorithms.Geometry.DelaunayTriangulation.DTSpec                  Algorithms.Geometry.WellSeparatedPairDecomposition.WSPDSpec                  Algorithms.Geometry.LineSegmentIntersection.BentleyOttmannSpec+                 Algorithms.Geometry.PolygonTriangulation.MakeMonotoneSpec+                 Algorithms.Geometry.PolygonTriangulation.TriangulateMonotoneSpec+                 Algorithms.Geometry.LowerEnvelope.LowerEnvSpec+                 Algorithms.Geometry.ConvexHull.ConvexHullSpec -                 QuickCheck.Instances                  Util     build-depends:        base-                      , hspec        >= 2.1-                      , QuickCheck   >= 2.5+                      , hspec                >= 2.1+                      , QuickCheck           >= 2.5+                      , quickcheck-instances    >= 0.3+                      , approximate-equality >= 1.1.0.2                       , hgeometry-                      , Frames                       , lens                       , data-clist                       , linear@@ -346,6 +473,8 @@                       , vector                       , containers                       , random+                      , singletons+                      , colour    default-extensions: TypeFamilies                     , GADTs@@ -405,23 +534,38 @@  benchmark benchmarks -  hs-source-dirs: benchmark test+  hs-source-dirs: benchmark test examples    main-is: Benchmarks.hs   type: exitcode-stdio-1.0 -  other-modules:+  other-modules: Benchmark.Util+                 Algorithms.Geometry.ConvexHull.Bench+                 Algorithms.Geometry.ConvexHull.GrahamV2+                 Algorithms.Geometry.ConvexHull.GrahamFam+                 Algorithms.Geometry.ConvexHull.GrahamFamPeano+                 Algorithms.Geometry.ConvexHull.GrahamFixed+                 Data.Geometry.Vector.VectorFamily6+                 Algorithms.Geometry.ConvexHull.GrahamFam6+                 Data.Geometry.IntervalTreeBench+                 Demo.ExpectedPairwiseDistance+                 Demo.TriangulateWorld+                 WSPDBench    build-depends:                 base-              , criterion >= 1.1.4.0 && < 1.2+              , criterion >= 1.1.4.0+              , fixed-vector+              , linear               , semigroups               , deepseq               , deepseq-generics               , hgeometry-              , Frames               , lens               , QuickCheck+              , bytestring+              , containers+              , optparse-applicative    ghc-options: -Wall -O2 -rtsopts -fno-warn-unticked-promoted-constructors 

@@ -0,0 +1,177 @@+{-# LANGUAGE ScopedTypeVariables #-}+module RenderCanvas where++import           Control.Lens+import           Data.Colour.SRGB(RGB(..), toSRGB24)+import           Data.Colour.Names(readColourName)+import           Data.Ext+import           Data.Geometry+import           Data.Geometry.Box+import           Data.Geometry.Ipe.Attributes+import           Data.Geometry.Ipe.Types hiding (ipeObject', width)+import qualified Data.Geometry.Ipe.Attributes as A+import           Data.Proxy+import           Data.Vinyl+import           Linear.V4 (V4(..))+import           Graphics.Rendering.Cairo.Canvas (Canvas)+import qualified Graphics.Rendering.Cairo.Canvas as Canvas+import qualified Data.Text as T+++rectangle    :: (Real r, Ord r, Num r) => Rectangle p r -> Canvas ()+rectangle r' = let r                       = bimap id realToFrac r'+                   (Point2 x y :+ _,_,_,_) = corners r+               in Canvas.rect $ Canvas.D x y (width r) (height r)++polygon     :: Real r => SimplePolygon p r -> Canvas ()+polygon pg' = let pg = bimap id realToFrac pg'+              in Canvas.polygon $ pg^..outerBoundary.traverse.core.vector.to toV2++lineSegment    :: Real r => LineSegment 2 p r -> Canvas ()+lineSegment s' = let s = bimap id realToFrac s'+                 in Canvas.line (s^.start.core.vector.to toV2) (s^.end.core.vector.to toV2)++polyLine    :: Real r => PolyLine 2 p r -> Canvas ()+polyLine p' = let p = bimap id realToFrac p'+              in Canvas.shape Canvas.ShapeLines $ p^..points.traverse.core.vector.to toV2++-- | draw a point as a small disk+point   :: Real r => Point 2 r -> Canvas ()+point p = Canvas.circle' (realToFrac <$> (p^.vector.to toV2)) 5++-- | draw as a point+point'   :: Real r => Point 2 r -> Canvas ()+point' p = Canvas.point . fmap realToFrac $ p^.vector.to toV2++pathSegment                     :: Real r => PathSegment r -> Canvas ()+pathSegment (PolyLineSegment p) = polyLine p+pathSegment (PolygonPath p)     = polygon p+pathSegment _                   = error "pathSegment: Not implemented yet"+++ipeUse              :: Real r => IpeSymbol r -> Canvas ()+ipeUse (Symbol p _) = Canvas.circle' (realToFrac <$> p^.vector.to toV2) 10++ipePath          :: Real r => Path r -> Canvas ()+ipePath (Path p) = mapM_ pathSegment p++ipeGroup :: RealFrac r => Group r -> Canvas ()+ipeGroup = mapM_ ipeObject . _groupItems++ipeObject'              :: forall g r. (RealFrac r, AllSatisfy ApplyAttr (AttributesOf g))+                        => (g r -> Canvas ())+                        -> g r :+ IpeAttributes g r+                        -> Canvas ()+ipeObject' f (i :+ ats) = do+                            Canvas.pushMatrix+                            applyAttributes (Proxy :: Proxy g) ats+                            f i+                            Canvas.popMatrix++ipeObject                  :: RealFrac r+                           => IpeObject r -> Canvas ()+ipeObject (IpeGroup g)     = ipeObject' ipeGroup g+ipeObject (IpeImage _)     = undefined+ipeObject (IpeTextLabel _) = undefined+ipeObject (IpeMiniPage _)  = undefined+ipeObject (IpeUse p)       = ipeObject' ipeUse  p+ipeObject (IpePath p)      = ipeObject' ipePath p+++applyAttributes               :: (RealFrac r, AllSatisfy ApplyAttr (AttributesOf g))+                              => proxy g -> IpeAttributes g r -> Canvas ()+applyAttributes _ (Attrs ats) = applyAttributes' ats++applyAttributes'            :: (RealFrac r, AllSatisfy ApplyAttr rs)+                            => Rec (Attr (AttrMapSym1 r)) rs+                            -> Canvas ()+applyAttributes' RNil       = pure ()+applyAttributes' (a :& ats) = applyAttribute a >> applyAttributes' ats+++newtype CanvasM = CanvasM { unCanvasM :: Canvas () }+instance Monoid CanvasM where+  mempty = CanvasM $ pure ()+  (CanvasM a) `mappend` (CanvasM b) = CanvasM $ a >> b+++applyAttribute' :: (RealFrac r, ApplyAttr label)+                => Attr (AttrMapSym1 r) label -> CanvasM+applyAttribute' = CanvasM . applyAttribute+++class ApplyAttr (label :: AttributeUniverse) where+  applyAttribute :: RealFrac r => Attr (AttrMapSym1 r) label -> Canvas ()+++instance ApplyAttr Stroke where+  applyAttribute NoAttr   = pure ()+  applyAttribute (Attr c) = maybe (pure ()) Canvas.stroke $ toCanvasColor c++instance ApplyAttr Fill where+  applyAttribute NoAttr   = pure ()+  applyAttribute (Attr c) = maybe (pure ()) Canvas.fill $ toCanvasColor c++instance ApplyAttr Pen where+  applyAttribute NoAttr            = pure ()+  applyAttribute (Attr (IpePen p)) = case p of+      Named _  -> pure () -- TODO+      Valued v -> Canvas.strokeWeight (realToFrac v)++instance ApplyAttr Clip where+  applyAttribute _ = pure ()++instance ApplyAttr Size where+  applyAttribute _ = pure ()++instance ApplyAttr Dash where+  applyAttribute _ = pure ()+++instance ApplyAttr Layer where+  applyAttribute _ = pure ()++instance ApplyAttr LineCap where+  applyAttribute _ = pure ()++instance ApplyAttr LineJoin where+  applyAttribute _ = pure ()++instance ApplyAttr A.Matrix where+  applyAttribute _ = pure ()++instance ApplyAttr Pin where+  applyAttribute _ = pure ()++instance ApplyAttr FillRule where+  applyAttribute _ = pure ()++instance ApplyAttr Arrow where+  applyAttribute _ = pure ()++instance ApplyAttr RArrow where+  applyAttribute _ = pure ()++instance ApplyAttr Opacity where+  applyAttribute _ = pure ()++instance ApplyAttr Tiling where+  applyAttribute _ = pure ()++instance ApplyAttr Gradient where+  applyAttribute _ = pure ()++instance ApplyAttr Transformations where+  applyAttribute _ = pure ()+++++-- | Looks up the colorname in the SVG colors if it is a name.+toCanvasColor :: RealFrac r => IpeColor r -> Maybe Canvas.Color+toCanvasColor (IpeColor c) = case c of+    Named t            -> h . toSRGB24 <$> readColourName (T.unpack $ T.toLower t)+    Valued v           -> Just $ f v+  where+    f (RGB r g b) = floor <$> V4 (255 *r) (255*g) (255*b) 255+    h (RGB r g b) = V4 r g b 255

@@ -0,0 +1,81 @@+{-# LANGUAGE OverloadedLabels #-}+module RenderUtil where++import           Control.Monad.IO.Class (MonadIO)+import           Control.Monad.Trans.Reader (runReaderT)+import           Data.GI.Base+import           Data.GI.Base.Signals (SignalInfo, HaskellCallbackType)+import           Data.IORef+import qualified Data.Text as T+import           Foreign.Ptr (castPtr)+import qualified GI.Cairo as GI.Cairo+import qualified GI.Gtk as Gtk+import           Graphics.Rendering.Cairo+import           Graphics.Rendering.Cairo.Canvas (Canvas)+import qualified Graphics.Rendering.Cairo.Canvas as Canvas+import qualified Graphics.Rendering.Cairo.Internal as Cairo.Internal+import           Graphics.Rendering.Cairo.Types (Cairo(Cairo))+import           Linear.V2 (V2(..))+import           Reactive.Banana+import           Reactive.Banana.Frameworks+import           Reactive.Banana.GI.Gtk+++++-- | This function bridges gi-cairo with the hand-written cairo+-- package. It takes a `GI.Cairo.Context` (as it appears in gi-cairo),+-- and a `Render` action (as in the cairo lib), and renders the+-- `Render` action into the given context.+renderWithContext      :: MonadIO m => GI.Cairo.Context -> Render a -> m a+renderWithContext ct r = liftIO $ withManagedPtr ct $ \p ->+                         runReaderT (Cairo.Internal.runRender r) (Cairo (castPtr p))+++renderCanvas           :: MonadIO m => GI.Cairo.Context -> V2 Double -> Canvas a -> m a+renderCanvas ct size c = liftIO $ Canvas.withRenderer (renderWithContext ct) size c+++showT :: Show a => a -> T.Text+showT = T.pack . show++-- | Get an 'Reactive.Banana.Event' from+-- a 'Data.GI.Base.Signals.SignalProxy' that produces one argument.+signalE1'+    ::+        ( HaskellCallbackType info ~ (a -> IO Bool)+        , SignalInfo info+        , Gtk.GObject self+        )+    => self+    -> SignalProxy self info+    -> (a -> IO b) -- ^ function to transform the Event with+    -> MomentIO (Event b)+signalE1' self signal h = signalEN self signal f >>= mapEventIO h+  where+    f g = \a -> g a >> return True -- we return True because the event has been+                                   -- handled, don't want to propagate it+                                   -- further+++draw                :: Gtk.DrawingArea -> Behavior (Canvas ()) -> MomentIO ()+draw drawingArea bc = do+    canvasRef <- liftIO . newIORef =<< valueB bc -- gets the initial canvas++    -- set up reactive-banana to update the canvasRef on changes, and triger a+    -- redraw+    c <- valueBLater bc+    liftIOLater $ writeIORef canvasRef c+    e <- changes bc+    reactimate' $ (fmap $ \c' -> do liftIO $ writeIORef canvasRef c'+                                    #queueDraw drawingArea+                  ) <$> e++    -- registers drawing event handler+    _ <- on drawingArea #draw $ \context -> do+        w <- realToFrac . fromIntegral <$> #getAllocatedWidth  drawingArea+        h <- realToFrac . fromIntegral <$> #getAllocatedHeight drawingArea+        canvas <- readIORef canvasRef+        renderCanvas context (V2 w h) canvas+        pure True+    pure ()

@@ -0,0 +1,172 @@+{-# LANGUAGE OverloadedLabels #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE ScopedTypeVariables #-}+module Main where++import           Control.Exception (catch)+import           Control.Lens+import           Control.Monad (forM_)+import           Data.Ext+import           Data.GI.Base+import qualified Data.List.NonEmpty as NonEmpty+import qualified Data.Text as T+import qualified GI.Gdk as Gdk+import qualified GI.Gtk as Gtk+import           Graphics.Rendering.Cairo.Canvas (Canvas)+import qualified Graphics.Rendering.Cairo.Canvas as Canvas+import           Linear.Affine ((.-.),(.+^))+import           Linear.V2+import           Linear.Vector ((*^))+import           Reactive.Banana+import           Reactive.Banana.Frameworks+import           Reactive.Banana.GI.Gtk+import qualified RenderCanvas as Render+import           RenderUtil++main :: IO ()+main = runGtk `catch` (\(e::Gtk.GError) -> Gtk.gerrorMessage e >>= putStrLn . T.unpack)+  where+    runGtk = do+      _ <- Gtk.init Nothing+      compile networkDescription >>= actuate+      Gtk.main++data ArrowKey = UpKey | DownKey | LeftKey | RightKey deriving (Show,Read,Eq,Bounded,Enum)++toArrowKey         :: T.Text -> Maybe ArrowKey+toArrowKey "Up"    = Just UpKey+toArrowKey "Down"  = Just DownKey+toArrowKey "Left"  = Just LeftKey+toArrowKey "Right" = Just RightKey+toArrowKey _       = Nothing+++++toDirection          :: Num a => ArrowKey -> V2 a+toDirection UpKey    = V2 0    1+toDirection DownKey  = V2 0    (-1)+toDirection LeftKey  = V2 (-1) 0+toDirection RightKey = V2 1    0+++networkDescription :: MomentIO ()+networkDescription = do+    b <- Gtk.builderNew+    _ <- Gtk.builderAddFromFile b "interactive/viewport.glade"++    window   <- castB b "window" Gtk.Window+    destroyE <- signalE0 window #destroy+    reactimate $ Gtk.mainQuit <$ destroyE++    -- mouseLabel <- castB b "mouseLabel" Gtk.Label++    drawingArea  <- castB b "canvas" Gtk.DrawingArea+    drawingAreaH <- realToFrac . fromIntegral . snd <$> #getPreferredHeight drawingArea+    drawingAreaW <- realToFrac . fromIntegral . snd <$> #getPreferredWidth  drawingArea++    Gtk.widgetAddEvents drawingArea (gflagsToWord [ Gdk.EventMaskPointerMotionMask+                                                  , Gdk.EventMaskButtonPressMask+                                                  , Gdk.EventMaskSmoothScrollMask+                                                  , Gdk.EventMaskKeyPressMask+                                                  ])++    -- scroll Events+    scrollE <- signalE1' drawingArea #scrollEvent $ \e ->+                     Gdk.getEventScrollDeltaY e++    -- events when we press a key+    keyPressedE <- signalE1' drawingArea #keyPressEvent $ \e -> do+                     v  <- Gdk.getEventKeyKeyval e+                     Gdk.keyvalName v+    -- events where we press an arrow key+    let arrowKeyE = filterJust . fmap (>>= toArrowKey) $ keyPressedE++    -- handle mouse clicks+    mousePressedE <- signalE1' drawingArea #buttonPressEvent $ \e -> do+                      x <- Gdk.getEventButtonX e+                      y <- Gdk.getEventButtonY e+                      return $! V2 x ((-1*y) + drawingAreaH)++    lastMousePressB <- stepper undefined mousePressedE++    -- mouse release+    -- mouseReleasedE <- signalE1' drawingArea #buttonReleaseEvent $ \e -> do+    --                     x <- Gdk.getEventButtonX e+    --                     y <- Gdk.getEventButtonY e+    --                     return $! V2 x ((-1*y) + drawingAreaH)++    -- mouse coordinates+    mouseMotionE <- signalE1' drawingArea #motionNotifyEvent $ \e -> do+                      x  <- Gdk.getEventMotionX e+                      y  <- Gdk.getEventMotionY e+                      st <- Gdk.getEventMotionState e+                      let !p = V2 x ((-1*y) + drawingAreaH)+                      return (p,st)+    mouseMotionB  <- stepper undefined mouseMotionE++    -- difference between the current mouse position and where we clicked last+    let dragOffsetB = (\p (q,_) -> p .-. q) <$> lastMousePressB <*> mouseMotionB++        -- sample the displacement vector whenever we are have a move event+        -- and the moude button is still on+        dragOffsetE = dragOffsetB+                   <@ filterE ((Gdk.ModifierTypeButton1Mask `elem`) . snd) mouseMotionE+++    zoomLevelB <- accumB 1 $ (\dy -> (+0.1*dy)) <$> scrollE+++    let lastPosE = unions [ (\k -> (.+^ 2 *^ toDirection k)) <$> arrowKeyE   -- key event+                          , (\v -> (.+^ v))                  <$> dragOffsetE -- drag event+                          ]+    viewPortPosB <- accumB (V2 (drawingAreaW/2) (drawingAreaH/2)) $ lastPosE++    let viewPortB = ViewPort <$> pure drawWorld'+                             <*> pure (V2 drawingAreaW drawingAreaH)+                             <*> viewPortPosB+                             <*> zoomLevelB+                             <*> pure 0++    -- draw everything+    draw drawingArea (mirrored drawingAreaH render <$> viewPortB)+    #showAll window++data ViewPort a = ViewPort { drawWorld             :: Canvas a+                           , screenSize            :: V2 Double+                           , clippwingWindowCenter :: V2 Double+                           , zoomLevel             :: Double+                           , rotation              :: Double+                           }++clippingWindow :: ViewPort a -> Canvas.Dim+clippingWindow (ViewPort _ (V2 w h) (V2 cx cy) z _) = let x = cx - z*w/2+                                                          y = cy - z*h/2+                                                      in Canvas.D x y (z*w) (z*h)++-- some drawing+drawWorld' :: Canvas ()+drawWorld' = do+    Canvas.background $ Canvas.gray 255+    Canvas.stroke $ Canvas.gray 0+    forM_ [1..20] $ \i ->+      forM_ [1..16] $ \j -> do+        Canvas.rect (Canvas.D (100*i) (100*j) 20 20)++-- | Mirror the canvas s.t. the bottom-left corner is the origin+mirrored       :: Double -> (a -> Canvas ()) -> a -> Canvas ()+mirrored h d x = do Canvas.scale     $ V2 1 (-1)+                    Canvas.translate $ V2 0 (-1*h)+                    d x++-- | Render the view+render    :: ViewPort a -> Canvas a+render vp = do+    let (Canvas.D x y _ _) = clippingWindow vp++    Canvas.scale $ V2 (1/zoomLevel vp) (1/zoomLevel vp) -- scale everything s.t. the+                                                        -- cippingWindow equals+                                                        -- the window size+    Canvas.translate $ V2 (-1*x) (-1*y) -- move screen to the origin+    drawWorld vp

@@ -12,7 +12,8 @@ import           Data.Semigroup import           Data.Semigroup.Foldable --- | \(O(n \log n)\) time ConvexHull using divide and conqueror.+-- | \(O(n \log n)\) time ConvexHull using divide and conqueror. The resulting polygon is+-- given in clockwise order. convexHull :: (Ord r, Num r)            => NonEmpty.NonEmpty (Point 2 r :+ p) -> ConvexPolygon p r convexHull = unMerge

@@ -49,8 +49,8 @@ hull'          :: (Ord r, Num r) => [Point 2 r :+ p] -> NonEmpty (Point 2 r :+ p) hull' (a:b:ps) = NonEmpty.fromList $ hull'' [b,a] ps   where-    hull'' h  []    = h-    hull'' h (p:ps) = hull'' (cleanMiddle (p:h)) ps+    hull'' h []      = h+    hull'' h (p:ps') = hull'' (cleanMiddle (p:h)) ps'      cleanMiddle h@[_,_]                         = h     cleanMiddle h@(z:y:x:rest)

@@ -11,9 +11,11 @@ import qualified Data.IntMap.Strict as IM import qualified Data.Map as M import qualified Data.Map.Strict as SM-import           Data.PlaneGraph+import qualified Data.PlaneGraph  as PG+import qualified Data.PlanarGraph as PPG import qualified Data.Vector as V + --------------------------------------------------------------------------------  -- We store all adjacency lists in clockwise order@@ -80,47 +82,18 @@  -- | convert the triangulation into a planarsubdivision ----- running time: \(O(n\log n)\).-toPlanarSubdivision :: proxy s -> Triangulation p r -> PlanarSubdivision s p () () r-toPlanarSubdivision px tr = PlanarSubdivision g-  where-    g = toPlaneGraph px tr & vertexData.traverse  %~ (\(v :+ e) -> VertexData v e)-                           & dartData.traverse._2 %~ EdgeData Visible-                           & faceData.traverse    %~ FaceData []+-- running time: \(O(n)\).+toPlanarSubdivision    :: (Ord r, Fractional r)+                       => proxy s -> Triangulation p r -> PlanarSubdivision s p () () r+toPlanarSubdivision px = fromPlaneGraph . toPlaneGraph px  -- | convert the triangulation into a plane graph ----- running time: \(O(n\log n)\).+-- running time: \(O(n)\). toPlaneGraph    :: forall proxy s p r.-                   proxy s -> Triangulation p r -> PlaneGraph s Primal_ p () () r-toPlaneGraph _ tr = g & vertexData .~ tr^.positions+                   proxy s -> Triangulation p r -> PG.PlaneGraph s p () () r+toPlaneGraph _ tr = PG.PlaneGraph $ g&PPG.vertexData .~ vtxData   where-    g       = fromAdjacencyLists . V.toList . V.imap f $ tr^.neighbours+    g       = PG.fromAdjacencyLists . V.toList . V.imap f $ tr^.neighbours     f i adj = (VertexId i, VertexId <$> adj)---  -- (planarGraph' . P.toCycleRep n $ perm)&vertexData .~ tr^.positions-  -- where-  --   neighs    = C.rightElements <$> tr^.neighbours-  --   n         = sum . fmap length $ neighs--  --   vtxIDs = [0..]-  --   perm = trd' . foldr toOrbit (ST mempty 0 mempty) $ zip vtxIDs (V.toList neighs)--  --   -- | Given a vertex with its adjacent vertices (u,vs) (in CCW order) convert this-  --   -- vertex with its adjacent vertices into an Orbit-  --   toOrbit                     :: (VertexID,[VertexID]) -> ST' [[Dart s]]-  --                               -> ST' [[Dart s]]-  --   toOrbit (u,vs) (ST m a dss) =-  --     let (ST m' a' ds') = foldr (toDart . (u,)) (ST m a mempty) vs-  --     in ST m' a' (ds':dss)---  --   -- | Given an edge (u,v) and a triplet (m,a,ds) we construct a new dart-  --   -- representing this edge.-  --   toDart                   :: (VertexID,VertexID) -> ST' [Dart s] -> ST' [Dart s]-  --   toDart (u,v) (ST m a ds) = let dir = if u < v then Positive else Negative-  --                                  t'  = (min u v, max u v)-  --                              in case M.lookup t' m of-  --     Just a' -> ST m                  a     (Dart (Arc a') dir : ds)-  --     Nothing -> ST (SM.insert t' a m) (a+1) (Dart (Arc a)  dir : ds)+    vtxData = (\(loc :+ p) -> VertexData loc p) <$> tr^.positions

@@ -25,14 +25,14 @@ -- running time: \(O(n \log n)\) euclideanMST     :: (Ord r, Fractional r)                  => NonEmpty.NonEmpty (Point 2 r :+ p) -> Tree (Point 2 r :+ p)-euclideanMST pts = (\v -> g^.vDataOf v) <$> t+euclideanMST pts = (\v -> g^.locationOf v :+ g^.dataOf v) <$> t   where     -- since we care only about the relative order of the edges we can use the     -- squared Euclidean distance rather than the Euclidean distance, thus     -- avoiding the Floating constraint     g = withEdgeDistances squaredEuclideanDist . toPlaneGraph (Proxy :: Proxy MSTW)       . delaunayTriangulation $ pts-    t = mst g+    t = mst $ g^.graph   data MSTW

@@ -0,0 +1,16 @@+module Algorithms.Geometry.LineSegmentIntersection where++import qualified Algorithms.Geometry.LineSegmentIntersection.BentleyOttmann as BO+import           Data.Geometry.LineSegment+import           Data.Geometry.Polygon++-- Tests if there are any interior intersections.+--+-- | \(O(n \log n)\)+hasInteriorIntersections :: (Ord r, Fractional r)+                         => [LineSegment 2 p r] -> Bool+hasInteriorIntersections = not . null . BO.interiorIntersections++-- | \(O(n \log n)\)+hasSelfIntersections :: (Ord r, Fractional r) => Polygon t p r -> Bool+hasSelfIntersections = hasInteriorIntersections . listEdges

@@ -1,12 +1,11 @@+{-# LANGUAGE ScopedTypeVariables #-} module Algorithms.Geometry.LineSegmentIntersection.BentleyOttmann where  import           Algorithms.Geometry.LineSegmentIntersection.Types import           Control.Lens hiding (contains)-import qualified Data.BalBST as SS -- status struct import           Data.Ext-import           Data.Function (on)+import qualified Data.Foldable as F import           Data.Geometry.Interval-import           Data.Geometry.Line import           Data.Geometry.LineSegment import           Data.Geometry.Point import           Data.Geometry.Properties@@ -16,12 +15,12 @@ import qualified Data.Map as M import           Data.Maybe import           Data.Ord (Down(..), comparing)+import           Data.OrdSeq (Compare)+import qualified Data.OrdSeq as SS -- status struct import           Data.Semigroup import qualified Data.Set as EQ -- event queue import           Data.Vinyl-import           Frames.CoRec--import           Debug.Trace+import           Data.Vinyl.CoRec  -------------------------------------------------------------------------------- @@ -30,12 +29,23 @@ -- \(O((n+k)\log n)\), where \(k\) is the number of intersections. intersections    :: (Ord r, Fractional r)                  => [LineSegment 2 p r] -> Intersections p r-intersections ss = merge $ sweep pts (SS.empty $ ordAtNav undefined)+intersections ss = merge $ sweep pts mempty   where-    pts = EQ.fromAscList . groupStarts . L.sort . concatMap f $ ss-    f s = let [p,q] = L.sortBy ordPoints [s^.start.core,s^.end.core]-          in [Event p (Start $ s :| []), Event q (End s)]+    pts = EQ.fromAscList . groupStarts . L.sort . concatMap asEventPts $ ss +-- | Computes all intersection points p s.t. p lies in the interior of at least+-- one of the segments.+--+--  \(O((n+k)\log n)\), where \(k\) is the number of intersections.+interiorIntersections :: (Ord r, Fractional r)+                       => [LineSegment 2 p r] -> Intersections p r+interiorIntersections = M.filter (not . isEndPointIntersection) . intersections++-- | Computes the event points for a given line segment+asEventPts   :: Ord r => LineSegment 2 p r -> [Event p r]+asEventPts s = let [p,q] = L.sortBy ordPoints [s^.start.core,s^.end.core]+               in [Event p (Start $ s :| []), Event q (End s)]+ -- | Group the segments with the intersection points merge :: Ord r =>  [IntersectionPoint p r] -> Intersections p r merge = foldr (\(IntersectionPoint p a) -> M.insertWith (<>) p a) M.empty@@ -47,14 +57,12 @@ groupStarts (Event p (Start s) : es) = Event p (Start ss) : groupStarts rest   where     (ss',rest) = L.span sameStart es-     -- sort the segs on lower endpoint     ss         = let (x:|xs) = s in x :| (xs ++ concatMap startSegs ss')      sameStart (Event q (Start _)) = p == q     sameStart _                   = False-groupStarts (e : es)           = e : groupStarts es-+groupStarts (e : es)                 = e : groupStarts es  -------------------------------------------------------------------------------- -- * Data type for Events@@ -89,7 +97,6 @@ ordPoints     :: Ord r => Point 2 r -> Point 2 r -> Ordering ordPoints a b = let f p = (Down $ p^.yCoord, p^.xCoord) in comparing f a b - -- | Get the segments that start at the given event point startSegs   :: Event p r -> [LineSegment 2 p r] startSegs e = case eventType e of@@ -98,24 +105,30 @@  -------------------------------------------------------------------------------- --- | The navigator that we use that orders the segments that intersect at a--- horizontal line (from left to right)-ordAtNav   :: (Ord r, Fractional r) => r -> SS.TreeNavigator r (LineSegment 2 p r)-ordAtNav y = SS.Nav (\s x -> h s <= x) (min `on` h)-  where-    h s = match (s `intersect` horizontalLine y) $-         (H $ \NoIntersection -> error "ordAtNav: No intersection")-      :& (H $ \p              -> p^.xCoord)-      :& (H $ \_              -> rightEndpoint s) -- the intersection is s itself-      :& RNil+-- | Compare based on the x-coordinate of the intersection with the horizontal+-- line through y+ordAt   :: (Fractional r, Ord r) => r -> Compare (LineSegment 2 p r)+ordAt y = comparing (xCoordAt y) +-- | Given a y coord and a line segment that intersects the horizontal line+-- through y, compute the x-coordinate of this intersection point.+--+-- note that we will pretend that the line segment is closed, even if it is not+xCoordAt             :: (Fractional r, Ord r) => r -> LineSegment 2 p r -> r+xCoordAt y (LineSegment' (Point2 px py :+ _) (Point2 qx qy :+ _))+      | py == qy     = px `max` qx  -- s is horizontal, and since it by the+                                    -- precondition it intersects the sweep+                                    -- line, we return the x-coord of the+                                    -- rightmost endpoint.+      | otherwise    = px + alpha * (qx - px)+  where+    alpha = (y - py) / (qy - py)  -------------------------------------------------------------------------------- -- * The Main Sweep -type EventQueue p r = EQ.Set (Event p r)--type StatusStructure p r = SS.BalBST r (LineSegment 2 p r)+type EventQueue      p r = EQ.Set (Event p r)+type StatusStructure p r = SS.OrdSeq (LineSegment 2 p r)  -- | Run the sweep handling all events sweep       :: (Ord r, Fractional r)@@ -124,8 +137,13 @@     Nothing      -> []     Just (e,eq') -> handle e eq' ss +isClosedStart                     :: Eq r => Point 2 r -> LineSegment 2 p r -> Bool+isClosedStart p (LineSegment s e)+  | p == s^.unEndPoint.core       = isClosed s+  | otherwise                     = isClosed e+ -- | Handle an event point-handle                           :: (Ord r, Fractional r)+handle                           :: forall r p. (Ord r, Fractional r)                                  => Event p r -> EventQueue p r -> StatusStructure p r                                  -> [IntersectionPoint p r] handle e@(eventPoint -> p) eq ss = toReport <> sweep eq' ss'@@ -133,22 +151,20 @@     starts                   = startSegs e     (before,contains',after) = extractContains p ss     (ends,contains)          = L.partition (endsAt p) contains'---    toReport = case starts ++ contains' of-                 (_:_:_) -> [IntersectionPoint p $ associated (starts <> ends) contains]+    -- starting segments, exluding those that have an open starting point+    starts'  = filter (isClosedStart p) starts+    toReport = case starts' ++ contains' of+                 (_:_:_) -> [IntersectionPoint p $ associated (starts' <> ends) contains]                  _       -> []      -- new status structure-    ss' = before `SS.join` newSegs `SS.join` after-+    ss' = before <> newSegs <> after     newSegs = toStatusStruct p $ starts ++ contains      -- the new eeventqueue     eq' = foldr EQ.insert eq es-     -- the new events:-    es | SS.null newSegs = maybeToList $ app (findNewEvent p) sl sr+    es | F.null newSegs  = maybeToList $ app (findNewEvent p) sl sr        | otherwise       = let s'  = fst <$> SS.minView newSegs                                s'' = fst <$> SS.maxView newSegs                            in catMaybes [ app (findNewEvent p) sl  s'@@ -164,24 +180,22 @@ extractContains      :: (Fractional r, Ord r)                      => Point 2 r -> StatusStructure p r                      -> (StatusStructure p r, [LineSegment 2 p r], StatusStructure p r)-extractContains p ss = (before, mid1 ++ mid2, after)+extractContains p ss = (before, F.toList $ mid1 <> mid2, after)   where-    n = ordAtNav (p^.yCoord)-    SS.Split before (mid1,mid2) after = SS.splitExtract pred' sel $ ss { SS.nav = n}--    pred' s = not $ SS.goLeft n s (p^.xCoord)-    sel   s = p `onSegment` s-+    (before, mid1, after') = SS.splitOn (xCoordAt $ p^.yCoord) (p^.xCoord) ss+    -- Make sure to also select the horizontal segments containing p+    (mid2, after) = SS.splitMonotonic (\s -> not $ p `onSegment` s) after'  -- | Given a point and the linesegements that contain it. Create a piece of -- status structure for it. toStatusStruct      :: (Fractional r, Ord r)                     => Point 2 r -> [LineSegment 2 p r] -> StatusStructure p r-toStatusStruct p xs = ss { SS.nav = ordAtNav $ p^.yCoord } `SS.join` hors+toStatusStruct p xs = ss <> hors+  -- ss { SS.nav = ordAtNav $ p^.yCoord } `SS.join` hors   where     (hors',rest) = L.partition isHorizontal xs-    ss           = SS.fromList (ordAtNav $ maxY xs) rest-    hors         = SS.fromList (SS.ordNavBy rightEndpoint) hors'+    ss           = SS.fromListBy (ordAt $ maxY xs) rest+    hors         = SS.fromListBy (comparing rightEndpoint) hors'      isHorizontal s  = s^.start.core.yCoord == s^.end.core.yCoord 

@@ -11,7 +11,7 @@ import qualified Data.Map as M import           Data.Semigroup import           Data.Vinyl-import           Frames.CoRec+import           Data.Vinyl.CoRec   -- | Compute all intersections (naively)

@@ -26,6 +26,7 @@                                  , _interiorTo        :: Set' (LineSegment 2 p r)                                  } deriving (Show) + instance (Eq p, Eq r) => Eq (Associated p r) where   (Associated es is) == (Associated es' is') = f es es' && f is is'     where@@ -65,6 +66,13 @@                     } deriving (Show,Eq) makeLenses ''IntersectionPoint ++-- | reports true if there is at least one segment for which this intersection+-- point is interior.+--+-- \(O(1)\)+isEndPointIntersection :: Associated p r -> Bool+isEndPointIntersection = Map.null . _interiorTo   -- newtype E a b = E (a -> b)

@@ -0,0 +1,49 @@+{-# LANGUAGE ScopedTypeVariables #-}+module Algorithms.Geometry.LowerEnvelope.DualCH where++import Data.Maybe(fromJust)+import Control.Lens((^.))+import Data.Ext+import Data.Geometry+import Algorithms.Geometry.ConvexHull.GrahamScan+import Data.List.NonEmpty(NonEmpty(..))+import qualified Data.List.NonEmpty as NonEmpty+import Data.Geometry.Duality+import Data.Proxy+import Data.Vinyl.CoRec++type Envelope a r = NonEmpty (Line 2 r :+ a)++-- | Given a list of non-vertical lines, computes the lower envelope using+-- duality.+--+-- \(O(n\log n)\)+lowerEnvelope :: (Ord r, Fractional r) => NonEmpty (Line 2 r :+ a) -> Envelope a r+lowerEnvelope = lowerEnvelopeWith upperHull+++type UpperHullAlgorithm a r = NonEmpty (Point 2 r :+ a) -> NonEmpty (Point 2 r :+ a)++-- | Given a list of non-vertical lines, computes the lower envelope by computing+-- the upper convex hull. It uses the given algorithm to do so+--+-- running time: O(time required by the given upper hull algorithm)+lowerEnvelopeWith        :: (Fractional r, Eq r)+                         => UpperHullAlgorithm (Line 2 r :+ a) r+                         -> NonEmpty (Line 2 r :+ a) -> Envelope a r+lowerEnvelopeWith chAlgo = fromPts . chAlgo . toPts+  where+    toPts   = fmap (\l -> dualPoint' (l^.core) :+ l)+    fromPts = fmap (^.extra)++-- | Computes the vertices of the envelope, in left to right order+vertices   :: (Ord r, Fractional r) => Envelope a r -> [Point 2 r :+ (a,a)]+vertices e = zipWith intersect' (NonEmpty.toList e) (NonEmpty.tail e)+++-- | Given two non-parallel lines, compute the intersection point and+-- return the pair of a's associated with the lines+intersect'                     :: forall r a. (Ord r, Fractional r)+                               => Line 2 r :+ a -> Line 2 r :+ a -> Point 2 r :+ (a,a)+intersect' (l :+ le) (r :+ re) = (:+ (le,re)) . fromJust+                               . asA (Proxy :: Proxy (Point 2 r)) $ l `intersect` r

@@ -0,0 +1,310 @@+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE ScopedTypeVariables #-}+module Algorithms.Geometry.PolygonTriangulation.MakeMonotone where++import           Algorithms.Geometry.LineSegmentIntersection.BentleyOttmann ( xCoordAt+                                                                            , ordAt)+import           Algorithms.Geometry.PolygonTriangulation.Types+import           Control.Lens+import           Control.Monad (forM_, when)+import           Control.Monad.Reader+import           Control.Monad.State.Strict+import           Control.Monad.Writer (WriterT, execWriterT,tell)+import           Data.Bifunctor+import           Data.CircularSeq (rotateL, rotateR, zip3LWith)+import qualified Data.DList as DList+import           Data.Ext+import qualified Data.Foldable as F+import           Data.Geometry.LineSegment+import           Data.Geometry.PlanarSubdivision.Basic+import           Data.Geometry.Point+import           Data.Geometry.Polygon+import qualified Data.IntMap as IntMap+import qualified Data.List.NonEmpty as NonEmpty+import           Data.Ord (comparing, Down(..))+import           Data.OrdSeq (OrdSeq)+import qualified Data.OrdSeq as SS+import           Data.Semigroup+import           Data.Util+import qualified Data.Vector as V+import qualified Data.Vector.Mutable as MV+++----------------------------------------------------------------------------------++data VertexType = Start | Merge | Split | End | Regular deriving (Show,Read,Eq)+++-- How about the hole vertices?++-- | assigns a vertex type to each vertex+--+-- pre: the polygon is given in CCW order+--+-- running time: \(O(n)\).+classifyVertices                     :: (Num r, Ord r)+                                     => Polygon t p r+                                     -> Polygon t (p :+ VertexType) r+classifyVertices p@(SimplePolygon _) = classifyVertices' p+classifyVertices (MultiPolygon vs h) = MultiPolygon vs' h'+  where+    (SimplePolygon vs') = classifyVertices' $ SimplePolygon vs+    h' = map (first (&extra %~ onHole) . classifyVertices') h++    -- the roles on hole vertices are slightly different+    onHole Start   = Split+    onHole Merge   = End+    onHole Split   = Start+    onHole End     = Merge+    onHole Regular = Regular++-- | assigns a vertex type to each vertex+--+-- pre: the polygon is given in CCW order+--+-- running time: \(O(n)\).+classifyVertices'                    :: (Num r, Ord r)+                                     => SimplePolygon p r+                                     -> SimplePolygon (p :+ VertexType) r+classifyVertices' (SimplePolygon vs) =+    SimplePolygon $ zip3LWith f (rotateL vs) vs (rotateR vs)+  where+    -- is the angle larger than > 180 degrees+    largeInteriorAngle p c n = case ccw (p^.core) (c^.core) (n^.core) of+           CCW -> False+           CW  -> True+           _   -> error "classifyVertices -> largeInteriorAngle: colinear points"++    f p c n = c&extra %~ (:+ vt)+      where+        vt = case (p `cmpSweep` c, n `cmpSweep` c, largeInteriorAngle p c n) of+               (LT, LT, False) -> Start+               (LT, LT, True)  -> Split+               (GT, GT, False) -> End+               (GT, GT, True)  -> Merge+               _               -> Regular++++-- | p < q = p.y < q.y || p.y == q.y && p.x > q.y+cmpSweep :: Ord r => Point 2 r :+ e -> Point 2 r :+ e -> Ordering+p `cmpSweep` q =+  comparing (^.core.yCoord) p q <> comparing (Down . (^.core.xCoord)) p q+++--------------------------------------------------------------------------------++type Event r = Point 2 r :+ (Two (LineSegment 2 Int r))++data StatusStruct r = SS { _statusStruct :: !(SS.OrdSeq (LineSegment 2 Int r))+                         , _helper       :: !(IntMap.IntMap Int)+                         -- ^ for every e_i, the id of the helper vertex+                         } deriving (Show)+makeLenses ''StatusStruct++ix'   :: Int -> Lens' (V.Vector a) a+ix' i = singular (ix i)++-- | Given a polygon, find a set of non-intersecting diagonals that partition+-- the polygon into y-monotone pieces.+--+-- running time: \(O(n\log n)\)+computeDiagonals    :: forall t r p. (Fractional r, Ord r)+                    => Polygon t p r -> [LineSegment 2 p r]+computeDiagonals p' = map f . sweep+                    . NonEmpty.sortBy (flip cmpSweep)+                    . polygonVertices . withIncidentEdges+                    . first (^._1) $ pg+  where+    -- remaps to get the p value rather than the vertexId+    f = first (\i -> vertexInfo^.ix' i._2)++    pg :: Polygon t (SP Int (p :+ VertexType)) r+    pg = numberVertices . classifyVertices . toCounterClockWiseOrder $ p'+    vertexInfo :: V.Vector (STR (Point 2 r) p VertexType)+    vertexInfo = let vs = polygonVertices pg+                     n  = F.length vs+                 in V.create $ do+                   v <- MV.new n+                   forM_ vs $ \(pt :+ SP i (p :+ vt)) ->+                     MV.write v i (STR pt p vt)+                   return v++    initialSS = SS mempty mempty++    sweep  es = flip runReader vertexInfo $ evalStateT (sweep' es) initialSS+    sweep' es = DList.toList <$> execWriterT (sweep'' es)++    sweep'' :: NonEmpty.NonEmpty (Event r) -> Sweep p r ()+    sweep'' = mapM_ handle++-- | Computes a set of diagionals that decompose the polygon into y-monotone+-- pieces.+--+-- running time: \(O(n\log n)\)+makeMonotone      :: (Fractional r, Ord r)+                  => proxy s -> Polygon t p r+                  -> PlanarSubdivision s p PolygonEdgeType PolygonFaceData r+makeMonotone px pg = let (e:es) = listEdges pg+                     in constructSubdivision px e es (computeDiagonals pg)++type Sweep p r = WriterT (DList.DList (LineSegment 2 Int r))+                   (StateT (StatusStruct r)+                     (Reader (V.Vector (VertexInfo p r))))++type VertexInfo p r = STR (Point 2 r) p VertexType+++tell' :: LineSegment 2 Int r -> Sweep p r ()+tell' = tell . DList.singleton++getIdx :: Event r -> Int+getIdx = view (extra._1.end.extra)++getVertexType   :: Int -> Sweep p r VertexType+getVertexType v = asks (^.ix' v._3)++getEventType :: Event r -> Sweep p r VertexType+getEventType = getVertexType . getIdx++handle   :: (Fractional r, Ord r) => Event r -> Sweep p r ()+handle e = let i = getIdx e in getEventType e >>= \case+    Start   -> handleStart   i e+    End     -> handleEnd     i e+    Split   -> handleSplit   i e+    Merge   -> handleMerge   i e+    Regular | isLeftVertex i e -> handleRegularL i e+            | otherwise        -> handleRegularR i e+++insertAt   :: (Ord r, Fractional r) => Point 2 r -> LineSegment 2 q r+           -> OrdSeq (LineSegment 2 q r) -> OrdSeq (LineSegment 2 q r)+insertAt v = SS.insertBy (ordAt $ v^.yCoord)++deleteAt   :: (Fractional r, Ord r) => Point 2 r -> LineSegment 2 p r+           -> OrdSeq (LineSegment 2 p r) -> OrdSeq (LineSegment 2 p r)+deleteAt v = SS.deleteAllBy (ordAt $ v^.yCoord)+++handleStart              :: (Fractional r, Ord r)+                         => Int -> Event r -> Sweep p r ()+handleStart i (v :+ adj) = modify $ \(SS t h) ->+                                SS (insertAt v (adj^._2) t)+                                   (IntMap.insert i i h)++handleEnd              :: (Fractional r, Ord r)+                       => Int -> Event r -> Sweep p r ()+handleEnd i (v :+ adj) = do let iPred = adj^._1.start.extra  -- i-1+                            -- lookup p's helper; if it is a merge vertex+                            -- we insert a new segment+                            tellIfMerge i v iPred+                            -- delete e_{i-1} from the status struct+                            modify $ \ss ->+                              ss&statusStruct %~ deleteAt v (adj^._1)++-- | Adds edge (i,j) if e_j's helper is a merge vertex+tellIfMerge       :: Int -> Point 2 r -> Int -> Sweep p r ()+tellIfMerge i v j = do SP u ut <- getHelper j+                       when (ut == Merge) (tell' $ ClosedLineSegment (v :+ i) u)++-- | Get the helper of edge i, and its vertex type+getHelper   :: Int -> Sweep p r (SP (Point 2 r :+ Int) VertexType)+getHelper i = do Just ui    <- gets (^.helper.at i)+                 STR u _ ut <- asks (^.ix' ui)+                 pure $ SP (u :+ ui) ut+++lookupLE     :: (Ord r, Fractional r)+             => Point 2 r -> OrdSeq (LineSegment 2 Int r)+             -> Maybe (LineSegment 2 Int r)+lookupLE v s = let (l,m,_) = SS.splitOn (xCoordAt $ v^.yCoord) (v^.xCoord) s+               in SS.lookupMax (l <> m)+++handleSplit              :: (Fractional r, Ord r) => Int -> Event r -> Sweep p r ()+handleSplit i (v :+ adj) = do Just ej <- gets $ \ss -> ss^.statusStruct.to (lookupLE v)+                              let j = ej^.start.extra+                              SP u _ <- getHelper j+                              -- update the status struct:+                              -- insert the new edge into the status Struct and+                              -- set the helper of e_j to be v_i+                              modify $ \(SS t h) ->+                                SS (insertAt v (adj^._2) t)+                                   (IntMap.insert i i . IntMap.insert j i $ h)+                              -- return the diagonal+                              tell' $ ClosedLineSegment (v :+ i) u++handleMerge              :: (Fractional r, Ord r) => Int -> Event r -> Sweep p r ()+handleMerge i (v :+ adj) = do let ePred = adj^._1.start.extra -- i-1+                              tellIfMerge i v ePred+                              -- delete e_{i-1} from the status struct+                              modify $ \ss -> ss&statusStruct %~ deleteAt v (adj^._1)+                              connectToLeft i v++-- | finds the edge j to the left of v_i, and connect v_i to it if the helper+-- of j is a merge vertex+connectToLeft     :: (Fractional r, Ord r) => Int -> Point 2 r -> Sweep p r ()+connectToLeft i v = do Just ej <- gets $ \ss -> ss^.statusStruct.to (lookupLE v)+                       let j = ej^.start.extra+                       tellIfMerge i v j+                       modify $ \ss -> ss&helper %~ IntMap.insert j i++-- | returns True if v the interior of the polygon is to the right of v+isLeftVertex              :: Ord r => Int -> Event r -> Bool+isLeftVertex i (v :+ adj) = case (adj^._1.start) `cmpSweep` (v :+ i) of+                              GT -> True+                              _  -> False+  -- if the predecessor occurs before the sweep, this must be a left vertex++handleRegularL              :: (Fractional r, Ord r) => Int -> Event r -> Sweep p r ()+handleRegularL i (v :+ adj) = do let ePred = adj^._1.start.extra -- i-1+                                 tellIfMerge i v ePred+                                 -- delete e_{i-1} from the status struct+                                 modify $ \ss ->+                                   ss&statusStruct %~ deleteAt v (adj^._1)+                                 -- insert a e_i in the status struct, and set its helper+                                 -- to be v_i+                                 modify $ \(SS t h) ->+                                     SS (insertAt v (adj^._2) t)+                                        (IntMap.insert i i h)++handleRegularR            :: (Fractional r, Ord r) => Int -> Event r -> Sweep p r ()+handleRegularR i (v :+ _) = connectToLeft i v+++++--------------------------------------------------------------------------------+++-- testPolygon :: SimplePolygon Int Rational+-- testPolygon = fromPoints [ point2 20 20 :+ 1+--                          , point2 18 19 :+ 2+--                          , point2 16 25 :+ 3+--                          , point2 13 23 :+ 4+--                          , point2 10 24 :+ 5+--                          , point2 6  22 :+ 6+--                          , point2 8  21 :+ 7+--                          , point2 7  18 :+ 8+--                          , point2 2  19 :+ 9+--                          , point2 1  10 :+ 10+--                          , point2 3  5  :+ 11+--                          , point2 11 7  :+ 12+--                          , point2 15 1  :+ 13+--                          , point2 12 15 :+ 14+--                          , point2 15 12 :+ 15+--                          ]++-- vertexTypes = [Start,Merge,Start,Merge,Start,Regular,Regular,Merge,Start,Regular,End,Split,End,Split,End]+++-- loadT = do pgs <- readAllFrom "/Users/frank/tmp/testPoly.ipe"+--                         :: IO [SimplePolygon () Rational :+ IpeAttributes Path Rational]+--            mapM_ print pgs+--            let diags = map (computeDiagonals . (^.core)) pgs+--                f = asIpeGroup . map (asIpeObject' mempty)+--                out = [ asIpeGroup $ map (\(pg :+ a) -> asIpeObject pg a) pgs+--                      , asIpeGroup $ map f diags+--                      ]+--                outFile = "/Users/frank/tmp/out.ipe"+--            writeIpeFile outFile . singlePageFromContent $ out

@@ -0,0 +1,70 @@+module Algorithms.Geometry.PolygonTriangulation.Triangulate where+++import qualified Algorithms.Geometry.PolygonTriangulation.MakeMonotone as MM+import qualified Algorithms.Geometry.PolygonTriangulation.TriangulateMonotone as TM+import           Algorithms.Geometry.PolygonTriangulation.Types+import           Control.Lens+import           Data.Either (lefts)+import           Data.Ext+import qualified Data.Foldable as F+import           Data.Geometry.LineSegment+import           Data.Geometry.PlanarSubdivision.Basic+import           Data.Geometry.Polygon+import           Data.PlaneGraph (PlaneGraph)+import           Data.Semigroup++-- | Triangulates a polygon of \(n\) vertices+--+-- running time: \(O(n \log n)\)+triangulate        :: (Ord r, Fractional r)+                   => proxy s -> Polygon t p r+                   -> PlanarSubdivision s p PolygonEdgeType PolygonFaceData r+triangulate px pg' = constructSubdivision px e es diags+  where+    (pg, diags)   = computeDiagonals' pg'+    (e:es)        = listEdges pg+++-- | Triangulates a polygon of \(n\) vertices+--+-- running time: \(O(n \log n)\)+triangulate'        :: (Ord r, Fractional r)+                   => proxy s -> Polygon t p r+                   -> PlaneGraph s p PolygonEdgeType PolygonFaceData r+triangulate' px pg' = constructGraph px e es diags+  where+    (pg, diags)   = computeDiagonals' pg'+    (e:es)        = listEdges pg+++-- | Computes a set of diagaonals that together triangulate the input polygon+-- of \(n\) vertices.+--+-- running time: \(O(n \log n)\)+computeDiagonals :: (Ord r, Fractional r) => Polygon t p r -> [LineSegment 2 p r]+computeDiagonals = snd . computeDiagonals'++-- | Computes a set of diagaonals that together triangulate the input polygon+-- of \(n\) vertices. Returns a copy of the input polygon, whose boundaries are+-- oriented in counter clockwise order, as well.+--+-- running time: \(O(n \log n)\)+computeDiagonals'     :: (Ord r, Fractional r)+                      => Polygon t p r -> (Polygon t p r, [LineSegment 2 p r])+computeDiagonals' pg' = (pg, monotoneDiags <> extraDiags)+  where+    pg            = toCounterClockWiseOrder pg'+    monotoneP     = MM.makeMonotone (Identity pg') pg -- use some arbitrary proxy type+    -- outerFaceId'  = outerFaceId monotoneP++    monotoneDiags = map (^._2.core) . filter (\e' -> e'^._2.extra == Diagonal)+                  . F.toList . edgeSegments $ monotoneP+    extraDiags    = concatMap (TM.computeDiagonals . toCounterClockWiseOrder')+                  . lefts . map (^._2.core)+                  -- . 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

@@ -0,0 +1,178 @@+module Algorithms.Geometry.PolygonTriangulation.TriangulateMonotone where++import           Control.Lens+import           Data.Bifunctor+import qualified Data.CircularSeq as C+import           Data.Ext+import qualified Data.Foldable as F+import           Data.Geometry.LineSegment+import           Data.Geometry.Point+import           Data.Geometry.Polygon+import qualified Data.List as L+import           Data.Ord (comparing, Down(..))+import           Data.Semigroup+import           Data.Util+import           Algorithms.Geometry.PolygonTriangulation.Types+import           Data.PlaneGraph (PlaneGraph)+import           Data.Geometry.PlanarSubdivision.Basic(PolygonFaceData, PlanarSubdivision)++--------------------------------------------------------------------------------++--+type MonotonePolygon p r = SimplePolygon p r++data LR = L | R deriving (Show,Eq)++-- | Triangulates a polygon of \(n\) vertices+--+-- running time: \(O(n \log n)\)+triangulate        :: (Ord r, Fractional r)+                   => proxy s -> MonotonePolygon p r+                   -> PlanarSubdivision s p PolygonEdgeType PolygonFaceData r+triangulate px pg' = constructSubdivision px e es (computeDiagonals pg)+  where+    pg     = toCounterClockWiseOrder pg'+    (e:es) = listEdges pg+  -- TODO: Find a way to construct the graph in O(n) time.++-- | Triangulates a polygon of \(n\) vertices+--+-- running time: \(O(n \log n)\)+triangulate'        :: (Ord r, Fractional r)+                    => proxy s -> MonotonePolygon p r+                    -> PlaneGraph s p PolygonEdgeType PolygonFaceData r+triangulate' px pg' = constructGraph px e es (computeDiagonals pg)+  where+    pg     = toCounterClockWiseOrder pg'+    (e:es) = listEdges pg+  -- TODO: Find a way to construct the graph in O(n) time.+++-- | Given a y-monotone polygon in counter clockwise order computes the diagonals+-- to add to triangulate the polygon+--+-- pre: the input polygon is y-monotone and has \(n \geq 3\) vertices+--+-- running time: \(O(n)\)+computeDiagonals    :: (Ord r, Num r)+                    => MonotonePolygon p r -> [LineSegment 2 p r]+computeDiagonals pg = diags'' <> diags'+  where+    -- | run the stack computation+    SP (_:stack') diags' = L.foldl' (\(SP stack acc) v' -> (<> acc) <$> process v' stack)+                                    (SP [v,u] []) vs'+    -- add vertices from the last guy w to all 'middle' guys of the final stack+    diags'' = map (seg w) $ init stack'+    -- extract the last vertex+    Just (vs',w) = unsnoc vs+    -- merge the two lists into one list for procerssing+    (u:v:vs) = uncurry (mergeBy $ comparing (\(Point2 x y :+ _) -> (Down y, x)))+             $ splitPolygon pg+++type P p r = Point 2 r :+ (LR :+ p)++type Stack a = [a]+++++-- type Scan p r = State (Stack (P p r))++chainOf :: P p r -> LR+chainOf = (^.extra.core)++toVtx :: P p r -> Point 2 r :+ p+toVtx = (&extra %~ (^.extra))++seg     :: P p r -> P p r -> LineSegment 2 p r+seg u v = ClosedLineSegment (toVtx u) (toVtx v)++process                    :: (Ord r, Num r)+                           => P p r -> Stack (P p r)+                           -> SP (Stack (P p r)) [LineSegment 2 p r]+process _ []               = error "TriangulateMonotone.process: absurd. empty stack"+process v stack@(u:ws)+  | chainOf v /= chainOf u = SP [v,u]      (map (seg v) . init $ stack)+  | otherwise              = SP (v:w:rest) (map (seg v) popped)+      where+        (popped,rest) = bimap (map fst) (map fst) . L.span (isInside v) $ zip ws stack+        w             = last $ u:popped+++-- | test if m does not block the line segment from v to u+isInside          :: (Ord r, Num r) => P p r -> (P p r, P p r) -> Bool+isInside v (u, m) = case ccw' v m u of+                     CoLinear -> False+                     CCW      -> chainOf v == R+                     CW       -> chainOf v == L++-- | given a comparison function, merge the two ordered lists+mergeBy     :: (a -> a -> Ordering) -> [a] -> [a] -> [a]+mergeBy cmp = go+  where+    go []     ys     = ys+    go xs     []     = xs+    go (x:xs) (y:ys) = case x `cmp` y of+                         GT -> y : go (x:xs) ys+                         _  -> x : go xs     (y:ys)+++-- | When the polygon is in counter clockwise order we return (leftChain,rightChain)+-- ordered from the top-down.+--+-- if there are multiple points with the maximum yCoord we pick the rightmost one,+-- if there are multiple point with the minimum yCoord we pick the leftmost one.+--+-- running time: \(O(n)\)+splitPolygon    :: Ord r => MonotonePolygon p r+                -> ([Point 2 r :+ (LR :+ p)], [Point 2 r :+ (LR :+ p)])+splitPolygon pg = bimap (f L) (f R)+                . second reverse+                . L.break (\v -> v^.core == vMinY)+                . F.toList . C.rightElements $ vs'+  where+    f x = map (&extra %~ (x :+))+    -- rotates the list to the vtx with max ycoord+    Just vs' = C.findRotateTo (\v -> v^.core == vMaxY)+             $ pg^.outerBoundary+    vMaxY = getY F.maximumBy+    vMinY = getY F.minimumBy+    swap' (Point2 x y) = Point2 y x+    getY ff = let p = ff (comparing (^.core.to swap')) $ pg^.outerBoundary+              in p^.core++++--------------------------------------------------------------------------------++-- testPolygon = fromPoints . map ext $ [ point2 10 10+--                                      , point2 5 20+--                                      , point2 3 14+--                                      , point2 1 1+--                                      , point2 8 8 ]+++++++testPoly5 :: SimplePolygon () Rational+testPoly5 = toCounterClockWiseOrder . fromPoints $ map ext $ [ Point2 176 736+                                                             , Point2 240 688+                                                             , Point2 240 608+                                                             , Point2 128 576+                                                             , Point2 64 640+                                                             , Point2 80 720+                                                             , Point2 128 752+                                                             ]+++-- testPoly5 :: SimplePolygon () Rational+-- testPoly5 = toCounterClockWiseOrder . fromPoints $ map ext $ [ Point2 320 320+--                                                              , Point2 256 320+--                                                              , Point2 224 320+--                                                              , Point2 128 240+--                                                              , Point2 64 224+--                                                              , Point2 256 192+--                                                              ]

@@ -0,0 +1,138 @@+{-# LANGUAGE ScopedTypeVariables #-}+module Algorithms.Geometry.PolygonTriangulation.Types where++import           Control.Lens+import           Control.Monad (forM_)+import           Data.Ext+import qualified Data.Foldable as F+import           Data.Geometry.LineSegment+import           Data.Geometry.PlanarSubdivision.Basic+import qualified Data.List.NonEmpty as NonEmpty+import           Data.List.NonEmpty (NonEmpty(..))+import qualified Data.PlaneGraph as PG+import           Data.Semigroup+import qualified Data.Vector as V+import qualified Data.Vector.Mutable as MV+--------------------------------------------------------------------------------++data PolygonEdgeType = Original | Diagonal+                     deriving (Show,Read,Eq)++-- | Given a list of original edges and a list of diagonals, creates a+-- planar-subdivision+--+--+-- running time: \(O(n\log n)\)+constructSubdivision                  :: forall proxy r s p. (Fractional r, Ord r)+                                      => proxy s+                                      -> LineSegment 2 p r -- ^ A counter-clockwise+                                                         -- edge along the outer+                                                         -- boundary+                                      -> [LineSegment 2 p r] -- ^ remaining original edges+                                      -> [LineSegment 2 p r] -- ^ diagonals+                                      -> PlanarSubdivision s+                                            p PolygonEdgeType PolygonFaceData r+constructSubdivision px e origs diags = fromPlaneGraph $ constructGraph px e origs diags++-- constructSubdivision px e origs diags =+--     subdiv & rawVertexData.traverse.dataVal  %~ NonEmpty.head+--            & rawFaceData                     %~ V.zipWith zipF faceData'+--            & rawDartData.traverse.dataVal    %~ snd+--   where+--     subdiv :: PlanarSubdivision s (NonEmpty p) (Bool,PolygonEdgeType) () r+--     subdiv = fromConnectedSegments px $ e' : origs' <> diags'++--     diags' = (:+ (True, Diagonal)) <$> diags+--     origs' = (:+ (False,Original)) <$> origs+--     e'     = e :+ (True, Original)++--     -- the darts incident to internal faces+--     queryDarts = concatMap shouldQuery . F.toList . edges' $ subdiv+--     shouldQuery d = case subdiv^.dataOf d of+--                       (True, Original) -> [d]+--                       (True, Diagonal) -> [d, twin d]+--                       _                -> []++--     -- the interior faces+--     intFaces = flip leftFace subdiv <$> queryDarts+--     faceData' = V.create $ do+--                   v' <- MV.replicate (numFaces subdiv) Outside+--                   forM_ intFaces $ \(PG.FaceId (PG.VertexId f)) ->+--                     MV.write v' f Inside+--                   pure v'++--     -- set the inside/outside data value+--     zipF x rfd = rfd&dataVal .~ x+-- -- TODO: Idea: generalize the face data assignment into a function+-- -- that does something like: [(Dart, fLeft, fRight] -> FaceData+++-- | Given a list of original edges and a list of diagonals, creates a+-- planar-subdivision+--+--+-- running time: \(O(n\log n)\)+constructGraph                  :: forall proxy r s p. (Fractional r, Ord r)+                                      => proxy s+                                      -> LineSegment 2 p r -- ^ A counter-clockwise+                                                         -- edge along the outer+                                                         -- boundary+                                      -> [LineSegment 2 p r] -- ^ remaining original edges+                                      -> [LineSegment 2 p r] -- ^ diagonals+                                      -> PG.PlaneGraph s+                                            p PolygonEdgeType PolygonFaceData r+constructGraph px e origs diags =+    subdiv & PG.vertexData.traverse  %~ NonEmpty.head+           & PG.faceData             .~ faceData'+           & PG.rawDartData.traverse  %~ snd+  where+    subdiv :: PG.PlaneGraph s (NonEmpty p) (Bool,PolygonEdgeType) () r+    subdiv = PG.fromConnectedSegments px $ e' : origs' <> diags'++    diags' = (:+ (True, Diagonal)) <$> diags+    origs' = (:+ (False,Original)) <$> origs+    e'     = e :+ (True, Original)++    -- the darts incident to internal faces+    queryDarts = concatMap shouldQuery . F.toList . PG.edges' $ subdiv+    shouldQuery d = case subdiv^.dataOf d of+                      (True, Original) -> [d]+                      (True, Diagonal) -> [d, twin d]+                      _                -> []++    -- the interior faces+    intFaces = flip PG.leftFace subdiv <$> queryDarts+    faceData' :: V.Vector PolygonFaceData+    faceData' = V.create $ do+                  v' <- MV.replicate (PG.numFaces subdiv) Outside+                  forM_ intFaces $ \(PG.FaceId (PG.VertexId f)) ->+                    MV.write v' f Inside+                  pure v'++-- -constructSubdivision px e origs diags =+-- -    subdiv & planeGraph.PG.vertexData.traverse        %~ NonEmpty.head+-- -           & planeGraph.PG.faceData                   .~ faceData'+-- -           & planeGraph.PG.rawDartData.traverse.eData %~ snd+-- -  where+-- -    subdiv = fromConnectedSegments px $ e' : origs' <> diags'+-- -+-- -    diags' = (:+ EdgeData Visible (True, Diagonal)) <$> diags+-- -    origs' = (:+ EdgeData Visible (False,Original)) <$> origs+-- -    e'     = e :+ EdgeData Visible (True, Original)+-- -+-- -    g = subdiv^.planeGraph+-- -+-- -    -- the darts incident to internal faces+-- -    queryDarts = concatMap shouldQuery . F.toList . PG.edges' $ g+-- -    shouldQuery d = case g^.dataOf d.eData of+-- -                      (True, Original) -> [d]+-- -                      (True, Diagonal) -> [d, twin d]+-- -                      _                -> []+-- -+-- -    -- the interior faces+-- -    intFaces = flip PG.leftFace g <$> queryDarts+-- -    faceData' = V.create $ do+-- -                  v' <- MV.replicate (PG.numFaces g) (FaceData [] Outside)+-- -                  forM_ intFaces $ \(PG.FaceId (PG.VertexId f)) ->+-- -                    MV.write v' f (FaceData [] Inside)+-- -                  pure v'

@@ -0,0 +1,151 @@+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE UndecidableInstances #-}+module Algorithms.Geometry.Sweep where++import qualified Data.Map as Map+import           Data.Map (Map)+import           Data.Proxy+import           Data.Reflection+import           Unsafe.Coerce++++newtype Tagged (s :: *) a = Tagged { unTag :: a} deriving (Show,Eq,Ord)++tag   :: proxy s -> a -> Tagged s a+tag _ = Tagged++newtype Timed s t a = Timed {atTime :: (Tagged s t) -> a }+++instance (Reifies s t, Ord k) => Ord (Timed s t k) where+  compare = compare_++instance (Reifies s t, Ord k) => Eq (Timed s t k) where+  a == b = a `compare` b == EQ++compare_                       :: forall s t k. (Ord k, Reifies s t)+                               => Timed s t k -> Timed s t k+                               -> Ordering+(Timed f) `compare_` (Timed g) = let t = reflect (Proxy :: Proxy s)+                                 in f (Tagged t) `compare` g (Tagged t)+++coerceTo :: proxy s -> f (Timed s' t k) v -> f (Timed s t k) v+coerceTo _ = unsafeCoerce++unTagged :: f (Timed s t k) v -> f (Timed () t k) v+unTagged = coerceTo (Proxy :: Proxy ())+++-- | Runs a computation at a given time.+runAt       :: forall s0 t k r f v. Ord k+            => t+            -> f (Timed s0 t k) v+            -> (forall s. Reifies s t => f (Timed s t k) v -> r)+            -> r+runAt t m f = reify t $ \prx -> f (coerceTo prx m)++getTime :: Timed s Int Int+getTime = Timed unTag++constT   :: proxy s -> Int -> Timed s Int Int+constT _ i = Timed (const i)+++test1 i = reify 5 $ \prx -> getTime < constT prx i++++++test2M   :: Reifies s Int => proxy s -> Map (Timed s Int Int) String+test2M p = Map.fromList [ (constT p 10, "ten")+                        , (getTime, "timed")+                        ]+++query :: forall s v. Ord (Timed s Int Int)+      => Map (Timed s Int Int) v -> Maybe v+query = fmap snd . Map.lookupGE (constT (Proxy :: Proxy s) 4)+++test2   :: Int -> Maybe String+test2 t = runAt t m query+  where+    m :: Map (Timed () Int Int) String+    m = reify 0 $ \p -> unTagged $ test2M p++++++-- test2 = reify 0 $ \p0 ->+--                     let m = unTagged $ test2M p0+--                     in runAt 10 m Map.lookup++++-- newtype Key s a b = Key { getKey :: a -> b }++-- instance (Eq b, Reifies s a) => Eq (Key s a b) where+--   (Key f) == (Key g) = let x = reflect (Proxy :: Proxy s)+--                        in f x == g x++-- instance (Ord b, Reifies s a) => Ord (Key s a b) where+--   Key f `compare` Key g = let x = reflect (Proxy :: Proxy s)+--                           in f x `compare` g x+++-- -- | Query the sweep+-- queryAt       :: a+--               -> (forall (s :: *). Reifies s a => Map (Key s a b) v -> res)+--               -> Map (a -> b) v -> res+-- queryAt x f m = reify x (\p -> f . coerceKeys p $ m)++-- updateAt      :: a+--               -> (forall (s :: *). Reifies s a =>+--                    Map (Key s a b) v -> Map (Key s a b) v')+--               -> Map (a -> b) v+--               -> Map (a -> b) v'+-- updateAt x f m = reify x (\p -> uncoerceKeys . f . coerceKeys p $ m)+++-- combineAt            :: a+--                      -> (forall (s :: *). Reifies s a =>+--                            Map (Key s a b) v -> Map (Key s a b) v+--                            -> Map (Key s a b) v)+--                      -> Map (a -> b) v+--                      -> Map (a -> b) v+--                      -> Map (a -> b) v+-- combineAt x uF m1 m2 = reify x (\p -> uncoerceKeys $+--                                         coerceKeys p m1 `uF` coerceKeys p m2)+++-- splitLookupAt       :: Ord b+--                     => a+--                     -> (a -> b)+--                     -> Map (a -> b) v+--                     -> (Map (a -> b) v, Maybe v, Map (a -> b) v)+-- splitLookupAt x k m = reify x (\p -> let (l,mv,r) = Map.splitLookup (Key k)+--                                                   $ coerceKeys p m+--                                    in (uncoerceKeys l, mv, uncoerceKeys r))+++-- --------------------------------------------------------------------------------++-- coerceKeys   :: proxy s -> Map (a -> b) v -> Map (Key s a b) v+-- coerceKeys _ = unsafeCoerce++-- uncoerceKeys :: Map (Key s a b) v -> Map (a -> b) v+-- uncoerceKeys = unsafeCoerce+++-- --------------------------------------------------------------------------------+++-- data Node a = Node2 a a+--             | Node3 a a a++-- data FT a = Single a+--           | Deep (FT (Node a)) a (FT (Node a))

@@ -1,4 +1,5 @@ {-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE UndecidableInstances  #-} module Algorithms.Geometry.WellSeparatedPairDecomposition.Types where  import           Control.Lens hiding (Level)

@@ -1,4 +1,3 @@-{-# LANGUAGE TemplateHaskell  #-} {-# LANGUAGE LambdaCase  #-} module Algorithms.Geometry.WellSeparatedPairDecomposition.WSPD where @@ -36,10 +35,8 @@ -- | Construct a split tree -- -- running time: \(O(n \log n)\)-fairSplitTree     :: (Fractional r, Ord r, Arity d, Index' 0 d,-                      KnownNat d-                        , Show r, Show p-+fairSplitTree     :: (Fractional r, Ord r, Arity d, 1 <= d+                     , Show r, Show p                      )                   => NonEmpty.NonEmpty (Point d r :+ p) -> SplitTree d p r () fairSplitTree pts = foldUp node' Leaf $ fairSplitTree' n pts'@@ -62,7 +59,7 @@ -- | Given a split tree, generate the Well separated pairs -- -- running time: \(O(s^d n)\)-wellSeparatedPairs   :: (Floating r, Ord r, AlwaysTrueWSPD d)+wellSeparatedPairs   :: (Floating r, Ord r, Arity d, Arity (d + 1))                      => r -> SplitTree d p r a -> [WSP d p r a] wellSeparatedPairs s = f   where@@ -107,7 +104,7 @@ -- (PointSeq))'s. Since we have the level assignment, we can compute these -- lists by traversing each original input list (i.e. one for every dimension) -- once, and partition the points based on their level assignment.-fairSplitTree'       :: (Fractional r, Ord r, Arity d, Index' 0 d, KnownNat d+fairSplitTree'       :: (Fractional r, Ord r, Arity d, 1 <= d                         , Show r, Show p                         )                      => Int -> GV.Vector d (PointSeq d (Idx :+ p) r)@@ -187,7 +184,7 @@ -- time. -- -- so, basically, run reIndex points in ST as well.-reIndexPoints      :: (Arity d, Index' 0 d)+reIndexPoints      :: (Arity d, 1 <= d)                    => GV.Vector d (PointSeq d (Idx :+ p) r)                    -> GV.Vector d (PointSeq d (Idx :+ p) r) reIndexPoints ptsV = fmap reIndex ptsV@@ -357,10 +354,10 @@ -------------------------------------------------------------------------------- -- * Finding Well Separated Pairs -type AlwaysTrueWSPD d = ( Arity d, KnownNat d-                        , AlwaysTruePFT d, AlwaysTrueTransformation d)+-- type AlwaysTrueWSPD d = ( Arity d, KnownNat d+--                         , AlwaysTruePFT d, AlwaysTrueTransformation d) -findPairs                     :: (Floating r, Ord r, AlwaysTrueWSPD d)+findPairs                     :: (Floating r, Ord r, Arity d, Arity (d + 1))                               => r -> SplitTree d p r a -> SplitTree d p r a                               -> [WSP d p r a] findPairs s l r@@ -370,7 +367,7 @@   -- | Test if the two sets are well separated with param s-areWellSeparated                     :: ( AlwaysTrueWSPD d, Fractional r, Ord r)+areWellSeparated                     :: (Arity d, Arity (d + 1), Fractional r, Ord r)                                      => r -- ^ separation factor                                      -> SplitTree d p r a                                      -> SplitTree d p r a -> Bool@@ -392,7 +389,7 @@ --     b' = translateBy v . scaleUniformlyBy s . translateBy ((-1) *^ v) $ b  -- | Test if the two boxes are sufficiently far appart-boxBox         :: (Fractional r, Ord r, AlwaysTruePFT d, AlwaysTrueTransformation d)+boxBox         :: (Fractional r, Ord r, Arity d, Arity (d + 1))                => r -> Box d p r -> Box d p r -> Bool boxBox s lb rb = boxBox' lb rb && boxBox' rb lb   where

@@ -13,7 +13,7 @@  -- | A State monad that can store earlier versions of the state. newtype PersistentStateT s m a =-  PersistentStateT { runPersistentStateT' :: StateT (NonEmpty s) m a }+  PersistentStateT (StateT (NonEmpty s) m a)   deriving (Functor,Applicative,Monad)            -- We store all the versions in reverse order 

@@ -1,6 +1,7 @@ {-# LANGUAGE RecordWildCards #-} module Data.BalBST where +import           Control.Applicative((<|>)) import           Data.Bifunctor import           Data.Function (on) import           Data.Functor.Contravariant@@ -93,8 +94,16 @@ member   :: Eq a => a -> BalBST k a -> Bool member x = isJust . lookup x --+-- | Search for the Predecessor+-- \(O(\log n)\)+lookupLE :: Ord k => k -> BalBST k a -> Maybe a+lookupLE kx (BalBST n@Nav{..} t) = lookup' t+  where+    lookup' Empty            = Nothing+    lookup' (Leaf y)         = if goLeft y kx then Just y else Nothing+    lookup' (Node _ _ l k r)+      | kx <= k              = lookup' l+      | otherwise            = lookup' r <|> lookupMax (BalBST n l)   -- | Insert an element in the BST.@@ -115,8 +124,15 @@  -- delete = undefined --- delete                        :: Eq a => a -> BalBST k a -> BalBST k a--- delete x (BalBST n@Nav{..} t) = delete' t+-- | Delete (one occurance of) an element.+-- \(O(\log n)\)+delete                        :: Eq a => a -> BalBST k a -> BalBST k a+delete x t = let Split l _ r = split x t+                 n           = nav t+             in BalBST n $ joinWith n l r+++-- (BalBST n@Nav{..} t) = delete' t --   where --     delete' Empty      = Empty --     delete' l@(Leaf y) = if x == y then Empty else l@@ -136,6 +152,9 @@     minView' (Leaf x)         = Just (x,Empty)     minView' (Node _ _ l _ r) = fmap (flip (joinWith n) r) <$> minView' l +lookupMin :: BalBST k b -> Maybe b+lookupMin = fmap fst . maxView+ -- | Extract the maximum from the tree -- \(O(\log n)\) maxView              :: BalBST k a -> Maybe (a, Tree k a)@@ -144,6 +163,10 @@     maxView' Empty            = Nothing     maxView' (Leaf x)         = Just (x,Empty)     maxView' (Node _ _ l _ r) = fmap (joinWith n l) <$> maxView' r++lookupMax :: BalBST k b -> Maybe b+lookupMax = fmap fst . maxView+  -- | Joins two BSTs. Assumes that the ranges are disjoint. It takes the left Tree nav --

@@ -32,17 +32,19 @@                        , isShiftOf                        ) where +import           Control.DeepSeq import           Control.Lens (lens, Lens', bimap) import qualified Data.Foldable as F import qualified Data.List as L import qualified Data.List.NonEmpty as NonEmpty import           Data.Maybe (listToMaybe) import           Data.Semigroup-import           Data.Semigroup.Foldable+import           Data.Semigroup.Foldable hiding (toNonEmpty) import           Data.Sequence ((|>),(<|),ViewL(..),ViewR(..),Seq) import qualified Data.Sequence as S import qualified Data.Traversable as T import           Data.Tuple (swap)+import           GHC.Generics (Generic)  -------------------------------------------------------------------------------- @@ -52,7 +54,10 @@  -- | Nonempty circular sequence data CSeq a = CSeq !(Seq a) !a !(Seq a)+  deriving (Generic)                      -- we keep the seq balanced, i.e. size left >= size right++instance NFData a => NFData (CSeq a)  instance Eq a => Eq (CSeq a) where   a == b = asSeq a == asSeq b

@@ -1,7 +1,14 @@ {-# LANGUAGE DeriveAnyClass  #-}+{-# LANGUAGE OverloadedStrings  #-}+{-|+Module    : Data.Ext+Description: A pair-like data type to represent a 'core' type that has extra information as well.+Copyright : (c) Frank Staals+License : See LICENCE file+-} module Data.Ext where -import Control.Lens+import Control.Lens hiding ((.=)) import Data.Biapplicative import Data.Bifoldable import Data.Bifunctor.Apply@@ -12,9 +19,12 @@ import Data.Semigroup.Bitraversable import GHC.Generics (Generic) import Control.DeepSeq-+import Data.Aeson+import Data.Aeson.Types(typeMismatch) -------------------------------------------------------------------------------- +-- | Our Ext type that represents the core datatype core extended with extra+-- information of type 'extra'. data core :+ extra = core :+ extra deriving (Show,Read,Eq,Ord,Bounded,Generic,NFData) infixr 1 :+ @@ -43,6 +53,17 @@ instance (Semigroup core, Semigroup extra) => Semigroup (core :+ extra) where   (c :+ e) <> (c' :+ e') = c <> c' :+ e <> e' ++instance (ToJSON core, ToJSON extra) => ToJSON (core :+ extra) where+  -- toJSON     (c :+ e) = toJSON     (c,e)+  -- toEncoding (c :+ e) = toEncoding (c,e)+  toJSON     (c :+ e) = object ["core" .= c, "extra" .= e]+  toEncoding (c :+ e) = pairs  ("core" .= c <> "extra" .= e)++instance (FromJSON core, FromJSON extra) => FromJSON (core :+ extra) where+  -- parseJSON = fmap (\(c,e) -> c :+ e) . parseJSON+  parseJSON (Object v) = (:+) <$> v .: "core" <*> v .: "extra"+  parseJSON invalid    = typeMismatch "Ext (:+)" invalid  _core :: (core :+ extra) -> core _core (c :+ _) = c

@@ -7,16 +7,16 @@ module Data.Geometry( module Data.Geometry.Properties                     , module Data.Geometry.Transformation                     , module Data.Geometry.Point-                    , module Data.Geometry.Vector+                    , module V                     , module Data.Geometry.Line                     , module Data.Geometry.LineSegment                     , module Data.Geometry.PolyLine                     , module Data.Geometry.Polygon-                    , module Linear.Affine-                    , module Linear.Vector+                    -- , module Linear.Affine+                    -- , module Linear.Vector                     ) where -+import Data.Geometry.Vector as V hiding (last) import Data.Geometry.Line import Data.Geometry.LineSegment import Data.Geometry.Point@@ -24,6 +24,5 @@ import Data.Geometry.Polygon hiding (fromPoints) import Data.Geometry.Properties import Data.Geometry.Transformation-import Data.Geometry.Vector-import Linear.Affine hiding (Point, Vector, origin)-import Linear.Vector+-- import Linear.Affine hiding (Point, Vector, origin)+-- import Linear.Vector

@@ -1,5 +1,11 @@ {-# LANGUAGE TemplateHaskell  #-} {-# LANGUAGE UndecidableInstances #-}+{-|+Module    : Data.Geometry.Ball+Description: \(d\)-dimensional Balls and Spheres+Copyright : (c) Frank Staals+License : See LICENCE file+-} module Data.Geometry.Ball where  import           Control.DeepSeq@@ -16,7 +22,7 @@ import qualified Data.List as L import qualified Data.Traversable as T import           Data.Vinyl-import           Frames.CoRec+import           Data.Vinyl.CoRec import           GHC.Generics (Generic) import           Linear.Matrix import           Linear.V3 (V3(..))@@ -36,7 +42,7 @@   deriving instance (Show r, Show p, Arity d)     => Show (Ball d p r)--- deriving instance (NFData p, NFData r, Arity d) => NFData (Ball d p r)+instance (NFData p, NFData r, Arity d) => NFData (Ball d p r) deriving instance (Eq r, Eq p, Arity d)         => Eq (Ball d p r)  type instance NumType   (Ball d p r) = r@@ -111,6 +117,7 @@  pattern Sphere     :: Point d r :+ p -> r -> Sphere d p r pattern Sphere c r = Boundary (Ball c r)+{-# COMPLETE Sphere #-}   @@ -122,12 +129,14 @@  pattern Disk     :: Point 2 r :+ p -> r -> Disk p r pattern Disk c r = Ball c r+{-# COMPLETE Disk #-}   type Circle p r = Sphere 2 p r  pattern Circle     :: Point 2 r :+ p ->  r -> Circle p r pattern Circle c r = Sphere c r+{-# COMPLETE Circle #-}  -- | Given three points, get the disk through the three points. If the three -- input points are colinear we return Nothing

@@ -7,8 +7,8 @@  -- | The boundary of a geometric object. newtype Boundary g = Boundary g-                   deriving (Show,Eq,Ord,Read,IsTransformable)-+                   deriving (Show,Eq,Ord,Read,IsTransformable+                            ,Functor,Foldable,Traversable)  type instance NumType (Boundary g)   = NumType g type instance Dimension (Boundary g) = Dimension g

@@ -1,7 +1,14 @@ {-# LANGUAGE TemplateHaskell  #-} {-# LANGUAGE ScopedTypeVariables  #-}+{-# LANGUAGE UndecidableInstances  #-} {-# LANGUAGE DeriveAnyClass  #-} {-# OPTIONS_GHC -fno-warn-orphans #-}+{-|+Module    : Data.Geometry.Box+Description: Orthogonal \(d\)-dimensiontal boxes (e.g. rectangles)+Copyright : (c) Frank Staals+License : See LICENCE file+-} module Data.Geometry.Box( module Data.Geometry.Box.Internal                         , topSide, leftSide, bottomSide, rightSide                         , sides, sides'

@@ -1,16 +1,17 @@ {-# LANGUAGE TemplateHaskell  #-} {-# LANGUAGE ScopedTypeVariables  #-}+{-# LANGUAGE UndecidableInstances  #-}+{-# LANGUAGE InstanceSigs  #-} module Data.Geometry.Box.Internal where  import           Control.DeepSeq import           Control.Lens import           Data.Bifunctor import           Data.Ext-import           Frames.CoRec (asA) import           Data.Geometry.Point import           Data.Geometry.Properties import           Data.Geometry.Transformation-import           Data.Geometry.Vector (Vector, Arity, Index',C(..))+import           Data.Geometry.Vector (Vector, Arity, C(..)) import qualified Data.Geometry.Vector as V import qualified Data.List.NonEmpty as NE import           Data.Proxy@@ -18,8 +19,9 @@ import           Data.Semigroup import qualified Data.Semigroup.Foldable as F import qualified Data.Vector.Fixed as FV-import           GHC.TypeLits+import           Data.Vinyl.CoRec (asA) import           GHC.Generics (Generic)+import           GHC.TypeLits   --------------------------------------------------------------------------------@@ -61,6 +63,12 @@                     (CWMax (Point $ fmap (^.R.upper.R.unEndPoint) rs) :+ mempty)  +-- | Given a center point and a vector specifying the box width's, construct a box.+fromCenter      :: (Arity d, Fractional r) => Point d r -> Vector d r -> Box d () r+fromCenter c ws = let f x r = R.ClosedRange (x-r) (x+r)+                  in fromExtent $ FV.zipWith f (toVec c) ((/2) <$> ws)++ -- | Center of the box centerPoint   :: (Arity d, Fractional r) => Box d p r -> Point d r centerPoint b = Point $ w V.^/ 2@@ -85,7 +93,12 @@       f = maybe (coRec NoIntersection) (coRec . fromExtent)       r `intersect'` s = asA (Proxy :: Proxy (R.Range r)) $ r `intersect` s -+instance Arity d => Bifunctor (Box d) where+  bimap :: forall p q r s. (p -> q) -> (r -> s) -> Box d p r -> Box d q s+  bimap f g (Box mi ma) = Box (bimap g' f mi) (bimap g' f ma)+    where+      g' :: Functor g => g (Point d r) -> g (Point d s)+      g' = fmap (fmap g)  -- -- In principle this should also just work for Boxes in higher dimensions. It is just -- -- that we need a better way to compute their corners@@ -119,7 +132,8 @@   pmap f (Box mi ma) = Box (first (fmap f) mi) (first (fmap f) ma)  -instance (Num r, AlwaysTruePFT d) => IsTransformable (Box d p r) where+instance (Fractional r, Arity d, Arity (d + 1))+         => IsTransformable (Box d p r) where   -- Note that this does not guarantee the box is still a proper box Only use   -- this to do translations and scalings. Other transformations may produce   -- unexpected results.@@ -171,7 +185,8 @@ -- 1 -- >>> widthIn (C :: C 3) (boundingBoxList' [origin, point3 1 2 3] :: Box 3 () Int) -- 3-widthIn   :: forall proxy p i d r. (Arity d, Num r, Index' (i-1) d) => proxy i -> Box d p r -> r+widthIn   :: forall proxy p i d r. (Arity d, Arity (i - 1), Num r, ((i-1)+1) <= d)+          => proxy i -> Box d p r -> r widthIn _ = view (V.element (C :: C (i - 1))) . size  

@@ -4,6 +4,7 @@ module Data.Geometry.HalfLine where  +import           Control.DeepSeq import           Control.Lens import           Data.Ext import qualified Data.Foldable as F@@ -18,7 +19,7 @@ import qualified Data.Traversable as T import           Data.UnBounded import           GHC.Generics (Generic)-import           Control.DeepSeq+import           GHC.TypeLits  -------------------------------------------------------------------------------- -- * d-dimensional Half-Lines@@ -50,7 +51,7 @@   supportingLine ~(HalfLine p v) = Line p v  -- Half-Lines are transformable-instance (Num r, AlwaysTruePFT d) => IsTransformable (HalfLine d r) where+instance (Fractional r, Arity d, Arity (d + 1)) => IsTransformable (HalfLine d r) where   transformBy t = toHalfLine . transformPointFunctor t . toLineSegment'     where       toLineSegment' :: (Num r, Arity d) => HalfLine d r -> LineSegment d () r

@@ -27,7 +27,7 @@ import           Data.Semigroup import qualified Data.Traversable as T import           Data.Vinyl-import           Frames.CoRec+import           Data.Vinyl.CoRec import           GHC.Generics (Generic)  --------------------------------------------------------------------------------@@ -74,7 +74,7 @@  pattern Interval     :: EndPoint (r :+ a) -> EndPoint (r :+ a) -> Interval a r pattern Interval l u = GInterval (Range l u)-+{-# COMPLETE Interval #-}  -------------------------------------------------------------------------------- 

@@ -1,3 +1,9 @@+{-|+Module    : Data.Geometry.Ipe+Description: Reexports the functionality for reading and writing Ipe files.+Copyright : (c) Frank Staals+License : See LICENCE file+-} module Data.Geometry.Ipe( module Data.Geometry.Ipe.Types                         , module Data.Geometry.Ipe.Writer                         , module Data.Geometry.Ipe.Reader

@@ -5,15 +5,16 @@ {-# LANGUAGE UndecidableInstances #-} module Data.Geometry.Ipe.Attributes where -import           Control.Lens hiding (rmap, Const)-import           Data.Semigroup-import           Data.Singletons-import           Data.Singletons.TH-import           Data.Text(Text)-import           Data.Vinyl-import           Data.Vinyl.Functor-import           Data.Vinyl.TypeLevel-import           GHC.Exts+import Control.Lens hiding (rmap, Const)+import Data.Colour.SRGB+import Data.Semigroup+import Data.Singletons+import Data.Singletons.TH+import Data.Text (Text)+import Data.Vinyl+import Data.Vinyl.Functor+import Data.Vinyl.TypeLevel+import GHC.Exts  -------------------------------------------------------------------------------- @@ -68,15 +69,20 @@  makeLenses ''Attr +pattern Attr   :: Apply f label -> Attr f label pattern Attr x = GAttr (Just x)++pattern NoAttr :: Attr f label pattern NoAttr = GAttr Nothing  -- | Give pref. to the *RIGHT*-instance Monoid (Attr f l) where-  mempty                 = NoAttr-  _ `mappend` b@(Attr _) = b-  a `mappend` _          = a+instance Semigroup (Attr f l) where+  _ <> b@(Attr _) = b+  a <> _          = a +instance Monoid (Attr f l) where+  mempty  = NoAttr+  mappend = (<>)  newtype Attributes (f :: TyFun u * -> *) (ats :: [u]) =   Attrs { _unAttrs :: Rec (Attr f) ats }@@ -187,16 +193,19 @@   -- | Many types either consist of a symbolc value, or a value of type v-data IpeValue v = Named Text | Valued v deriving (Show,Eq,Ord)+data IpeValue v = Named Text | Valued v deriving (Show,Eq,Ord,Functor,Foldable,Traversable)  instance IsString (IpeValue v) where   fromString = Named . fromString -type Colour = Text -- TODO: Make this a Colour.Colour+newtype IpeSize  r = IpeSize  (IpeValue r)          deriving (Show,Eq,Ord)+newtype IpePen   r = IpePen   (IpeValue r)          deriving (Show,Eq,Ord)+newtype IpeColor r = IpeColor (IpeValue (RGB r))    deriving (Show,Eq) -newtype IpeSize r = IpeSize  (IpeValue r)      deriving (Show,Eq,Ord)-newtype IpePen  r = IpePen   (IpeValue r)      deriving (Show,Eq,Ord)-newtype IpeColor  = IpeColor (IpeValue Colour) deriving (Show,Eq,Ord)+instance Ord r => Ord (IpeColor r) where+  (IpeColor c) `compare` (IpeColor c') = fmap f c `compare` fmap f c'+    where+      f (RGB r g b) = (r,g,b)   -- -- | And the corresponding types@@ -317,26 +326,6 @@ -- GroupAttributeUniverse instance IpeAttrName Clip     where attrName _ = "clip" ---- | Wrap up a value with a capability given by its type-data GDict (c :: k -> Constraint) (a :: k) where-  GDict :: c a => Proxy a -> GDict c a---- -- | Sometimes we may know something for /all/ fields of a record, but when--- -- you expect to be able to /each/ of the fields, you are then out of luck.--- -- Surely given @∀x:u.φ(x)@ we should be able to recover @x:u ⊢ φ(x)@! Sadly,--- -- the constraint solver is not quite smart enough to realize this and we must--- -- make it patently obvious by reifying the constraint pointwise with proof.--- gReifyConstraint---   :: RecAll f rs c---   => proxy c---   -> Rec f rs---   -> Rec (GDict c :. f) rs--- gReifyConstraint prx RNil      = RNil--- gReifyConstraint prx (x :& xs) = Compose (mkDict x) :& reifyConstraint prx xs---   where---     mkDict   :: f l -> (GDict c l)---     mkDict _ = GDict (Proxy :: Proxy l)  -- | Function that states that all elements in xs satisfy a given constraint c type family AllSatisfy (c :: k -> Constraint) (xs :: [k]) :: Constraint where

@@ -1,49 +1,146 @@+{-# LANGUAGE OverloadedStrings #-} module Data.Geometry.Ipe.FromIpe where -import           Control.Lens+import           Control.Lens hiding (Simple) import           Data.Ext+import           Data.Geometry.Ipe.Reader import           Data.Geometry.Ipe.Types-import           Data.Geometry.Line import           Data.Geometry.LineSegment import qualified Data.Geometry.PolyLine as PolyLine import           Data.Geometry.Polygon+import           Data.Geometry.Properties import qualified Data.Seq2 as S2-import qualified Data.Traversable as Tr+import qualified Data.List.NonEmpty as NonEmpty +--------------------------------------------------------------------------------+-- $setup+-- >>> :{+-- import           Data.Geometry.Ipe.Attributes+--+-- let testPath :: Path Int+--     testPath = Path . S2.l1Singleton  . PolyLineSegment . PolyLine.fromPoints . map ext+--              $ [ origin, point2 10 10, point2 200 100 ]+--+--     testPathAttrs :: IpeAttributes Path Int+--     testPathAttrs = attr SStroke (IpeColor (Named "red")) +--     testObject :: IpeObject Int+--     testObject = IpePath (testPath :+ testPathAttrs)+-- :}++ -- | Try to convert a path into a line segment, fails if the path is not a line -- segment or a polyline with more than two points.+--+-- _asLineSegment :: Prism' (Path r) (LineSegment 2 () r) _asLineSegment = prism' seg2path path2seg   where     seg2path   = review _asPolyLine . PolyLine.fromLineSegment     path2seg p = PolyLine.asLineSegment' =<< preview _asPolyLine p - -- | Convert to a polyline. Ignores all non-polyline parts+--+-- >>> testPath ^? _asPolyLine+-- Just (PolyLine {_points = Seq2 (Point2 [0,0] :+ ()) (fromList [Point2 [10,10] :+ ()]) (Point2 [200,100] :+ ())}) _asPolyLine :: Prism' (Path r) (PolyLine.PolyLine 2 () r) _asPolyLine = prism' poly2path path2poly   where     poly2path = Path . S2.l1Singleton  . PolyLineSegment-    path2poly = preview (pathSegments.Tr.traverse._PolyLineSegment)+    path2poly = preview (pathSegments.traverse._PolyLineSegment)     -- TODO: Check that the path actually is a polyline, rather     -- than ignoring everything that does not fit --- | Convert to a simple polygon: simply takes the first closed path-_asSimplePolygon :: Prism' (Path r) (SimplePolygon () r)-_asSimplePolygon = prism' poly2path path2poly+-- | Convert to a simple polygon+_asSimplePolygon :: Prism' (Path r) (Polygon Simple () r)+_asSimplePolygon = prism' polygonToPath path2poly   where-    poly2path = Path . S2.l1Singleton . PolygonPath-    path2poly = preview (pathSegments.Tr.traverse._PolygonPath)-    -- TODO: Check that the path actually is a simple polygon, rather-    -- than ignoring everything that does not fit+    path2poly p = pathToPolygon p >>= either pure (const Nothing) +-- | Convert to a multipolygon+_asMultiPolygon :: Prism' (Path r) (MultiPolygon () r)+_asMultiPolygon = prism' polygonToPath path2poly+  where+    path2poly p = pathToPolygon p >>= either (const Nothing) pure++polygonToPath                      :: Polygon t () r -> Path r+polygonToPath pg@(SimplePolygon _) = Path . S2.l1Singleton . PolygonPath $ pg+polygonToPath (MultiPolygon vs hs) = Path . S2.viewL1FromNonEmpty . fmap PolygonPath+                                   $ SimplePolygon vs NonEmpty.:| hs+++pathToPolygon   :: Path r -> Maybe (Either (SimplePolygon () r) (MultiPolygon () r))+pathToPolygon p = case p^..pathSegments.traverse._PolygonPath of+                    []                   -> Nothing+                    [pg]                 -> Just . Left  $ pg+                    SimplePolygon vs: hs -> Just . Right $ MultiPolygon vs hs+++ -- | use the first prism to select the ipe object to depicle with, and the second -- how to select the geometry object from there on. Then we can select the geometry -- object, directly with its attributes here.+--+-- >>> testObject ^? _withAttrs _IpePath _asPolyLine+-- Just (PolyLine {_points = Seq2 (Point2 [0,0] :+ ()) (fromList [Point2 [10,10] :+ ()]) (Point2 [200,100] :+ ())} :+ Attrs {_unAttrs = {GAttr {_getAttr = Nothing}, GAttr {_getAttr = Nothing}, GAttr {_getAttr = Nothing}, GAttr {_getAttr = Nothing}, GAttr {_getAttr = Just (IpeColor (Named "red"))}, GAttr {_getAttr = Nothing}, GAttr {_getAttr = Nothing}, GAttr {_getAttr = Nothing}, GAttr {_getAttr = Nothing}, GAttr {_getAttr = Nothing}, GAttr {_getAttr = Nothing}, GAttr {_getAttr = Nothing}, GAttr {_getAttr = Nothing}, GAttr {_getAttr = Nothing}, GAttr {_getAttr = Nothing}, GAttr {_getAttr = Nothing}}}) _withAttrs       :: Prism' (IpeObject r) (i r :+ IpeAttributes i r) -> Prism' (i r) g                  -> Prism' (IpeObject r) (g :+ IpeAttributes i r) _withAttrs po pg = prism' g2o o2g   where     g2o    = review po . over core (review pg)     o2g o  = preview po o >>= \(i :+ ats) -> (:+ ats) <$> preview pg i++++++-- instance HasDefaultIpeObject Path where+--   defaultIpeObject' = _IpePath+++-- class HasDefaultFromIpe g where+--   type DefaultFromIpe g :: * -> *+--   defaultIpeObject :: proxy g -> Prism' (IpeObject r) (DefaultFromIpe g r :+ IpeAttributes (DefaultFromIpe g) r)+--   defaultFromIpe   :: proxy g -> Prism' (DefaultFromIpe g (NumType g)) g+++class HasDefaultFromIpe g where+  type DefaultFromIpe g :: * -> *+  defaultFromIpe :: (r ~ NumType g)+                 => Prism' (IpeObject r) (g :+ IpeAttributes (DefaultFromIpe g) r)++-- instance HasDefaultFromIpe (Point 2 r) where+--   type DefaultFromIpe (Point 2 r) = IpeSymbol+--   defaultFromIpe = _withAttrs _IpeUse symbolPoint+++instance HasDefaultFromIpe (LineSegment 2 () r) where+  type DefaultFromIpe (LineSegment 2 () r) = Path+  defaultFromIpe = _withAttrs _IpePath _asLineSegment++instance HasDefaultFromIpe (PolyLine.PolyLine 2 () r) where+  type DefaultFromIpe (PolyLine.PolyLine 2 () r) = Path+  defaultFromIpe = _withAttrs _IpePath _asPolyLine+++instance HasDefaultFromIpe (SimplePolygon () r) where+  type DefaultFromIpe (SimplePolygon () r) = Path+  defaultFromIpe = _withAttrs _IpePath _asSimplePolygon++instance HasDefaultFromIpe (MultiPolygon () r) where+  type DefaultFromIpe (MultiPolygon () r) = Path+  defaultFromIpe = _withAttrs _IpePath _asMultiPolygon+++-- | Read all g's from some ipe page(s).+readAll :: (HasDefaultFromIpe g, r ~ NumType g, Foldable f)+        => f (IpePage r) -> [g :+ IpeAttributes (DefaultFromIpe g) r]+readAll = foldMap (^..content.traverse.defaultFromIpe)+++-- | Convenience function from reading all g's from an ipe file. If there+-- is an error reading or parsing the file the error is "thrown away".+readAllFrom    :: (HasDefaultFromIpe g, r ~ NumType g, Coordinate r, Eq r)+               => FilePath -> IO [g :+ IpeAttributes (DefaultFromIpe g) r]+readAllFrom fp = readAll <$> readSinglePageFile fp+

@@ -1,14 +1,17 @@ {-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE ScopedTypeVariables #-} module Data.Geometry.Ipe.IpeOut where -import           Control.Lens+import           Control.Lens hiding (Simple) import           Data.Bifunctor import           Data.Ext import           Data.Geometry.Ball import           Data.Geometry.Boundary import           Data.Geometry.Box import           Data.Geometry.Ipe.Attributes+import           Data.Geometry.Ipe.FromIpe import           Data.Geometry.Ipe.Types+import           Data.Geometry.Line import           Data.Geometry.LineSegment import           Data.Geometry.Point import           Data.Geometry.PolyLine@@ -16,9 +19,12 @@ import           Data.Geometry.Polygon.Convex import           Data.Geometry.Properties import           Data.Geometry.Transformation+import           Data.Maybe (fromMaybe)+import           Data.Proxy import           Data.Semigroup import qualified Data.Seq2 as S2 import           Data.Text (Text)+import           Data.Vinyl.CoRec  -------------------------------------------------------------------------------- @@ -44,7 +50,8 @@ -- ipe geometry object, the geometry object, and a record with its attributes, -- construct an ipe Object representing it. asIpeObjectWith          :: (ToObject i, NumType g ~ r)-                      => IpeOut g (IpeObject' i r) -> g -> IpeAttributes i r -> IpeObject r+                         => IpeOut g (IpeObject' i r) -> g -> IpeAttributes i r+                         -> IpeObject r asIpeObjectWith io g ats = asIpe (ipeObject io ats) g  @@ -83,6 +90,10 @@   --                  => IpeOut g (IpeObject (NumType g))   -- defaultIpeObject = IpeOut $ flip asIpeObject mempty +-- instance HasDefaultIpeOut g => HasDefaultIpeOut [g] where+--   type DefaultIpeOut [g] = Group+--   defaultIpeOut = IpeOut $ asIpeGroup . map (asIpeObject' mempty)+ instance HasDefaultIpeOut (Point 2 r) where   type DefaultIpeOut (Point 2 r) = IpeSymbol   defaultIpeOut = ipeDiskMark@@ -99,11 +110,15 @@   type DefaultIpeOut (PolyLine 2 p r) = Path   defaultIpeOut = noAttrs ipePolyLine -instance HasDefaultIpeOut (SimplePolygon p r) where-  type DefaultIpeOut (SimplePolygon p r) = Path-  defaultIpeOut = flip addAttributes ipeSimplePolygon $-                    mempty <> attr SFill (IpeColor "red")+instance HasDefaultIpeOut (Polygon t p r) where+  type DefaultIpeOut (Polygon t p r) = Path+  defaultIpeOut = flip addAttributes ipePolygon $+                    mempty <> attr SFill (IpeColor "0.722 0.145 0.137") +instance HasDefaultIpeOut (SomePolygon p r) where+  type DefaultIpeOut (SomePolygon p r) = Path+  defaultIpeOut = IpeOut $ either (asIpe defaultIpeOut) (asIpe defaultIpeOut)+ instance HasDefaultIpeOut (ConvexPolygon p r) where   type DefaultIpeOut (ConvexPolygon p r) = Path   defaultIpeOut = IpeOut $ asIpe defaultIpeOut . view simplePolygon@@ -132,15 +147,22 @@ defaultClipRectangle = boundingBox (point2 (-200) (-200)) <>                        boundingBox (point2 1000 1000) --- -- | An ipe out to draw a line, by clipping it to stay within a rectangle of--- -- default size.--- line :: IpeOut (Line 2 r) (IpeObject' Path r)--- line = line' defaultClipRectangle---- -- | An ipe out to draw a line, by clipping it to stay within the rectangle--- line'   :: Rectangle p r -> IpeOut (Line 2 r) (IpeObject' Path r)--- line' r = IpeOut $ \l -> error "not implemented yet"+-- | An ipe out to draw a line, by clipping it to stay within a rectangle of+-- default size.+line :: (Fractional r, Ord r) => IpeOut (Line 2 r) (IpeObject' Path r)+line = lineWith defaultClipRectangle +-- | An ipe out to draw a line, by clipping it to stay within the rectangle.+--+-- pre: intersection of the line and the rectangle is a line segment+-- (otherwise it arbitrarily inserts the bottom of the rectangle as the path)+lineWith   :: forall p r. (Ord r, Fractional r)+              => Rectangle p r -> IpeOut (Line 2 r) (IpeObject' Path r)+lineWith r = IpeOut (asIpe defaultIpeOut . clip)+  where+    def    = bimap (const ()) id $ bottomSide r+    clip l = fromMaybe def . asA (Proxy :: Proxy (LineSegment 2 () r))+           $ l `intersect` r  ipeLineSegment :: IpeOut (LineSegment 2 p r) (IpeObject' Path r) ipeLineSegment = noAttrs $ fromPathSegment ipeLineSegment'@@ -175,11 +197,17 @@ fromPathSegment    :: IpeOut g (PathSegment r) -> IpeOut g (Path r) fromPathSegment io = IpeOut $ Path . S2.l1Singleton . asIpe io -ipeSimplePolygon :: IpeOut (SimplePolygon p r) (Path r)-ipeSimplePolygon = fromPathSegment . IpeOut $ PolygonPath . dropExt++ipePolygon :: IpeOut (Polygon t p r) (Path r)+ipePolygon = IpeOut $ io . first (const ())   where-    dropExt                    :: SimplePolygon p r -> SimplePolygon () r-    dropExt (SimplePolygon vs) = SimplePolygon $ fmap (&extra .~ ()) vs+    io                       :: forall t r. Polygon t () r -> Path r+    io pg@(SimplePolygon _)  = pg^.re _asSimplePolygon+    io pg@(MultiPolygon _ _) = pg^.re _asMultiPolygon++++   

@@ -2,7 +2,8 @@ {-# Language OverloadedStrings  #-} module Data.Geometry.Ipe.ParserPrimitives( runP, runP'                                          , pMany, pMany1, pChoice-                                         , pChar, pSpace, pWhiteSpace, pInteger, pNatural+                                         , pChar, pSpace, pWhiteSpace, pInteger+                                         , pNatural, pPaddedNatural                                          , (<*><>) , (<*><)                                          , (<***>) , (<***) , (***>)                                          , pMaybe , pCount , pSepBy@@ -43,6 +44,11 @@  pNatural :: Parser Integer pNatural = read <$> pMany1 digit++-- | parses an integer with a prefix of zeros. Returns the total length of the+-- string parced (i.e. number of digits) and the resulting antural number.+pPaddedNatural :: Parser (Int, Integer)+pPaddedNatural = (\s -> (length s, read s)) <$> pMany1 digit  pInteger :: Parser Integer pInteger = pNatural

@@ -1,4 +1,5 @@ {-# Language OverloadedStrings #-}+{-# Language DefaultSignatures #-} module Data.Geometry.Ipe.PathParser where  import           Data.Bifunctor@@ -14,29 +15,31 @@ import           Data.Semigroup import           Data.Text (Text) import qualified Data.Text as T-import           Numeric import           Text.Parsec.Error (messageString, errorMessages)   ----------------------------------------------------------------------- -- | Represent stuff that can be used as a coordinate in ipe. (similar to show/read) -class Num r => Coordinate r where-    fromSeq :: Integer -> Maybe Integer -> r--defaultFromSeq :: (Ord r, Fractional r) => Integer -> Maybe Integer -> r-defaultFromSeq x Nothing  = fromInteger x-defaultFromSeq x (Just y) = let x'        = fromInteger x-                                y'        = fromInteger y-                                asDecimal = head . dropWhile (>= 1) . iterate (* 0.1)-                            in signum x' * (abs x' + asDecimal y')+class Fractional r => Coordinate r where+    -- reads a coordinate. The input is an integer representing the+    -- part before the decimal point, and a length and an integer+    -- representing the part after the decimal point+    fromSeq :: Integer -> Maybe (Int, Integer) -> r+    default fromSeq :: (Ord r, Fractional r) => Integer -> Maybe (Int, Integer) -> r+    fromSeq = defaultFromSeq -instance Coordinate Double where-  fromSeq = defaultFromSeq+defaultFromSeq                :: (Ord r, Fractional r)+                              => Integer -> Maybe (Int, Integer) -> r+defaultFromSeq x Nothing      = fromInteger x+defaultFromSeq x (Just (l,y)) = let x'          = fromInteger x+                                    y'          = fromInteger y+                                    asDecimal a =  a * (0.1 ^ l)+                                    z           = if x' < 0 then (-1) else 1+                                in z * (abs x' + asDecimal y') -instance Coordinate (Ratio Integer) where-    fromSeq x  Nothing = fromInteger x-    fromSeq x (Just y) = fst . head $ readSigned readFloat (show x ++ "." ++ show y)+instance Coordinate Double+instance Coordinate (Ratio Integer)  ----------------------------------------------------------------------- -- | Running the parsers@@ -52,13 +55,16 @@  -- Collect errors data Either' l r = Left' l | Right' r deriving (Show,Eq)++instance (Semigroup l, Semigroup r) => Semigroup (Either' l r) where+  (Left' l)  <> (Left' l')  = Left' $ l <> l'+  (Left' l)  <> _           = Left' l+  _          <> (Left' l')  = Left' l'+  (Right' r) <> (Right' r') = Right' $ r <> r'+ instance (Semigroup l, Semigroup r, Monoid r) => Monoid (Either' l r) where   mempty = Right' mempty-  (Left' l)  `mappend` (Left' l')  = Left' $ l <> l'-  (Left' l)  `mappend` _           = Left' l-  _          `mappend` (Left' l')  = Left' l'-  (Right' r) `mappend` (Right' r') = Right' $ r <> r'-+  mappend = (<>) either' :: (l -> a) -> (r -> a) -> Either' l r -> a either' lf _  (Left' l)  = lf l either' _  rf (Right' r) = rf r@@ -130,8 +136,9 @@ pCoordinate :: Coordinate r => Parser r pCoordinate = fromSeq <$> pInteger <*> pDecimal               where-                pDecimal = pMaybe (pChar '.' *> pInteger)+                pDecimal  = pMaybe (pChar '.' *> pPaddedNatural) + pRectangle :: Coordinate r => Parser (Rectangle () r) pRectangle = (\p q -> box (ext p) (ext q)) <$> pPoint                                            <*  pWhiteSpace@@ -144,9 +151,9 @@  -- | Generate a matrix from a list of 6 coordinates. mkMatrix               :: Coordinate r => [r] -> Matrix 3 3 r-mkMatrix [a,b,c,d,e,f] = Matrix $ v3 (v3 a c e)-                                     (v3 b d f)-                                     (v3 0 0 1)+mkMatrix [a,b,c,d,e,f] = Matrix $ Vector3 (Vector3 a c e)+                                          (Vector3 b d f)+                                          (Vector3 0 0 1)                            -- We need the matrix in the following order:                          -- 012                          -- 345

@@ -22,15 +22,19 @@                                , ipeReadObject                                , ipeReadAttrs                                , ipeReadRec++                               , Coordinate(..)                                ) where +import           Control.Applicative((<|>)) import           Control.Lens hiding (Const, rmap) import qualified Data.ByteString as B+import           Data.Colour.SRGB (RGB(..)) import           Data.Either (rights) import           Data.Ext import           Data.Geometry.Box import           Data.Geometry.Ipe.Attributes-import           Data.Geometry.Ipe.ParserPrimitives (pInteger)+import           Data.Geometry.Ipe.ParserPrimitives (pInteger, pWhiteSpace) import           Data.Geometry.Ipe.PathParser import           Data.Geometry.Ipe.Types import           Data.Geometry.Point@@ -40,14 +44,14 @@ import qualified Data.List as L import qualified Data.List.NonEmpty as NE import           Data.Maybe (fromMaybe, mapMaybe)-import           Data.Monoid import           Data.Proxy import qualified Data.Seq2 as S2+import           Data.Semigroup ((<>)) import           Data.Singletons-import qualified Data.Text as T import           Data.Text (Text)+import qualified Data.Text as T import qualified Data.Traversable as Tr-import           Data.Vinyl+import           Data.Vinyl hiding (Label) import           Data.Vinyl.Functor import           Data.Vinyl.TypeLevel import           Text.XML.Expat.Tree@@ -59,19 +63,22 @@   -- | Given a file path, tries to read an ipe file-readRawIpeFile :: Coordinate r => FilePath -> IO (Either ConversionError (IpeFile r))+readRawIpeFile :: (Coordinate r, Eq r)+               => FilePath -> IO (Either ConversionError (IpeFile r)) readRawIpeFile = fmap fromIpeXML . B.readFile   -- | Given a file path, tries to read an ipe file. This function applies all -- matrices to objects.-readIpeFile :: Coordinate r => FilePath -> IO (Either ConversionError (IpeFile r))+readIpeFile :: (Coordinate r, Eq r)+            => FilePath -> IO (Either ConversionError (IpeFile r)) readIpeFile = fmap (bimap id applyMatrices) . readRawIpeFile   -- | Since most Ipe file contain only one page, we provide a shortcut for that -- as well. This function applies all matrices.-readSinglePageFile :: Coordinate r => FilePath -> IO (Either ConversionError (IpePage r))+readSinglePageFile :: (Coordinate r, Eq r)+                   => FilePath -> IO (Either ConversionError (IpePage r)) readSinglePageFile = fmap f . readIpeFile   where     f (Left e)  = Left e@@ -151,9 +158,17 @@ instance Coordinate r => IpeReadText (Rectangle () r) where   ipeReadText = readRectangle -instance IpeReadText IpeColor where-  ipeReadText = fmap IpeColor . ipeReadTextWith Right+instance Coordinate r => IpeReadText (RGB r) where+  ipeReadText = runParser (pRGB <|> pGrey)+    where+      pGrey = (\c -> RGB c c c) <$> pCoordinate+      pRGB  = RGB <$> pCoordinate <* pWhiteSpace+                  <*> pCoordinate <* pWhiteSpace+                  <*> pCoordinate +instance Coordinate r => IpeReadText (IpeColor r) where+  ipeReadText = fmap IpeColor . ipeReadTextWith ipeReadText+ instance Coordinate r => IpeReadText (IpePen r) where   ipeReadText = fmap IpePen . ipeReadTextWith readCoordinate @@ -164,7 +179,7 @@ instance Coordinate r => IpeReadText [Operation r] where   ipeReadText = readPathOperations -instance Coordinate r => IpeReadText (NE.NonEmpty (PathSegment r)) where+instance (Coordinate r, Eq r) => IpeReadText (NE.NonEmpty (PathSegment r)) where   ipeReadText t = ipeReadText t >>= fromOpsN     where       fromOpsN xs = case fromOps xs of@@ -179,7 +194,7 @@       fromOps' _ []             = Left "Found only a MoveTo operation"       fromOps' s (LineTo q:ops) = let (ls,xs) = span' _LineTo ops                                       pts  = map ext $ s:q:mapMaybe (^?_LineTo) ls-                                      poly = Polygon.fromPoints pts+                                      poly = Polygon.fromPoints . dropRepeats $ pts                                       pl   = fromPoints pts                                   in case xs of                                        (ClosePath : xs') -> PolygonPath poly   <<| xs'@@ -190,7 +205,11 @@        x <<| xs = (x:) <$> fromOps xs -instance Coordinate r => IpeReadText (Path r) where++dropRepeats :: Eq a => [a] -> [a]+dropRepeats = map head . L.group++instance (Coordinate r, Eq r) => IpeReadText (Path r) where   ipeReadText = fmap (Path . S2.viewL1FromNonEmpty) . ipeReadText  --------------------------------------------------------------------------------@@ -253,17 +272,17 @@ ipeReadAttrs _ _ = fmap Attrs . ipeReadRec (Proxy :: Proxy f) (Proxy :: Proxy ats)  -testSym :: B.ByteString-testSym = "<use name=\"mark/disk(sx)\" pos=\"320 736\" size=\"normal\" stroke=\"black\"/>"+-- testSym :: B.ByteString+-- testSym = "<use name=\"mark/disk(sx)\" pos=\"320 736\" size=\"normal\" stroke=\"black\"/>"     -- readAttrsFromXML :: B.ByteString -> Either -readSymAttrs :: Either ConversionError (IpeAttributes IpeSymbol Double)-readSymAttrs = readXML testSym-               >>= ipeReadAttrs (Proxy :: Proxy IpeSymbol) (Proxy :: Proxy Double)+-- readSymAttrs :: Either ConversionError (IpeAttributes IpeSymbol Double)+-- readSymAttrs = readXML testSym+--                >>= ipeReadAttrs (Proxy :: Proxy IpeSymbol) (Proxy :: Proxy Double)   @@ -301,7 +320,7 @@     unT (Text t) = Right t     unT _        = Left "allText: Expected Text, found an Element" -instance Coordinate r => IpeRead (Path r) where+instance (Coordinate r, Eq r) => IpeRead (Path r) where   ipeRead (Element "path" _ chs) = allText chs >>= ipeReadText   ipeRead _                      = Left "path: expected element, found text" @@ -333,7 +352,7 @@   ipeRead (Element "image" ats _) = Image () <$> (lookup' "rect" ats >>= ipeReadText)   ipeRead _                       = Left "Image: Element expected, text found" -instance Coordinate r => IpeRead (IpeObject r) where+instance (Coordinate r, Eq r) => IpeRead (IpeObject r) where   ipeRead x = firstRight [ IpeUse       <$> ipeReadObject (Proxy :: Proxy IpeSymbol) r x                          , IpePath      <$> ipeReadObject (Proxy :: Proxy Path)      r x                          , IpeGroup     <$> ipeReadObject (Proxy :: Proxy Group)     r x@@ -348,8 +367,8 @@ firstRight = maybe (Left "No matching object") Right . firstOf (traverse._Right)  -instance Coordinate r => IpeRead (Group r) where-  ipeRead (Element "group" _ chs) = Group <$> mapM ipeRead chs+instance (Coordinate r, Eq r) => IpeRead (Group r) where+  ipeRead (Element "group" _ chs) = Right . Group . rights . map ipeRead $ chs   ipeRead _                       = Left "ipeRead Group: expected Element, found Text"  @@ -366,7 +385,7 @@  -- TODO: this instance throws away all of our error collecting (and is pretty -- slow/stupid since it tries parsing all children with all parsers)-instance Coordinate r => IpeRead (IpePage r) where+instance (Coordinate r, Eq r) => IpeRead (IpePage r) where   ipeRead (Element "page" _ chs) = Right $ IpePage (readAll chs) (readAll chs) (readAll chs)   ipeRead _                      = Left "page: Element expected, text found"       -- withDef   :: b -> Either a b -> Either c b@@ -381,7 +400,7 @@ readAll   = rights . map ipeRead  -instance Coordinate r => IpeRead (IpeFile r) where+instance (Coordinate r, Eq r) => IpeRead (IpeFile r) where   ipeRead (Element "ipe" _ chs) = case readAll chs of                                     []  -> Left "Ipe: no pages found"                                     pgs -> Right $ IpeFile Nothing [] (NE.fromList pgs)@@ -389,59 +408,6 @@   -testz :: Either ConversionError (IpeObject Double)-testz = (bimap (T.pack . show) id $ parse' defaultParseOptions testSym)-               >>= ipeRead -- Object (Proxy :: Proxy IpeSymbol) (Proxy :: Proxy Double)---------------  -- ipeReadAttrs (Element _ ats _) =--instance Coordinate r => IpeReadText (PolyLine 2 () r) where-  ipeReadText t = readPathOperations t >>= fromOps-    where-      fromOps (MoveTo p:LineTo q:ops) = (\ps -> fromPoints' $ [p,q] ++ ps)-                                     <$> validateAll "Expected LineTo p" _LineTo ops-      fromOps _                       = Left "Expected MoveTo p:LineTo q:... "--validateAll         :: ConversionError -> Prism' (Operation r) (Point 2 r) -> [Operation r]-                    -> Either ConversionError [Point 2 r]-validateAll err fld = bimap T.unlines id . validateAll' err fld---validateAll' :: err -> Prism' (Operation r) (Point 2 r) -> [Operation r]-               -> Either [err] [Point 2 r]-validateAll' err field = toEither . foldr (\op' res -> f op' <> res) (Right' [])-  where-    f op' = maybe (Left' [err]) (\p -> Right' [p]) $ op' ^? field-    toEither = either' Left Right---- This is a bit of a hack-instance Coordinate r => IpeRead (PolyLine 2 () r) where-  ipeRead (Element "path" _ ts) = ipeReadText . T.unlines . map unText $ ts-                                    -- apparently hexpat already splits the text into lines-  ipeRead _                     = Left "iperead: no polyline."--unText          :: Node t t1 -> t1-unText (Text t) = t-unText _        = error "unText: element found, text expected"--instance Coordinate r => IpeRead (PathSegment r) where-  ipeRead = fmap PolyLineSegment . ipeRead--testP :: B.ByteString-testP = "<path stroke=\"black\">\n128 656 m\n224 768 l\n304 624 l\n432 752 l\n</path>"   --------------------------------------------------------------------------------

@@ -8,7 +8,7 @@  import           Control.Lens import           Data.Proxy-import           Data.Vinyl+import           Data.Vinyl hiding (Label)  import           Data.Ext import           Data.Geometry.Box(Rectangle)@@ -50,7 +50,7 @@ type instance NumType   (Image r) = r type instance Dimension (Image r) = 2 -instance Num r => IsTransformable (Image r) where+instance Fractional r => IsTransformable (Image r) where   transformBy t = over rect (transformBy t)  --------------------------------------------------------------------------------@@ -68,10 +68,10 @@ type instance NumType   (MiniPage r) = r type instance Dimension (MiniPage r) = 2 -instance Num r => IsTransformable (TextLabel r) where+instance Fractional r => IsTransformable (TextLabel r) where   transformBy t (Label txt p) = Label txt (transformBy t p) -instance Num r => IsTransformable (MiniPage r) where+instance Fractional r => IsTransformable (MiniPage r) where   transformBy t (MiniPage txt p w) = MiniPage txt (transformBy t p) w  width                  :: MiniPage t -> t@@ -90,7 +90,7 @@ type instance NumType   (IpeSymbol r) = r type instance Dimension (IpeSymbol r) = 2 -instance Num r => IsTransformable (IpeSymbol r) where+instance Fractional r => IsTransformable (IpeSymbol r) where   transformBy t = over symbolPoint (transformBy t)  @@ -118,7 +118,7 @@ type instance NumType   (PathSegment r) = r type instance Dimension (PathSegment r) = 2 -instance Num r => IsTransformable (PathSegment r) where+instance Fractional r => IsTransformable (PathSegment r) where   transformBy t (PolyLineSegment p) = PolyLineSegment $ transformBy t p   transformBy t (PolygonPath p)     = PolygonPath $ transformBy t p   transformBy _ _                   = error "transformBy: not implemented yet"@@ -132,7 +132,7 @@ type instance NumType   (Path r) = r type instance Dimension (Path r) = 2 -instance Num r => IsTransformable (Path r) where+instance Fractional r => IsTransformable (Path r) where   transformBy t (Path s) = Path $ fmap (transformBy t) s  -- | type that represents a path in ipe.@@ -161,9 +161,9 @@   AttrMap r Pin             = PinType   AttrMap r Transformations = TransformationTypes -  AttrMap r Stroke = IpeColor+  AttrMap r Stroke = IpeColor r   AttrMap r Pen    = IpePen r-  AttrMap r Fill   = IpeColor+  AttrMap r Fill   = IpeColor r   AttrMap r Size   = IpeSize r    AttrMap r Dash     = IpeDash r@@ -203,7 +203,7 @@ type instance NumType   (Group r) = r type instance Dimension (Group r) = 2 -instance Num r => IsTransformable (Group r) where+instance Fractional r => IsTransformable (Group r) where   transformBy t (Group s) = Group $ fmap (transformBy t) s  @@ -256,7 +256,7 @@ instance ToObject IpeSymbol  where ipeObject' s a = IpeUse       (s :+ a) instance ToObject Path       where ipeObject' p a = IpePath      (p :+ a) -instance Num r => IsTransformable (IpeObject r) where+instance Fractional r => IsTransformable (IpeObject r) where   transformBy t (IpeGroup i)     = IpeGroup     $ i&core %~ transformBy t   transformBy t (IpeImage i)     = IpeImage     $ i&core %~ transformBy t   transformBy t (IpeTextLabel i) = IpeTextLabel $ i&core %~ transformBy t@@ -288,6 +288,17 @@     s (IpeUse i)       a = IpeUse       $ i&select .~ a     s (IpePath i)      a = IpePath      $ i&select .~ a +-- | collect all non-group objects+flattenGroups :: [IpeObject r] -> [IpeObject r]+flattenGroups = concatMap flattenGroups'+  where+    flattenGroups'                              :: IpeObject r -> [IpeObject r]+    flattenGroups' (IpeGroup (Group gs :+ ats)) =+      map (applyAts ats) . concatMap flattenGroups' $ gs+        where+          applyAts _ = id+    flattenGroups' o                            = [o]+ --------------------------------------------------------------------------------  @@ -372,7 +383,7 @@     (mm,ats') = takeAttr (Proxy :: Proxy AT.Matrix) ats  -- | Applies the matrix to an ipe object if it has one.-applyMatrix                  :: Num r => IpeObject r -> IpeObject r+applyMatrix                  :: Fractional r => IpeObject r -> IpeObject r applyMatrix (IpeGroup i)     = IpeGroup . applyMatrix'                              $ i&core.groupItems.traverse %~ applyMatrix                              -- note that for a group we first (recursively)@@ -384,10 +395,10 @@ applyMatrix (IpeUse i)       = IpeUse       $ applyMatrix' i applyMatrix (IpePath i)      = IpePath      $ applyMatrix' i -applyMatrices   :: Num r => IpeFile r -> IpeFile r+applyMatrices   :: Fractional r => IpeFile r -> IpeFile r applyMatrices f = f&pages.traverse %~ applyMatricesPage -applyMatricesPage   :: Num r => IpePage r -> IpePage r+applyMatricesPage   :: Fractional r => IpePage r -> IpePage r applyMatricesPage p = p&content.traverse %~ applyMatrix  

@@ -6,6 +6,7 @@ import           Control.Lens ((^.),(^..),(.~),(&), to) import qualified Data.ByteString as B import qualified Data.ByteString.Char8 as C+import           Data.Colour.SRGB (RGB(..)) import           Data.Ext import           Data.Fixed import qualified Data.Foldable as F@@ -13,11 +14,10 @@ import           Data.Geometry.Ipe.Attributes import qualified Data.Geometry.Ipe.Attributes as IA import           Data.Geometry.Ipe.Types-import qualified Data.Geometry.Ipe.Types as IT import           Data.Geometry.LineSegment import           Data.Geometry.Point import           Data.Geometry.PolyLine-import           Data.Geometry.Polygon (SimplePolygon, outerBoundary)+import           Data.Geometry.Polygon (Polygon, outerBoundary, holeList, asSimplePolygon) import qualified Data.Geometry.Transformation as GT import           Data.Geometry.Vector import           Data.Maybe (catMaybes, mapMaybe, fromMaybe)@@ -26,17 +26,15 @@ import           Data.Semigroup import qualified Data.Seq2 as S2 import           Data.Singletons+import           Data.Text (Text) import qualified Data.Text as T import qualified Data.Text as Text-import           Data.Text (Text)-import qualified Data.Traversable as Tr-import           Data.Vinyl+import           Data.Vinyl hiding (Label) import           Data.Vinyl.Functor import           Data.Vinyl.TypeLevel import           System.IO (hPutStrLn,stderr) import           Text.XML.Expat.Format (format') import           Text.XML.Expat.Tree- --------------------------------------------------------------------------------  -- | Given a prism to convert something of type g into an ipe file, a file path,@@ -92,10 +90,16 @@ class IpeWrite t where   ipeWrite :: t -> Maybe (Node Text Text) +instance (IpeWrite l, IpeWrite r) => IpeWrite (Either l r) where+  ipeWrite = either ipeWrite ipeWrite  instance IpeWriteText (Apply f at) => IpeWriteText (Attr f at) where   ipeWriteText att = _getAttr att >>= ipeWriteText +instance (IpeWriteText l, IpeWriteText r) => IpeWriteText (Either l r) where+  ipeWriteText = either ipeWriteText ipeWriteText++ -- | Functon to write all attributes in a Rec ipeWriteAttrs           :: ( AllSatisfy IpeAttrName rs                            , RecAll (Attr f) rs IpeWriteText@@ -175,16 +179,19 @@   ipeWriteText Horizontal = Just "h"   ipeWriteText Vertical   = Just "v" +instance IpeWriteText r => IpeWriteText (RGB r) where+  ipeWriteText (RGB r g b) = unwords' . map ipeWriteText $ [r,g,b]+ deriving instance IpeWriteText r => IpeWriteText (IpeSize  r) deriving instance IpeWriteText r => IpeWriteText (IpePen   r)-deriving instance IpeWriteText IpeColor+deriving instance IpeWriteText r => IpeWriteText (IpeColor r)  instance IpeWriteText r => IpeWriteText (IpeDash r) where   ipeWriteText (DashNamed t) = Just t   ipeWriteText (DashPattern xs x) = (\ts t -> mconcat [ "["                                                       , Text.intercalate " " ts                                                       , "] ", t ])-                                    <$> Tr.mapM ipeWriteText xs+                                    <$> mapM ipeWriteText xs                                     <*> ipeWriteText x  instance IpeWriteText FillType where@@ -196,7 +203,7 @@                                                             <*> ipeWriteText s  instance IpeWriteText r => IpeWriteText (Path r) where-  ipeWriteText = fmap concat' . Tr.sequence . fmap ipeWriteText . _pathSegments+  ipeWriteText = fmap concat' . sequence . fmap ipeWriteText . _pathSegments     where       concat' = F.foldr1 (\t t' -> t <> "\n" <> t') @@ -227,28 +234,37 @@   instance IpeWriteText r => IpeWriteText (Operation r) where-  ipeWriteText (MoveTo p)      = unwords' [ ipeWriteText p, Just "m"]-  ipeWriteText (LineTo p)      = unwords' [ ipeWriteText p, Just "l"]-  ipeWriteText (CurveTo p q r) = unwords' [ ipeWriteText p-                                          , ipeWriteText q-                                          , ipeWriteText r, Just "m"]-  ipeWriteText (Ellipse m)     = unwords' [ ipeWriteText m, Just "e"]-  -- TODO: The rest-  ipeWriteText ClosePath       = Just "h"+  ipeWriteText (MoveTo p)         = unwords' [ ipeWriteText p, Just "m"]+  ipeWriteText (LineTo p)         = unwords' [ ipeWriteText p, Just "l"]+  ipeWriteText (CurveTo p q r)    = unwords' [ ipeWriteText p+                                             , ipeWriteText q+                                             , ipeWriteText r, Just "c"]+  ipeWriteText (QCurveTo p q)     = unwords' [ ipeWriteText p+                                             , ipeWriteText q, Just "q"]+  ipeWriteText (Ellipse m)        = unwords' [ ipeWriteText m, Just "e"]+  ipeWriteText (ArcTo m p)        = unwords' [ ipeWriteText m+                                             , ipeWriteText p, Just "a"]+  ipeWriteText (Spline pts)       = unlines' $ map ipeWriteText pts <> [Just "s"]+  ipeWriteText (ClosedSpline pts) = unlines' $ map ipeWriteText pts <> [Just "u"]+  ipeWriteText ClosePath          = Just "h"   instance IpeWriteText r => IpeWriteText (PolyLine 2 () r) where-  ipeWriteText pl = case pl^..points.Tr.traverse.core of+  ipeWriteText pl = case pl^..points.traverse.core of     (p : rest) -> unlines' . map ipeWriteText $ MoveTo p : map LineTo rest     -- the polyline type guarantees that there is at least one point -instance IpeWriteText r => IpeWriteText (SimplePolygon () r) where-  ipeWriteText pg = case pg^..outerBoundary.to F.toList.Tr.traverse.core of-    (p : rest) -> unlines' . map ipeWriteText $ MoveTo p : map LineTo rest ++ [ClosePath]-    _          -> Nothing+instance IpeWriteText r => IpeWriteText (Polygon t () r) where+  ipeWriteText pg = fmap mconcat . traverse f $ asSimplePolygon pg : holeList pg+    where+      f pg' = case pg'^..outerBoundary.traverse.core of+        (p : rest) -> unlines' . map ipeWriteText+                    $ MoveTo p : map LineTo rest ++ [ClosePath]+        _          -> Nothing     -- TODO: We are not really guaranteed that there is at least one point, it would     -- be nice if the type could guarantee that. + instance IpeWriteText r => IpeWriteText (PathSegment r) where   ipeWriteText (PolyLineSegment p) = ipeWriteText p   ipeWriteText (PolygonPath     p) = ipeWriteText p@@ -369,7 +385,7 @@ instance (IpeWriteText r, IpeWrite p) => IpeWrite (PolyLine 2 p r) where   ipeWrite p = ipeWrite path     where-      path = fromPolyLine $ p & points.Tr.traverse.extra .~ ()+      path = fromPolyLine $ p & points.traverse.extra .~ ()       -- TODO: Do something with the p's  fromPolyLine :: PolyLine 2 () r -> Path r

@@ -1,9 +1,10 @@+{-# LANGUAGE UndecidableInstances  #-} {-# LANGUAGE ScopedTypeVariables #-} module Data.Geometry.KDTree where  import           Control.Lens hiding (imap, element, Empty, (:<)) import           Data.BinaryTree-import           Data.Coerce+import           Unsafe.Coerce(unsafeCoerce) import           Data.Ext import qualified Data.Foldable as F import           Data.Geometry.Box@@ -70,11 +71,11 @@ -- | Expects the input to be a set, i.e. no duplicates -- -- running time: \(O(n \log n)\)-buildKDTree :: (Arity d, KnownNat d, Index' 0 d, Ord r)+buildKDTree :: (Arity d, 1 <= d, Ord r)             => [Point d r :+ p] -> KDTree d p r buildKDTree = maybe Empty (Tree . buildKDTree') . NonEmpty.nonEmpty -buildKDTree' :: (Arity d, KnownNat d, Index' 0 d, Ord r)+buildKDTree' :: (Arity d, 1 <= d, Ord r)              => NonEmpty.NonEmpty (Point d r :+ p) -> KDTree' d p r buildKDTree' = KDT . addBoxes . build (Coord 1) . toPointSet . Seq.fromNonEmpty   where     -- compute one tree with bounding boxes, then merge them together@@ -102,7 +103,7 @@                   in (f p, p^.core) `compare` (f q, q^.core)  -build      :: (Index' 0 d, Arity d, KnownNat d, Ord r)+build      :: (1 <= d, Arity d, Ord r)            => Coord d            -> PointSet (LSeq 1) d p r            -> BinLeafTree (Split' d r) (Point d r :+ p)@@ -179,8 +180,8 @@     unzip' = bimap vectorFromListUnsafe vectorFromListUnsafe . unzip . F.toList  -asSingleton   :: (Index' 0 d, Arity d) => PointSet (LSeq 1) d p r+asSingleton   :: (1 <= d, Arity d) => PointSet (LSeq 1) d p r               -> Either (Point d r :+ p) (PointSet (LSeq 2) d p r) asSingleton v = case Seq.viewl $ v^.element (C :: C 0) of-                  _ :< _ Seq.:<< _ -> Right $ coerce v+                  _ :< _ Seq.:<< _ -> Right $ unsafeCoerce v                   p :< _           -> Left p -- only one element

@@ -1,18 +1,34 @@ {-# OPTIONS_GHC -fno-warn-orphans #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE UnicodeSyntax #-} module Data.Geometry.Line( module Data.Geometry.Line.Internal                          ) where -import Data.Geometry.Line.Internal-import Data.Geometry.LineSegment-import Data.Geometry.Box-import Data.Geometry.Properties-import Data.Geometry.Point-import Data.Geometry.Boundary-import Data.Geometry.Transformation-import Data.Geometry.Vector+import           Control.Lens ((^.), re, bimap)+import           Data.Ext+import           Data.Geometry.Boundary+import           Data.Geometry.Box+import           Data.Geometry.Line.Internal+import           Data.Geometry.LineSegment+import           Data.Geometry.Point+import           Data.Geometry.Properties+import           Data.Geometry.SubLine+import           Data.Geometry.Transformation+import           Data.Geometry.Vector+import qualified Data.List as L+import           Data.Maybe (mapMaybe)+import           Data.Proxy+import           Data.UnBounded+import           Data.Vinyl.CoRec+import           Data.Vinyl.Core+import           Data.Vinyl.Lens+import           GHC.TypeLits +--------------------------------------------------------------------------------+ -- | Lines are transformable, via line segments-instance (Num r, AlwaysTruePFT d) => IsTransformable (Line d r) where+instance (Fractional r, Arity d, Arity (d + 1)) => IsTransformable (Line d r) where   transformBy t = supportingLine . transformPointFunctor t . toLineSegment'     where       toLineSegment' :: (Num r, Arity d) => Line d r -> LineSegment d () r@@ -23,17 +39,50 @@   [ NoIntersection, Point 2 r, (Point 2 r, Point 2 r) , LineSegment 2 () r]  --- instance (Eq r, Fractional r)---          => (Line 2 r) `IsIntersectableWith` (Boundary (Rectangle p r)) where---   nonEmptyIntersection = defaultNonEmptyIntersection+instance (Ord r, Fractional r)+         => (Line 2 r) `IsIntersectableWith` (Boundary (Rectangle p r)) where+  nonEmptyIntersection = defaultNonEmptyIntersection ---   _    `intersect` (Boundary Empty) = coRec NoIntersection---   line `intersect` (Boundary rect)  = error "TODO"---     where---       (t,r,b,l) = sides' rect+  line' `intersect` (Boundary rect)  = case asA' segP of+      [sl'] -> coRec . bimap id unVal $ sl'^.re _SubLine+      []    -> case nub' . map (fmap unVal) $ asA' pointP of+        [p]   -> coRec p+        [p,q] -> coRec (p,q)+        _     -> coRec NoIntersection+      _     -> error "intersect; ine x boundary rect; absurd"+    where+      (t,r,b,l) = sides' rect+      ints = map (\s -> sl `intersect` toSL s) [t,r,b,l] ---       ints = map (line `intersect`) [t,r,b,l]+      nub' = map head . L.group . L.sort +      sl = fromLine line'+      -- wrap a segment into an potentially unbounded subline+      toSL s = bimap (const ()) Val $ s^._SubLine +      unVal (Val x) = x+      unVal _       = error "intersect; line x boundary rect: unVal Unbounded"++      asA'    :: (t ∈ IntersectionOf (SubLine 2 () (UnBounded r))+                                     (SubLine 2 () (UnBounded r)))+              => proxy t -> [t]+      asA' px = mapMaybe (asA px) ints++      segP   = Proxy :: Proxy (SubLine 2 () (UnBounded r))+      pointP = Proxy :: Proxy (Point 2      (UnBounded r))++ type instance IntersectionOf (Line 2 r) (Rectangle p r) =   [ NoIntersection, Point 2 r, LineSegment 2 () r]+++instance (Ord r, Fractional r)+         => (Line 2 r) `IsIntersectableWith` (Rectangle p r) where+  nonEmptyIntersection = defaultNonEmptyIntersection++  line' `intersect` rect  = match (line' `intersect` (Boundary rect)) $+       (H $ \NoIntersection -> coRec NoIntersection)+    :& (H $ \p@(Point2 _ _) -> coRec p)+    :& (H $ \(p,q)          -> coRec $ ClosedLineSegment (ext p) (ext q))+    :& (H $ \s              -> coRec s)+    :& RNil

@@ -13,7 +13,7 @@ import           Data.Ord (comparing) import qualified Data.Traversable as T import           Data.Vinyl-import           Frames.CoRec+import           Data.Vinyl.CoRec import           GHC.Generics (Generic)  @@ -46,14 +46,14 @@ lineThrough p q = Line p (q .-. p)  verticalLine   :: Num r => r -> Line 2 r-verticalLine x = Line (point2 x 0) (v2 0 1)+verticalLine x = Line (point2 x 0) (Vector2 0 1)  horizontalLine   :: Num r => r -> Line 2 r-horizontalLine y = Line (point2 0 y) (v2 1 0)+horizontalLine y = Line (point2 0 y) (Vector2 1 0)  -- | Given a line l with anchor point p, get the line perpendicular to l that also goes through p. perpendicularTo                           :: Num r => Line 2 r -> Line 2 r-perpendicularTo (Line p ~(Vector2 vx vy)) = Line p (v2 (-vy) vx)+perpendicularTo (Line p ~(Vector2 vx vy)) = Line p (Vector2 (-vy) vx)   @@ -89,6 +89,11 @@ p `onLine` (Line q v) = p == q || (p .-. q) `isScalarMultipleOf` v  +-- | Specific 2d version of testing if apoint lies on a line.+onLine2 :: (Ord r, Num r) => Point 2 r -> Line 2 r -> Bool+p `onLine2` (Line q v) = ccw p q (q .+^ v) == CoLinear++ -- | The intersection of two lines is either: NoIntersection, a point or a line. type instance IntersectionOf (Line 2 r) (Line 2 r) = [ NoIntersection                                                      , Point 2 r@@ -146,7 +151,7 @@  -- | Create a line from the linear function ax + b fromLinearFunction     :: Num r => r -> r -> Line 2 r-fromLinearFunction a b = Line (point2 0 b) (v2 1 a)+fromLinearFunction a b = Line (point2 0 b) (Vector2 1 a)  -- | get values a,b s.t. the input line is described by y = ax + b. -- returns Nothing if the line is vertical@@ -194,3 +199,25 @@ bisector p q = let v = q .-. p                    h = p .+^ (v ^/ 2)                in perpendicularTo (Line h v)+++-- | Compares the lines on slope. Vertical lines are considered larger than+-- anything else.+--+-- >>> (Line origin (Vector2 5 1)) `cmpSlope` (Line origin (Vector2 3 3))+-- LT+-- >>> (Line origin (Vector2 5 1)) `cmpSlope` (Line origin (Vector2 (-3) 3))+-- GT+-- >>> (Line origin (Vector2 5 1)) `cmpSlope` (Line origin (Vector2 0 1))+-- LT+cmpSlope :: (Num r, Ord r) => Line 2 r -> Line 2 r -> Ordering+(Line _ u) `cmpSlope` (Line _ v) = case ccw origin (f u) (f v) of+                                     CCW      -> LT+                                     CW       -> GT+                                     CoLinear -> EQ+  where+    f w@(Vector2 x y) = Point $ case (x `compare` 0, y >= 0) of+                                  (GT,_)    -> w+                                  (EQ,True) -> w+                                  _         -> (-1) *^ w+                                  -- x < 0, or (x==0 and y <0 ; i.e. a vertical line)

@@ -1,5 +1,11 @@ {-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE UndecidableInstances #-}+{-|+Module    : Data.Geometry.LineSegment+Description: Line segment data type and some basic functions on line segments+Copyright : (c) Frank Staals+License : See LICENCE file+-} module Data.Geometry.LineSegment( LineSegment                                 , pattern LineSegment                                 , pattern LineSegment'@@ -17,24 +23,25 @@                                 , flipSegment                                 ) where -import           Data.Ord(comparing)-import           Control.Arrow((&&&))+import           Control.Arrow ((&&&)) import           Control.Lens import           Data.Bifunctor-import           Data.Semigroup import           Data.Ext+import qualified Data.Foldable as F import           Data.Geometry.Box.Internal-import           Data.Geometry.Interval+import           Data.Geometry.Interval hiding (width, midPoint) import           Data.Geometry.Line.Internal import           Data.Geometry.Point import           Data.Geometry.Properties import           Data.Geometry.SubLine import           Data.Geometry.Transformation import           Data.Geometry.Vector-import           Data.Vinyl+import           Data.Ord (comparing)+import           Data.Semigroup import           Data.UnBounded-import           Frames.CoRec-import qualified Data.Foldable as F+import           Data.Vinyl+import           Data.Vinyl.CoRec+import           GHC.TypeLits  -------------------------------------------------------------------------------- -- * d-dimensional LineSegments@@ -43,32 +50,33 @@ -- | Line segments. LineSegments have a start and end point, both of which may -- contain additional data of type p. We can think of a Line-Segment being defined as -----  data LineSegment d p r = LineSegment (EndPoint (Point d r :+ p))---                                       (EndPoint (Point d r :+ p))+--+-- >>>  data LineSegment d p r = LineSegment (EndPoint (Point d r :+ p)) (EndPoint (Point d r :+ p)) newtype LineSegment d p r = GLineSegment { _unLineSeg :: Interval p (Point d r)}  makeLenses ''LineSegment  -- | Pattern that essentially models the line segment as a: ----- data LineSegment d p r = LineSegment (EndPoint (Point d r :+ p))---                                      (EndPoint (Point d r :+ p))+-- >>> data LineSegment d p r = LineSegment (EndPoint (Point d r :+ p)) (EndPoint (Point d r :+ p)) pattern LineSegment           :: EndPoint (Point d r :+ p)                               -> EndPoint (Point d r :+ p)                               -> LineSegment d p r pattern LineSegment       s t = GLineSegment (Interval s t)+{-# COMPLETE LineSegment #-}  -- | Gets the start and end point, but forgetting if they are open or closed. pattern LineSegment'          :: Point d r :+ p                               -> Point d r :+ p                               -> LineSegment d p r pattern LineSegment'      s t <- ((^.start) &&& (^.end) -> (s,t))-+{-# COMPLETE LineSegment' #-}  pattern ClosedLineSegment     :: Point d r :+ p                               -> Point d r :+ p                               -> LineSegment d p r pattern ClosedLineSegment s t = GLineSegment (ClosedInterval s t)+{-# COMPLETE ClosedLineSegment #-}  type instance Dimension (LineSegment d p r) = d type instance NumType   (LineSegment d p r) = r@@ -86,6 +94,7 @@  _SubLine :: (Fractional r, Eq r, Arity d) => Iso' (LineSegment d p r) (SubLine d p r) _SubLine = iso segment2SubLine subLineToSegment+{-# INLINE _SubLine #-}  segment2SubLine    :: (Fractional r, Eq r, Arity d)                    => LineSegment d p r -> SubLine d p r@@ -107,6 +116,7 @@ instance (Num r, Arity d) => HasSupportingLine (LineSegment d p r) where   supportingLine s = lineThrough (s^.start.core) (s^.end.core) + instance (Show r, Show p, Arity d) => Show (LineSegment d p r) where   show ~(LineSegment p q) = concat ["LineSegment (", show p, ") (", show q, ")"] @@ -121,7 +131,7 @@ instance Arity d => IsBoxable (LineSegment d p r) where   boundingBox l = boundingBox (l^.start.core) <> boundingBox (l^.end.core) -instance (Num r, AlwaysTruePFT d) => IsTransformable (LineSegment d p r) where+instance (Fractional r, Arity d, Arity (d + 1)) => IsTransformable (LineSegment d p r) where   transformBy = transformPointFunctor  instance Arity d => Bifunctor (LineSegment d) where@@ -207,7 +217,10 @@ p `onSegment` l = let s          = l^.start.core                       t          = l^.end.core                       inRange' x = 0 <= x && x <= 1-                  in maybe False inRange' $ scalarMultiple (p .-. s) (t .-. s)+                  in+                  if s == t -- zero length segment+                  then p == s+                  else maybe False inRange' $ scalarMultiple (p .-. s) (t .-. s)   -- | The left and right end point (or left below right if they have equal x-coords)

@@ -1,132 +1,107 @@ {-# LANGUAGE TemplateHaskell #-}-module Data.Geometry.PlanarSubdivision where+{-# LANGUAGE PartialTypeSignatures #-}+{-# LANGUAGE ScopedTypeVariables #-}+module Data.Geometry.PlanarSubdivision( module Data.Geometry.PlanarSubdivision.Basic+                                      , fromPolygon, fromPolygons+                                      ) where -import           Control.Lens-import qualified Data.BalBST as SS-import           Data.Bifunctor.Apply-import qualified Data.CircularSeq as C+-- import           Algorithms.Geometry.PolygonTriangulation.Triangulate+import           Control.Lens hiding (holes, holesOf, (.=))+import qualified Data.CircularSeq as CSeq import           Data.Ext import qualified Data.Foldable as F-import           Data.Geometry.Interval-import           Data.Geometry.LineSegment-import           Data.Geometry.Point+import           Data.Geometry.Box+import           Data.Geometry.PlanarSubdivision.Basic import           Data.Geometry.Polygon-import           Data.Geometry.Properties-import qualified Data.Map as M-import           Data.PlanarGraph-import           Data.PlaneGraph-import           Data.Semigroup-import           Data.Util+import           Data.Geometry.Vector+import qualified Data.Geometry.Vector as Vec+import           Data.List.NonEmpty (NonEmpty(..)) import qualified Data.Vector as V+import qualified Data.PlaneGraph as PG+import Data.Proxy  --- | Note that the functor instance is in v-data VertexData r v = VertexData { _location :: !(Point 2 r)-                                 , _vData    :: !v-                                 } deriving (Show,Eq,Ord,Functor,Foldable,Traversable)-makeLenses ''VertexData+-- | Construct a planar subdivision from a polygon. Since our PlanarSubdivision+-- models only connected planar subdivisions, this may add dummy/invisible+-- edges.+--+-- running time: \(O(n)\) for a simple polygon, \(O(n\log n)\) for a polygon+-- with holes.+fromPolygon                              :: forall proxy t p f r s.+                                            (Ord r, Fractional r) => proxy s+                                         -> Polygon t p r+                                         -> f -- ^ data inside+                                         -> f -- ^ data outside the polygon+                                         -> PlanarSubdivision s p () f r+fromPolygon p pg@(SimplePolygon _) iD oD = fromSimplePolygon p pg iD oD+fromPolygon _ (MultiPolygon vs hs) iD oD = PlanarSubdivision cs vd dd fd+  where+    wp = Proxy :: Proxy (Wrap s) -instance Bifunctor VertexData where-  bimap f g (VertexData p v) = VertexData (fmap f p) (g v)+    -- the components+    cs = undefined+    cs' = PG.fromSimplePolygon wp (SimplePolygon vs) iD oD+        : map (\h -> PG.fromSimplePolygon wp h oD iD) hs +    vd = undefined+    dd = undefined+    fd = undefined -data EdgeType = Visible | Invisible deriving (Show,Read,Eq,Ord) -data EdgeData e = EdgeData { _edgeType :: !EdgeType-                           , _eData    :: !e-                           } deriving (Show,Eq,Ord,Functor,Foldable,Traversable)-makeLenses ''EdgeData  --- | The Face data consists of the data itself and a list of holes-data FaceData h f = FaceData { _holes :: [h]-                             , _fData :: !f-                             } deriving (Show,Eq,Ord,Functor,Foldable,Traversable)-makeLenses ''FaceData---newtype PlanarSubdivision s v e f r = PlanarSubdivision { _graph ::-    PlanarGraph s Primal_ (VertexData r v) (EdgeData e) (FaceData (Dart s) f) }-      deriving (Show,Eq)-makeLenses ''PlanarSubdivision--instance Functor (PlanarSubdivision s v e f) where-  fmap f s = s&graph.vertexData.traverse.location %~ fmap f+-- | Given a list of *disjoint* polygons, construct a planarsubdivsion+-- representing them. This may create dummy vertices which have no vertex data,+-- hence the 'Maybe p' data type for the vertices.+--+-- running time: \(O(n\log n)\)+fromPolygons           :: (Ord r, Fractional r)+                       => proxy s+                       -> NonEmpty (SimplePolygon p r :+ f)+                       -> f -- ^ data outside the polygons+                       -> PlanarSubdivision s (Maybe p) () f r+fromPolygons px pgs oD = undefined+  -- subd&planeGraph.faceData .~ faceData'+  --                            &planeGraph.vertexData.traverse %~ getP+  -- where+  --   faceData' = fmap (\(fi, FaceData hs _) -> FaceData hs (getFData fi)) . faces $ subd +  --   -- given a faceId lookup the+  --   getFData fi = let v = boundaryVertices fi subd V.! 0+  --                 in subd^.dataOf v.to holeData ---------------------------------------------------------------------------------+  --   -- note that we intentionally reverse the order of iDd and oD in the call below,+  --   -- as our holes are now outside+  --   subd = fromPolygon px (MultiPolygon (CSeq.fromList [a,b,c,d]) holes') (Just oD) Nothing --- | Construct a planar subdivision from a polygon------ running time: \(O(n)\).-fromPolygon                            :: proxy s-                                       -> SimplePolygon p r-                                       -> f -- ^ data inside-                                       -> f -- ^ data outside the polygon-                                       -> PlanarSubdivision s p () f r-fromPolygon p (SimplePolygon vs) iD oD = PlanarSubdivision g'-  where-    g      = fromVertices p vs-    fData' = V.fromList [FaceData [] iD, FaceData [] oD]+  --   -- for every polygon, construct a hole.+  --   holes' = map withF . F.toList $ pgs+  --   -- add the facedata to the vertex data+  --   withF (pg :+ f) = bimap (\p -> Hole f p) id pg -    g'     = g & faceData .~ fData'-               & dartData.traverse._2 .~ EdgeData Visible ()--- The following does not really work anymore--- frompolygon p (MultiPolygon vs hs) iD oD = PlanarSubdivision g'---   where---     g      = fromVertices p vs---     hs'    = map (\h -> fromPolygon p h oD iD) hs---            -- note that oD and iD are exchanged---     fData' = V.fromList [FaceData iD hs', FaceData oD []]+  --   -- corners of the slightly enlarged boundingbox+  --   (a,b,c,d) = corners . bimap (const $ Outer oD) id+  --             . grow 1 . boundingBoxList . fmap (^.core) $ pgs +    --TODO: We need to mark the edges of the outer square as invisible. -fromVertices      :: proxy s-                  -> C.CSeq (Point 2 r :+ p)-                  -> PlanarGraph s Primal_ (VertexData r p) () ()-fromVertices _ vs = g&vertexData .~ vData'-  where-    n = length vs-    g = planarGraph [ [ (Dart (Arc i)               Positive, ())-                      , (Dart (Arc $ (i+1) `mod` n) Negative, ())-                      ]-                    | i <- [0..(n-1)]]-    vData' = V.fromList . map (\(p :+ e) -> VertexData p e) . F.toList $ vs+    -- Main Idea: Assign the vertices the hole-number on which they occur. For+    -- each face we then find an incident vertex to find the data corresponding+    -- to that face. +data HoleData f p = Outer !f | Hole !f !p deriving (Show,Eq) --- | Constructs a connected planar subdivision.------ pre: the segments form a single connected component--- running time: \(O(n\log n)\)-fromConnectedSegments       :: (Foldable f, Ord r, Num r)-                            => proxy s-                            -> f (LineSegment 2 p r :+ EdgeData e)-                            -> PlanarSubdivision s [p] e () r-fromConnectedSegments px ss = PlanarSubdivision $-    fromConnectedSegments' px ss & faceData.traverse %~ FaceData []+holeData            :: HoleData f p -> f+holeData (Outer f)  = f+holeData (Hole f _) = f --- | Constructs a planar graph------ pre: The segments form a single connected component------ running time: \(O(n\log n)\)-fromConnectedSegments'      :: (Foldable f, Ord r, Num r)-                            => proxy s-                            -> f (LineSegment 2 p r :+ e)-                            -> PlanarGraph s Primal_ (VertexData r [p]) e ()-fromConnectedSegments' _ ss = planarGraph dts & vertexData .~ vxData-  where-    pts         = M.fromListWith (<>) . concatMap f . zipWith g [0..] . F.toList $ ss-    f (s :+ e)  = [ ( s^.start.core-                    , SP [s^.start.extra] [(s^.end.core)   :+ h Positive e])-                  , ( s^.end.core-                    , SP [s^.end.extra]   [(s^.start.core) :+ h Negative e])-                  ]-    g i (s :+ e) = s :+ (Arc i :+ e)-    h d (a :+ e) = (Dart a d, e)+getP            :: HoleData f p -> Maybe p+getP (Outer _)  = Nothing+getP (Hole _ p) = Just p -    vts    = map (\(p,sp) -> (p,map (^.extra) . sortArround (ext p) <$> sp))-           . M.assocs $ pts-    -- vertex Data-    vxData = V.fromList . map (\(p,sp) -> VertexData p (sp^._1)) $ vts-    -- The darts-    dts    = map (^._2._2) vts+-- | grows the box by x on all sides+grow     :: (Num r, Arity d) => r -> Box d p r -> Box d p r+grow x b = let mi = minPoint b+               ma = maxPoint b+               v  = Vec.replicate x+           in box (mi&core %~ (.+^ ((-1) *^ v))) (ma&core %~ (.+^ v))

@@ -0,0 +1,598 @@+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE PartialTypeSignatures #-}+{-# LANGUAGE ScopedTypeVariables #-}+module Data.Geometry.PlanarSubdivision.Basic( VertexId', FaceId'+                                           , VertexData(VertexData), PG.vData, PG.location++                                           , FaceData(FaceData), holes, fData++                                           , PlanarSubdivision(PlanarSubdivision)+                                           , Wrap++                                           , Component, ComponentId++                                           , PolygonFaceData(..)+                                           , PlanarGraph+                                           , PlaneGraph+                                           , fromSimplePolygon+                                           , fromConnectedSegments+                                           , fromPlaneGraph, fromPlaneGraph'++                                           , numVertices, numEdges, numFaces, numDarts+                                           , dual++                                           , components, component+                                           , vertices', vertices+                                           , edges', edges+                                           , faces', faces, internalFaces+                                           , darts'++                                           , headOf, tailOf, twin, endPoints++                                           , incidentEdges, incomingEdges, outgoingEdges+                                           , nextIncidentEdge+                                           , neighboursOf++                                           , leftFace, rightFace+                                           , outerBoundaryDarts, boundaryVertices, holesOf+                                           , outerFaceId+                                           , boundary'++                                           , locationOf+                                           , HasDataOf(..)++                                           , endPointsOf, endPointData++                                           , edgeSegment, edgeSegments+                                           , rawFacePolygon, rawFaceBoundary+                                           , rawFacePolygons++                                           , VertexId(..), FaceId(..), Dart, World(..)+++                                           , rawVertexData, rawDartData, rawFaceData+                                           , dataVal++                                           , dartMapping, Raw(..)+                                           ) where++import           Control.Lens hiding (holes, holesOf, (.=))+import           Data.Aeson+import           Data.Coerce+import           Data.Ext+import qualified Data.Foldable as F+import           Data.Geometry.Box+import           Data.Geometry.LineSegment+import           Data.Geometry.Point+import           Data.Geometry.Polygon+import           Data.Geometry.Properties+import qualified Data.List as L+import           Data.List.NonEmpty (NonEmpty(..))+import qualified Data.List.NonEmpty as NonEmpty+import           Data.Permutation (ix')+import           Data.PlanarGraph (toAdjacencyLists,buildFromJSON, isPositive, allDarts)+import qualified Data.PlaneGraph as PG+import           Data.PlaneGraph( PlaneGraph, PlanarGraph, dual+                                , Dart, VertexId(..), FaceId(..), twin+                                , World(..)+                                , VertexId', FaceId'+                                , VertexData, location, vData+                                , HasDataOf(..)+                                )+import qualified Data.Sequence as Seq+import qualified Data.Vector as V+import qualified Data.Vector.Mutable as MV+import           GHC.Generics (Generic)++--------------------------------------------------------------------------------++-- | The Face data consists of the data itself and a list of holes+data FaceData h f = FaceData { _holes :: (Seq.Seq h)+                             , _fData :: !f+                             } deriving (Show,Eq,Ord,Functor,Foldable,Traversable,Generic)+makeLenses ''FaceData++instance Bifunctor FaceData where+  bimap f g (FaceData hs x) = FaceData (fmap f hs) (g x)+++instance (FromJSON h, FromJSON f) => FromJSON (FaceData h f)+instance (ToJSON h, ToJSON f)     => ToJSON (FaceData h f) where+  toEncoding = genericToEncoding defaultOptions+++--------------------------------------------------------------------------------+++data Wrap' s+type family Wrap (s :: k) :: k where+  Wrap s = Wrap' s++newtype ComponentId s = ComponentId { unCI :: Int }+  deriving (Show,Eq,Ord,Generic,Bounded,Enum)+++data Raw s ia a = Raw { _compId  :: {-# UNPACK #-} !(ComponentId s)+                      , _idxVal  :: {-# UNPACK #-} !ia+                      , _dataVal :: !a+                      } deriving (Eq,Show,Functor,Foldable,Traversable)+makeLenses ''Raw++++--------------------------------------------------------------------------------++-- | A connected component.+--+-- For every face f, and every hole in this face, the facedata points to a dart+-- d on the hole s.t. this dart has the face f on its left. i.e.+-- leftFace d = f+type Component s r = PlaneGraph (Wrap s)+                                (VertexId' s) (Dart s) (FaceData (Dart s) (FaceId' s))+                                r+++-- | A planarsubdivision is essentially a bunch of plane-graphs; one for every+-- connected component. These graphs store the global ID's (darts, vertexId's, faceId's)+-- in their data values. This essentially gives us a mapping between the two.+--+-- note that a face may actually occur in multiple graphs, hence when we store+-- the edges to the the holes, we store the global edgeId's rather than the+-- 'local' edgeId (dart)'s.+--+-- invariant: the outerface has faceId 0+data PlanarSubdivision s v e f r =+  PlanarSubdivision { _components    :: V.Vector (Component s r)+                    , _rawVertexData :: V.Vector (Raw s (VertexId' (Wrap s)) v)+                    , _rawDartData   :: V.Vector (Raw s (Dart      (Wrap s)) e)+                    , _rawFaceData   :: V.Vector (Raw s (FaceId'   (Wrap s)) f)+                    } deriving (Show,Eq,Functor)+makeLenses ''PlanarSubdivision+++type instance NumType   (PlanarSubdivision s v e f r) = r+type instance Dimension (PlanarSubdivision s v e f r) = 2++instance IsBoxable (PlanarSubdivision s v e f r) where+  boundingBox = boundingBoxList' . V.toList . _components+++component    :: ComponentId s -> Lens' (PlanarSubdivision s v e f r)+                                       (Component s r)+component ci = components.ix' (unCI ci)++--------------------------------------------------------------------------------++-- | Constructs a planarsubdivision from a PlaneGraph+--+-- runningTime: \(O(n)\)+fromPlaneGraph   :: forall s v e f r. (Ord r, Fractional r)+                      => PlaneGraph s v e f r -> PlanarSubdivision s v e f r+fromPlaneGraph g = fromPlaneGraph' g (PG.outerFaceDart g)++-- | Given a (connected) PlaneGraph and a dart that has the outerface on its left+-- | Constructs a planarsubdivision+--+-- runningTime: \(O(n)\)+fromPlaneGraph'        :: forall s v e f r. PlaneGraph s v e f r -> Dart s+                       -> PlanarSubdivision s v e f r+fromPlaneGraph' g ofD = PlanarSubdivision (V.singleton . coerce $ g') vd ed fd+  where+    c = ComponentId 0+    vd = V.imap    (\i v   -> Raw c (VertexId i) v)             $ g^.PG.vertexData+    ed = V.zipWith (\d dd  -> Raw c d            dd) allDarts'' $ g^.PG.rawDartData+    fd = V.imap (\i f      -> Raw c (mkFaceId i) f) . swapOf    $ g^.PG.faceData++    g' :: PlaneGraph s (VertexId' s) (Dart s) (FaceData (Dart s) (FaceId' s)) r+    g' = g&PG.faceData    %~ V.imap (\i _ -> mkFaceData i)+          &PG.vertexData  %~ V.imap (\i _ -> VertexId i)+          &PG.rawDartData .~ allDarts''++    allDarts'' :: forall s'. V.Vector (Dart s')+    allDarts'' = allDarts' (PG.numDarts g)++    -- make sure the outerFaceId is 0+    (FaceId (VertexId of')) = PG.leftFace ofD g++    -- at index i we are storing the outerface+    mkFaceData i | i == of'  = faceData (Seq.singleton ofD) 0+                 | i == 0    = faceData mempty of'+                 | otherwise = faceData mempty i+    faceData xs i = FaceData xs (FaceId . VertexId $ i)+++    mkFaceId :: forall s'. Int -> FaceId' s'+    mkFaceId = FaceId . VertexId . mkFaceId'+    mkFaceId' i | i == 0    = of'+                | i == of'  = 0+                | otherwise = i+    swapOf = V.modify (\v -> MV.swap v 0 of')++++++-- | Construct a planar subdivision from a simple polygon+--+-- running time: \(O(n)\).+fromSimplePolygon            :: (Ord r, Fractional r)+                             => proxy s+                             -> SimplePolygon p r+                             -> f -- ^ data inside+                             -> f -- ^ data outside the polygon+                             -> PlanarSubdivision s p () f r+fromSimplePolygon p pg iD oD =+  fromPlaneGraph (PG.fromSimplePolygon p pg iD oD)++-- | Constructs a connected planar subdivision.+--+-- pre: the segments form a single connected component+-- running time: \(O(n\log n)\)+fromConnectedSegments    :: (Foldable f, Ord r, Fractional r)+                         => proxy s+                         -> f (LineSegment 2 p r :+ e)+                         -> PlanarSubdivision s (NonEmpty p) e () r+fromConnectedSegments px = fromPlaneGraph . PG.fromConnectedSegments px++--------------------------------------------------------------------------------++-- | Data type that expresses whether or not we are inside or outside the+-- polygon.+data PolygonFaceData = Inside | Outside deriving (Show,Read,Eq)+++--------------------------------------------------------------------------------+-- * Basic Graph information++-- | Get the number of vertices+--+-- >>> numVertices myGraph+-- 4+numVertices :: PlanarSubdivision s v e f r  -> Int+numVertices = V.length . _rawVertexData++-- | Get the number of Darts+--+-- >>> numDarts myGraph+-- 12+numDarts :: PlanarSubdivision s v e f r  -> Int+numDarts = V.length . _rawDartData++-- | Get the number of Edges+--+-- >>> numEdges myGraph+-- 6+numEdges :: PlanarSubdivision s v e f r  -> Int+numEdges = (`div` 2) . V.length . _rawDartData++-- | Get the number of faces+--+-- >>> numFaces myGraph+-- 4+numFaces :: PlanarSubdivision s v e f r  -> Int+numFaces = V.length . _rawFaceData++-- | Enumerate all vertices+--+-- >>> vertices' myGraph+-- [VertexId 0,VertexId 1,VertexId 2,VertexId 3]+vertices'   :: PlanarSubdivision s v e f r  -> V.Vector (VertexId' s)+vertices' ps = let n = numVertices ps+               in V.fromList $ map VertexId [0..n-1]++-- | Enumerate all vertices, together with their vertex data++-- >>> vertices myGraph+-- [(VertexId 0,()),(VertexId 1,()),(VertexId 2,()),(VertexId 3,())]+vertices    :: PlanarSubdivision s v e f r  -> V.Vector (VertexId' s, VertexData r v)+vertices ps = (\vi -> (vi,ps^.vertexDataOf vi)) <$> vertices' ps++-- | Enumerate all darts+darts' :: PlanarSubdivision s v e f r  -> V.Vector (Dart s)+darts' = allDarts' . numDarts++allDarts'   :: forall s'. Int -> V.Vector (Dart s')+allDarts' n = V.fromList $ take n allDarts+++-- | Enumerate all edges. We report only the Positive darts+edges' :: PlanarSubdivision s v e f r  -> V.Vector (Dart s)+edges' = V.filter isPositive . darts'++-- | Enumerate all edges with their edge data. We report only the Positive+-- darts.+--+-- >>> mapM_ print $ edges myGraph+-- (Dart (Arc 2) +1,"c+")+-- (Dart (Arc 1) +1,"b+")+-- (Dart (Arc 0) +1,"a+")+-- (Dart (Arc 5) +1,"g+")+-- (Dart (Arc 4) +1,"e+")+-- (Dart (Arc 3) +1,"d+")+edges    :: PlanarSubdivision s v e f r  -> V.Vector (Dart s, e)+edges ps = (\e -> (e,ps^.dataOf e)) <$> edges' ps+++faces'    :: PlanarSubdivision s v e f r -> V.Vector (FaceId' s)+faces' ps = let n = numFaces ps+            in V.fromList $ map (FaceId . VertexId) [0..n-1]++faces    :: PlanarSubdivision s v e f r -> V.Vector (FaceId' s, FaceData (Dart s) f)+faces ps = (\fi -> (fi,ps^.faceDataOf fi)) <$> faces' ps++-- | Enumerates all faces with their face data exlcluding  the outer face+internalFaces    :: (Ord r, Fractional r)+                 => PlanarSubdivision s v e f r+                 -> V.Vector (FaceId' s, FaceData (Dart s) f)+internalFaces ps = let i = outerFaceId ps+                   in V.filter (\(j,_) -> i /= j) $ faces ps+++-- -- | lens to access the Dart Data+-- dartData :: Lens (PlanarSubdivision s v e f r) (PlanarSubdivision s v e' f r)+--                  (V.Vector (Dart s, e)) (V.Vector (Dart s, e'))+-- dartData = graph.PG.dartData++++++-- | The tail of a dart, i.e. the vertex this dart is leaving from+--+-- running time: \(O(1)\)+tailOf      :: Dart s -> PlanarSubdivision s v e f r  -> VertexId' s+tailOf d ps = let (_,d',g) = asLocalD d ps+              in g^.dataOf (PG.tailOf d' g)+++-- | The vertex this dart is heading in to+--+-- running time: \(O(1)\)+headOf       :: Dart s -> PlanarSubdivision s v e f r  -> VertexId' s+headOf d ps = let (_,d',g) = asLocalD d ps+              in g^.dataOf (PG.headOf d' g)+++-- | endPoints d g = (tailOf d g, headOf d g)+--+-- running time: \(O(1)\)+endPoints      :: Dart s -> PlanarSubdivision s v e f r+               -> (VertexId' s, VertexId' s)+endPoints d ps = (tailOf d ps, headOf d ps)+++-- | All edges incident to vertex v, in counterclockwise order around v.+--+-- running time: \(O(k)\), where \(k\) is the number of edges reported.+incidentEdges                 :: VertexId' s -> PlanarSubdivision s v e f r+                              -> V.Vector (Dart s)+incidentEdges v ps=  let (_,v',g) = asLocalV v ps+                         ds       = PG.incidentEdges v' g+                     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.+--+-- running time: \(O(1)\)+nextIncidentEdge      :: Dart s -> PlanarSubdivision s v e f r -> Dart s+nextIncidentEdge d ps = let (_,d',g) = asLocalD d ps+                            d''      = PG.nextIncidentEdge d' g+                        in g^.dataOf d''++-- | All incoming edges incident to vertex v, in counterclockwise order around v.+incomingEdges      :: VertexId' s -> PlanarSubdivision s v e f r -> V.Vector (Dart s)+incomingEdges v ps = V.filter (not . isPositive) $ incidentEdges v ps++-- | All outgoing edges incident to vertex v, in counterclockwise order around v.+outgoingEdges      :: VertexId' s -> PlanarSubdivision s v e f r  -> V.Vector (Dart s)+outgoingEdges v ps = V.filter isPositive $ incidentEdges v ps+++-- | Gets the neighbours of a particular vertex, in counterclockwise order+-- around the vertex.+--+-- running time: \(O(k)\), where \(k\) is the output size+neighboursOf      :: VertexId' s -> PlanarSubdivision s v e f r -> V.Vector (VertexId' s)+neighboursOf v ps = otherVtx <$> incidentEdges v ps+  where+    otherVtx d = let u = tailOf d ps in if u == v then headOf d ps else u+++-- | The face to the left of the dart+--+-- running time: \(O(1)\).+leftFace      :: Dart s -> PlanarSubdivision s v e f r  -> FaceId' s+leftFace d ps = let (_,d',g) = asLocalD d ps+                    fi       = PG.leftFace d' g+                in g^.dataOf fi.fData++-- | The face to the right of the dart+--+-- running time: \(O(1)\).+rightFace      :: Dart s -> PlanarSubdivision s v e f r  -> FaceId' s+rightFace d ps = let (_,d',g) = asLocalD d ps+                     fi       = PG.rightFace d' g+                in g^.dataOf fi.fData++-- | The darts on 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.+outerBoundaryDarts      :: FaceId' s -> PlanarSubdivision s v e f r  -> V.Vector (Dart s)+outerBoundaryDarts f ps = let (_,f',g) = asLocalF f ps+                              ds       = PG.boundary f' g+                          in (\d -> g^.dataOf d) <$> ds++++-- | 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.+boundaryVertices      :: FaceId' s -> PlanarSubdivision s v e f r+                      -> V.Vector (VertexId' s)+boundaryVertices f ps = (\d -> headOf d ps) <$> outerBoundaryDarts f ps+++-- | Lists the holes in this face, given as a list of darts to arbitrary darts+-- on those faces.+--+-- running time: \(O(k)\), where \(k\) is the number of darts returned.+holesOf      :: FaceId' s -> PlanarSubdivision s v e f r -> Seq.Seq (Dart s)+holesOf f ps = ps^.faceDataOf f.holes+++--------------------------------------------------------------------------------+-- * Access data++++asLocalD      :: Dart s -> PlanarSubdivision s v e f r+              -> (ComponentId s, Dart (Wrap s), Component s r)+asLocalD d ps = let (Raw ci d' _) = ps^.rawDartData.ix' (fromEnum d)+                in (ci,d',ps^.component ci)+++++asLocalV                 :: VertexId' s -> PlanarSubdivision s v e f r+                         -> (ComponentId s, VertexId' (Wrap s), Component s r)+asLocalV (VertexId v) ps = let (Raw ci v' _) = ps^.rawVertexData.ix' v+                           in (ci,v',ps^.component ci)++asLocalF                          :: FaceId' s -> PlanarSubdivision s v e f r+                                  -> (ComponentId s, FaceId' (Wrap s), Component s r)+asLocalF (FaceId (VertexId f)) ps = let (Raw ci f' _) = ps^.rawFaceData.ix' f+                                    in (ci,f',ps^.component ci)++++-- | Note that using the setting part of this lens may be very expensive!!+vertexDataOf               :: VertexId' s+                           -> Lens' (PlanarSubdivision s v e f r ) (VertexData r v)+vertexDataOf (VertexId vi) = lens get' set''+  where+    get' ps = let (Raw ci wvdi x) = ps^.rawVertexData.ix' vi+                  vd              = ps^.component ci.PG.vertexDataOf wvdi+              in vd&vData .~ x+    set'' ps x = let (Raw ci wvdi _)  = ps^.rawVertexData.ix' vi+                 in ps&rawVertexData.ix' vi.dataVal               .~ (x^.vData)+                      &component ci.PG.vertexDataOf wvdi.location .~ (x^.location)++locationOf   :: VertexId' s -> Lens' (PlanarSubdivision s v e f r ) (Point 2 r)+locationOf v = vertexDataOf v.location+++faceDataOf    :: FaceId' s -> Lens' (PlanarSubdivision s v e f r)+                                    (FaceData (Dart s) f)+faceDataOf fi = lens getF setF+  where+    (FaceId (VertexId i)) = fi+    getF ps = let (Raw ci wfi x) = ps^.rawFaceData.ix' i+                  fd             = ps^.component ci.dataOf wfi+              in fd&fData .~ x++    setF ps fd = let (Raw ci wfi _) = ps^.rawFaceData.ix' i+                     fd'            = fd&fData .~ fi+                     x              = fd^.fData+                 in ps&component ci.dataOf wfi   .~ fd'+                      &rawFaceData.ix' i.dataVal .~ x++instance HasDataOf (PlanarSubdivision s v e f r) (VertexId' s) where+  type DataOf (PlanarSubdivision s v e f r) (VertexId' s) = v+  dataOf v = vertexDataOf v.vData++instance HasDataOf (PlanarSubdivision s v e f r) (Dart s) where+  type DataOf (PlanarSubdivision s v e f r) (Dart s) = e+  dataOf d = rawDartData.ix' (fromEnum d).dataVal++instance HasDataOf (PlanarSubdivision s v e f r) (FaceId' s) where+  type DataOf (PlanarSubdivision s v e f r) (FaceId' s) = f+  dataOf f = faceDataOf f.fData++++-- | Getter for the data at the endpoints of a dart+--+-- running time: \(O(1)\)+endPointsOf   :: Dart s -> Getter (PlanarSubdivision s v e f r )+                                  (VertexData r v, VertexData r v)+endPointsOf d = to (endPointData d)++-- | data corresponding to the endpoints of the dart+--+-- running time: \(O(1)\)+endPointData      :: Dart s -> PlanarSubdivision s v e f r+                  ->  (VertexData r v, VertexData r v)+endPointData d ps = let (u,v) = endPoints d ps+                    in (ps^.vertexDataOf u, ps^.vertexDataOf v)+++--------------------------------------------------------------------------------++-- | gets the id of the outer face+--+-- running time: \(O(1)\)+outerFaceId :: PlanarSubdivision s v e f r -> FaceId' s+outerFaceId = const . FaceId . VertexId $ 0+  -- our invariant tells us the outerface is always at faceId 0++--------------------------------------------------------------------------------++-- | Reports all edges as line segments+edgeSegments    :: PlanarSubdivision s v e f r -> V.Vector (Dart s, LineSegment 2 v r :+ e)+edgeSegments ps = (\d -> (d,edgeSegment d ps)) <$> edges' ps+++-- | Given a dart and the subdivision constructs the line segment representing it+--+-- \(O(1)\)+edgeSegment      :: Dart s -> PlanarSubdivision s v e f r -> LineSegment 2 v r :+ e+edgeSegment d ps = let (p,q) = bimap PG.vtxDataToExt PG.vtxDataToExt $ ps^.endPointsOf d+                   in ClosedLineSegment p q :+ ps^.dataOf d+++-- | Generates the darts incident to a face, starting with the given dart.+--+--+-- \(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+--+-- \(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 = fromPoints 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)+++-- | 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+                        [] -> Left  res                               :+ x+                        hs -> Right (MultiPolygon vs $ map toHole hs) :+ x+  where+    res@(SimplePolygon vs) :+ x = rawFaceBoundary i ps+    toHole d = (rawFaceBoundary (leftFace d ps) ps)^.core++-- | Lists all 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)) . faces' $ ps++++dartMapping ps = ps^.component (ComponentId 0).PG.dartData