hgeometry 0.7.0.0 → 0.8.0.0

raw patch · 122 files changed

+5468/−1873 lines, 122 filesdep +filepathdep +profunctorsdep ~basedep ~yamlbinary-addedPVP ok

version bump matches the API change (PVP)

Dependencies added: filepath, profunctors

Dependency ranges changed: base, yaml

API changes (from Hackage documentation)

Files

• README.md +6/−4
• benchmark/Algorithms/Geometry/ClosestPair/Bench.hs +48/−0
• benchmark/Algorithms/Geometry/ConvexHull/Bench.hs +2/−2
• benchmark/Benchmarks.hs +1/−1
• benchmark/Data/Geometry/Vector/VectorFamily6.hs +0/−1
• changelog.md +10/−0
• docs/Data/PlanarGraph/testG.png binary
• doctests.hs +33/−20
• examples/Demo/Delaunay.hs +2/−2
• examples/Demo/DrawGPX.hs +1/−1
• examples/Demo/GPXParser.hs +1/−2
• examples/Demo/MinDisk.hs +5/−14
• examples/Demo/TriangulateWorld.hs +16/−22
• examples/Demo/WriteEnsemble.hs +1/−4
• hgeometry.cabal +80/−38
• src/Algorithms/Geometry/ClosestPair/DivideAndConquer.hs +129/−0
• src/Algorithms/Geometry/ClosestPair/Naive.hs +48/−0
• src/Algorithms/Geometry/ConvexHull/DivideAndConquer.hs +38/−0
• src/Algorithms/Geometry/ConvexHull/DivideAndConqueror.hs +0/−27
• src/Algorithms/Geometry/ConvexHull/GrahamScan.hs +0/−1
• src/Algorithms/Geometry/DelaunayTriangulation/DivideAndConquer.hs +296/−0
• src/Algorithms/Geometry/DelaunayTriangulation/DivideAndConqueror.hs +0/−296
• src/Algorithms/Geometry/DelaunayTriangulation/Types.hs +8/−5
• src/Algorithms/Geometry/EuclideanMST/EuclideanMST.hs +14/−5
• src/Algorithms/Geometry/LineSegmentIntersection/BentleyOttmann.hs +0/−1
• src/Algorithms/Geometry/LineSegmentIntersection/Naive.hs +0/−1
• src/Algorithms/Geometry/LineSegmentIntersection/Types.hs +1/−1
• src/Algorithms/Geometry/PolyLineSimplification/DouglasPeucker.hs +25/−19
• src/Algorithms/Geometry/PolygonTriangulation/MakeMonotone.hs +3/−4
• src/Algorithms/Geometry/PolygonTriangulation/Triangulate.hs +2/−1
• src/Algorithms/Geometry/PolygonTriangulation/TriangulateMonotone.hs +0/−1
• src/Algorithms/Geometry/PolygonTriangulation/Types.hs +0/−1
• src/Algorithms/Geometry/SmallestEnclosingBall/Naive.hs +17/−6
• src/Algorithms/Geometry/SmallestEnclosingBall/RandomizedIncrementalConstruction.hs +12/−1
• src/Algorithms/Geometry/SmallestEnclosingBall/Types.hs +11/−1
• src/Algorithms/Geometry/Sweep.hs +23/−3
• src/Algorithms/Geometry/WellSeparatedPairDecomposition/Types.hs +12/−3
• src/Algorithms/Geometry/WellSeparatedPairDecomposition/WSPD.hs +26/−23
• src/Algorithms/Util.hs +0/−20
• src/Data/BinaryTree.hs +18/−5
• src/Data/BinaryTree/Zipper.hs +0/−1
• src/Data/CircularSeq.hs +0/−1
• src/Data/Ext.hs +18/−10
• src/Data/Geometry.hs +13/−6
• src/Data/Geometry/Arrangement.hs +25/−0
• src/Data/Geometry/Arrangement/Draw.hs +24/−0
• src/Data/Geometry/Arrangement/Internal.hs +328/−0
• src/Data/Geometry/Ball.hs +10/−6
• src/Data/Geometry/Box.hs +10/−6
• src/Data/Geometry/Box/Internal.hs +17/−3
• src/Data/Geometry/HalfLine.hs +8/−13
• src/Data/Geometry/HyperPlane.hs +94/−0
• src/Data/Geometry/Interval.hs +1/−1
• src/Data/Geometry/Ipe.hs +14/−6
• src/Data/Geometry/Ipe/Attributes.hs +44/−22
• src/Data/Geometry/Ipe/Color.hs +110/−0
• src/Data/Geometry/Ipe/FromIpe.hs +42/−15
• src/Data/Geometry/Ipe/IpeOut.hs +126/−152
• src/Data/Geometry/Ipe/PathParser.hs +0/−1
• src/Data/Geometry/Ipe/Reader.hs +4/−3
• src/Data/Geometry/Ipe/Types.hs +28/−14
• src/Data/Geometry/Ipe/Value.hs +24/−0
• src/Data/Geometry/Ipe/Writer.hs +24/−7
• src/Data/Geometry/KDTree.hs +12/−11
• src/Data/Geometry/Line.hs +22/−12
• src/Data/Geometry/Line/Internal.hs +43/−1
• src/Data/Geometry/LineSegment.hs +28/−29
• src/Data/Geometry/PlanarSubdivision.hs +21/−27
• src/Data/Geometry/PlanarSubdivision/Basic.hs +165/−68
• src/Data/Geometry/PlanarSubdivision/Draw.hs +48/−11
• src/Data/Geometry/Point.hs +25/−7
• src/Data/Geometry/PolyLine.hs +12/−19
• src/Data/Geometry/Polygon.hs +11/−7
• src/Data/Geometry/Polygon/Convex.hs +12/−6
• src/Data/Geometry/Properties.hs +13/−6
• src/Data/Geometry/SegmentTree.hs +0/−14
• src/Data/Geometry/SegmentTree/Generic.hs +30/−20
• src/Data/Geometry/Slab.hs +6/−7
• src/Data/Geometry/SubLine.hs +94/−44
• src/Data/Geometry/Transformation.hs +43/−0
• src/Data/Geometry/Triangle.hs +169/−9
• src/Data/Geometry/Vector.hs +31/−3
• src/Data/Geometry/Vector/VectorFamily.hs +37/−3
• src/Data/Geometry/Vector/VectorFamilyPeano.hs +1/−5
• src/Data/Geometry/Vector/VectorFixed.hs +0/−1
• src/Data/LSeq.hs +310/−0
• src/Data/OrdSeq.hs +0/−1
• src/Data/Permutation.hs +18/−14
• src/Data/PlanarGraph.hs +195/−167
• src/Data/PlaneGraph.hs +87/−15
• src/Data/PlaneGraph/Draw.hs +7/−12
• src/Data/Range.hs +27/−12
• src/Data/Seq.hs +0/−236
• src/Data/Seq2.hs +0/−221
• src/Data/SlowSeq.hs +0/−1
• src/Data/Tree/Util.hs +154/−0
• src/Data/UnBounded.hs +15/−7
• src/Data/Util.hs +59/−3
• src/Data/Yaml/Util.hs +71/−4
• src/Graphics/Camera.hs +165/−0
• src/Test/QuickCheck/HGeometryInstances.hs +28/−10
• test/Algorithms/Geometry/ClosestPair/ClosestPairSpec.hs +45/−0
• test/Algorithms/Geometry/ConvexHull/ConvexHullSpec.hs +3/−4
• test/Algorithms/Geometry/DelaunayTriangulation/DTSpec.hs +3/−3
• test/Algorithms/Geometry/WellSeparatedPairDecomposition/WSPDSpec.hs +3/−3
• test/Data/Geometry/ArrangementSpec.hs +102/−0
• test/Data/Geometry/IntervalSpec.hs +0/−1
• test/Data/Geometry/Ipe/ReaderSpec.hs +1/−2
• test/Data/Geometry/KDTreeSpec.hs +5/−5
• test/Data/Geometry/LineSegmentSpec.hs +13/−9
• test/Data/Geometry/PlanarSubdivisionSpec.hs +76/−3
• test/Data/Geometry/PointSpec.hs +6/−0
• test/Data/Geometry/SubLineSpec.hs +13/−7
• test/Data/Geometry/TriangleSpec.hs +107/−0
• test/Data/Geometry/arrangement.ipe +296/−0
• test/Data/Geometry/arrangement.ipe.out.ipe +274/−0
• test/Data/Geometry/pointInTriangle.ipe +310/−0
• test/Data/PlanarGraphSpec.hs +33/−14
• test/Data/PlaneGraphSpec.hs +72/−0
• test/Data/myGraph.yaml +98/−0
• test/Data/myPlaneGraph.yaml +131/−0
• test/Util.hs +65/−6

@@ -17,7 +17,7 @@ ```haskell newtype Point (d :: Nat) (r :: *) = Point { toVec :: Vector d r } ```-- the vertices of a `PolyLine d p r` are stored in a `Data.Seq2` which enforces+- the vertices of a `PolyLine d p r` are stored in a `Data.LSeq` which enforces that a polyline is a proper polyline, and thus has at least two vertices.  Please note that aspect (2), implementing good algorithms, is much work in@@ -25,11 +25,11 @@ some improvements. Currently, HGeometry provides the following algorithms:  * two \(O(n \log n)\) time algorithms for convex hull in-  $\mathbb{R}^2$: the typical Graham scan, and a divide and conqueror algorithm,+  $\mathbb{R}^2$: the typical Graham scan, and a divide and conquer algorithm, * an \(O(n)\) expected time algorithm for smallest enclosing disk in $\mathbb{R}^$2, * the well-known Douglas Peucker polyline line simplification algorithm, * an \(O(n \log n)\) time algorithm for computing the Delaunay triangulation-(using divide and conqueror).+(using divide and conquer). * an \(O(n \log n)\) time algorithm for computing the Euclidean Minimum Spanning Tree (EMST), based on computing the Delaunay Triangulation. * an \(O(\log^2 n)\) time algorithm to find extremal points and tangents on/to a@@ -38,6 +38,8 @@ polygons. * An \(O(1/\varepsilon^dn\log n)\) time algorithm for constructing a Well-Separated pair   decomposition.+* The classic (optimal) \(O(n\log n)\) time divide and conquer algorithm to+  compute the closest pair among a set of \(n\) points in \(\mathbb{R}^2\).  It also has some geometric data structures. In particular, HGeometry contans an implementation of@@ -108,7 +110,7 @@            let pts  = syms&traverse.core %~ (^.symbolPoint)                pts' = NonEmpty.fromList pts                dt   = delaunayTriangulation $ pts'-               out  = [asIpe drawTriangulation dt]+               out  = [iO $ drawTriangulation dt]            writeIpeFile outFile . singlePageFromContent $ out ``` 

@@ -0,0 +1,48 @@+module Algorithms.Geometry.ClosestPair.Bench where++import qualified Algorithms.Geometry.ClosestPair.DivideAndConquer as DivideAndConquer+import qualified Algorithms.Geometry.ClosestPair.Naive as Naive+import           Benchmark.Util+import           Control.DeepSeq+import           Criterion.Main+import           Criterion.Types+import           Data.Ext+import           Data.Geometry.Point+import           Data.Proxy+import           Test.QuickCheck+import           Test.QuickCheck.HGeometryInstances ()+import           Data.LSeq (LSeq)+import qualified Data.LSeq as LSeq++--------------------------------------------------------------------------------++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 (LSeq 2 (Point 2 r :+ ()))+genPts _ n | n >= 2    = generate (LSeq.promise . LSeq.fromList <$> vectorOf n arbitrary)+           | otherwise = error "genPts: Need at least 2 points"++-- | Benchmark computing the closest pair+benchBuild    :: (Ord r, Num r, NFData r) => LSeq 2 (Point 2 r :+ ()) -> Benchmark+benchBuild ps = bgroup "closestPair" [ bgroup (show n) (build $ take' n ps)+                                     | n <- sizes' ps+                                     ]+  where+    take' n = LSeq.promise . LSeq.take n+    sizes' pts = let n = length pts in [ n*i `div` 100 | i <- [10,20,25,50,75,100]]++    build pts = [ bench "sort"     $ nf LSeq.unstableSort pts+                , bench "Div&Conq" $ nf DivideAndConquer.closestPair pts+                , bench "Naive"    $ nf Naive.closestPair pts+                ]

@@ -1,6 +1,6 @@ module Algorithms.Geometry.ConvexHull.Bench where -import qualified Algorithms.Geometry.ConvexHull.DivideAndConqueror as DivideAndConqueror+import qualified Algorithms.Geometry.ConvexHull.DivideAndConquer as DivideAndConquer import qualified Algorithms.Geometry.ConvexHull.GrahamScan as GrahamScan  -- | copies of the convex hull algo with different point types@@ -62,7 +62,7 @@                 , bench "grahamScan_Family"    $ nf GFam.convexHull       ptsFam                 , bench "grahamScan_Fixed"     $ nf GFix.convexHull       ptsFix -                , bench "Div&Conq"             $ nf DivideAndConqueror.convexHull pts+                , bench "Div&Conq"             $ nf DivideAndConquer.convexHull pts                 ]       where         ptsV2       = fmap (GV.fromP) pts

@@ -1,6 +1,6 @@ module Main where -import qualified Algorithms.Geometry.ConvexHull.Bench as M+import qualified Algorithms.Geometry.ClosestPair.Bench as M  main :: IO () main = M.main

@@ -10,7 +10,6 @@ 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(..))

@@ -1,3 +1,13 @@+#### 0.8 ###++- Compatibility with GHC 8.6+- Added \(O(n\log n)\) time closest pair algorithm.+- Added arrangement data type+- Various Bugfixes+- Added Camera data type with some world to screen transformations.+- Additional read/show instances+- Updated some of the show instances for Ipe related types.+ #### 0.7 ###  - Compatibility with GHC 8.0-8.4

binary file changed (absent → 38819 bytes)

@@ -33,27 +33,40 @@           , "FlexibleContexts"           ] -files = mconcat [ geomModules-                , dataModules-                ]+files = map toFile modules -prefixWith s = map (\s' -> "src/" <> s <> s')+toFile = (\s -> "src/" <> s <> ".hs") . replace '.' '/' +replace     :: Eq a => a -> a -> [a] -> [a]+replace a b = go+  where+    go []                 = []+    go (c:cs) | c == a    = b:go cs+              | otherwise = c:go cs -dataModules = prefixWith "Data/" [ "Range.hs"-                                 , "CircularList/Util.hs"-                                 , "Permutation.hs"-                                 , "CircularSeq.hs"-                                 , "PlanarGraph.hs"-                                 ]+modules =+  [ "Data.Range"+  , "Data.CircularList.Util"+  , "Data.Permutation"+  , "Data.CircularSeq"+  , "Data.LSeq"+  , "Data.PlanarGraph"+  , "Data.Tree.Util" -geomModules = prefixWith "Data/Geometry/" [ "Point.hs"-                                          , "Vector.hs"-                                          , "Line.hs"-                                          , "Line/Internal.hs"                                                                        , "Interval.hs"-                                          , "LineSegment.hs"-                                          , "PolyLine.hs"-                                          , "Polygon.hs"-                                          , "Ball.hs"-                                          , "Box.hs"-                                          ]+  , "Data.Geometry.Point"+  , "Data.Geometry.Vector"+  , "Data.Geometry.Transformation"+  , "Data.Geometry.Line"+  , "Data.Geometry.Line.Internal"+  , "Data.Geometry.Interval"+  , "Data.Geometry.LineSegment"+  , "Data.Geometry.PolyLine"+  , "Data.Geometry.Polygon"+  , "Data.Geometry.Ball"+  , "Data.Geometry.Box"++  , "Data.Geometry.Ipe.IpeOut"+  , "Data.Geometry.Ipe.FromIpe"++  -- , "Algorithms.Geometry.HiddenSurfaceRemoval.HiddenSurfaceRemoval"+  ]

@@ -1,7 +1,7 @@ {-# LANGUAGE ScopedTypeVariables #-} module Demo.Delaunay where -import           Algorithms.Geometry.DelaunayTriangulation.DivideAndConqueror+import           Algorithms.Geometry.DelaunayTriangulation.DivideAndConquer import           Algorithms.Geometry.DelaunayTriangulation.Types import           Algorithms.Geometry.EuclideanMST.EuclideanMST import           Control.Lens@@ -45,7 +45,7 @@                pts' = NonEmpty.fromList pts                dt   = delaunayTriangulation $ pts'                emst = euclideanMST pts'-               out  = [asIpe drawTriangulation dt, asIpe drawTree' emst]+               out  = [iO $ drawTriangulation dt, iO $ drawTree' emst]            -- print $ length $ edges' dt            -- print $ toPlaneGraph (Proxy :: Proxy DT) dt            writeIpeFile outFile . singlePageFromContent $ out

@@ -61,7 +61,7 @@     tks   <- concatMap (_tracks . combineTracks) <$> mapM readGPXFile files     let polies  = mapMaybe asPolyLine tks         polies' = map (douglasPeucker 0.01 . scaleUniformlyBy 100) polies-        pg = singlePageFromContent $ map (asIpeObject' mempty) polies'+        pg = singlePageFromContent $ map (iO . defIO) polies'     -- print pg     writeIpeFile outPath pg 

@@ -6,9 +6,8 @@ module Demo.GPXParser where  -import           Control.Applicative+ import           Control.Lens-import           Control.Monad import qualified Data.ByteString.Lazy as B import           Data.Ext import           Data.Geometry.Point

@@ -3,7 +3,6 @@  import           Algorithms.Geometry.SmallestEnclosingBall.RandomizedIncrementalConstruction import           Algorithms.Geometry.SmallestEnclosingBall.Types-import           Control.Applicative import           Control.Lens import           Data.Data import           Data.Ext@@ -11,16 +10,10 @@ import           Data.Geometry import           Data.Geometry.Ball import           Data.Geometry.Ipe-import           Data.Geometry.Ipe.Types import           Data.Geometry.Line-import           Data.Geometry.PolyLine-import           Data.Maybe-import           Data.Semigroup-import           Data.Seq2+import qualified Data.LSeq as LSeq import qualified Data.Traversable as Tr-import           Data.Vinyl import           Options.Applicative-import           System.Environment (getArgs) import           System.Random  --------------------------------------------------------------------------------@@ -39,12 +32,10 @@  -------------------------------------------------------------------------------- -diskResult :: Floating r => IpeOut (DiskResult p r) (IpeObject r)-diskResult = IpeOut f+diskResult :: Floating r => IpeOut (DiskResult p r) Group r+diskResult (DiskResult d pts) = ipeGroup (iO' d  : (F.toList . fmap g $ pts))   where-    f (DiskResult d pts) = asIpeGroup (asIpeObject d mempty : (F.toList . fmap g $ pts))-    g p = asIpeObject (p^.core) mempty-+    g p = iO' (p^.core)  mainWith              :: Options -> IO () mainWith (Options fp) = do@@ -56,7 +47,7 @@         case map ext $ ipeP^..content.Tr.traverse._IpeUse.core.symbolPoint of           pts@(_:_:_) -> do                            let res = smallestEnclosingDisk gen pts-                           printAsIpeSelection . asIpe diskResult $ res+                           printAsIpeSelection . iO . diskResult $ res           _           -> putStrLn "Not enough points!"  

@@ -37,20 +37,6 @@                          <> 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 ()@@ -68,14 +54,10 @@                      . 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 ]+          out     = mconcat [ [ iO' pg | pg <- polies ]+                            , [ iO' s  | s  <- segs ]+                            , [ iO' pg | pg <- yMonotones ]+                            ]       mapM_ print . map (\pg -> pg^.core.to polygonVertices.to length) $ polies'       writeIpeFile outFile . singlePageFromContent $ out @@ -104,3 +86,15 @@ --     where --       applyAts ats = id -- flattenGroups' o                            = [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

@@ -3,17 +3,14 @@ import Control.Lens import Data.Data import Data.Ext-import Control.Applicative import Data.Fixed import qualified Data.Text as T import qualified Data.Text.IO as TIO import Data.Geometry.Ipe import Data.Geometry import Data.Geometry.PolyLine(fromPoints)-import System.Environment import System.Directory import Data.List(isSuffixOf)-import Data.Monoid import Data.Time.Calendar import Options.Applicative @@ -54,7 +51,7 @@           _        -> asTempPt     polies <- mapM (fmap (asPts f) . readFile' . ((inPath ++ "/") ++)) inFiles     let polies' = map (fromPoints . take 100) . trim $ polies-    writeIpeFile outPath . singlePageFromContent . map (flip asIpeObject mempty) $ polies'+    writeIpeFile outPath . singlePageFromContent . map iO' $ polies'  readFile'    :: String -> IO T.Text readFile' fp = putStrLn fp >> TIO.readFile fp

@@ -2,7 +2,7 @@ -- documentation, see http://haskell.org/cabal/users-guide/  name:                hgeometry-version:             0.7.0.0+version:             0.8.0.0 synopsis:            Geometric Algorithms, Data structures, and Data types. description:   HGeometry provides some basic geometry types, and geometric algorithms and@@ -30,7 +30,13 @@                      test/Algorithms/Geometry/PolygonTriangulation/simplepolygon6.ipe                      test/Algorithms/Geometry/SmallestEnclosingDisk/manual.ipe                      test/Data/Geometry/pointInPolygon.ipe+                     test/Data/Geometry/pointInTriangle.ipe                      test/Data/Geometry/Polygon/Convex/convexTests.ipe+                     test/Data/Geometry/arrangement.ipe+                     test/Data/Geometry/arrangement.ipe.out.ipe+                     test/Data/myGraph.yaml+                     test/Data/myPlaneGraph.yaml+                      examples/BAPC2014/sample.in                      examples/BAPC2014/sample.out                      examples/BAPC2014/testdata.in@@ -40,11 +46,16 @@                      examples/BAPC2012/sampleG.in                      examples/BAPC2012/sampleG.out -+                     -- in the future (cabal >=2.4) we can use+                     -- examples/**/*.in+                     -- examples/**/*.out  extra-source-files:  README.md                      changelog.md +Extra-doc-files:     docs/Data/PlanarGraph/testG.png+                     -- docs/**/*.png+ cabal-version:       2.0 source-repository head   type:     git@@ -93,6 +104,7 @@                     Data.Geometry.SubLine                     Data.Geometry.HalfLine                     Data.Geometry.PolyLine+                    Data.Geometry.HyperPlane                     Data.Geometry.Triangle                     -- Data.Geometry.Plane                     Data.Geometry.Slab@@ -108,16 +120,20 @@                     Data.Geometry.SegmentTree.Generic                      Data.Geometry.KDTree+                     Data.Geometry.PlanarSubdivision                     Data.Geometry.PlanarSubdivision.Basic                     Data.Geometry.PlanarSubdivision.Draw +                    Data.Geometry.Arrangement+                    Data.Geometry.Arrangement.Internal+                    Data.Geometry.Arrangement.Draw+                     -- * Algorithms-                    Algorithms.Util                      -- * Geometric Algorithms                     Algorithms.Geometry.ConvexHull.GrahamScan-                    Algorithms.Geometry.ConvexHull.DivideAndConqueror+                    Algorithms.Geometry.ConvexHull.DivideAndConquer                      Algorithms.Geometry.LowerEnvelope.DualCH @@ -126,7 +142,7 @@                     Algorithms.Geometry.SmallestEnclosingBall.Naive                      Algorithms.Geometry.DelaunayTriangulation.Types-                    Algorithms.Geometry.DelaunayTriangulation.DivideAndConqueror+                    Algorithms.Geometry.DelaunayTriangulation.DivideAndConquer                     Algorithms.Geometry.DelaunayTriangulation.Naive                      Algorithms.Geometry.PolyLineSimplification.DouglasPeucker@@ -150,6 +166,12 @@                     Algorithms.Geometry.LineSegmentIntersection.BentleyOttmann                     Algorithms.Geometry.LineSegmentIntersection.Types +                    -- Algorithms.Geometry.HiddenSurfaceRemoval.HiddenSurfaceRemoval++                    Algorithms.Geometry.ClosestPair.Naive+                    Algorithms.Geometry.ClosestPair.DivideAndConquer++                     -- * Graph Algorithms                     Algorithms.Graph.DFS                     Algorithms.Graph.MST@@ -157,6 +179,8 @@                     -- * Ipe Types                     Data.Geometry.Ipe                     Data.Geometry.Ipe.Literal+                    Data.Geometry.Ipe.Value+                    Data.Geometry.Ipe.Color                     Data.Geometry.Ipe.Attributes                     Data.Geometry.Ipe.Types                     Data.Geometry.Ipe.Writer@@ -169,8 +193,7 @@                     Data.UnBounded                     Data.Range                     Data.Ext-                    Data.Seq2-                    Data.Seq+                    Data.LSeq                     Data.CircularSeq                     Data.Sequence.Util                     Data.BinaryTree@@ -180,6 +203,7 @@                     Data.BalBST                     Data.OrdSeq                     Data.SlowSeq+                    Data.Tree.Util                     Data.Util                      -- * Planar Graphs@@ -188,6 +212,9 @@                     Data.PlaneGraph                     Data.PlaneGraph.Draw +                    -- * Graphics stuff+                    Graphics.Camera+                     -- * Other                     System.Random.Shuffle                     Control.Monad.State.Persistent@@ -203,42 +230,44 @@    -- other-extensions:   build-depends:-                  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-                , semigroupoids    >= 5-                , semigroups       >= 0.18-                , singletons       >= 2.0-                , text             >= 1.1.1.0-                , vinyl            >= 0.6-                , deepseq          >= 1.1-                , fingertree       >= 0.1-                , colour           >= 2.3.3-                , reflection       >= 2.1+                base             >= 4.10      &&     < 5+              , bifunctors       >= 4.1+              , bytestring       >= 0.10+              , containers       >= 0.5.5+              , dlist            >= 0.7+              , contravariant    >= 1.4+              , profunctors      >= 5.2.1+              , lens             >= 4.2+              , linear           >= 1.10+              , semigroupoids    >= 5+              , semigroups       >= 0.18+              , singletons       >= 2.0+              , text             >= 1.1.1.0+              , vinyl            >= 0.6        && < 0.9+              , deepseq          >= 1.1+              , fingertree       >= 0.1+              , colour           >= 2.3.3+              , reflection       >= 2.1 -               -- , validation       >= 0.4 -               , parsec           >= 3-               -- , tranformers      > 0.3+              -- , validation       >= 0.4 -               , vector           >= 0.11-               , fixed-vector     >= 1.0-               , data-clist       >= 0.0.7.2+              , parsec           >= 3+                -- , tranformers      > 0.3 -               , hexpat           >= 0.20.9-               , aeson            >= 1.0-               , yaml             >= 0.8+              , vector           >= 0.11+              , fixed-vector     >= 1.0+              , data-clist       >= 0.0.7.2 +              , hexpat           >= 0.20.9+              , aeson            >= 1.0+              , yaml             >= 0.8 -               , mtl-               , random-               , template-haskell +              , mtl+              , random+              , template-haskell+   if flag(with-quickcheck)     build-depends: QuickCheck              >= 2.5                  , quickcheck-instances    >= 0.3@@ -416,7 +445,8 @@   ghc-options:   -threaded   main-is:       doctests.hs   build-depends: base-               , doctest >= 0.8+               , doctest             >= 0.8+--               , doctest-discover    default-language:    Haskell2010 @@ -425,7 +455,9 @@   default-language:     Haskell2010   hs-source-dirs:       test   main-is:              Spec.hs-  ghc-options:   -O2 -fno-warn-unticked-promoted-constructors+  ghc-options:   -fno-warn-unticked-promoted-constructors+                 -fno-warn-partial-type-signatures+                 -fno-warn-missing-signatures    build-tool-depends: hspec-discover:hspec-discover @@ -433,6 +465,7 @@   other-modules: Data.RangeSpec                  Data.EdgeOracleSpec                  Data.PlanarGraphSpec+                 Data.PlaneGraphSpec                  Data.OrdSeqSpec                  Data.Geometry.Ipe.ReaderSpec                  Data.Geometry.PolygonSpec@@ -445,6 +478,8 @@                  Data.Geometry.LineSpec                  Data.Geometry.SubLineSpec                  Data.Geometry.PlanarSubdivisionSpec+                 Data.Geometry.TriangleSpec+                 Data.Geometry.ArrangementSpec                   Algorithms.Geometry.SmallestEnclosingDisk.RISpec                  Algorithms.Geometry.DelaunayTriangulation.DTSpec@@ -454,6 +489,8 @@                  Algorithms.Geometry.PolygonTriangulation.TriangulateMonotoneSpec                  Algorithms.Geometry.LowerEnvelope.LowerEnvSpec                  Algorithms.Geometry.ConvexHull.ConvexHullSpec+                 -- Algorithms.Geometry.HiddenSurfaceRemoval.HiddenSurfaceRemovalSpec+                 Algorithms.Geometry.ClosestPair.ClosestPairSpec                   Util @@ -475,7 +512,11 @@                       , random                       , singletons                       , colour+                      , filepath+                      , directory+                      , yaml +   default-extensions: TypeFamilies                     , GADTs                     , KindSignatures@@ -551,6 +592,7 @@                  Demo.ExpectedPairwiseDistance                  Demo.TriangulateWorld                  WSPDBench+                 Algorithms.Geometry.ClosestPair.Bench    build-depends:                 base

@@ -0,0 +1,129 @@+{-# LANGUAGE ScopedTypeVariables #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Algorithms.Geometry.ClosestPair.DivideAndConquer+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Classical \(O(n\log n)\) time divide and conquer algorithm to compute the+-- closest pair among a set of \(n\) points in \(\mathbb{R}^2\).+--+--------------------------------------------------------------------------------+module Algorithms.Geometry.ClosestPair.DivideAndConquer where++import           Control.Lens+import           Data.BinaryTree+import           Data.Ext+import qualified Data.Foldable as F+import           Data.Geometry.Point+import           Data.LSeq (LSeq)+import qualified Data.LSeq as LSeq+import qualified Data.List as List+import           Data.List.NonEmpty (NonEmpty(..))+import qualified Data.List.NonEmpty as NonEmpty+import           Data.Ord (comparing)+import           Data.Semigroup.Foldable(foldMap1)+import           Data.UnBounded+import           Data.Util++--------------------------------------------------------------------------------++-- | Classical divide and conquer algorithm to compute the closest pair among+-- \(n\) points.+--+-- running time: \(O(n)\)+closestPair :: ( Ord r, Num r) => LSeq 2 (Point 2 r :+ p) -> Two (Point 2 r :+ p)+closestPair = f . foldMap1 mkCCP . asBalancedBinLeafTree . LSeq.toNonEmpty+            . LSeq.unstableSortBy (comparing (^.core))+  where+    mkCCP (Elem p) = CCP (p :| []) Top+    f = \case+          CCP _ (ValT (SP cp _)) -> cp+          CCP _ Top              -> error "closestPair: absurd."+++-- | the closest pair and its (squared) distance+type CP q r = Top (SP (Two q) r)++-- | Type used in the closest pair computation. The fields represent the points+-- ordered on increasing y-order and the closest pair (if we know it)+data CCP p r = CCP (NonEmpty (Point 2 r :+ p))+                   !(CP (Point 2 r :+ p) r)+            deriving (Show,Eq)++instance (Num r, Ord r) => Semigroup (CCP p r) where+  (CCP ptsl cpl) <> (CCP ptsr cpr) = CCP (mergeSortedBy cmp ptsl ptsr)+                                         (mergePairs (minBy getDist cpl cpr) ptsl ptsr)+    where+      -- compare on y first then on x+      cmp     :: Point 2 r :+ p -> Point 2 r :+ p -> Ordering+      cmp p q = comparing (^.core.yCoord) p q <> comparing (^.core.xCoord) p q+++-- | Given an ordering and two nonempty sequences ordered according to that+-- ordering, merge them+mergeSortedBy           :: (a -> a -> Ordering) -> NonEmpty a -> NonEmpty a -> NonEmpty a+mergeSortedBy cmp ls rs = NonEmpty.fromList $ go (F.toList ls) (F.toList rs)+  where+    go []         ys     = ys+    go xs         []     = xs+    go xs@(x:xs') ys@(y:ys') = case x `cmp` y of+                                 LT -> x : go xs' ys+                                 EQ -> x : go xs' ys+                                 GT -> y : go xs  ys'++-- | Function that does the actual merging work+mergePairs            :: forall p r. (Ord r, Num r)+                      => CP (Point 2 r :+ p) r -- ^ current closest pair and its dist+                      -> NonEmpty (Point 2 r :+ p) -- ^ pts on the left+                      -> NonEmpty (Point 2 r :+ p) -- ^ pts on the right+                      -> CP (Point 2 r :+ p) r+mergePairs cp' ls' rs' = go cp' (NonEmpty.toList ls') (NonEmpty.toList rs')+  where++    -- scan through the points on the right in increasing order.+    go              :: CP (Point 2 r :+ p) r -> [Point 2 r :+ p] -> [Point 2 r :+ p]+                    -> CP (Point 2 r :+ p) r+    go cp _  []     = cp+    go cp ls (r:rs) = let ls'' = trim (getDist cp) ls r+                          cp'' = run cp r ls'' -- try to find a new closer pair with r.+                      in go cp'' ls'' rs   -- and then process the remaining points++    -- ditch the points on the left that are too low anyway+    trim               :: Top r -> [Point 2 r :+ q] -> Point 2 r :+ a+                       -> [Point 2 r :+ q]+    trim (ValT d) ls r = List.dropWhile (\l -> sqVertDist l r > d) ls+    trim _        ls _ = ls++    -- the squared vertical distance (in case r lies above l) or 0 otherwise+    sqVertDist l r = let d = 0 `max` (r^.core.yCoord - l^.core.yCoord) in d*d++    -- try and find a new closest pair with r. If we get to points that are too far above+    -- r we stop (since none of those points will be closer to r anyway)+    run          :: CP (Point 2 r :+ p) r -> Point 2 r :+ p -> [Point 2 r :+ p]+                 -> CP (Point 2 r :+ p) r+    run cp'' r ls =+      runWhile cp'' ls+               (\cp l -> ValT (l^.core.yCoord - r^.core.yCoord) < (getDist cp))+               (\cp l -> minBy getDist cp (ValT $ SP (Two l r) (dist l r)))++    dist (p :+ _) (q :+ _) = squaredEuclideanDist p q+++-- | Given some function that decides when to keep things while maintaining some state.+runWhile           :: s -> [a] -> (s -> a -> Bool) -> (s -> a -> s) -> s+runWhile s' ys p f = go s' ys+  where+    go s []                 = s+    go s (x:xs) | p s x     = go (f s x) xs  -- continue with new state+                | otherwise = s -- stop, return the current state++-- | returns the minimum element according to some function.+minBy                   :: Ord b => (a -> b) -> a -> a -> a+minBy f a b | f a < f b = a+            | otherwise = b++-- | Get the distance of a (candidate) closest pair+getDist :: CP a r -> Top r+getDist = fmap (view _2)

@@ -0,0 +1,48 @@+--------------------------------------------------------------------------------+-- |+-- Module      :  Algorithms.Geometry.ClosestPair.Naive+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Naive \O(n\^2)\) time algorithm to compute the closest pair of points among+-- \(n\) points in \(\mathbb{R}^d\).+--+--------------------------------------------------------------------------------+module Algorithms.Geometry.ClosestPair.Naive where++import           Data.Ext+import qualified Data.Foldable as F+import           Data.Geometry (qdA)+import           Data.Geometry.Point+import           Data.Geometry.Vector (Arity)+import           Data.LSeq (LSeq)+import qualified Data.List.NonEmpty as NonEmpty+import           Data.Semigroup+import           Data.Util++--------------------------------------------------------------------------------++-- | Naive algorithm to compute the closest pair in \(d\) dimensions. Runs in+-- \(O(n^2)\) time (for any constant \(d\)). Note that we need at least two elements+-- for there to be a closest pair.+closestPair :: ( Ord r, Arity d, Num r+               ) => LSeq 2 (Point d r :+ p) -> Two (Point d r :+ p)+closestPair = getVal . getMin . sconcat . fmap (uncurry' mkPair) . pairs+  where+    uncurry' f (Two a b) = f a b+    getVal (Arg _ x) = x++-- | A pair of points+type PP d p r = ArgMin r (Two (Point d r :+ p))++-- | Create a pair of points+mkPair                         :: (Arity d, Num r)+                               => Point d r :+ p -> Point d r :+ p -> PP d p r+mkPair pp@(p :+ _) qq@(q :+ _) = let dst = qdA p q+                                 in Min (Arg dst (Two pp qq))++-- | Produce all lists from a vec of elements. Since the Vec contains at least two+-- elements, the resulting list is non-empty+pairs :: LSeq 2 a -> NonEmpty.NonEmpty (Two a)+pairs = NonEmpty.fromList . uniquePairs . F.toList

@@ -0,0 +1,38 @@+--------------------------------------------------------------------------------+-- |+-- Module      :  Algorithms.Geometry.ConvexHull.DivideAndConquer+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- \(O(n\log n)\) time divide and conquer algorithm to compute the convex hull+-- of a set of \(n\) points in \(\mathbb{R}^2\).+--+--------------------------------------------------------------------------------+module Algorithms.Geometry.ConvexHull.DivideAndConquer(convexHull) where++import           Control.Lens ((^.))+import           Data.BinaryTree+import           Data.Ext+import           Data.Function (on)+import           Data.Geometry.Point+import           Data.Geometry.Polygon+import           Data.Geometry.Polygon.Convex+import qualified Data.List.NonEmpty as NonEmpty+import           Data.Semigroup.Foldable++--------------------------------------------------------------------------------++-- | \(O(n \log n)\) time ConvexHull using divide and conquer. The resulting polygon is+-- given in clockwise order.+convexHull :: (Ord r, Num r)+           => NonEmpty.NonEmpty (Point 2 r :+ p) -> ConvexPolygon p r+convexHull = unMerge+           . foldMap1 (Merge . ConvexPolygon . fromPoints . (:[]) . _unElem)+           . asBalancedBinLeafTree+           . NonEmpty.sortBy (compare `on` (^.core))++newtype Merge r p = Merge { unMerge :: ConvexPolygon p r }++instance (Num r, Ord r) => Semigroup (Merge r p) where+  (Merge lp) <> (Merge rp) = let (ch,_,_) = merge lp rp in Merge ch

@@ -1,27 +0,0 @@-module Algorithms.Geometry.ConvexHull.DivideAndConqueror( convexHull-                                                        ) where--import           Control.Lens ((^.))-import           Data.BinaryTree-import           Data.Ext-import           Data.Function (on)-import           Data.Geometry.Point-import           Data.Geometry.Polygon-import           Data.Geometry.Polygon.Convex-import qualified Data.List.NonEmpty as NonEmpty-import           Data.Semigroup-import           Data.Semigroup.Foldable---- | \(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-           . foldMap1 (Merge . ConvexPolygon . fromPoints . (:[]) . _unElem)-           . asBalancedBinLeafTree-           . NonEmpty.sortBy (compare `on` (^.core))--newtype Merge r p = Merge { unMerge :: ConvexPolygon p r }--instance (Num r, Ord r) => Semigroup (Merge r p) where-  (Merge lp) <> (Merge rp) = let (ch,_,_) = merge lp rp in Merge ch

@@ -10,7 +10,6 @@ import           Data.Geometry.Polygon.Convex (ConvexPolygon(..)) import qualified Data.List.NonEmpty as NonEmpty import           Data.List.NonEmpty (NonEmpty(..))-import           Data.Monoid   -- | \(O(n \log n)\) time ConvexHull using Graham-Scan. The resulting polygon is

@@ -0,0 +1,296 @@+{-# LANGUAGE ScopedTypeVariables #-}+module Algorithms.Geometry.DelaunayTriangulation.DivideAndConquer where++import           Algorithms.Geometry.ConvexHull.GrahamScan as GS+import           Algorithms.Geometry.DelaunayTriangulation.Types+import           Control.Lens+import           Control.Monad.Reader+import           Control.Monad.State+import           Data.BinaryTree+import qualified Data.CircularList as CL+import qualified Data.CircularSeq as CS+import qualified Data.CircularList.Util as CU+import           Data.Ext+import qualified Data.Foldable as F+import           Data.Function (on)+import           Data.Geometry hiding (rotateTo)+import           Data.Geometry.Ball (disk, insideBall)+import           Data.Geometry.Polygon+import qualified Data.Geometry.Polygon.Convex as Convex+import           Data.Geometry.Polygon.Convex (ConvexPolygon(..), simplePolygon)+import qualified Data.IntMap.Strict as IM+import qualified Data.List as L+import qualified Data.List.NonEmpty as NonEmpty+import qualified Data.Map as M+import           Data.Maybe (fromJust, fromMaybe)+import qualified Data.Vector as V++-------------------------------------------------------------------------------+-- * Divide & Conqueror Delaunay Triangulation+--+-- Implementation of the Divide & Conqueror algorithm as described in:+--+-- Two Algorithms for Constructing a Delaunay Triangulation+-- Lee and Schachter+-- International Journal of Computer and Information Sciences, Vol 9, No. 3, 1980+--+-- We store all adjacency lists in clockwise order+--+-- : If v on the convex hull, then its first entry in the adj. lists is its CCW+-- successor (i.e. its predecessor) on the convex hull+--+-- Rotating Right <-> rotate clockwise+++-- | Computes the delaunay triangulation of a set of points.+--+-- Running time: \(O(n \log n)\)+-- (note: We use an IntMap in the implementation. So maybe actually \(O(n \log^2 n)\))+--+-- pre: the input is a *SET*, i.e. contains no duplicate points. (If the+-- input does contain duplicate points, the implementation throws them away)+delaunayTriangulation      :: (Ord r, Fractional r)+                           => NonEmpty.NonEmpty (Point 2 r :+ p) -> Triangulation p r+delaunayTriangulation pts' = Triangulation vtxMap ptsV adjV+  where+    pts    = nub' . NonEmpty.sortBy (compare `on` (^.core)) $ pts'+    ptsV   = V.fromList . F.toList $ pts+    vtxMap = M.fromList $ zip (map (^.core) . V.toList $ ptsV) [0..]++    tr     = _unElem <$> asBalancedBinLeafTree pts++    (adj,_) = delaunayTriangulation' tr (vtxMap,ptsV)+    adjV    = V.fromList . IM.elems $ adj++++-- : pre: - Input points are sorted lexicographically+delaunayTriangulation' :: (Ord r, Fractional r)+                       => BinLeafTree Size (Point 2 r :+ p)+                       -> Mapping p r+                       -> (Adj, ConvexPolygon (p :+ VertexID) r)+delaunayTriangulation' pts mapping'@(vtxMap,_)+  | size' pts == 1 = let (Leaf p) = pts+                         i        = lookup' vtxMap (p^.core)+                     in (IM.singleton i CL.empty, ConvexPolygon $ fromPoints [withID p i])+  | size' pts <= 3 = let pts'  = NonEmpty.fromList+                               . map (\p -> withID p (lookup' vtxMap (p^.core)))+                               . F.toList $ pts+                         ch    = GS.convexHull pts'+                     in (fromHull mapping' ch, ch)+  | otherwise      = let (Node lt _ rt) = pts+                         (ld,lch)       = delaunayTriangulation' lt mapping'+                         (rd,rch)       = delaunayTriangulation' rt mapping'+                         (ch, bt, ut)   = Convex.merge lch rch+                     in (merge ld rd bt ut mapping' (firsts ch), ch)++--------------------------------------------------------------------------------+-- * Implementation++-- | Mapping that says for each vtx in the convex hull what the first entry in+-- the adj. list should be. The input polygon is given in Clockwise order+firsts :: ConvexPolygon (p :+ VertexID) r -> IM.IntMap VertexID+firsts = IM.fromList . map (\s -> (s^.end.extra.extra, s^.start.extra.extra))+       . F.toList . outerBoundaryEdges . _simplePolygon+++-- | Given a polygon; construct the adjacency list representation+-- pre: at least two elements+fromHull              :: Ord r => Mapping p r -> ConvexPolygon (p :+ q) r -> Adj+fromHull (vtxMap,_) p = let vs@(u:v:vs') = map (lookup' vtxMap . (^.core))+                                         . F.toList . CS.rightElements+                                         $ p^.simplePolygon.outerBoundary+                            es           = zipWith3 f vs (tail vs ++ [u]) (vs' ++ [u,v])+                            f prv c nxt  = (c,CL.fromList . L.nub $ [prv, nxt])+                        in IM.fromList es+++-- | Merge the two delaunay triangulations.+--+-- running time: \(O(n)\) (although we cheat a bit by using a IntMap)+merge                            :: (Ord r, Fractional r)+                                 => Adj+                                 -> Adj+                                 -> LineSegment 2 (p :+ VertexID) r -- ^ lower tangent+                                 -> LineSegment 2 (p :+ VertexID) r -- ^ upper tangent+                                 -> Mapping p r+                                 -> Firsts+                                 -> Adj+merge ld rd bt ut mapping'@(vtxMap,_) fsts =+    flip runReader (mapping', fsts) . flip execStateT adj $ moveUp (tl,tr) l r+  where+    l   = lookup' vtxMap (bt^.start.core)+    r   = lookup' vtxMap (bt^.end.core)+    tl  = lookup' vtxMap (ut^.start.core)+    tr  = lookup' vtxMap (ut^.end.core)+    adj = ld `IM.union` rd++type Merge p r = StateT Adj (Reader (Mapping p r, Firsts))++type Firsts = IM.IntMap VertexID++-- | Merges the two delaunay traingulations.+moveUp          :: (Ord r, Fractional r)+                => (VertexID,VertexID) -> VertexID -> VertexID -> Merge p r ()+moveUp ut l r+  | (l,r) == ut = insert l r+  | otherwise   = do+                     insert l r+                     -- Get the neighbours of r and l along the convex hull+                     r1 <- pred' . rotateTo l . lookup'' r <$> get+                     l1 <- succ' . rotateTo r . lookup'' l <$> get++                     (r1',a) <- rotateR l r r1+                     (l1',b) <- rotateL l r l1+                     c       <- qTest l r r1' l1'+                     let (l',r') = case (a,b,c) of+                                     (True,_,_)          -> (focus' l1', r)+                                     (False,True,_)      -> (l,          focus' r1')+                                     (False,False,True)  -> (l,          focus' r1')+                                     (False,False,False) -> (focus' l1', r)+                     moveUp ut l' r'+++-- | ''rotates'' around r and removes all neighbours of r that violate the+-- delaunay condition. Returns the first vertex (as a Neighbour of r) that+-- should remain in the Delaunay Triangulation, as well as a boolean A that+-- helps deciding if we merge up by rotating left or rotating right (See+-- description in the paper for more info)+rotateR        :: (Ord r, Fractional r)+               => VertexID -> VertexID -> Vertex -> Merge p r (Vertex, Bool)+rotateR l r r1 = focus' r1 `isLeftOf` (l, r) >>= \case+                   True  -> (,False) <$> rotateR' l r r1 (pred' r1)+                   False -> pure (r1,True)++-- | The code that does the actual rotating+rotateR'     :: (Ord r, Fractional r)+             => VertexID -> VertexID -> Vertex -> Vertex -> Merge p r Vertex+rotateR' l r = go+  where+    go r1 r2 = qTest l r r1 r2 >>= \case+                 True  -> pure r1+                 False -> do modify $ delete r (focus' r1)+                             go r2 (pred' r2)+++-- | Symmetric to rotateR+rotateL     :: (Ord r, Fractional r)+                     => VertexID -> VertexID -> Vertex -> Merge p r (Vertex, Bool)+rotateL l r l1 = focus' l1 `isRightOf` (r, l) >>= \case+                   True  -> (,False) <$> rotateL' l r l1 (succ' l1)+                   False -> pure (l1,True)++-- | The code that does the actual rotating. Symmetric to rotateR'+rotateL'     :: (Ord r, Fractional r)+             => VertexID -> VertexID -> Vertex -> Vertex -> Merge p r Vertex+rotateL' l r = go+  where+    go l1 l2 = qTest l r l1 l2 >>= \case+                 True  -> pure l1+                 False -> do modify $ delete l (focus' l1)+                             go l2 (succ' l2)++--------------------------------------------------------------------------------+-- * Primitives used by the Algorithm++-- | returns True if the forth point (vertex) does not lie in the disk defined+-- by the first three points.+qTest         :: (Ord r, Fractional r)+              => VertexID -> VertexID -> Vertex -> Vertex -> Merge p r Bool+qTest h i j k = withPtMap . snd . fst <$> ask+  where+    withPtMap ptMap = let h' = ptMap V.! h+                          i' = ptMap V.! i+                          j' = ptMap V.! (focus' j)+                          k' = ptMap V.! (focus' k)+                      in not . maybe True ((k'^.core) `insideBall`) $ disk' h' i' j'+    disk' p q r = disk (p^.core) (q^.core) (r^.core)++-- | Inserts an edge into the right position.+insert     :: (Num r, Ord r) => VertexID -> VertexID -> Merge p r ()+insert u v = do+               (mapping',fsts) <- ask+               modify $ insert' u v mapping'+               rotateToFirst u fsts+               rotateToFirst v fsts+++-- | make sure that the first vtx in the adj list of v is its predecessor on the CH+rotateToFirst        :: VertexID -> Firsts -> Merge p r ()+rotateToFirst v fsts = modify $ IM.adjust f v+  where+    mfst   = IM.lookup v fsts+    f  cl  = fromMaybe cl $ mfst >>= flip CL.rotateTo cl+++-- | Inserts an edge (and makes sure that the vertex is inserted in the+-- correct. pos in the adjacency lists)+insert'               :: (Num r, Ord r)+                      => VertexID -> VertexID -> Mapping p r -> Adj -> Adj+insert' u v (_,ptMap) = IM.adjustWithKey (insert'' v) u+                      . IM.adjustWithKey (insert'' u) v+  where+    -- inserts b into the adjacency list of a+    insert'' bi ai = CU.insertOrdBy (cwCmpAround (ptMap V.! ai) `on` (ptMap V.!)) bi+++-- | Deletes an edge+delete     :: VertexID -> VertexID -> Adj -> Adj+delete u v = IM.adjust (delete' v) u . IM.adjust (delete' u) v+  where+    delete' x = CL.filterL (/= x) -- should we rotate left or right if it is the focus?+++++-- | Lifted version of Convex.IsLeftOf+isLeftOf           :: (Ord r, Num r)+                   => VertexID -> (VertexID, VertexID) -> Merge p r Bool+p `isLeftOf` (l,r) = withPtMap . snd . fst <$> ask+  where+    withPtMap ptMap = (ptMap V.! p) `Convex.isLeftOf` (ptMap V.! l, ptMap V.! r)++-- | Lifted version of Convex.IsRightOf+isRightOf           :: (Ord r, Num r)+                    => VertexID -> (VertexID, VertexID) -> Merge p r Bool+p `isRightOf` (l,r) = withPtMap . snd . fst <$> ask+  where+    withPtMap ptMap = (ptMap V.! p) `Convex.isRightOf` (ptMap V.! l, ptMap V.! r)++--------------------------------------------------------------------------------+-- * Some Helper functions+++lookup'     :: Ord k => M.Map k a -> k -> a+lookup' m x = fromJust $ M.lookup x m++size'              :: BinLeafTree Size a -> Size+size' (Leaf _)     = 1+size' (Node _ s _) = s++-- | an 'unsafe' version of rotateTo that assumes the element to rotate to+-- occurs in the list.+rotateTo   :: Eq a => a -> CL.CList a -> CL.CList a+rotateTo x = fromJust . CL.rotateTo x++-- | Adjacency lists are stored in clockwise order, so pred means rotate right+pred' :: CL.CList a -> CL.CList a+pred' = CL.rotR++-- | Adjacency lists are stored in clockwise order, so pred and succ rotate left+succ' :: CL.CList a -> CL.CList a+succ' = CL.rotL++focus' :: CL.CList a -> a+focus' = fromJust . CL.focus++-- | Removes duplicates from a sorted list+nub' :: Eq a => NonEmpty.NonEmpty (a :+ b) -> NonEmpty.NonEmpty (a :+ b)+nub' = fmap NonEmpty.head . NonEmpty.groupBy1 ((==) `on` (^.core))+++withID     :: c :+ e -> e' -> c :+ (e :+ e')+withID p i = p&extra %~ (:+i)++lookup'' :: Int -> IM.IntMap a -> a+lookup'' k m = fromJust . IM.lookup k $ m

@@ -1,296 +0,0 @@-{-# LANGUAGE ScopedTypeVariables #-}-module Algorithms.Geometry.DelaunayTriangulation.DivideAndConqueror where--import           Algorithms.Geometry.ConvexHull.GrahamScan as GS-import           Algorithms.Geometry.DelaunayTriangulation.Types-import           Control.Lens-import           Control.Monad.Reader-import           Control.Monad.State-import           Data.BinaryTree-import qualified Data.CircularList as CL-import qualified Data.CircularSeq as CS-import qualified Data.CircularList.Util as CU-import           Data.Ext-import qualified Data.Foldable as F-import           Data.Function (on)-import           Data.Geometry-import           Data.Geometry.Ball (disk, insideBall)-import           Data.Geometry.Polygon-import qualified Data.Geometry.Polygon.Convex as Convex-import           Data.Geometry.Polygon.Convex (ConvexPolygon(..), simplePolygon)-import qualified Data.IntMap.Strict as IM-import qualified Data.List as L-import qualified Data.List.NonEmpty as NonEmpty-import qualified Data.Map as M-import           Data.Maybe (fromJust, fromMaybe)-import qualified Data.Vector as V------------------------------------------------------------------------------------ * Divide & Conqueror Delaunay Triangulation------ Implementation of the Divide & Conqueror algorithm as described in:------ Two Algorithms for Constructing a Delaunay Triangulation--- Lee and Schachter--- International Journal of Computer and Information Sciences, Vol 9, No. 3, 1980------ We store all adjacency lists in clockwise order------ : If v on the convex hull, then its first entry in the adj. lists is its CCW--- successor (i.e. its predecessor) on the convex hull------ Rotating Right <-> rotate clockwise----- | Computes the delaunay triangulation of a set of points.------ Running time: \(O(n \log n)\)--- (note: We use an IntMap in the implementation. So maybe actually \(O(n \log^2 n)\))------ pre: the input is a *SET*, i.e. contains no duplicate points. (If the--- input does contain duplicate points, the implementation throws them away)-delaunayTriangulation      :: (Ord r, Fractional r)-                           => NonEmpty.NonEmpty (Point 2 r :+ p) -> Triangulation p r-delaunayTriangulation pts' = Triangulation vtxMap ptsV adjV-  where-    pts    = nub' . NonEmpty.sortBy (compare `on` (^.core)) $ pts'-    ptsV   = V.fromList . F.toList $ pts-    vtxMap = M.fromList $ zip (map (^.core) . V.toList $ ptsV) [0..]--    tr     = _unElem <$> asBalancedBinLeafTree pts--    (adj,_) = delaunayTriangulation' tr (vtxMap,ptsV)-    adjV    = V.fromList . IM.elems $ adj------ : pre: - Input points are sorted lexicographically-delaunayTriangulation' :: (Ord r, Fractional r)-                       => BinLeafTree Size (Point 2 r :+ p)-                       -> Mapping p r-                       -> (Adj, ConvexPolygon (p :+ VertexID) r)-delaunayTriangulation' pts mapping'@(vtxMap,_)-  | size' pts == 1 = let (Leaf p) = pts-                         i        = lookup' vtxMap (p^.core)-                     in (IM.singleton i CL.empty, ConvexPolygon $ fromPoints [withID p i])-  | size' pts <= 3 = let pts'  = NonEmpty.fromList-                               . map (\p -> withID p (lookup' vtxMap (p^.core)))-                               . F.toList $ pts-                         ch    = GS.convexHull pts'-                     in (fromHull mapping' ch, ch)-  | otherwise      = let (Node lt _ rt) = pts-                         (ld,lch)       = delaunayTriangulation' lt mapping'-                         (rd,rch)       = delaunayTriangulation' rt mapping'-                         (ch, bt, ut)   = Convex.merge lch rch-                     in (merge ld rd bt ut mapping' (firsts ch), ch)------------------------------------------------------------------------------------- * Implementation---- | Mapping that says for each vtx in the convex hull what the first entry in--- the adj. list should be. The input polygon is given in Clockwise order-firsts :: ConvexPolygon (p :+ VertexID) r -> IM.IntMap VertexID-firsts = IM.fromList . map (\s -> (s^.end.extra.extra, s^.start.extra.extra))-       . F.toList . outerBoundaryEdges . _simplePolygon----- | Given a polygon; construct the adjacency list representation--- pre: at least two elements-fromHull              :: Ord r => Mapping p r -> ConvexPolygon (p :+ q) r -> Adj-fromHull (vtxMap,_) p = let vs@(u:v:vs') = map (lookup' vtxMap . (^.core))-                                         . F.toList . CS.rightElements-                                         $ p^.simplePolygon.outerBoundary-                            es           = zipWith3 f vs (tail vs ++ [u]) (vs' ++ [u,v])-                            f prv c nxt  = (c,CL.fromList . L.nub $ [prv, nxt])-                        in IM.fromList es----- | Merge the two delaunay triangulations.------ running time: \(O(n)\) (although we cheat a bit by using a IntMap)-merge                            :: (Ord r, Fractional r)-                                 => Adj-                                 -> Adj-                                 -> LineSegment 2 (p :+ VertexID) r -- ^ lower tangent-                                 -> LineSegment 2 (p :+ VertexID) r -- ^ upper tangent-                                 -> Mapping p r-                                 -> Firsts-                                 -> Adj-merge ld rd bt ut mapping'@(vtxMap,_) fsts =-    flip runReader (mapping', fsts) . flip execStateT adj $ moveUp (tl,tr) l r-  where-    l   = lookup' vtxMap (bt^.start.core)-    r   = lookup' vtxMap (bt^.end.core)-    tl  = lookup' vtxMap (ut^.start.core)-    tr  = lookup' vtxMap (ut^.end.core)-    adj = ld `IM.union` rd--type Merge p r = StateT Adj (Reader (Mapping p r, Firsts))--type Firsts = IM.IntMap VertexID---- | Merges the two delaunay traingulations.-moveUp          :: (Ord r, Fractional r)-                => (VertexID,VertexID) -> VertexID -> VertexID -> Merge p r ()-moveUp ut l r-  | (l,r) == ut = insert l r-  | otherwise   = do-                     insert l r-                     -- Get the neighbours of r and l along the convex hull-                     r1 <- pred' . rotateTo l . lookup'' r <$> get-                     l1 <- succ' . rotateTo r . lookup'' l <$> get--                     (r1',a) <- rotateR l r r1-                     (l1',b) <- rotateL l r l1-                     c       <- qTest l r r1' l1'-                     let (l',r') = case (a,b,c) of-                                     (True,_,_)          -> (focus' l1', r)-                                     (False,True,_)      -> (l,          focus' r1')-                                     (False,False,True)  -> (l,          focus' r1')-                                     (False,False,False) -> (focus' l1', r)-                     moveUp ut l' r'----- | ''rotates'' around r and removes all neighbours of r that violate the--- delaunay condition. Returns the first vertex (as a Neighbour of r) that--- should remain in the Delaunay Triangulation, as well as a boolean A that--- helps deciding if we merge up by rotating left or rotating right (See--- description in the paper for more info)-rotateR        :: (Ord r, Fractional r)-               => VertexID -> VertexID -> Vertex -> Merge p r (Vertex, Bool)-rotateR l r r1 = focus' r1 `isLeftOf` (l, r) >>= \case-                   True  -> (,False) <$> rotateR' l r r1 (pred' r1)-                   False -> pure (r1,True)---- | The code that does the actual rotating-rotateR'     :: (Ord r, Fractional r)-             => VertexID -> VertexID -> Vertex -> Vertex -> Merge p r Vertex-rotateR' l r = go-  where-    go r1 r2 = qTest l r r1 r2 >>= \case-                 True  -> pure r1-                 False -> do modify $ delete r (focus' r1)-                             go r2 (pred' r2)----- | Symmetric to rotateR-rotateL     :: (Ord r, Fractional r)-                     => VertexID -> VertexID -> Vertex -> Merge p r (Vertex, Bool)-rotateL l r l1 = focus' l1 `isRightOf` (r, l) >>= \case-                   True  -> (,False) <$> rotateL' l r l1 (succ' l1)-                   False -> pure (l1,True)---- | The code that does the actual rotating. Symmetric to rotateR'-rotateL'     :: (Ord r, Fractional r)-             => VertexID -> VertexID -> Vertex -> Vertex -> Merge p r Vertex-rotateL' l r = go-  where-    go l1 l2 = qTest l r l1 l2 >>= \case-                 True  -> pure l1-                 False -> do modify $ delete l (focus' l1)-                             go l2 (succ' l2)------------------------------------------------------------------------------------- * Primitives used by the Algorithm---- | returns True if the forth point (vertex) does not lie in the disk defined--- by the first three points.-qTest         :: (Ord r, Fractional r)-              => VertexID -> VertexID -> Vertex -> Vertex -> Merge p r Bool-qTest h i j k = withPtMap . snd . fst <$> ask-  where-    withPtMap ptMap = let h' = ptMap V.! h-                          i' = ptMap V.! i-                          j' = ptMap V.! (focus' j)-                          k' = ptMap V.! (focus' k)-                      in not . maybe True ((k'^.core) `insideBall`) $ disk' h' i' j'-    disk' p q r = disk (p^.core) (q^.core) (r^.core)---- | Inserts an edge into the right position.-insert     :: (Num r, Ord r) => VertexID -> VertexID -> Merge p r ()-insert u v = do-               (mapping',fsts) <- ask-               modify $ insert' u v mapping'-               rotateToFirst u fsts-               rotateToFirst v fsts----- | make sure that the first vtx in the adj list of v is its predecessor on the CH-rotateToFirst        :: VertexID -> Firsts -> Merge p r ()-rotateToFirst v fsts = modify $ IM.adjust f v-  where-    mfst   = IM.lookup v fsts-    f  cl  = fromMaybe cl $ mfst >>= flip CL.rotateTo cl----- | Inserts an edge (and makes sure that the vertex is inserted in the--- correct. pos in the adjacency lists)-insert'               :: (Num r, Ord r)-                      => VertexID -> VertexID -> Mapping p r -> Adj -> Adj-insert' u v (_,ptMap) = IM.adjustWithKey (insert'' v) u-                      . IM.adjustWithKey (insert'' u) v-  where-    -- inserts b into the adjacency list of a-    insert'' bi ai = CU.insertOrdBy (cwCmpAround (ptMap V.! ai) `on` (ptMap V.!)) bi----- | Deletes an edge-delete     :: VertexID -> VertexID -> Adj -> Adj-delete u v = IM.adjust (delete' v) u . IM.adjust (delete' u) v-  where-    delete' x = CL.filterL (/= x) -- should we rotate left or right if it is the focus?------- | Lifted version of Convex.IsLeftOf-isLeftOf           :: (Ord r, Num r)-                   => VertexID -> (VertexID, VertexID) -> Merge p r Bool-p `isLeftOf` (l,r) = withPtMap . snd . fst <$> ask-  where-    withPtMap ptMap = (ptMap V.! p) `Convex.isLeftOf` (ptMap V.! l, ptMap V.! r)---- | Lifted version of Convex.IsRightOf-isRightOf           :: (Ord r, Num r)-                    => VertexID -> (VertexID, VertexID) -> Merge p r Bool-p `isRightOf` (l,r) = withPtMap . snd . fst <$> ask-  where-    withPtMap ptMap = (ptMap V.! p) `Convex.isRightOf` (ptMap V.! l, ptMap V.! r)------------------------------------------------------------------------------------- * Some Helper functions---lookup'     :: Ord k => M.Map k a -> k -> a-lookup' m x = fromJust $ M.lookup x m--size'              :: BinLeafTree Size a -> Size-size' (Leaf _)     = 1-size' (Node _ s _) = s---- | an 'unsafe' version of rotateTo that assumes the element to rotate to--- occurs in the list.-rotateTo   :: Eq a => a -> CL.CList a -> CL.CList a-rotateTo x = fromJust . CL.rotateTo x---- | Adjacency lists are stored in clockwise order, so pred means rotate right-pred' :: CL.CList a -> CL.CList a-pred' = CL.rotR---- | Adjacency lists are stored in clockwise order, so pred and succ rotate left-succ' :: CL.CList a -> CL.CList a-succ' = CL.rotL--focus' :: CL.CList a -> a-focus' = fromJust . CL.focus---- | Removes duplicates from a sorted list-nub' :: Eq a => NonEmpty.NonEmpty (a :+ b) -> NonEmpty.NonEmpty (a :+ b)-nub' = fmap NonEmpty.head . NonEmpty.groupBy1 ((==) `on` (^.core))---withID     :: c :+ e -> e' -> c :+ (e :+ e')-withID p i = p&extra %~ (:+i)--lookup'' :: Int -> IM.IntMap a -> a-lookup'' k m = fromJust . IM.lookup k $ m

@@ -40,7 +40,10 @@                          deriving (Show,Eq) makeLenses ''Triangulation +type instance NumType   (Triangulation p r) = r+type instance Dimension (Triangulation p r) = 2 + type Mapping p r = (M.Map (Point 2 r) VertexID, V.Vector (Point 2 r :+ p))  @@ -59,11 +62,11 @@ tEdges = concatMap (\(i,ns) -> map (i,) . filter (> i) . C.toList $ ns)        . zip [0..] . V.toList . _neighbours -drawTriangulation :: IpeOut (Triangulation p r) (IpeObject r)-drawTriangulation = IpeOut $ \tr ->-    let es = map (uncurry ClosedLineSegment) . triangulationEdges $ tr-    in asIpeGroup $ map (\e -> asIpeObjectWith ipeLineSegment e mempty) es-+drawTriangulation :: IpeOut (Triangulation p r) Group r+drawTriangulation tr =+  ipeGroup [ iO $ ipeLineSegment e+           | e <- map (uncurry ClosedLineSegment) . triangulationEdges $ tr+           ]  -------------------------------------------------------------------------------- 

@@ -1,6 +1,17 @@+--------------------------------------------------------------------------------+-- |+-- Module      :  Algorithms.Geometry.EuclideanMST.EuclideanMST+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- \(O(n\log n)\) time algorithm algorithm to compute the Euclidean minimum+-- spanning tree of a set of \(n\) points in \(\mathbb{R}^2\).+--+-------------------------------------------------------------------------------- module Algorithms.Geometry.EuclideanMST.EuclideanMST where -import           Algorithms.Geometry.DelaunayTriangulation.DivideAndConqueror+import           Algorithms.Geometry.DelaunayTriangulation.DivideAndConquer import           Algorithms.Geometry.DelaunayTriangulation.Types import           Algorithms.Graph.MST import           Control.Lens@@ -12,7 +23,6 @@ import           Data.Proxy import           Data.Tree - --------------------------------------------------------------------------------  -- | Computes the Euclidean Minimum Spanning Tree. We compute the Delaunay@@ -38,9 +48,8 @@ data MSTW  -drawTree' :: IpeOut (Tree (Point 2 r :+ p)) (IpeObject r)-drawTree' = IpeOut $-  asIpeGroup . map (asIpeObject' mempty . uncurry ClosedLineSegment) . treeEdges+drawTree' :: IpeOut (Tree (Point 2 r :+ p)) Group r+drawTree' = ipeGroup . map (iO . defIO  . uncurry ClosedLineSegment) . treeEdges   treeEdges              :: Tree a -> [(a,a)]

@@ -17,7 +17,6 @@ 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           Data.Vinyl.CoRec

@@ -9,7 +9,6 @@ import           Data.Geometry.Point import           Data.Geometry.Properties import qualified Data.Map as M-import           Data.Semigroup import           Data.Vinyl import           Data.Vinyl.CoRec 

@@ -11,8 +11,8 @@ import qualified Data.List.NonEmpty as NonEmpty import qualified Data.List as L import qualified Data.Map as Map-import           Data.Semigroup +--------------------------------------------------------------------------------  -- get the endpoints of a line segment endPoints'   :: (HasEnd s, HasStart s) => s -> (StartCore s, EndCore s)

@@ -1,17 +1,18 @@ module Algorithms.Geometry.PolyLineSimplification.DouglasPeucker where -import Data.Semigroup-import Data.Ord(comparing)-import Control.Lens hiding (only)-import Data.Ext-import Data.Geometry.PolyLine-import Data.Geometry.Point-import Data.Geometry.Vector-import Data.Geometry.LineSegment-import qualified Data.Seq2 as S2-import qualified Data.Sequence as S+import           Control.Lens hiding (only)+import           Data.Ext import qualified Data.Foldable as F+import           Data.Geometry.LineSegment+import           Data.Geometry.Point+import           Data.Geometry.PolyLine+import           Data.Geometry.Vector+import qualified Data.LSeq as LSeq+import           Data.LSeq (LSeq, pattern (:|>))+import           Data.Ord (comparing) +--------------------------------------------------------------------------------+ -- | Line simplification with the well-known Douglas Peucker alogrithm. Given a distance -- value eps adn a polyline pl, constructs a simplification of pl (i.e. with -- vertices from pl) s.t. all other vertices are within dist eps to the@@ -24,7 +25,9 @@     | dst <= (eps*eps) = fromPoints [a,b]     | otherwise        = douglasPeucker eps pref `merge` douglasPeucker eps subf   where-    pts@(S2.Seq2 a _ b) = pl^.points+    pts     = pl^.points+    a       = LSeq.head pts+    b       = LSeq.last pts     (i,dst)             = maxDist pts (ClosedLineSegment a b)      (pref,subf)         = split i pl@@ -34,24 +37,27 @@  -- | Concatenate the two polylines, dropping their shared vertex merge          :: PolyLine d p r -> PolyLine d p r -> PolyLine d p r-merge pref sub = PolyLine $ pref' >+< (sub^.points)+merge pref sub = PolyLine $ pref' `append` (sub^.points)   where-    (pref' S2.:>> _) = S2.viewr $ pref^.points-    ~(a S2.:< as) >+< bs = S2.fromSeqUnsafe $ a S.<| as <> S2.toSeq bs+    (pref' :|> _) = pref^.points+    append     :: LSeq.LSeq n a ->  LSeq.LSeq m a -> LSeq.LSeq m a+    append a b = LSeq.promise $ LSeq.append a b+    -- our seq is actually even longer + -- | Split the polyline at the given vertex. Both polylines contain this vertex split                  :: Int -> PolyLine d p r                        -> (PolyLine d p r, PolyLine d p r) split i (PolyLine pts) = bimap f f (as,bs)   where-    f = PolyLine . S2.fromSeqUnsafe-    as = S2.take (i+1) pts-    bs = S2.drop i     pts+    f = PolyLine . LSeq.forceLSeq (C  :: C 2)+    as = LSeq.take (i+1) pts+    bs = LSeq.drop i     pts  -- | Given a sequence of points, find the index of the point that has the -- Furthest distance to the LineSegment. The result is the index of the point -- and this distance. maxDist       :: (Ord r, Fractional r, Arity d)-              => S2.Seq2 (Point d r :+ p) -> LineSegment d p r -> (Int,r)-maxDist pts s = F.maximumBy (comparing snd) . S2.mapWithIndex (\i (p :+ _) ->+              => LSeq n (Point d r :+ p) -> LineSegment d p r -> (Int,r)+maxDist pts s = F.maximumBy (comparing snd) . LSeq.mapWithIndex (\i (p :+ _) ->                                                      (i,sqDistanceToSeg p s)) $ pts

@@ -24,7 +24,6 @@ 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@@ -209,7 +208,7 @@  -- | 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)+getHelper i = do ui         <- gets (^?!helper.ix i)                  STR u _ ut <- asks (^.ix' ui)                  pure $ SP (u :+ ui) ut @@ -222,7 +221,7 @@   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)+handleSplit i (v :+ adj) = do ej <- gets $ \ss -> ss^?!statusStruct.to (lookupLE v)._Just                               let j = ej^.start.extra                               SP u _ <- getHelper j                               -- update the status struct:@@ -244,7 +243,7 @@ -- | 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)+connectToLeft i v = do ej <- gets $ \ss -> ss^?!statusStruct.to (lookupLE v)._Just                        let j = ej^.start.extra                        tellIfMerge i v j                        modify $ \ss -> ss&helper %~ IntMap.insert j i

@@ -12,7 +12,8 @@ import           Data.Geometry.PlanarSubdivision.Basic import           Data.Geometry.Polygon import           Data.PlaneGraph (PlaneGraph)-import           Data.Semigroup++--------------------------------------------------------------------------------  -- | Triangulates a polygon of \(n\) vertices --

@@ -10,7 +10,6 @@ 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)

@@ -10,7 +10,6 @@ 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 --------------------------------------------------------------------------------

@@ -1,3 +1,14 @@+--------------------------------------------------------------------------------+-- |+-- Module      :  Algorithms.Geometry.SmallestEnclosingBall.Naive+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Naive implementation to compute the smallest enclosing disk of a set of+-- points in \(\mathbb{R}^2\)+--+-------------------------------------------------------------------------------- module Algorithms.Geometry.SmallestEnclosingBall.Naive where  -- just for the types@@ -6,10 +17,10 @@ import Algorithms.Geometry.SmallestEnclosingBall.Types import Data.Geometry.Ball import Data.Geometry.Point-import Data.List(minimumBy)-import Data.Function(on)-import Data.Maybe(fromMaybe)-import Algorithms.Util+import Data.List (minimumBy)+import Data.Function (on)+import Data.Maybe (fromMaybe)+import Data.Util(STR(..),SP(..), uniquePairs, uniqueTriplets)  -------------------------------------------------------------------------------- @@ -24,12 +35,12 @@ smallestEnclosingDisk _           = error "smallestEnclosingDisk: Too few points"  pairs     :: Fractional r => [Point 2 r :+ p] -> [DiskResult p r]-pairs pts = [DiskResult (fromDiameter (a^.core) (b^.core)) (Two a b)+pairs pts = [ DiskResult (fromDiameter (a^.core) (b^.core)) (Two a b)             | SP a b <- uniquePairs pts]  triplets     :: (Ord r, Fractional r) => [Point 2 r :+ p] -> [DiskResult p r] triplets pts = [DiskResult (disk' a b c) (Three a b c)-               | ST a b c <- uniqueTriplets pts]+               | STR a b c <- uniqueTriplets pts]  disk'       :: (Ord r, Fractional r)             => Point 2 r :+ p -> Point 2 r :+ p -> Point 2 r :+ p -> Disk () r

@@ -1,5 +1,16 @@ {-# LANGUAGE DeriveFunctor  #-} {-# LANGUAGE TemplateHaskell  #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Algorithms.Geometry.SmallestEnclosingBall.RandomizedIncrementalConstruction+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- An randomized algorithm to compute the smallest enclosing disk of a set of+-- \(n\) points in \(\mathbb{R}^2\). The expected running time is \(O(n)\).+--+-------------------------------------------------------------------------------- module Algorithms.Geometry.SmallestEnclosingBall.RandomizedIncrementalConstruction where  import           Algorithms.Geometry.SmallestEnclosingBall.Types@@ -13,7 +24,7 @@ import           System.Random import           System.Random.Shuffle (shuffle) -+--------------------------------------------------------------------------------  -- | O(n) expected time algorithm to compute the smallest enclosing disk of a -- set of points. we need at least two points.

@@ -1,8 +1,18 @@ {-# LANGUAGE DeriveFunctor  #-} {-# LANGUAGE TemplateHaskell  #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Algorithms.Geometry.SmallestEnclosingBall.Types+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Types to represent the smallest enclosing disk of a set of points in+-- \(\mathbb{R}^2\)+--+-------------------------------------------------------------------------------- module Algorithms.Geometry.SmallestEnclosingBall.Types where -import           Data.Monoid import qualified Data.Foldable as F import           Data.Geometry import           Data.Geometry.Ball

@@ -1,5 +1,15 @@ {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE UndecidableInstances #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Algorithms.Geometry.Sweep+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Helper types and functions for implementing Sweep line algorithms.+--+-------------------------------------------------------------------------------- module Algorithms.Geometry.Sweep where  import qualified Data.Map as Map@@ -8,13 +18,15 @@ 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 ++-- | Represent a computation that needs a particular time as input. newtype Timed s t a = Timed {atTime :: (Tagged s t) -> a }  @@ -24,6 +36,7 @@ instance (Reifies s t, Ord k) => Eq (Timed s t k) where   a == b = a `compare` b == EQ +-- | Comparison function for timed values compare_                       :: forall s t k. (Ord k, Reifies s t)                                => Timed s t k -> Timed s t k                                -> Ordering@@ -31,9 +44,11 @@                                  in f (Tagged t) `compare` g (Tagged t)  -coerceTo :: proxy s -> f (Timed s' t k) v -> f (Timed s t k) v+-- | Coerce timed values+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 ()) @@ -46,10 +61,15 @@             -> 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     :: proxy s -> Int -> Timed s Int Int constT _ i = Timed (const i)  

@@ -1,5 +1,15 @@ {-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE UndecidableInstances  #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Algorithms.Geometry.WellSeparatedPairDecomposition.Types+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Data types that can represent a well separated pair decomposition (wspd).+--+-------------------------------------------------------------------------------- module Algorithms.Geometry.WellSeparatedPairDecomposition.Types where  import           Control.Lens hiding (Level)@@ -8,8 +18,7 @@ import           Data.Geometry.Box import           Data.Geometry.Point import           Data.Geometry.Vector-import           Data.Semigroup-import qualified Data.Seq2 as S2+import qualified Data.LSeq as LSeq import qualified Data.Sequence as S import qualified Data.Traversable as Tr @@ -46,7 +55,7 @@ -------------------------------------------------------------------------------- -- * Implementation types -type PointSeq d p r = S2.ViewL1 (Point d r :+ p)+type PointSeq d p r = LSeq.LSeq 1 (Point d r :+ p)   data Level = Level { _unLevel   :: Int

@@ -1,4 +1,13 @@-{-# LANGUAGE LambdaCase  #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Algorithms.Geometry.WellSeparatedPairDecomposition.WSPD+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Algorithm to construct a well separated pair decomposition (wspd).+--+-------------------------------------------------------------------------------- module Algorithms.Geometry.WellSeparatedPairDecomposition.WSPD where  import           Algorithms.Geometry.WellSeparatedPairDecomposition.Types@@ -9,26 +18,26 @@ import           Data.Ext import qualified Data.Foldable as F import           Data.Geometry.Box-import           Data.Geometry.Transformation-import           Data.Geometry.Properties import           Data.Geometry.Point+import           Data.Geometry.Properties+import           Data.Geometry.Transformation import           Data.Geometry.Vector import qualified Data.Geometry.Vector as GV+import qualified Data.IntMap.Strict as IntMap+import qualified Data.LSeq as LSeq+import           Data.LSeq (LSeq(..),toSeq,ViewL(..),ViewR(..),pattern (:<|),pattern (:|>)) import qualified Data.List as L import qualified Data.List.NonEmpty as NonEmpty import           Data.Maybe import           Data.Ord (comparing) import           Data.Range import qualified Data.Range as Range-import           Data.Semigroup-import qualified Data.Seq2 as S2 import qualified Data.Sequence as S import qualified Data.Vector as V import qualified Data.Vector.Mutable as MV import           GHC.TypeLits-import qualified Data.IntMap.Strict as IntMap -import Debug.Trace+import           Debug.Trace  -------------------------------------------------------------------------------- @@ -45,7 +54,7 @@     n    = length $ pts'^.GV.element (C :: C 0)      sortOn' i = NonEmpty.sortWith (^.core.unsafeCoord i)-    sortOn  i = S2.viewL1FromNonEmpty . sortOn' (i + 1)+    sortOn  i = LSeq.fromNonEmpty . sortOn' (i + 1)     -- sorts the points on the first coordinate, and then associates each point     -- with an index,; its rank in terms of this first coordinate.     g = NonEmpty.zipWith (\i (p :+ e) -> p :+ (i :+ e)) (NonEmpty.fromList [0..])@@ -110,7 +119,7 @@                      => Int -> GV.Vector d (PointSeq d (Idx :+ p) r)                      -> BinLeafTree Int (Point d r :+ p) fairSplitTree' n pts-    | n <= 1    = let (p S2.:< _) = pts^.GV.element (C :: C 0) in Leaf (dropIdx p)+    | n <= 1    = let p = LSeq.head $ pts^.GV.element (C :: C 0) in Leaf (dropIdx p)     | otherwise = foldr node' (V.last path) $ V.zip nodeLevels (V.init path)   where     -- note that points may also be assigned level 'Nothing'.@@ -168,7 +177,7 @@     level p = maybe (k-1) _unLevel $ levels V.! (p^.extra.core)     append v i p = MV.read v i >>= MV.write v i . (S.|> p) -+fromSeqUnsafe = LSeq.promise . LSeq.fromSeq  -- | Given a sequence of points, whose index is increasing in the first -- dimension, i.e. if idx p < idx q, then p[0] < q[0].@@ -250,9 +259,9 @@ -- | Remove allready assigned points from the sequence -- -- pre: there are points remaining-compactEnds'               :: PointSeq d (Idx :+ p) r-                           -> RST s (PointSeq d (Idx :+ p) r)-compactEnds' (l0 S2.:< s0) = fmap fromSeqUnsafe . goL $ l0 S.<| s0+compactEnds'              :: PointSeq d (Idx :+ p) r+                          -> RST s (PointSeq d (Idx :+ p) r)+compactEnds' (l0 :<| s0) = fmap fromSeqUnsafe . goL $ l0 S.<| toSeq s0   where     goL s@(S.viewl -> l S.:< s') = hasLevel l >>= \case                                      False -> goR s@@ -288,7 +297,7 @@                                  -> RST s ( PointSeq d (Idx :+ p) r                                           , PointSeq d (Idx :+ p) r                                           )-findAndCompact j (l0 S2.:< s0) m = fmap select . stepL $ l0 S.<| s0+findAndCompact j (l0 :<| s0) m = fmap select . stepL $ l0 S.<| toSeq s0   where     -- stepL and stepR together build a data structure (FAC l r S) that     -- contains the left part of the list, i.e. the points before midpoint, and@@ -345,10 +354,9 @@ -- -- pre: points are sorted according to their dimension extends :: Arity d => GV.Vector d (PointSeq d p r) -> GV.Vector d (Range r)-extends = GV.imap (\i pts@(l S2.:< _) ->-                     let (_ S2.:> r) = S2.viewL1toR1 pts-                     in ClosedRange (l^.core.unsafeCoord (i + 1))-                                    (r^.core.unsafeCoord (i + 1)))+extends = GV.imap (\i pts ->+                     ClosedRange ((LSeq.head pts)^.core.unsafeCoord (i + 1))+                                 ((LSeq.last pts)^.core.unsafeCoord (i + 1)))   --------------------------------------------------------------------------------@@ -442,11 +450,6 @@ children'              :: BinLeafTree v a -> [BinLeafTree v a] children' (Leaf _)     = [] children' (Node l _ r) = [l,r]---fromSeqUnsafe                         :: S.Seq a -> S2.ViewL1 a-fromSeqUnsafe (S.viewl -> (l S.:< s)) = l S2.:< s-fromSeqUnsafe _                       = error "fromSeqUnsafe: Empty seq"   -- | Turn a traversal into lens

@@ -1,20 +0,0 @@-module Algorithms.Util where--import qualified Data.List as L---data SP a b = SP !a !b deriving (Eq,Ord,Show)---- | Given a list xs, generate all unique (unordered) pairs.-uniquePairs    :: [a] -> [SP a a]-uniquePairs xs = [ SP x y | (x:ys) <- nonEmptyTails xs, y <- ys ]--nonEmptyTails :: [a] -> [[a]]-nonEmptyTails = L.init . L.tails---data ST a b c = ST !a !b !c deriving (Eq,Ord,Show)---- | All unieuqe unordered triplets.-uniqueTriplets    :: [a] -> [ST a a a]-uniqueTriplets xs = [ ST x y z | (x:ys) <- nonEmptyTails xs, SP y z <- uniquePairs ys]

@@ -1,12 +1,21 @@ {-# Language DeriveFunctor#-} {-# Language FunctionalDependencies #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.BinaryTree+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Several types of Binary trees.+--+-------------------------------------------------------------------------------- module Data.BinaryTree where  import           Control.DeepSeq import           Data.List.NonEmpty (NonEmpty(..),(<|)) import qualified Data.List.NonEmpty as NonEmpty import           Data.Maybe (mapMaybe)-import           Data.Semigroup import           Data.Semigroup.Foldable import qualified Data.Tree as Tree import qualified Data.Vector as V@@ -14,6 +23,8 @@  -------------------------------------------------------------------------------- +-- | Binary tree that stores its values (of type a) in the leaves. Internal+-- nodes store something of type v. data BinLeafTree v a = Leaf !a                      | Node (BinLeafTree v a) !v (BinLeafTree v a)                      deriving (Show,Read,Eq,Ord,Functor,Generic)@@ -146,7 +157,7 @@ -------------------------------------------------------------------------------- -- * Internal Node Tree -+-- | Binary tree in which we store the values of type a in internal nodes. data BinaryTree a = Nil                   | Internal (BinaryTree a) !a (BinaryTree a)                   deriving (Show,Read,Eq,Ord,Functor,Foldable,Traversable,Generic)@@ -157,9 +168,9 @@ access Nil              = Nothing access (Internal _ x _) = Just x --- | Create a balanced binary tree+-- | Create a balanced binary tree. ----- \(O(n)\)+-- running time: \(O(n)\) asBalancedBinTree :: [a] -> BinaryTree a asBalancedBinTree = mkTree . V.fromList   where@@ -170,7 +181,7 @@                             else Internal (mkTree $ V.slice 0 h v) x                                           (mkTree $ V.slice (h+1) (n - h -1) v) -+-- | Fold function for folding over a binary tree. foldBinaryUp                      :: b -> (a -> b -> b -> b)                                   -> BinaryTree a -> BinaryTree (a,b) foldBinaryUp _ _ Nil              = Nil@@ -180,9 +191,11 @@                                         b  = f x (g l') (g r')                                     in Internal l' (x,b) r' +-- | Convert a @BinaryTree@ into a RoseTree toRoseTree'                  :: BinaryTree a -> Maybe (Tree.Tree a) toRoseTree' Nil              = Nothing toRoseTree' (Internal l v r) = Just $ Tree.Node v $ mapMaybe toRoseTree' [l,r] +-- | Draw a binary tree. drawTree' :: Show a => BinaryTree a -> String drawTree' = maybe "Nil" (Tree.drawTree . fmap show) . toRoseTree'

@@ -1,7 +1,6 @@ module Data.BinaryTree.Zipper where  import Data.BinaryTree-import Data.Semigroup  -------------------------------------------------------------------------------- 

@@ -38,7 +38,6 @@ import qualified Data.List as L import qualified Data.List.NonEmpty as NonEmpty import           Data.Maybe (listToMaybe)-import           Data.Semigroup import           Data.Semigroup.Foldable hiding (toNonEmpty) import           Data.Sequence ((|>),(<|),ViewL(..),ViewR(..),Seq) import qualified Data.Sequence as S

@@ -1,26 +1,32 @@ {-# 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+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- A pair-like data type to represent a 'core' type that has extra information+-- as well.+--+-------------------------------------------------------------------------------- module Data.Ext where +import Control.DeepSeq import Control.Lens hiding ((.=))+import Data.Aeson+import Data.Aeson.Types (typeMismatch) import Data.Biapplicative import Data.Bifoldable import Data.Bifunctor.Apply import Data.Bitraversable import Data.Functor.Apply (liftF2)-import Data.Semigroup+import Data.Geometry.Properties import Data.Semigroup.Bifoldable 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@@ -28,6 +34,8 @@ data core :+ extra = core :+ extra deriving (Show,Read,Eq,Ord,Bounded,Generic,NFData) infixr 1 :+ +type instance NumType   (core :+ ext) = NumType   core+type instance Dimension (core :+ ext) = Dimension core  instance Bifunctor (:+) where   bimap f g (c :+ e) = f c :+ g e

@@ -1,9 +1,14 @@-{-|-Module    : Data.Geometry-Description: Basic Geometry types-Copyright : (c) Frank Staals-License : See LICENCE file--}++--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Basic Geometry Types+--+-------------------------------------------------------------------------------- module Data.Geometry( module Data.Geometry.Properties                     , module Data.Geometry.Transformation                     , module Data.Geometry.Point@@ -26,3 +31,5 @@ import Data.Geometry.Transformation -- import Linear.Affine hiding (Point, Vector, origin) -- import Linear.Vector++--------------------------------------------------------------------------------

@@ -0,0 +1,25 @@+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.Arrangement+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Data type for representing an Arrangement of lines in \(\mathbb{R}^2\).+--+--------------------------------------------------------------------------------+module Data.Geometry.Arrangement( Arrangement(..)+                                , inputLines, subdivision, boundedArea, unboundedIntersections+                                , ArrangementBoundary++                                , constructArrangement+                                , constructArrangementInBox+                                , constructArrangementInBox'++                                , traverseLine+                                , findStart, findStartVertex, findStartDart+                                , follow+                                ) where+++import Data.Geometry.Arrangement.Internal

@@ -0,0 +1,24 @@+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.Arrangement.Draw+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Functions for Drawing arrangements++--------------------------------------------------------------------------------+module Data.Geometry.Arrangement.Draw where++import Control.Lens+import Data.Geometry.Arrangement+import Data.Geometry.Ipe+import Data.Geometry.PlanarSubdivision.Draw++-- | Draws an arrangement+drawArrangement :: IpeOut (Arrangement s l v e f r) Group r+drawArrangement = drawPlanarSubdivision' . view subdivision++-- | Draws an arrangement+drawColoredArrangement :: IpeOut (Arrangement s l v e (Maybe (IpeColor r)) r) Group r+drawColoredArrangement = drawColoredPlanarSubdivision . view subdivision

@@ -0,0 +1,328 @@+{-# LANGUAGE TemplateHaskell  #-}+{-# LANGUAGE ScopedTypeVariables  #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.Arrangement.Internal+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Data type for representing an Arrangement of lines in \(\mathbb{R}^2\).+--+--------------------------------------------------------------------------------+module Data.Geometry.Arrangement.Internal where++import           Control.Lens+import qualified Data.CircularSeq as CSeq+import           Data.Ext+import qualified Data.Foldable as F+import           Data.Geometry.Boundary+import           Data.Geometry.Box+import           Data.Geometry.Line+import           Data.Geometry.LineSegment+import           Data.Geometry.PlanarSubdivision+import           Data.Geometry.Point+import           Data.Geometry.Properties+import qualified Data.List as List+import           Data.Maybe+import           Data.Ord (Down(..))+import           Data.Proxy+import           Data.Sequence.Util+import qualified Data.Vector as V+import           Data.Vinyl.CoRec++--------------------------------------------------------------------------------++type ArrangementBoundary s e r = V.Vector (Point 2 r, VertexId' s, Maybe (Line 2 r :+ e))++-- | Data type representing a two dimensional planar arrangement+data Arrangement s l v e f r = Arrangement {+    _inputLines             :: V.Vector (Line 2 r :+ l)+  , _subdivision            :: PlanarSubdivision s v e f r+  , _boundedArea            :: Rectangle () r+  , _unboundedIntersections :: ArrangementBoundary s l r+  } deriving (Show,Eq)+  -- unboundedIntersections also stores the corners of the box. They are not+  -- associated with any line+makeLenses ''Arrangement++type instance NumType   (Arrangement s l v e f r) = r+type instance Dimension (Arrangement s l v e f r) = 2++--------------------------------------------------------------------------------++-- | Builds an arrangement of \(n\) lines+--+-- running time: \(O(n^2\log n\)+constructArrangement       :: (Ord r, Fractional r)+                           => proxy s+                           -> [Line 2 r :+ l]+                           -> Arrangement s l () (Maybe l) () r+constructArrangement px ls = let b  = makeBoundingBox ls+                             in constructArrangementInBox' px b ls++-- | Constructs the arrangemnet inside the box.  note that the resulting box+-- may be larger than the given box to make sure that all vertices of the+-- arrangement actually fit.+--+-- running time: \(O(n^2\log n\)+constructArrangementInBox            :: (Ord r, Fractional r)+                                     => proxy s+                                     -> Rectangle () r+                                     -> [Line 2 r :+ l]+                                     -> Arrangement s l () (Maybe l) () r+constructArrangementInBox px rect ls = let b  = makeBoundingBox ls+                                       in constructArrangementInBox' px (b <> rect) ls+++-- | Constructs the arrangemnet inside the box. (for parts to be useful, it is+-- assumed this boxfits at least the boundingbox of the intersections in the+-- Arrangement)+constructArrangementInBox'            :: (Ord r, Fractional r)+                                      => proxy s+                                      -> Rectangle () r+                                      -> [Line 2 r :+ l]+                                      -> Arrangement s l () (Maybe l) () r+constructArrangementInBox' px rect ls =+    Arrangement (V.fromList ls) subdiv rect (link parts' subdiv)+  where+    subdiv = fromConnectedSegments px segs+                & rawVertexData.traverse.dataVal .~ ()+    (segs,parts') = computeSegsAndParts rect ls++computeSegsAndParts         :: (Ord r, Fractional r)+                            => Rectangle () r+                            -> [Line 2 r :+ l]+                            -> ( [LineSegment 2 () r :+ Maybe l]+                               , [(Point 2 r, Maybe (Line 2 r :+ l))]+                               )+computeSegsAndParts rect ls = ( segs <> boundarySegs, parts')+  where+    segs         = map (&extra %~ Just)+                 . concatMap (\(l,ls') -> perLine rect l ls') $ makePairs ls+    boundarySegs = map (:+ Nothing) . toSegments . dupFirst $ map fst parts'+    dupFirst = \case []       -> []+                     xs@(x:_) -> xs ++ [x]+    parts'       = unBoundedParts rect ls+++perLine       :: forall r l. (Ord r, Fractional r)+              => Rectangle () r -> Line 2 r :+ l -> [Line 2 r :+ l]+              -> [LineSegment 2 () r :+ l]+perLine b m ls = map (:+ m^.extra) . toSegments . rmDuplicates . List.sort $ vs <> vs'+  where+    rmDuplicates = map head . List.group+    vs  = mapMaybe (m `intersectionPoint`) ls+    vs' = maybe [] (\(p,q) -> [p,q]) . asA (Proxy :: Proxy (Point 2 r, Point 2 r))+        $ (m^.core) `intersect` (Boundary b)+++intersectionPoint                   :: (Ord r, Fractional r)+                                    => Line 2 r :+ l -> Line 2 r :+ l -> Maybe (Point 2 r)+intersectionPoint (l :+ _) (m :+ _) = asA (Proxy :: Proxy (Point 2 r)) $ l `intersect` m+++toSegments      :: Ord r => [Point 2 r] -> [LineSegment 2 () r]+toSegments ps = let pts = map ext $ ps in+  zipWith ClosedLineSegment pts (tail pts)+++-- | Constructs a boundingbox containing all intersections+--+-- running time: \(O(n^2)\), where \(n\) is the number of input lines+makeBoundingBox :: (Ord r, Fractional r) => [Line 2 r :+ l] -> Rectangle () r+makeBoundingBox = grow 1 . boundingBoxList' . intersections++-- | Computes all intersections+intersections :: (Ord r, Fractional r) => [Line 2 r :+ l] -> [Point 2 r]+intersections = mapMaybe (uncurry intersectionPoint) . allPairs+++-- intersections :: forall p r. (Ord r, Fractional r)+--               => [Line 2 r :+ p] -> Map.Map (Point 2 r) (NonEmpty (Line 2 r :+ p))+-- intersections = Map.map sortNub . collect+--               . mapMaybe (\(l,m) -> (l, m,) <$> f l m) . allPairs+--   where+--     f (l :+ _) (m :+ _) = asA (Proxy :: Proxy (Point 2 r)) $ l `intersect` m++++-- collect :: Ord k => [(v,v,k)] -> Map.Map k (NonEmpty v)+-- collect = foldr f mempty+--   where+--     f (l,m,p) = Map.insertWith (<>) p (NonEmpty.fromList [l,m])++-- sortNub :: Ord r => NonEmpty (Line 2 r :+ p) -> NonEmpty (Line 2 r :+ p)+-- sortNub = fmap (NonEmpty.head) .  groupLines++-- groupLines :: Ord r => NonEmpty (Line 2 r :+ p)+--            -> NonEmpty (NonEmpty (Line 2 r :+ p))+-- groupLines = NonEmpty.groupWith1 L2 . NonEmpty.sortWith L2+++-- -- | Newtype wrapper that allows us to sort lines+-- newtype L2 r p = L2 (Line 2 r :+ p) deriving (Show)++-- instance Eq r => Eq (L2 r p) where+--   (L2 (Line p u :+ _)) == (L2 (Line q v :+ _)) = (p,u) == (q,v)+-- instance Ord r => Ord (L2 r p) where+--   (L2 (Line p u :+ _)) `compare` (L2 (Line q v :+ _)) = p `compare` q <> u `compare` v++-- -- | Collect the intersection points per line+-- byLine :: Ord r+--        => Map.Map (Point 2 r) (NonEmpty (Line 2 r :+ p))+--        -> Map.Map (L2 r p)    (NonEmpty (Point 2 r))+-- byLine = foldr f mempty . flatten . Map.assocs+--   where+--     flatten = concatMap (\(p,ls) -> map (\l -> (L2 l,p)) $ NonEmpty.toList ls)+--     f (l,p) = Map.insertWith (<>) l $ NonEmpty.fromList [p]+++-- | Computes the intersections with a particular side+sideIntersections      :: (Ord r, Fractional r)+                       => [Line 2 r :+ l] -> LineSegment 2 q r+                       -> [(Point 2 r, Line 2 r :+ l)]+sideIntersections ls s = let l   = supportingLine s :+ undefined+                         in List.sortOn fst . filter (flip onSegment s . fst)+                          . mapMaybe (\m -> (,m) <$> l `intersectionPoint` m) $ ls++-- | Constructs the unbounded intersections. Reported in clockwise direction.+unBoundedParts         :: (Ord r, Fractional r)+                       => Rectangle () r+                       -> [Line 2 r :+ l]+                       -> [(Point 2 r, Maybe (Line 2 r :+ l))]+unBoundedParts rect ls = [tl] <> t <> [tr] <> reverse r <> [br] <> reverse b <> [bl] <> l+  where+    sideIntersections' = over (traverse._2) Just . sideIntersections ls+    (t,r,b,l)     = map4 sideIntersections'      $ sides   rect+    (tl,tr,br,bl) = map4 ((,Nothing) . (^.core)) $ corners rect+++map4              :: (a -> b) -> (a,a,a,a) -> (b,b,b,b)+map4 f (a,b',c,d) = (f a, f b', f c, f d)++-- | Links the vertices  of the outer boundary with those in the subdivision+link       :: Eq r => [(Point 2 r, a)] -> PlanarSubdivision s v (Maybe e) f r+           -> V.Vector (Point 2 r, VertexId' s, a)+link vs ps = V.fromList . map (\((p,x),(_,y)) -> (p,y,x)) . F.toList+           . fromJust' $ alignWith (\(p,_) (q,_) -> p == q) (CSeq.fromList vs) vs'+  where+    vs' = CSeq.fromList . map (\v -> (ps^.locationOf v,v) ) . V.toList+        $ boundaryVertices (outerFaceId ps) ps+    fromJust' = fromMaybe (error "Data.Geometry.Arrangement.link: fromJust")++--------------------------------------------------------------------------------++makePairs :: [a] -> [(a,[a])]+makePairs = go+  where+    go []     = []+    go (x:xs) = (x,xs) : map (\(y,ys) -> (y,x:ys)) (go xs)++allPairs    :: [a] -> [(a,a)]+allPairs ys = go ys+  where+    go []     = []+    go (x:xs) = map (x,) xs ++ go xs++-- | Given a predicate that tests if two elements of a CSeq match, find a+-- rotation of the seqs such at they match.+--+-- Running time: \(O(n)\)+alignWith         :: (a -> b -> Bool) -> CSeq.CSeq a -> CSeq.CSeq b+                  -> Maybe (CSeq.CSeq (a,b))+alignWith p xs ys = CSeq.zipL xs <$> CSeq.findRotateTo (p (CSeq.focus xs)) ys++--------------------------------------------------------------------------------++-- | Given an Arrangement and a line in the arrangement, follow the line+-- through he arrangement.+--+traverseLine       :: (Ord r, Fractional r)+                   => Line 2 r -> Arrangement s l v (Maybe e) f r -> [Dart s]+traverseLine l arr = let md    = findStart l arr+                         dup x = (x,x)+                     in maybe [] (List.unfoldr (fmap dup . follow arr)) md++-- | Find the starting point of the line  the arrangement+findStart       :: (Ord r, Fractional r)+                => Line 2 r -> Arrangement s l v (Maybe e) f r -> Maybe (Dart s)+findStart l arr = do+    (p,_)   <- asA (Proxy :: Proxy (Point 2 r, Point 2 r)) $+                 l `intersect` (Boundary $ arr^.boundedArea)+    (_,v,_) <- findStartVertex p arr+    findStartDart (arr^.subdivision) v++++-- | Given a point on the boundary of the boundedArea box; find the vertex+--  this point corresponds to.+--+-- running time: \(O(\log n)\)+--+-- basically; maps every point to a tuple of the point and the side the+-- point occurs on. We then binary search to find the point we are looking+-- for.+findStartVertex       :: (Ord r, Fractional r)+                      => Point 2 r+                      -> Arrangement s l v e f r+                      -> Maybe (Point 2 r, VertexId' s, Maybe (Line 2 r :+ l))+findStartVertex p arr = do+    ss <- findSide p+    i  <- binarySearchVec (pred' ss) (arr^.unboundedIntersections)+    pure $ arr^.unboundedIntersections.singular (ix i)+  where+    (t,r,b,l) = sides'' $ arr^.boundedArea+    sides'' = map4 (\(ClosedLineSegment a c) -> LineSegment (Closed a) (Open c)) . sides++    findSide q = fmap fst . List.find (onSegment q . snd) $ zip [1..] [t,r,b,l]++    pred' ss (q,_,_) = let Just j = findSide q+                           x      = before (ss,p) (j,q)+                       in  x == LT || x == EQ++    before (i,p') (j,q') = case i `compare` j of+                                LT -> LT+                                GT -> GT+                                EQ | i == 2 || i == 3 -> Down p' `compare` Down q'+                                   | otherwise        -> p' `compare` q'+++-- | Find the starting dart of the given vertex v. Reports a dart s.t.+-- tailOf d = v+--+-- running me: \(O(k)\) where \(k\) is the degree of the vertex+findStartDart      :: PlanarSubdivision s v (Maybe e) f r -> VertexId' s -> Maybe (Dart s)+findStartDart ps v = V.find (\d -> isJust $ ps^.dataOf d) $ incidentEdges v ps+    -- the "real" dart is the one that has ata associated to it.+++-- | Given a dart d that incoming to v (headOf d == v), find the outgoing dart+-- colinear with the incoming one. Again reports dart d' s.t. tailOf d' = v+--+-- running time: \(O(k)\), where k is the degree of the vertex d points to+follow       :: (Ord r, Num r) => Arrangement s l v e f r -> Dart s -> Maybe (Dart s)+follow arr d = V.find extends $ incidentEdges v ps+  where+    ps = arr^.subdivision+    v  = headOf d ps+    (up,vp) = over both (^.location) $ endPointData d ps++    extends d' = let wp = ps^.locationOf (headOf d' ps)+                 in d' /= twin d && ccw up vp wp == CoLinear++--------------------------------------------------------------------------------++-- TODO: we can skip the findStart by just traversing from all boundary points++-- computeFaceData :: (Arrangement s v e f r -> Dart s -> f')+--                -> Arrangement s v e f r -> V.Vertex f'+-- computeFaceData arr f = fmap fromJust . V.create $ do+--                           v <- MV.replicate (arr^.subdivision.to numFaces) Nothing+--                           mapM_ (computeFaceData' arr f v) $ arr^.inputLines+--                           pure v+++-- computeFaceData' arr f v l = mapM_ (assign ) traverseLine arr l++--------------------------------------------------------------------------------

@@ -1,11 +1,15 @@ {-# 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+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- \(d\)-dimensional Balls and Spheres+--+-------------------------------------------------------------------------------- module Data.Geometry.Ball where  import           Control.DeepSeq

@@ -3,12 +3,16 @@ {-# 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+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Orthogonal \(d\)-dimensiontal boxes (e.g. rectangles)+--+-------------------------------------------------------------------------------- module Data.Geometry.Box( module Data.Geometry.Box.Internal                         , topSide, leftSide, bottomSide, rightSide                         , sides, sides'

@@ -2,6 +2,16 @@ {-# LANGUAGE ScopedTypeVariables  #-} {-# LANGUAGE UndecidableInstances  #-} {-# LANGUAGE InstanceSigs  #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.Box.Internal+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Orthogonal \(d\)-dimensiontal boxes (e.g. rectangles)+--+-------------------------------------------------------------------------------- module Data.Geometry.Box.Internal where  import           Control.DeepSeq@@ -11,19 +21,17 @@ import           Data.Geometry.Point import           Data.Geometry.Properties import           Data.Geometry.Transformation-import           Data.Geometry.Vector (Vector, Arity, C(..))+import           Data.Geometry.Vector import qualified Data.Geometry.Vector as V import qualified Data.List.NonEmpty as NE import           Data.Proxy import qualified Data.Range as R-import           Data.Semigroup import qualified Data.Semigroup.Foldable as F import qualified Data.Vector.Fixed as FV import           Data.Vinyl.CoRec (asA) import           GHC.Generics (Generic) import           GHC.TypeLits - --------------------------------------------------------------------------------  -- | Coordinate wize minimum@@ -56,6 +64,12 @@ -- coordinates, create a box. box          :: Point d r :+ p -> Point d r :+ p -> Box d p r box low high = Box (low&core %~ CWMin) (high&core %~ CWMax)++-- | 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 v = V.replicate x+           in b&minP.core.cwMin %~ (.-^ v)+               &maxP.core.cwMax %~ (.+^ v)  -- | Build a d dimensional Box given d ranges. fromExtent    :: Arity d => Vector d (R.Range r) -> Box d () r

@@ -60,23 +60,18 @@ --------------------------------------------------------------------------------  halfLineToSubLine                :: (Arity d, Num r)-                                 => HalfLine d r -> SubLine d () (UnBounded r)-halfLineToSubLine (HalfLine p v) = let l = fmap Val $ Line p v+                                 => HalfLine d r -> SubLine d () (UnBounded r) r+halfLineToSubLine (HalfLine p v) = let l = Line p v                                    in SubLine l (Interval (Closed $ ext (Val 0))                                                           (Open   $ ext MaxInfinity))  -fromSubLine               :: (Num r, Arity d) => SubLine d p (UnBounded r) -> Maybe (HalfLine d r)-fromSubLine (SubLine l' i) = case (i^.start.core, i^.end.core) of-                               (Val x, MaxInfinity) -> f x-                               (MinInfinity, Val x) -> f x-                               _                    -> Nothing-  where-    f x = (\l@(Line _ v) -> HalfLine (pointAt x l) v)-       <$> T.mapM unBoundedToMaybe l'---+fromSubLine               :: (Num r, Arity d) => SubLine d p (UnBounded r) r+                          -> Maybe (HalfLine d r)+fromSubLine (SubLine l i) = case (i^.start.core, i^.end.core) of+   (Val x, MaxInfinity) -> Just $ HalfLine (pointAt x l) (l^.direction)+   (MinInfinity, Val x) -> Just $ HalfLine (pointAt x l) ((-1) *^ l^.direction)+   _                    -> Nothing  type instance IntersectionOf (HalfLine 2 r) (Line 2 r) = [ NoIntersection                                                          , Point 2 r

@@ -0,0 +1,94 @@+{-# LANGUAGE DeriveAnyClass  #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE TemplateHaskell  #-}+module Data.Geometry.HyperPlane where++import Control.DeepSeq+import Control.Lens+import Data.Geometry.Line+import Data.Geometry.Point+import Data.Geometry.Properties+import Data.Geometry.Transformation+import Data.Geometry.Vector+import GHC.Generics (Generic)+import GHC.TypeLits++--------------------------------------------------------------------------------++-- | Hyperplanes embedded in a \(d\) dimensional space.+data HyperPlane (d :: Nat) (r :: *) = HyperPlane { _inPlane   :: !(Point d r)+                                                 , _normalVec :: !(Vector d r)+                                                 } deriving Generic+makeLenses ''HyperPlane++type instance Dimension (HyperPlane d r) = d+type instance NumType   (HyperPlane d r) = r++deriving instance (Arity d, Show r)   => Show    (HyperPlane d r)+deriving instance (Arity d, Eq r)     => Eq      (HyperPlane d r)+deriving instance (NFData r, Arity d) => NFData  (HyperPlane d r)+deriving instance Arity d => Functor     (HyperPlane d)+deriving instance Arity d => Foldable    (HyperPlane d)+deriving instance Arity d => Traversable (HyperPlane d)++instance (Arity d, Arity (d + 1), Fractional r) => IsTransformable (HyperPlane d r) where+  transformBy t (HyperPlane p v) = HyperPlane (transformBy t p) (transformBy t v)++--------------------------------------------------------------------------------++-- | Test if a point lies on a hyperplane.+onHyperPlane                      :: (Num r, Eq r, Arity d)+                                  => Point d r -> HyperPlane d r -> Bool+q `onHyperPlane` (HyperPlane p n) = n `dot` (q .-. p) == 0+++-- -- | Compute a transformation that maps the last dimension (i.e. the d-axis) to+-- -- the normal vector of the plane. The origin of the coordinate system will+-- -- correspond to the inPlane point.+-- changeCoordinateSystem                  :: Floating r => HyperPlane 3 r -> Vector 3 r+--                                         -> Transformation 3 r+-- changeCoordinateSystem (HyperPlane p n) u = rotateTo (Vector3 u v w)+--                                        |.| translation (origin .-. p)+--   where+--     v = undefined+--     w = normalize n++-- toPlaneCoordinates :: HyperPlane d r ->++--------------------------------------------------------------------------------+-- * 3 Dimensional planes++type Plane = HyperPlane 3++pattern Plane     :: Point 3 r -> Vector 3 r -> Plane r+pattern Plane p n = HyperPlane p n++from3Points       :: Num r => Point 3 r -> Point 3 r -> Point 3 r -> HyperPlane 3 r+from3Points p q r = let u = q .-. p+                        v = r .-. p+                    in HyperPlane p (u `cross` v)+++type instance IntersectionOf (Line 3 r) (Plane r) = [NoIntersection, Point 3 r, Line 3 r]++instance (Eq r, Fractional r) => (Line 3 r) `IsIntersectableWith` (Plane r) where+  nonEmptyIntersection = defaultNonEmptyIntersection+  l@(Line p v) `intersect` (HyperPlane q n)+      | denum == 0 = if num == 0 then coRec l else coRec NoIntersection+      | otherwise  = coRec $ p .+^ (num / denum) *^ v+    where+      num   = (q .-. p) `dot` n+      denum = v `dot` n+    -- see https://en.wikipedia.org/wiki/Line%E2%80%93plane_intersection+++--------------------------------------------------------------------------------+-- * Supporting Planes++-- | Types for which we can compute a supporting hyperplane, i.e. a hyperplane+-- that contains the thing of type t.+class HasSupportingPlane t where+  supportingPlane :: t -> HyperPlane (Dimension t) (NumType t)++instance HasSupportingPlane (HyperPlane d r) where+  supportingPlane = id

@@ -24,7 +24,7 @@ import qualified Data.Foldable as F import           Data.Geometry.Properties import           Data.Range-import           Data.Semigroup+import           Data.Semigroup(Arg(..)) import qualified Data.Traversable as T import           Data.Vinyl import           Data.Vinyl.CoRec

@@ -1,15 +1,21 @@-{-|-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+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Reexports the functionality for reading and writing Ipe files.+--+-------------------------------------------------------------------------------- module Data.Geometry.Ipe( module Data.Geometry.Ipe.Types                         , module Data.Geometry.Ipe.Writer                         , module Data.Geometry.Ipe.Reader                         , module Data.Geometry.Ipe.IpeOut                         , module Data.Geometry.Ipe.FromIpe                         , module Data.Geometry.Ipe.Attributes+                        , module Data.Geometry.Ipe.Value+                        , module Data.Geometry.Ipe.Color                         ) where  import Data.Geometry.Ipe.Types@@ -18,3 +24,5 @@ import Data.Geometry.Ipe.IpeOut import Data.Geometry.Ipe.FromIpe import Data.Geometry.Ipe.Attributes+import Data.Geometry.Ipe.Value+import Data.Geometry.Ipe.Color(IpeColor(..))

@@ -3,11 +3,20 @@ {-# LANGUAGE UnicodeSyntax #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE UndecidableInstances #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.Ipe.Attributes+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Possible Attributes we can assign to items in an Ipe file+--+-------------------------------------------------------------------------------- module Data.Geometry.Ipe.Attributes where  import Control.Lens hiding (rmap, Const)-import Data.Colour.SRGB-import Data.Semigroup+import Data.Geometry.Ipe.Value import Data.Singletons import Data.Singletons.TH import Data.Text (Text)@@ -15,6 +24,8 @@ import Data.Vinyl.Functor import Data.Vinyl.TypeLevel import GHC.Exts+import Text.Read(lexP, step, parens, prec, (+++)+                , Lexeme(Ident), readPrec, readListPrec, readListPrecDefault)  -------------------------------------------------------------------------------- @@ -62,19 +73,44 @@                                    -- Labels in universe u to concrete types              (label :: u) = GAttr { _getAttr :: Maybe (Apply f label) } -deriving instance Show (Apply f label) => Show (Attr f label)-deriving instance Read (Apply f label) => Read (Attr f label)++ deriving instance Eq   (Apply f label) => Eq   (Attr f label) deriving instance Ord  (Apply f label) => Ord  (Attr f label)  makeLenses ''Attr +-- | Constructor for constructing an Attr given an actual value. pattern Attr   :: Apply f label -> Attr f label pattern Attr x = GAttr (Just x) +-- | An Attribute that is not set pattern NoAttr :: Attr f label pattern NoAttr = GAttr Nothing+{-# COMPLETE NoAttr, Attr #-} +instance Show (Apply f label) => Show (Attr f label) where+  showsPrec d NoAttr   = showParen (d > app_prec) $ showString "NoAttr"+    where app_prec = 10+  showsPrec d (Attr a) = showParen (d > up_prec) $+                           showString "Attr " . showsPrec (up_prec+1) a+    where up_prec  = 5++instance Read (Apply f label) => Read (Attr f label) where+  readPrec = parens $ (prec app_prec $ do+                                         Ident "NoAttr" <- lexP+                                         pure NoAttr)+                  +++ (prec up_prec $ do+                                         Ident "Attr" <- lexP+                                         a <- step readPrec+                                         pure $ Attr a)+    where+      app_prec = 10+      up_prec = 5+  readListPrec = readListPrecDefault+++ -- | Give pref. to the *RIGHT* instance Semigroup (Attr f l) where   _ <> b@(Attr _) = b@@ -84,17 +120,13 @@   mempty  = NoAttr   mappend = (<>) -newtype Attributes (f :: TyFun u * -> *) (ats :: [u]) =-  Attrs { _unAttrs :: Rec (Attr f) ats }--makeLenses ''Attributes----- type All' c i = RecAll (Attr (IpeObjectSymbolF i)) (IpeObjectAttrF i) c+newtype Attributes (f :: TyFun u * -> *) (ats :: [u]) = Attrs (Rec (Attr f) ats) --- deriving instance All' Show atsShow (Attributes f ats)+unAttrs :: Lens (Attributes f ats) (Attributes f' ats') (Rec (Attr f) ats) (Rec (Attr f') ats')+unAttrs = lens (\(Attrs r) -> r) (const Attrs)  deriving instance (RecAll (Attr f) ats Show) => Show (Attributes f ats)+-- deriving instance (RecAll (Attr f) ats Read) => Read (Attributes f ats)  instance (RecAll (Attr f) ats Eq)   => Eq   (Attributes f ats) where   (Attrs a) == (Attrs b) = and . recordToList@@ -192,20 +224,10 @@ --                              deriving (Show,Eq)  --- | 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,Functor,Foldable,Traversable) -instance IsString (IpeValue v) where-  fromString = Named . fromString- 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) -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

@@ -0,0 +1,110 @@+{-# LANGUAGE OverloadedStrings #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.Ipe.Color+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Data type for representing colors in ipe as well as the colors available in+-- the standard ipe stylesheet.+--+--------------------------------------------------------------------------------+module Data.Geometry.Ipe.Color where++import           Data.Colour.SRGB (RGB(..))+import           Data.Geometry.Ipe.Value+import           Data.Text+--------------------------------------------------------------------------------++newtype IpeColor r = IpeColor (IpeValue (RGB r))    deriving (Show,Read,Eq)++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)++-- | Creates a named color+named :: Text -> IpeColor r+named = IpeColor . Named++--------------------------------------------------------------------------------+-- * Basic Named colors++red :: IpeColor r+red = named "red"++green :: IpeColor r+green = named "green"++blue :: IpeColor r+blue = named "blue"++yellow :: IpeColor r+yellow = named "yellow"++orange :: IpeColor r+orange = named "orange"++gold :: IpeColor r+gold = named "gold"++purple :: IpeColor r+purple = named "purple"++gray :: IpeColor r+gray = named "gray"++brown :: IpeColor r+brown = named "brown"++navy :: IpeColor r+navy = named "navy"++pink :: IpeColor r+pink = named "pink"++seagreen :: IpeColor r+seagreen = named "seagreen"++turquoise :: IpeColor r+turquoise = named "turquoise"++violet :: IpeColor r+violet = named "violet"++darkblue :: IpeColor r+darkblue = named "darkblue"++darkcyan :: IpeColor r+darkcyan = named "darkcyan"++darkgray :: IpeColor r+darkgray = named "darkgray"++darkgreen :: IpeColor r+darkgreen = named "darkgreen"++darkmagenta :: IpeColor r+darkmagenta = named "darkmagenta"++darkorange :: IpeColor r+darkorange = named "darkorange"++darkred :: IpeColor r+darkred = named "darkred"++lightblue :: IpeColor r+lightblue = named "lightblue"++lightcyan :: IpeColor r+lightcyan = named "lightcyan"++lightgray :: IpeColor r+lightgray = named "lightgray"++lightgreen :: IpeColor r+lightgreen = named "lightgreen"++lightyellow :: IpeColor r+lightyellow = named "lightyellow"

@@ -1,4 +1,14 @@ {-# LANGUAGE OverloadedStrings #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.Ipe.FromIpe+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Functions that help reading geometric values from ipe images.+--+-------------------------------------------------------------------------------- module Data.Geometry.Ipe.FromIpe where  import           Control.Lens hiding (Simple)@@ -9,21 +19,23 @@ import qualified Data.Geometry.PolyLine as PolyLine import           Data.Geometry.Polygon import           Data.Geometry.Properties-import qualified Data.Seq2 as S2-import qualified Data.List.NonEmpty as NonEmpty+import           Data.Geometry.Triangle+import qualified Data.LSeq as LSeq+import           Data.List.NonEmpty (NonEmpty(..))  -------------------------------------------------------------------------------- -- $setup+-- >>> :set -XOverloadedStrings+-- >>> import Data.Geometry.Ipe.Attributes+-- >>> import Data.Geometry.Ipe.Color(IpeColor(..))+-- >>> import Data.Geometry.Point -- >>> :{--- import           Data.Geometry.Ipe.Attributes--- -- let testPath :: Path Int---     testPath = Path . S2.l1Singleton  . PolyLineSegment . PolyLine.fromPoints . map ext+--     testPath = Path . fromSingleton  . PolyLineSegment+--              . PolyLine.fromPoints . map ext --              $ [ origin, point2 10 10, point2 200 100 ]--- --     testPathAttrs :: IpeAttributes Path Int---     testPathAttrs = attr SStroke (IpeColor (Named "red"))-+--     testPathAttrs = attr SStroke (IpeColor "red") --     testObject :: IpeObject Int --     testObject = IpePath (testPath :+ testPathAttrs) -- :}@@ -42,11 +54,11 @@ -- | 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] :+ ())})+-- Just (PolyLine {_points = LSeq (fromList [Point2 [0,0] :+ (),Point2 [10,10] :+ (),Point2 [200,100] :+ ()])}) _asPolyLine :: Prism' (Path r) (PolyLine.PolyLine 2 () r) _asPolyLine = prism' poly2path path2poly   where-    poly2path = Path . S2.l1Singleton  . PolyLineSegment+    poly2path = Path . fromSingleton  . PolyLineSegment     path2poly = preview (pathSegments.traverse._PolyLineSegment)     -- TODO: Check that the path actually is a polyline, rather     -- than ignoring everything that does not fit@@ -57,6 +69,19 @@   where     path2poly p = pathToPolygon p >>= either pure (const Nothing) +-- | Convert to a triangle+_asTriangle :: Prism' (Path r) (Triangle 2 () r)+_asTriangle = prism' triToPath path2tri+  where+    triToPath (Triangle p q r) = polygonToPath . fromPoints . map (&extra .~ ()) $ [p,q,r]+    path2tri p = case p^..pathSegments.traverse._PolygonPath of+                    []   -> Nothing+                    [pg] -> case polygonVertices pg of+                              (a :| [b,c]) -> Just $ Triangle a b c+                              _            -> Nothing+                    _    -> Nothing++ -- | Convert to a multipolygon _asMultiPolygon :: Prism' (Path r) (MultiPolygon () r) _asMultiPolygon = prism' polygonToPath path2poly@@ -64,9 +89,9 @@     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+polygonToPath pg@(SimplePolygon _) = Path . fromSingleton . PolygonPath $ pg+polygonToPath (MultiPolygon vs hs) = Path . LSeq.fromNonEmpty . fmap PolygonPath+                                   $ (SimplePolygon vs) :| hs   pathToPolygon   :: Path r -> Maybe (Either (SimplePolygon () r) (MultiPolygon () r))@@ -77,12 +102,12 @@   --- | use the first prism to select the ipe object to depicle with, and the second+-- | 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}}})+-- Just (PolyLine {_points = LSeq (fromList [Point2 [0,0] :+ (),Point2 [10,10] :+ (),Point2 [200,100] :+ ()])} :+ Attrs {NoAttr, NoAttr, NoAttr, NoAttr, Attr IpeColor (Named "red"), NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr}) _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@@ -144,3 +169,5 @@                => FilePath -> IO [g :+ IpeAttributes (DefaultFromIpe g) r] readAllFrom fp = readAll <$> readSinglePageFile fp +fromSingleton :: a -> LSeq.LSeq 1 a+fromSingleton = LSeq.fromNonEmpty . (:| [])

@@ -1,7 +1,22 @@ {-# LANGUAGE OverloadedStrings #-} {-# LANGUAGE ScopedTypeVariables #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.Ipe.IpeOut+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Functions that help drawing geometric values in ipe. An "IpeOut" is+-- essenitally a function that converts a geometric type g into an IpeObject.+--+-- We also proivde a "HasDefaultIpeOut" typeclass that defines a default+-- conversion function from a geometry type g to an ipe type.+--+-------------------------------------------------------------------------------- module Data.Geometry.Ipe.IpeOut where + import           Control.Lens hiding (Simple) import           Data.Bifunctor import           Data.Ext@@ -9,72 +24,83 @@ import           Data.Geometry.Boundary import           Data.Geometry.Box import           Data.Geometry.Ipe.Attributes+import           Data.Geometry.Ipe.Color (IpeColor(..)) 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+import           Data.Geometry.PolyLine(PolyLine,fromLineSegment) import           Data.Geometry.Polygon import           Data.Geometry.Polygon.Convex import           Data.Geometry.Properties import           Data.Geometry.Transformation-import           Data.Maybe (fromMaybe)+import qualified Data.LSeq as LSeq+import           Data.List.NonEmpty (NonEmpty(..)) import           Data.Proxy-import           Data.Semigroup-import qualified Data.Seq2 as S2 import           Data.Text (Text)+import           Data.Maybe (fromMaybe)+import           Linear.Affine ((.+^)) import           Data.Vinyl.CoRec  -------------------------------------------------------------------------------- --- | An IpeOut is essentially a funciton to convert a geometry object of type--- 'g' into an ipe object of type 'i'.-newtype IpeOut g i = IpeOut { asIpe :: g -> i } deriving (Functor)----- | Given an geometry object, and a record with its attributes, construct an ipe--- Object representing it using the default conversion.-asIpeObject :: (HasDefaultIpeOut g, DefaultIpeOut g ~ i, NumType g ~ r)-            => g -> IpeAttributes i r -> IpeObject r-asIpeObject = asIpeObjectWith defaultIpeOut---- | asIpeObject with its arguments flipped. Convenient if you don't want to--- map asIpeObject over a list or so.-asIpeObject' :: (HasDefaultIpeOut g, DefaultIpeOut g ~ i, NumType g ~ r)-             => IpeAttributes i r -> g -> IpeObject r-asIpeObject' = flip asIpeObject+-- $setup+-- >>> :set -XOverloadedStrings+-- >>> :{+-- let myPolygon = fromPoints . map ext $ [origin, Point2 10 10, Point2 100 200]+-- :} +--------------------------------------------------------------------------------+-- * The IpeOut type and the default combinator to use it --- -- | Given a IpeOut that specifies how to convert a geometry object into an--- 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-asIpeObjectWith io g ats = asIpe (ipeObject io ats) g+type IpeOut g i r = g -> IpeObject' i r  --- | Create an ipe group without group attributes-asIpeGroup :: [IpeObject r] -> IpeObject r-asIpeGroup = flip asIpeGroup' mempty---- | Creates a group out of ipe-asIpeGroup'        :: [IpeObject r] -> IpeAttributes Group r -> IpeObject r-asIpeGroup' gs ats = IpeGroup $ Group gs :+ ats+-- | Add attributes to an IpeObject'+(!)       :: IpeObject' i r -> IpeAttributes i r -> IpeObject' i r+(!) i ats = i&extra %~ (<> ats) ---------------------------------------------------------------------------------+-- | Render an ipe object+--+--+-- >>> :{+--   iO $ defIO myPolygon ! attr SFill (IpeColor "blue")+--                        ! attr SLayer "alpha"+--                        ! attr SLayer "beta"+-- :}+-- IpePath (Path {_pathSegments = LSeq (fromList [PolygonPath SimplePolygon CSeq [Point2 [0,0] :+ (),Point2 [10,10] :+ (),Point2 [100,200] :+ ()]])} :+ Attrs {Attr LayerName {_layerName = "beta"}, NoAttr, NoAttr, NoAttr, NoAttr, Attr IpeColor (Named "blue"), NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr})+--+-- >>> :{+--   iO $ ipeGroup [ iO $ ipePolygon myPolygon ! attr SFill (IpeColor "red")+--                 ] ! attr SLayer "alpha"+-- :}+-- IpeGroup (Group [IpePath (Path {_pathSegments = LSeq (fromList [PolygonPath SimplePolygon CSeq [Point2 [0,0] :+ (),Point2 [10,10] :+ (),Point2 [100,200] :+ ()]])} :+ Attrs {NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, Attr IpeColor (Named "red"), NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr})] :+ Attrs {Attr LayerName {_layerName = "alpha"}, NoAttr, NoAttr, NoAttr, NoAttr})+--+iO :: ToObject i => IpeObject' i r -> IpeObject r+iO = mkIpeObject --- | Helper to construct an IpeOut g IpeObject , if we already know how to--- construct a specific Ipe type.-ipeObject        :: (ToObject i, NumType g ~ r)-                   => IpeOut g (IpeObject' i r) -> IpeAttributes i r -> IpeOut g (IpeObject r)-ipeObject io ats = IpeOut $ \g -> let (i :+ ats') = asIpe io g-                                    in ipeObject' i (ats' <> ats)+-- | Render to an ipe object using the defIO IpeOut+--+--+-- >>> :{+--   iO'' myPolygon $  attr SFill (IpeColor "red")+--                  <> attr SLayer "alpha"+--                  <> attr SLayer "beta"+-- :}+-- IpePath (Path {_pathSegments = LSeq (fromList [PolygonPath SimplePolygon CSeq [Point2 [0,0] :+ (),Point2 [10,10] :+ (),Point2 [100,200] :+ ()]])} :+ Attrs {Attr LayerName {_layerName = "beta"}, NoAttr, NoAttr, NoAttr, NoAttr, Attr IpeColor (Named "red"), NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr})+--+-- >>> iO'' [ myPolygon , myPolygon ] $ attr SLayer "alpha"+-- IpeGroup (Group [IpePath (Path {_pathSegments = LSeq (fromList [PolygonPath SimplePolygon CSeq [Point2 [0,0] :+ (),Point2 [10,10] :+ (),Point2 [100,200] :+ ()]])} :+ Attrs {NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr}),IpePath (Path {_pathSegments = LSeq (fromList [PolygonPath SimplePolygon CSeq [Point2 [0,0] :+ (),Point2 [10,10] :+ (),Point2 [100,200] :+ ()]])} :+ Attrs {NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr, NoAttr})] :+ Attrs {Attr LayerName {_layerName = "alpha"}, NoAttr, NoAttr, NoAttr, NoAttr})+iO''       :: ( HasDefaultIpeOut g, NumType g ~ r+             , DefaultIpeOut g ~ i, ToObject i+             ) => g -> IpeAttributes i r+           -> IpeObject r+iO'' g ats = iO $ defIO g ! ats --- | Construct an ipe object from the core of an Ext-coreOut    :: IpeOut g i -> IpeOut (g :+ a) i-coreOut io = IpeOut $ asIpe io . (^.core)+-- | generate an ipe object without any specific attributes+iO' :: HasDefaultIpeOut g => g -> IpeObject (NumType g)+iO' = iO . defIO  -------------------------------------------------------------------------------- -- * Default Conversions@@ -84,157 +110,105 @@ class ToObject (DefaultIpeOut g) => HasDefaultIpeOut g where   type DefaultIpeOut g :: * -> * -  defaultIpeOut :: IpeOut g (IpeObject' (DefaultIpeOut g) (NumType g))+  defIO :: IpeOut g (DefaultIpeOut g) (NumType g) -  -- defaultIpeObject :: RecApplicative (AttributesOf (DefaultIpeOut g))-  --                  => IpeOut g (IpeObject (NumType g))-  -- defaultIpeObject = IpeOut $ flip asIpeObject mempty+instance (HasDefaultIpeOut g, a ~ IpeAttributes (DefaultIpeOut g) (NumType g))+        => HasDefaultIpeOut (g :+ a) where+  type DefaultIpeOut (g :+ a) = DefaultIpeOut g+  defIO (g :+ ats) = defIO g ! ats --- instance HasDefaultIpeOut g => HasDefaultIpeOut [g] where---   type DefaultIpeOut [g] = Group---   defaultIpeOut = IpeOut $ asIpeGroup . map (asIpeObject' mempty)+instance HasDefaultIpeOut a => HasDefaultIpeOut [a] where+  type DefaultIpeOut [a] = Group+  defIO = ipeGroup . map (iO .  defIO)  instance HasDefaultIpeOut (Point 2 r) where   type DefaultIpeOut (Point 2 r) = IpeSymbol-  defaultIpeOut = ipeDiskMark+  defIO = ipeDiskMark  instance HasDefaultIpeOut (LineSegment 2 p r) where   type DefaultIpeOut (LineSegment 2 p r) = Path-  defaultIpeOut = ipeLineSegment--instance Floating r => HasDefaultIpeOut (Disk p r) where-  type DefaultIpeOut (Disk p r) = Path-  defaultIpeOut = ipeDisk+  defIO = ipeLineSegment  instance HasDefaultIpeOut (PolyLine 2 p r) where   type DefaultIpeOut (PolyLine 2 p r) = Path-  defaultIpeOut = noAttrs ipePolyLine+  defIO = ipePolyLine +instance (Fractional r, Ord r) => HasDefaultIpeOut (Line 2 r) where+  type DefaultIpeOut (Line 2 r) = Path+  defIO = ipeLineSegment . toSeg+    where+      b :: Rectangle () r+      b = box (ext $ Point2 (-200) (-200)) (ext $ Point2 1200 1200)+      naive (Line p v) = ClosedLineSegment (ext p) (ext $ p .+^ v)+      toSeg l = fromMaybe (naive l) . asA (Proxy :: Proxy (LineSegment 2 () r))+              $ l `intersect` b+ 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")+  defIO = ipePolygon  instance HasDefaultIpeOut (SomePolygon p r) where   type DefaultIpeOut (SomePolygon p r) = Path-  defaultIpeOut = IpeOut $ either (asIpe defaultIpeOut) (asIpe defaultIpeOut)+  defIO = either defIO defIO  instance HasDefaultIpeOut (ConvexPolygon p r) where   type DefaultIpeOut (ConvexPolygon p r) = Path-  defaultIpeOut = IpeOut $ asIpe defaultIpeOut . view simplePolygon+  defIO = defIO . view simplePolygon ++instance Floating r => HasDefaultIpeOut (Disk p r) where+  type DefaultIpeOut (Disk p r) = Path+  defIO = ipeDisk+ -------------------------------------------------------------------------------- -- * Point Converters -ipeMark   :: Text -> IpeOut (Point 2 r) (IpeObject' IpeSymbol r)-ipeMark n = noAttrs . IpeOut $ flip Symbol n+ipeMark     :: Text -> IpeOut (Point 2 r) IpeSymbol r+ipeMark n p = Symbol p n :+ mempty -ipeDiskMark :: IpeOut (Point 2 r) (IpeObject' IpeSymbol r)+ipeDiskMark :: IpeOut (Point 2 r) IpeSymbol r ipeDiskMark = ipeMark "mark/disk(sx)"  ----------------------------------------------------------------------------------noAttrs :: Monoid extra => IpeOut g core -> IpeOut g (core :+ extra)-noAttrs = addAttributes mempty--addAttributes :: extra -> IpeOut g core -> IpeOut g (core :+ extra)-addAttributes ats io = IpeOut $ \g -> asIpe io g :+ ats----- | Default size of the cliping rectangle used to clip lines. This is--- Rectangle is large enough to cover the normal page size in ipe.-defaultClipRectangle :: (Num r, Ord r) => Rectangle () r-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 :: (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'--ipeLineSegment' :: IpeOut (LineSegment 2 p r) (PathSegment r)-ipeLineSegment' = IpeOut $ PolyLineSegment . fromLineSegment . first (const ())---ipePolyLine :: IpeOut (PolyLine 2 p r) (Path r)-ipePolyLine = fromPathSegment ipePolyLine'+-- * Path Converters -ipePolyLine' :: IpeOut (PolyLine 2 a r) (PathSegment r)-ipePolyLine' = IpeOut $ PolyLineSegment . first (const ())+ipeLineSegment   :: IpeOut (LineSegment 2 p r) Path r+ipeLineSegment s = (path . pathSegment $ s) :+ mempty -ipeDisk :: Floating r => IpeOut (Disk p r) (IpeObject' Path r)-ipeDisk = noAttrs . IpeOut $ asIpe ipeCircle . Boundary+ipePolyLine   :: IpeOut (PolyLine 2 p r) Path r+ipePolyLine p = (path . PolyLineSegment . first (const ()) $ p) :+ mempty -ipeCircle :: Floating r => IpeOut (Circle p r) (Path r)-ipeCircle = fromPathSegment ipeCircle'+ipeDisk   :: Floating r => IpeOut (Disk p r) Path r+ipeDisk d = ipeCircle (Boundary d) ! attr SFill (IpeColor "0.722 0.145 0.137") -ipeCircle' :: Floating r => IpeOut (Circle p r) (PathSegment r)-ipeCircle' = IpeOut circle''-  where-    circle'' (Circle (c :+ _) r) = EllipseSegment m+ipeCircle                     :: Floating r => IpeOut (Circle p r) Path r+ipeCircle (Circle (c :+ _) r) = (path $ EllipseSegment m) :+ mempty       where         m = translation (toVec c) |.| uniformScaling (sqrt r) ^. transformationMatrix         -- m is the matrix s.t. if we apply m to the unit circle centered at the origin, we         -- get the input circle. ---- | Helper to construct a IpeOut g Path, for when we already have an IpeOut g PathSegment-fromPathSegment    :: IpeOut g (PathSegment r) -> IpeOut g (Path r)-fromPathSegment io = IpeOut $ Path . S2.l1Singleton . asIpe io---ipePolygon :: IpeOut (Polygon t p r) (Path r)-ipePolygon = IpeOut $ io . first (const ())-  where-    io                       :: forall t r. Polygon t () r -> Path r-    io pg@(SimplePolygon _)  = pg^.re _asSimplePolygon-    io pg@(MultiPolygon _ _) = pg^.re _asMultiPolygon------------ ls = (ClosedLineSegment (ext origin) (ext (point2 1 1)))----- testzz :: IpeObject Integer--- testzz = asIpeObjectWith ipeLineSegment ls $ mempty <> attr SStroke (IpeColor "red")---+-- | Helper to construct a path from a singleton item+path :: PathSegment r -> Path r+path = Path . LSeq.fromNonEmpty . (:| []) --- test' :: Attributes (PathAttrElfSym1 Integer) (AttributesOf (Path Integer) (PathAttrElfSym1 Integer))--- -- test' :: RecApplicative (AttributesOf (Path Integer) (IpeObjectSymbolF (Path Integer)))--- --       => IpeAttributes (Path Integer)--- test' = mempty+pathSegment :: LineSegment 2 p r -> PathSegment r+pathSegment = PolyLineSegment . fromLineSegment . first (const ()) +-- | Draw a polygon+ipePolygon                          :: IpeOut (Polygon t p r) Path r+ipePolygon (first (const ()) -> pg) = case pg of+               (SimplePolygon _)  -> pg^.re _asSimplePolygon :+ mempty+               (MultiPolygon _ _) -> pg^.re _asMultiPolygon  :+ mempty  +--------------------------------------------------------------------------------+-- * Group Converters --- -- test' :: IpeObject Integer ('IpePath '[])--- test' = asIpeObject' ls emptyPathAttributes+ipeGroup    :: IpeOut [IpeObject r] Group r+ipeGroup xs = Group xs :+ mempty    --- emptyPathAttributes :: Rec (PathAttribute r) '[]--- emptyPathAttributes = RNil+--------------------------------------------------------------------------------

@@ -12,7 +12,6 @@ import           Data.Geometry.Transformation import           Data.Geometry.Vector import           Data.Ratio-import           Data.Semigroup import           Data.Text (Text) import qualified Data.Text as T import           Text.Parsec.Error (messageString, errorMessages)

@@ -37,6 +37,8 @@ import           Data.Geometry.Ipe.ParserPrimitives (pInteger, pWhiteSpace) import           Data.Geometry.Ipe.PathParser import           Data.Geometry.Ipe.Types+import           Data.Geometry.Ipe.Value+import           Data.Geometry.Ipe.Color(IpeColor(..)) import           Data.Geometry.Point import           Data.Geometry.PolyLine import qualified Data.Geometry.Polygon as Polygon@@ -45,8 +47,7 @@ import qualified Data.List.NonEmpty as NE import           Data.Maybe (fromMaybe, mapMaybe) import           Data.Proxy-import qualified Data.Seq2 as S2-import           Data.Semigroup ((<>))+import qualified Data.LSeq as LSeq import           Data.Singletons import           Data.Text (Text) import qualified Data.Text as T@@ -210,7 +211,7 @@ dropRepeats = map head . L.group  instance (Coordinate r, Eq r) => IpeReadText (Path r) where-  ipeReadText = fmap (Path . S2.viewL1FromNonEmpty) . ipeReadText+  ipeReadText = fmap (Path . LSeq.fromNonEmpty) . ipeReadText  -------------------------------------------------------------------------------- -- Reading attributes

@@ -3,6 +3,16 @@ {-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE QuasiQuotes #-} {-# LANGUAGE OverloadedStrings #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.Ipe.Types+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Data type modeling the various elements in Ipe files.+--+-------------------------------------------------------------------------------- module Data.Geometry.Ipe.Types where  @@ -22,6 +32,7 @@ import           Data.Singletons.TH(genDefunSymbols)  import           Data.Geometry.Ipe.Literal+import           Data.Geometry.Ipe.Color import qualified Data.Geometry.Ipe.Attributes as AT import           Data.Geometry.Ipe.Attributes hiding (Matrix) import           Data.Text(Text)@@ -31,7 +42,7 @@   import qualified Data.List.NonEmpty as NE-import qualified Data.Seq2     as S2+import qualified Data.LSeq          as LSeq  -------------------------------------------------------------------------------- @@ -125,7 +136,7 @@   -- | A path is a non-empty sequence of PathSegments.-newtype Path r = Path { _pathSegments :: S2.ViewL1 (PathSegment r) }+newtype Path r = Path { _pathSegments :: LSeq.LSeq 1 (PathSegment r) }                  deriving (Show,Eq) makeLenses ''Path @@ -197,8 +208,7 @@   -- | A group is essentially a list of IpeObjects.-newtype Group r = Group { _groupItems :: [IpeObject r] }-                  deriving (Show,Eq)+newtype Group r = Group [IpeObject r] deriving (Show,Eq)  type instance NumType   (Group r) = r type instance Dimension (Group r) = 2@@ -238,23 +248,26 @@   deriving instance (Show r) => Show (IpeObject r)+-- deriving instance (Read r) => Read (IpeObject r) deriving instance (Eq r)   => Eq   (IpeObject r)  type instance NumType   (IpeObject r) = r type instance Dimension (IpeObject r) = 2  makePrisms ''IpeObject-makeLenses ''Group +groupItems :: Lens (Group r) (Group s) [IpeObject r] [IpeObject s]+groupItems = lens (\(Group xs) -> xs) (const Group)+ class ToObject i where-  ipeObject' :: i r -> IpeAttributes i r -> IpeObject r+  mkIpeObject :: IpeObject' i r -> IpeObject r -instance ToObject Group      where ipeObject' g a = IpeGroup     (g :+ a)-instance ToObject Image      where ipeObject' p a = IpeImage     (p :+ a)-instance ToObject TextLabel  where ipeObject' p a = IpeTextLabel (p :+ a)-instance ToObject MiniPage   where ipeObject' p a = IpeMiniPage  (p :+ a)-instance ToObject IpeSymbol  where ipeObject' s a = IpeUse       (s :+ a)-instance ToObject Path       where ipeObject' p a = IpePath      (p :+ a)+instance ToObject Group      where mkIpeObject = IpeGroup+instance ToObject Image      where mkIpeObject = IpeImage+instance ToObject TextLabel  where mkIpeObject = IpeTextLabel+instance ToObject MiniPage   where mkIpeObject = IpeMiniPage+instance ToObject IpeSymbol  where mkIpeObject = IpeUse+instance ToObject Path       where mkIpeObject = IpePath  instance Fractional r => IsTransformable (IpeObject r) where   transformBy t (IpeGroup i)     = IpeGroup     $ i&core %~ transformBy t@@ -264,8 +277,9 @@   transformBy t (IpeUse i)       = IpeUse       $ i&core %~ transformBy t   transformBy t (IpePath i)      = IpePath      $ i&core %~ transformBy t --+-- | Shorthand for constructing ipeObjects+ipeObject'     :: ToObject i => i r -> IpeAttributes i r -> IpeObject r+ipeObject' i a = mkIpeObject $ i :+ a  commonAttributes :: Lens' (IpeObject r) (Attributes (AttrMapSym1 r) CommonAttributes) commonAttributes = lens (Attrs . g) (\x (Attrs a) -> s x a)

@@ -0,0 +1,24 @@+{-# LANGUAGE OverloadedStrings #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.Ipe.Value+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Data type for representing values in ipe.+--+--------------------------------------------------------------------------------+module Data.Geometry.Ipe.Value where++import GHC.Exts+import Data.Text++--------------------------------------------------------------------------------++-- | Many types either consist of a symbolc value, or a value of type v+data IpeValue v = Named Text | Valued v+  deriving (Show,Read,Eq,Ord,Functor,Foldable,Traversable)++instance IsString (IpeValue v) where+  fromString = Named . fromString

@@ -1,9 +1,18 @@ {-# LANGUAGE OverloadedStrings #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE UndecidableInstances #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.Ipe.Writer+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+-- Description :  Converting data types into IpeTypes+--+-------------------------------------------------------------------------------- module Data.Geometry.Ipe.Writer where -import           Control.Lens ((^.),(^..),(.~),(&), to)+import           Control.Lens ((^.), (^..), (.~), (&)) import qualified Data.ByteString as B import qualified Data.ByteString.Char8 as C import           Data.Colour.SRGB (RGB(..))@@ -13,18 +22,20 @@ import           Data.Geometry.Box import           Data.Geometry.Ipe.Attributes import qualified Data.Geometry.Ipe.Attributes as IA+import           Data.Geometry.Ipe.Color (IpeColor(..)) import           Data.Geometry.Ipe.Types+import           Data.Geometry.Ipe.Value import           Data.Geometry.LineSegment import           Data.Geometry.Point import           Data.Geometry.PolyLine import           Data.Geometry.Polygon (Polygon, outerBoundary, holeList, asSimplePolygon) import qualified Data.Geometry.Transformation as GT import           Data.Geometry.Vector+import qualified Data.LSeq as LSeq+import           Data.List.NonEmpty (NonEmpty(..)) import           Data.Maybe (catMaybes, mapMaybe, fromMaybe) import           Data.Proxy import           Data.Ratio-import           Data.Semigroup-import qualified Data.Seq2 as S2 import           Data.Singletons import           Data.Text (Text) import qualified Data.Text as T@@ -35,6 +46,7 @@ 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,@@ -90,6 +102,11 @@ class IpeWrite t where   ipeWrite :: t -> Maybe (Node Text Text) +instance IpeWrite t => IpeWrite [t] where+  ipeWrite gs = case mapMaybe ipeWrite gs of+                  [] -> Nothing+                  ns -> (Just $ Element "group" [] ns)+ instance (IpeWrite l, IpeWrite r) => IpeWrite (Either l r) where   ipeWrite = either ipeWrite ipeWrite @@ -252,6 +269,7 @@ instance IpeWriteText r => IpeWriteText (PolyLine 2 () r) where   ipeWriteText pl = case pl^..points.traverse.core of     (p : rest) -> unlines' . map ipeWriteText $ MoveTo p : map LineTo rest+    _          -> error "ipeWriteText. absurd. no vertices polyline"     -- the polyline type guarantees that there is at least one point  instance IpeWriteText r => IpeWriteText (Polygon t () r) where@@ -269,6 +287,7 @@   ipeWriteText (PolyLineSegment p) = ipeWriteText p   ipeWriteText (PolygonPath     p) = ipeWriteText p   ipeWriteText (EllipseSegment  m) = ipeWriteText $ Ellipse m+  ipeWriteText _                   = error "ipeWriteText: PathSegment, not implemented yet."  instance IpeWriteText r => IpeWrite (Path r) where   ipeWrite p = (\t -> Element "path" [] [Text t]) <$> ipeWriteText p@@ -277,9 +296,7 @@   instance (IpeWriteText r) => IpeWrite (Group r) where-  ipeWrite (Group gs) = case mapMaybe ipeWrite gs of-                          [] -> Nothing-                          ns -> (Just $ Element "group" [] ns)+  ipeWrite (Group gs) = ipeWrite gs   instance ( AllSatisfy IpeAttrName rs@@ -389,7 +406,7 @@       -- TODO: Do something with the p's  fromPolyLine :: PolyLine 2 () r -> Path r-fromPolyLine = Path . S2.l1Singleton . PolyLineSegment+fromPolyLine = Path . LSeq.fromNonEmpty . (:| []) . PolyLineSegment   instance (IpeWriteText r) => IpeWrite (LineSegment 2 p r) where

@@ -14,8 +14,8 @@ import qualified Data.List.NonEmpty as NonEmpty import           Data.Maybe (fromJust) import           Data.Proxy-import           Data.Seq (LSeq(..), ViewL(..))-import qualified Data.Seq as Seq+import           Data.LSeq (LSeq(..), pattern (:<|))+import qualified Data.LSeq as LSeq import           Data.Util import qualified Data.Vector.Fixed as FV import           GHC.TypeLits@@ -77,7 +77,7 @@  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+buildKDTree' = KDT . addBoxes . build (Coord 1) . toPointSet . LSeq.fromNonEmpty   where     -- compute one tree with bounding boxes, then merge them together     addBoxes t = let bbt = foldUpData (\l _ r -> boundingBoxList' [l,r])                                       (boundingBox . (^.core)) t@@ -94,7 +94,7 @@            => LSeq n (Point d r :+ p) -> PointSet (LSeq n) d p r toPointSet = FV.imap sort . FV.replicate   where-    sort i = Seq.unstableSortBy (compareOn $ 1 + i)+    sort i = LSeq.unstableSortBy (compareOn $ 1 + i)   compareOn       :: (Ord r, Arity d)@@ -166,13 +166,13 @@     -- i = traceShow (c,j) j      m = let xs = fromJust $ pts^?element' (i-1)-        in xs `Seq.index` (F.length xs `div` 2)+        in xs `LSeq.index` (F.length xs `div` 2)      -- Since the input seq has >= 2 elems, F.length xs / 2 >= 1. It follows     -- that the both sets thus have at least one elemnt.     -- f :: LSeq 2 _ -> (LSeq 1 _, LSeq 1 _)-    f = bimap Seq.promise Seq.promise-      . Seq.partition (\p -> compareOn i p m == LT)+    f = bimap LSeq.promise LSeq.promise+      . LSeq.partition (\p -> compareOn i p m == LT)      (l,r) = unzip' . fmap f $ pts @@ -180,8 +180,9 @@     unzip' = bimap vectorFromListUnsafe vectorFromListUnsafe . unzip . F.toList  -asSingleton   :: (1 <= 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 $ unsafeCoerce v-                  p :< _           -> Left p -- only one element+asSingleton v = case v^.element (C :: C 0) of+                  (p :<| s) | null s -> Left p -- only one lement+                  _                  -> Right $ unsafeCoerce v

@@ -2,10 +2,20 @@ {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE UnicodeSyntax #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.Line+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- \(d\)-dimensional lines.+--+-------------------------------------------------------------------------------- module Data.Geometry.Line( module Data.Geometry.Line.Internal                          ) where -import           Control.Lens ((^.), re, bimap)+import           Control.Lens import           Data.Ext import           Data.Geometry.Boundary import           Data.Geometry.Box@@ -44,12 +54,14 @@   nonEmptyIntersection = defaultNonEmptyIntersection    line' `intersect` (Boundary rect)  = case asA' segP of-      [sl'] -> coRec . bimap id unVal $ sl'^.re _SubLine-      []    -> case nub' . map (fmap unVal) $ asA' pointP of+      [sl'] -> case fromUnbounded sl' of+        Nothing   -> error "intersect: line x boundary rect; unbounded line? absurd"+        Just sl'' -> coRec $ sl''^.re _SubLine+      []    -> case nub' $ asA' pointP of         [p]   -> coRec p         [p,q] -> coRec (p,q)         _     -> coRec NoIntersection-      _     -> error "intersect; ine x boundary rect; absurd"+      _     -> error "intersect; line x boundary rect; absurd"     where       (t,r,b,l) = sides' rect       ints = map (\s -> sl `intersect` toSL s) [t,r,b,l]@@ -58,18 +70,16 @@        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"+      toSL  :: LineSegment 2 p r -> SubLine 2 () (UnBounded r) r+      toSL s = s^._SubLine.re _unBounded.to dropExtra -      asA'    :: (t ∈ IntersectionOf (SubLine 2 () (UnBounded r))-                                     (SubLine 2 () (UnBounded r)))+      asA'    :: (t ∈ IntersectionOf (SubLine 2 () (UnBounded r) r)+                                     (SubLine 2 () (UnBounded r) r))               => proxy t -> [t]       asA' px = mapMaybe (asA px) ints -      segP   = Proxy :: Proxy (SubLine 2 () (UnBounded r))-      pointP = Proxy :: Proxy (Point 2      (UnBounded r))+      segP   = Proxy :: Proxy (SubLine 2 () (UnBounded r) r)+      pointP = Proxy :: Proxy (Point 2 r)   type instance IntersectionOf (Line 2 r) (Rectangle p r) =

@@ -2,6 +2,16 @@ {-# LANGUAGE DeriveAnyClass  #-} {-# LANGUAGE ScopedTypeVariables  #-} {-# LANGUAGE UndecidableInstances #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.Line.Internal+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- \(d\)-dimensional lines.+--+-------------------------------------------------------------------------------- module Data.Geometry.Line.Internal where  import           Control.DeepSeq@@ -29,13 +39,22 @@  instance (Show r, Arity d) => Show (Line d r) where   show (Line p v) = concat [ "Line (", show p, ") (", show v, ")" ]-deriving instance (Eq r,   Arity d)   => Eq            (Line d r)++-- -- TODO:+-- instance (Read r, Arity d)   => Read (Line d r) where++++ deriving instance (NFData r, Arity d) => NFData        (Line d r) deriving instance Arity d             => Functor       (Line d) deriving instance Arity d             => F.Foldable    (Line d) deriving instance Arity d             => T.Traversable (Line d) +instance (Arity d, Eq r, Fractional r) => Eq (Line d r) where+  l@(Line p _) == m = l `isParallelTo` m && p `onLine` m + type instance Dimension (Line d r) = d type instance NumType   (Line d r) = r @@ -92,6 +111,29 @@ -- | 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+++-- | Get the point at the given position along line, where 0 corresponds to the+-- anchorPoint of the line, and 1 to the point anchorPoint .+^ directionVector+pointAt              :: (Num r, Arity d) => r -> Line d r -> Point d r+pointAt a (Line p v) = p .+^ (a *^ v)+++-- | Given point p and a line (Line q v), Get the scalar lambda s.t.+-- p = q + lambda v. If p does not lie on the line this returns a Nothing.+toOffset              :: (Eq r, Fractional r, Arity d) => Point d r -> Line d r -> Maybe r+toOffset p (Line q v) = scalarMultiple (p .-. q) v+++-- | Given point p *on* a line (Line q v), Get the scalar lambda s.t.+-- p = q + lambda v. (So this is an unsafe version of 'toOffset')+--+-- pre: the input point p lies on the line l.+toOffset'             :: (Eq r, Fractional r, Arity d) => Point d r -> Line d r -> r+toOffset' p = fromJust' . toOffset p+  where+    fromJust' (Just x) = x+    fromJust' _        = error "toOffset: Nothing"   -- | The intersection of two lines is either: NoIntersection, a point or a line.

@@ -1,11 +1,15 @@ {-# 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+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Line segment data type and some basic functions on line segments+--+-------------------------------------------------------------------------------- module Data.Geometry.LineSegment( LineSegment                                 , pattern LineSegment                                 , pattern LineSegment'@@ -37,8 +41,6 @@ import           Data.Geometry.Transformation import           Data.Geometry.Vector import           Data.Ord (comparing)-import           Data.Semigroup-import           Data.UnBounded import           Data.Vinyl import           Data.Vinyl.CoRec import           GHC.TypeLits@@ -92,22 +94,21 @@   end = unLineSeg.end  -_SubLine :: (Fractional r, Eq r, Arity d) => Iso' (LineSegment d p r) (SubLine d p r)+_SubLine :: (Num r, Arity d) => Iso' (LineSegment d p r) (SubLine d p r r) _SubLine = iso segment2SubLine subLineToSegment {-# INLINE _SubLine #-} -segment2SubLine    :: (Fractional r, Eq r, Arity d)-                   => LineSegment d p r -> SubLine d p r-segment2SubLine ss = SubLine l (Interval s e)+segment2SubLine    :: (Num r, Arity d)+                   => LineSegment d p r -> SubLine d p r r+segment2SubLine ss = SubLine (Line p (q .-. p)) (Interval s e)   where-    l = supportingLine ss-    f = flip toOffset l-    (Interval p q)  = ss^.unLineSeg--    s = p&unEndPoint.core %~ f-    e = q&unEndPoint.core %~ f+    p = ss^.start.core+    q = ss^.end.core+    (Interval a b)  = ss^.unLineSeg+    s = a&unEndPoint.core .~ 0+    e = b&unEndPoint.core .~ 1 -subLineToSegment    :: (Num r, Arity d) => SubLine d p r -> LineSegment d p r+subLineToSegment    :: (Num r, Arity d) => SubLine d p r r -> LineSegment d p r subLineToSegment sl = let (Interval s' e') = (fixEndPoints sl)^.subRange                           s = s'&unEndPoint %~ (^.extra)                           e = e'&unEndPoint %~ (^.extra)@@ -125,8 +126,8 @@ deriving instance Arity d                   => Functor (LineSegment d p)  instance PointFunctor (LineSegment d p) where-  pmap f ~(LineSegment s e) = LineSegment (s&unEndPoint %~ first f)-                                          (e&unEndPoint %~ first f)+  pmap f ~(LineSegment s e) = LineSegment (s&unEndPoint.core %~ f)+                                          (e&unEndPoint.core %~ f)  instance Arity d => IsBoxable (LineSegment d p r) where   boundingBox l = boundingBox (l^.start.core) <> boundingBox (l^.end.core)@@ -174,14 +175,12 @@          (LineSegment 2 p r) `IsIntersectableWith` (Line 2 r) where   nonEmptyIntersection = defaultNonEmptyIntersection -  ~s@(LineSegment p q) `intersect` l = let f  = bimap (fmap Val) (const ())-                                           s' = LineSegment (p&unEndPoint %~ f)-                                                            (q&unEndPoint %~ f)-                                    in match ((s'^._SubLine) `intersect` (fromLine l)) $-         (H   coRec)-      :& (H $ coRec . fmap (_unUnBounded))-      :& (H $ const (coRec s))-      :& RNil+  s `intersect` l = let ubSL = s^._SubLine.re _unBounded.to dropExtra+                    in match (ubSL `intersect` (fromLine l)) $+                            (H   coRec)+                         :& (H $ coRec)+                         :& (H $ const (coRec s))+                         :& RNil  -- * Functions on LineSegments 

@@ -1,25 +1,36 @@ {-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE PartialTypeSignatures #-} {-# LANGUAGE ScopedTypeVariables #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.PlnarSubdivision+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Data type to represent a PlanarSubdivision+--+-------------------------------------------------------------------------------- module Data.Geometry.PlanarSubdivision( module Data.Geometry.PlanarSubdivision.Basic-                                      , fromPolygon, fromPolygons+                                      -- , module Data.Geometry.PlanarSubdivision.Build+                                      , fromPolygon                                       ) where  -- 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.Box import           Data.Geometry.PlanarSubdivision.Basic import           Data.Geometry.Polygon-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+import           Data.Proxy+import           Data.Util+import           Control.Lens+import           Data.Bitraversable+import qualified Data.Foldable as F+import qualified Data.Vector as V+import qualified Data.Vector.Mutable as MV +--------------------------------------------------------------------------------  -- | Construct a planar subdivision from a polygon. Since our PlanarSubdivision -- models only connected planar subdivisions, this may add dummy/invisible@@ -50,17 +61,7 @@   --- | 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@@ -98,10 +99,3 @@ getP            :: HoleData f p -> Maybe p getP (Outer _)  = Nothing getP (Hole _ p) = Just p---- | 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))

@@ -2,60 +2,71 @@ {-# LANGUAGE OverloadedStrings #-} {-# LANGUAGE PartialTypeSignatures #-} {-# LANGUAGE ScopedTypeVariables #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.PlnarSubdivision.Basic+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+-- Description :  Basic data types to represent a PlanarSubdivision+--+-------------------------------------------------------------------------------- module Data.Geometry.PlanarSubdivision.Basic( VertexId', FaceId'-                                           , VertexData(VertexData), PG.vData, PG.location+                                            , VertexData(VertexData), PG.vData, PG.location -                                           , FaceData(FaceData), holes, fData+                                            , FaceData(FaceData), holes, fData -                                           , PlanarSubdivision(PlanarSubdivision)-                                           , Wrap+                                            , PlanarSubdivision(PlanarSubdivision)+                                            , Wrap -                                           , Component, ComponentId+                                            , Component, ComponentId -                                           , PolygonFaceData(..)-                                           , PlanarGraph-                                           , PlaneGraph-                                           , fromSimplePolygon-                                           , fromConnectedSegments-                                           , fromPlaneGraph, fromPlaneGraph'+                                            , PolygonFaceData(..)+                                            , PlanarGraph+                                            , PlaneGraph+                                            , fromSimplePolygon+                                            , fromConnectedSegments+                                            , fromPlaneGraph, fromPlaneGraph' -                                           , numVertices, numEdges, numFaces, numDarts-                                           , dual+                                            , numVertices, numEdges, numFaces, numDarts+                                            , dual -                                           , components, component-                                           , vertices', vertices-                                           , edges', edges-                                           , faces', faces, internalFaces-                                           , darts'+                                            , components, component+                                            , vertices', vertices+                                            , edges', edges+                                            , faces', faces, internalFaces+                                            , darts'+                                            -- , traverseVertices, traverseDarts, traverseFaces -                                           , headOf, tailOf, twin, endPoints+                                            , headOf, tailOf, twin, endPoints -                                           , incidentEdges, incomingEdges, outgoingEdges-                                           , nextIncidentEdge-                                           , neighboursOf+                                            , incidentEdges, incomingEdges, outgoingEdges+                                            , nextIncidentEdge+                                            , neighboursOf -                                           , leftFace, rightFace-                                           , outerBoundaryDarts, boundaryVertices, holesOf-                                           , outerFaceId-                                           , boundary'+                                            , leftFace, rightFace+                                            , outerBoundaryDarts, boundaryVertices, holesOf+                                            , outerFaceId+                                            , boundary' -                                           , locationOf-                                           , HasDataOf(..)+                                            , locationOf+                                            , HasDataOf(..) -                                           , endPointsOf, endPointData+                                            , endPointsOf, endPointData -                                           , edgeSegment, edgeSegments-                                           , rawFacePolygon, rawFaceBoundary-                                           , rawFacePolygons+                                            , edgeSegment, edgeSegments+                                            , rawFacePolygon, rawFaceBoundary+                                            , rawFacePolygons -                                           , VertexId(..), FaceId(..), Dart, World(..)+                                            , VertexId(..), FaceId(..), Dart, World(..)  -                                           , rawVertexData, rawDartData, rawFaceData-                                           , dataVal+                                            , rawVertexData, rawDartData, rawFaceData+                                            , vertexData, dartData, faceData+                                            , dataVal -                                           , dartMapping, Raw(..)-                                           ) where+                                            , dartMapping, Raw(..)+                                            ) where  import           Control.Lens hiding (holes, holesOf, (.=)) import           Data.Aeson@@ -67,11 +78,8 @@ 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           Data.PlanarGraph (isPositive, allDarts) import qualified Data.PlaneGraph as PG import           Data.PlaneGraph( PlaneGraph, PlanarGraph, dual                                 , Dart, VertexId(..), FaceId(..), twin@@ -110,17 +118,23 @@   Wrap s = Wrap' s  newtype ComponentId s = ComponentId { unCI :: Int }-  deriving (Show,Eq,Ord,Generic,Bounded,Enum)+  deriving (Show,Eq,Ord,Generic,Bounded,Enum,ToJSON,FromJSON)  -data Raw s ia a = Raw { _compId  :: {-# UNPACK #-} !(ComponentId s)-                      , _idxVal  :: {-# UNPACK #-} !ia+data Raw s ia a = Raw { _compId  :: !(ComponentId s)+                      , _idxVal  :: !ia                       , _dataVal :: !a-                      } deriving (Eq,Show,Functor,Foldable,Traversable)-makeLenses ''Raw+                      } deriving (Eq,Show,Functor,Foldable,Traversable,Generic) +instance (FromJSON ia, FromJSON a) => FromJSON (Raw s ia a)+instance (ToJSON ia, ToJSON a) => ToJSON (Raw s ia a) where+  toEncoding = genericToEncoding defaultOptions +-- | get the dataVal of a Raw+dataVal :: Lens (Raw s ia a) (Raw s ia b) a b+dataVal = lens (\(Raw _ _ x) -> x) (\(Raw c i _) y -> Raw c i y) + --------------------------------------------------------------------------------  -- | A connected component.@@ -147,7 +161,7 @@                     , _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)+                    } deriving (Show,Eq,Functor,Generic) makeLenses ''PlanarSubdivision  @@ -160,8 +174,14 @@  component    :: ComponentId s -> Lens' (PlanarSubdivision s v e f r)                                        (Component s r)-component ci = components.ix' (unCI ci)+component ci = components.singular (ix $ unCI ci) +-- instance (ToJSON v, ToJSON v, ToJSON e, ToJSON f, ToJSON r)+--          => ToJSON (PlanarSubdivision s v e f r) where+--   toEncoding = genericToEncoding defaultOptions+++ --------------------------------------------------------------------------------  -- | Constructs a planarsubdivision from a PlaneGraph@@ -196,10 +216,10 @@     (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)+    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'@@ -235,6 +255,27 @@                          -> PlanarSubdivision s (NonEmpty p) e () r fromConnectedSegments px = fromPlaneGraph . PG.fromConnectedSegments px +-- g1 = PG.fromConnectedSegments (Identity Test1) testSegs+-- ps1 = fromConnectedSegments (Identity Test1) testSegs++-- data Test1 = Test1++-- draw = V.filter isEmpty . rawFacePolygons+--   where+--     isEmpty (_,Left  p :+ _) = (< 3) . length . polygonVertices $ p+--     isEmpty (_,Right p :+ _) = (< 3) . length . polygonVertices $ p++-- testSegs = map (\(p,q) -> ClosedLineSegment (ext p) (ext q) :+ ())+--                    [ (origin, Point2 10 10)+--                    , (origin, Point2 12 10)+--                    , (origin, Point2 20 5)+--                    , (origin, Point2 13 20)+--                    , (Point2 10 10, Point2 12 10)+--                    , (Point2 10 10, Point2 13 20)+--                    , (Point2 12 10, Point2 20 5)+--                    ]++ --------------------------------------------------------------------------------  -- | Data type that expresses whether or not we are inside or outside the@@ -321,6 +362,9 @@ 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@@ -329,15 +373,42 @@                    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+-- | 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 = lens getF setF+  where+    getF     = V.imap (\i x -> (toEnum i, x^.dataVal)) . _rawDartData+    setF ps ds' = ps&rawDartData %~ mkDS' ds' +    -- create a new dartData vector to assign the values to+    mkDS' ds' ds = V.create $ do+                     v <- MV.new (V.length ds)+                     mapM_ (assignDart ds v) ds'+                     pure v +    assignDart ds v (d,x) = let i = fromEnum d+                                y = ds V.! i+                            in MV.write v i (y&dataVal .~ x)  +-- | Lens to the facedata of the faces themselves. The indices correspond to the faceIds+faceData :: Lens (PlanarSubdivision s v e f r) (PlanarSubdivision s v e f' r)+                 (V.Vector f)                  (V.Vector f')+faceData = lens getF setF+  where+    getF = fmap (^.dataVal) . _rawFaceData+    setF ps v' = ps&rawFaceData %~ V.zipWith (\x' x -> x&dataVal .~ x') v' +-- | Lens to the facedata of the vertexdata themselves. The indices correspond to the vertexId's+vertexData :: Lens (PlanarSubdivision s v e f r) (PlanarSubdivision s v' e f r)+                   (V.Vector v)                  (V.Vector v')+vertexData = lens getF setF+  where+    getF = fmap (^.dataVal) . _rawVertexData+    setF ps v' = ps&rawVertexData %~ V.zipWith (\x' x -> x&dataVal .~ x') v'++ -- | The tail of a dart, i.e. the vertex this dart is leaving from -- -- running time: \(O(1)\)@@ -453,7 +524,7 @@  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)+asLocalD d ps = let (Raw ci d' _) = ps^?!rawDartData.ix (fromEnum d)                 in (ci,d',ps^.component ci)  @@ -461,12 +532,12 @@  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+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+asLocalF (FaceId (VertexId f)) ps = let (Raw ci f' _) = ps^?!rawFaceData.ix f                                     in (ci,f',ps^.component ci)  @@ -476,11 +547,11 @@                            -> 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+    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)+    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)@@ -492,15 +563,15 @@ faceDataOf fi = lens getF setF   where     (FaceId (VertexId i)) = fi-    getF ps = let (Raw ci wfi x) = ps^.rawFaceData.ix' i+    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+    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+                 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@@ -508,14 +579,39 @@  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+  dataOf d = rawDartData.singular (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  +-- -- | Traverse the vertices+-- --+-- traverseVertices   :: Applicative m+--                    => (VertexId' s -> v -> m v')+--                    -> PlanarSubdivision s v e f r+--                    -> m (PlanarSubdivision s v' e f r)+-- traverseVertices f = itraverseOf (vertexData.itraversed) (\i -> f (VertexId i)) +-- -- | Traverses the darts+-- --+-- traverseDarts   :: Applicative m+--                 => (Dart s -> e -> m e')+--                 -> PlanarSubdivision s v e f r+--                 -> m (PlaneGraph s v e' f r)+-- traverseDarts f = traverseOf (dart) (PG.traverseDarts f)+++-- -- | Traverses the faces+-- --+-- traverseFaces   :: Applicative m+--                 => (FaceId' s  -> f -> m f')+--                 -> PlaneGraph s v e f r+--                 -> m (PlaneGraph s v e f' r)+-- traverseFaces f = traverseOf graph (PG.traverseFaces f)++ -- | Getter for the data at the endpoints of a dart -- -- running time: \(O(1)\)@@ -579,7 +675,7 @@ -- | 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+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@@ -595,4 +691,5 @@   +dartMapping    :: PlanarSubdivision s v e f r -> V.Vector (Dart (Wrap s), Dart s) dartMapping ps = ps^.component (ComponentId 0).PG.dartData

@@ -1,17 +1,54 @@ module Data.Geometry.PlanarSubdivision.Draw where -import           Data.Geometry.Ipe-import           Data.Ext import           Control.Lens+import           Data.Ext+import           Data.Geometry.Ipe+import           Data.Geometry.LineSegment+import           Data.Geometry.PlanarSubdivision+import           Data.Geometry.Polygon+import           Data.Maybe (mapMaybe) import qualified Data.Vector as V-import Data.Geometry.PlanarSubdivision -drawPlanarSubdivision :: forall s v e f r. IpeOut (PlanarSubdivision s v e f r) (IpeObject r)-drawPlanarSubdivision = IpeOut draw++drawColoredPlanarSubdivision  ::  IpeOut (PlanarSubdivision s v e (Maybe (IpeColor r)) r)+                                          Group r+drawColoredPlanarSubdivision ps = drawPlanarSubdivision+    (ps&vertexData.traverse  .~ Just mempty+       &dartData.traverse._2 .~ Just mempty+       &faceData.traverse    %~ fmap (attr SFill)+    )++-- | Draws only the values for which we have a Just attribute+drawPlanarSubdivision :: forall s r.+                         IpeOut (PlanarSubdivision s (Maybe (IpeAttributes IpeSymbol r))+                                                     (Maybe (IpeAttributes Path      r))+                                                     (Maybe (IpeAttributes Path      r))+                                r) Group r+drawPlanarSubdivision = drawPlanarSubdivisionWith fv fe ff   where-    draw   :: PlanarSubdivision s v e f r -> IpeObject r-    draw g = asIpeGroup $ concatMap V.toList [vs, es, fs]-      where-        vs = (\(_,VertexData p _) -> asIpeObject p mempty) <$> vertices g-        es = (\(_,s :+ _) -> asIpeObject s mempty) <$> edgeSegments g-        fs = (\(_,f :+ _) -> asIpeObject f mempty) <$> rawFacePolygons g+    fv                     :: (VertexId' s, VertexData r (Maybe (IpeAttributes IpeSymbol r)))+                           -> Maybe (IpeObject' IpeSymbol r)+    fv (_,VertexData p ma) = (\a -> defIO p ! a) <$> ma -- draws a point+    fe (_,s :+ ma)         = (\a -> defIO s ! a) <$> ma -- draw segment+    ff (_,f :+ ma)         = (\a -> defIO f ! a) <$> ma -- draw a face+++-- | Draw everything using the defaults+drawPlanarSubdivision'    :: forall s v e f r. IpeOut (PlanarSubdivision s v e f r) Group r+drawPlanarSubdivision' ps = drawPlanarSubdivision+  (ps&vertexData.traverse   .~ Just (mempty :: IpeAttributes IpeSymbol r)+     &dartData.traverse._2  .~ Just (mempty :: IpeAttributes Path      r)+     &faceData.traverse     .~ Just (mempty :: IpeAttributes Path      r))++type MIO g i r = g -> Maybe (IpeObject' i r)++drawPlanarSubdivisionWith            :: (ToObject vi, ToObject ei, ToObject fi)+                                     => MIO (VertexId' s, VertexData r v)          vi r+                                     -> MIO (Dart s,      LineSegment 2 v r :+ e)  ei r+                                     -> MIO (FaceId' s,   SomePolygon v r :+ f)    fi r+                                     -> IpeOut (PlanarSubdivision s v e f r) Group r+drawPlanarSubdivisionWith fv fe ff g = ipeGroup . concat $ [vs, es, fs]+  where+    vs = mapMaybe (fmap iO . fv) . V.toList . vertices        $ g+    es = mapMaybe (fmap iO . fe) . V.toList . edgeSegments    $ g+    fs = mapMaybe (fmap iO . ff) . V.toList . rawFacePolygons $ g

@@ -1,11 +1,15 @@ {-# LANGUAGE ScopedTypeVariables  #-} {-# LANGUAGE UndecidableInstances #-}-{-|-Module    : Data.Geometry.Point-Description: \(d\)-dimensional points-Copyright : (c) Frank Staals-License : See LICENCE file--}+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.Point+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- \(d\)-dimensional points.+--+-------------------------------------------------------------------------------- module Data.Geometry.Point where  import           Control.DeepSeq@@ -22,8 +26,11 @@ import           Data.Proxy import           GHC.Generics (Generic) import           GHC.TypeLits---------------------------------------------------------------------------------+import           Text.ParserCombinators.ReadP (ReadP, string,pfail)+import           Text.ParserCombinators.ReadPrec (lift)+import           Text.Read (Read(..),readListPrecDefault, readPrec_to_P,minPrec) + -------------------------------------------------------------------------------- -- $setup -- >>> :{@@ -43,6 +50,17 @@   show (Point v) = mconcat [ "Point", show $ F.length v , " "                            , show $ F.toList v                            ]+instance (Read r, Arity d) => Read (Point d r) where+  readPrec     = lift readPt+  readListPrec = readListPrecDefault++readPt :: forall d r. (Arity d, Read r) => ReadP (Point d r)+readPt = do let d = natVal (Proxy :: Proxy d)+            _  <- string $ "Point" <> show d <> " "+            rs <- readPrec_to_P readPrec minPrec+            case pointFromList rs of+              Just p -> pure p+              _      -> pfail  deriving instance (Eq r, Arity d)     => Eq (Point d r) deriving instance (Ord r, Arity d)    => Ord (Point d r)

@@ -6,23 +6,23 @@ import           Control.Lens import           Data.Bifunctor import           Data.Ext+import qualified Data.Foldable as F import           Data.Geometry.Box import           Data.Geometry.LineSegment import           Data.Geometry.Point import           Data.Geometry.Properties import           Data.Geometry.Transformation import           Data.Geometry.Vector+import           Data.LSeq (LSeq, pattern (:<|))+import qualified Data.LSeq as LSeq import qualified Data.List.NonEmpty as NE-import           Data.Semigroup-import qualified Data.Seq2 as S2-import qualified Data.Sequence as Seq import           GHC.TypeLits  -------------------------------------------------------------------------------- -- * d-dimensional Polygonal Lines (PolyLines) --- | A Poly line in R^d-newtype PolyLine d p r = PolyLine { _points :: S2.Seq2 (Point d r :+ p) }+-- | A Poly line in R^d has at least 2 vertices+newtype PolyLine d p r = PolyLine { _points :: LSeq 2 (Point d r :+ p) } makeLenses ''PolyLine  deriving instance (Show r, Show p, Arity d) => Show    (PolyLine d p r)@@ -53,7 +53,7 @@  -- | pre: The input list contains at least two points fromPoints :: [Point d r :+ p] -> PolyLine d p r-fromPoints = PolyLine . S2.fromList+fromPoints = PolyLine . LSeq.forceLSeq (C  :: C 2) . LSeq.fromList  -- | pre: The input list contains at least two points. All extra vields are -- initialized with mempty.@@ -66,19 +66,12 @@ fromLineSegment ~(LineSegment' p q) = fromPoints [p,q]  -- | Convert to a closed line segment by taking the first two points.-asLineSegment                              :: PolyLine d p r -> LineSegment d p r-asLineSegment (PolyLine (S2.Seq2 p mid q)) = ClosedLineSegment p (f $ Seq.viewl mid)-  where-    f Seq.EmptyL    = q-    f (q' Seq.:< _) = q'+asLineSegment                            :: PolyLine d p r -> LineSegment d p r+asLineSegment (PolyLine (p :<| q :<| _)) = ClosedLineSegment p q  -- | Stricter version of asLineSegment that fails if the Polyline contains more -- than two points.-asLineSegment'                            :: PolyLine d p r -> Maybe (LineSegment d p r)-asLineSegment' (PolyLine (S2.Seq2 p m q))-  | Seq.null m                            = Just $ ClosedLineSegment p q-  | otherwise                             = Nothing----- polylineEdges :: Polyline d p r -> NonEmpty.NonEmpty (LineSegment d p r)--- polylineEdges (Polyline )+asLineSegment'                :: PolyLine d p r -> Maybe (LineSegment d p r)+asLineSegment' (PolyLine pts) = case F.toList pts of+                                  [p,q] -> Just $ ClosedLineSegment p q+                                  _     -> Nothing

@@ -1,10 +1,14 @@ {-# LANGUAGE ScopedTypeVariables  #-}-{-|-Module    : Data.Geometry.Polygon-Description: A Polygon data type and some basic functions to interact with them.-Copyright : (c) Frank Staals-License : See LICENCE file--}+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.Polygon+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- A Polygon data type and some basic functions to interact with them.+--+-------------------------------------------------------------------------------- module Data.Geometry.Polygon where  import           Control.DeepSeq@@ -26,7 +30,7 @@ import qualified Data.List.NonEmpty as NonEmpty import           Data.Maybe (mapMaybe) import           Data.Proxy-import           Data.Semigroup+import           Data.Semigroup(sconcat) import qualified Data.Sequence as Seq import           Data.Util import           Data.Vinyl.CoRec (asA)

@@ -1,11 +1,15 @@ {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TemplateHaskell #-}-{-|-Module    : Data.Geometry.Polygon.Convex-Description: Convex Polygons-Copyright : (c) Frank Staals-License : See LICENCE file--}+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.Polygon.Convex+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Convex Polygons+--+-------------------------------------------------------------------------------- module Data.Geometry.Polygon.Convex( ConvexPolygon(..), simplePolygon                                    , merge                                    , lowerTangent, upperTangent@@ -41,8 +45,10 @@  -- import           Data.Geometry.Ipe -- import           Debug.Trace+ -------------------------------------------------------------------------------- +-- | Data Type representing a convex polygon newtype ConvexPolygon p r = ConvexPolygon {_simplePolygon :: SimplePolygon p r }                           deriving (Show,Eq,NFData) makeLenses ''ConvexPolygon

@@ -2,12 +2,17 @@ {-# LANGUAGE ImpredicativeTypes #-} {-# LANGUAGE UnicodeSyntax #-} {-# LANGUAGE DefaultSignatures #-}-{-|-Module    : Data.Geometry.Properties-Description: Defines some generic geometric properties e.g. Dimensions, NumType, and Intersection types.-Copyright : (c) Frank Staals-License : See LICENCE file--}+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.Properties+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Defines some generic geometric properties e.g. Dimensions, NumType, and+-- Intersection types.+--+-------------------------------------------------------------------------------- module Data.Geometry.Properties where  import Data.Maybe (isNothing)@@ -26,6 +31,8 @@  -- | A type family for types that have an associated numeric type. type family NumType t :: *++type instance NumType [t] = NumType t  -- | A simple data type expressing that there are no intersections data NoIntersection = NoIntersection deriving (Show,Read,Eq,Ord)

@@ -1,18 +1,4 @@-{-# LANGUAGE TemplateHaskell #-} module Data.Geometry.SegmentTree( module Data.Geometry.SegmentTree.Generic                                 ) where --import           Data.Ext-import           Data.Semigroup-import           Control.Lens-import           Data.List.NonEmpty (NonEmpty)-import qualified Data.List.NonEmpty as NonEmpty-import qualified Data.List as List-import           Data.BinaryTree-import           Data.Range-import           Data.Geometry.Interval import           Data.Geometry.SegmentTree.Generic-import           Data.Geometry.Interval.Util-import           Data.Geometry.Properties-import           Data.Geometry.IntervalTree (IntervalLike(..))

@@ -1,4 +1,13 @@ {-# LANGUAGE TemplateHaskell #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.SegmentTree.Generic+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+-- Description :  Implementation of SegmentTrees+--+-------------------------------------------------------------------------------- module Data.Geometry.SegmentTree.Generic( NodeData(..), splitPoint, range, assoc                                         , LeafData(..), atomicRange, leafAssoc @@ -19,14 +28,12 @@ import           Control.DeepSeq import           Control.Lens import           Data.BinaryTree-import           Data.Ext import           Data.Geometry.Interval import           Data.Geometry.IntervalTree (IntervalLike(..)) import           Data.Geometry.Properties import qualified Data.List as List import           Data.List.NonEmpty (NonEmpty) import qualified Data.List.NonEmpty as NonEmpty-import           Data.Semigroup import           GHC.Generics (Generic)  --------------------------------------------------------------------------------@@ -258,13 +265,16 @@ -------------------------------------------------------------------------------- -- * Counting the number of segments intersected -newtype Count = Count { getCount :: Int}+newtype Count = Count { getCount :: Word }               deriving (Show,Eq,Ord,Num,Integral,Enum,Real,Generic,NFData)  newtype C a = C { _unC :: a} deriving (Show,Read,Eq,Ord,Generic,NFData)  instance Semigroup Count where   a <> b = Count $ getCount a + getCount b+instance Monoid Count where+  mempty = 0+  mappend = (<>)  instance Measured Count (C i) where   measure _ = 1@@ -277,29 +287,29 @@ -------------------------------------------------------------------------------- -- * Testing stuff -test'' = fromIntervals' . NonEmpty.fromList $ test-test = [Interval (Closed (238 :+ ())) (Open (309 :+ ())), Interval (Closed (175 :+ ())) (Closed (269 :+ ())),Interval (Closed (255 :+ ())) (Open (867 :+ ())),Interval (Open (236 :+ ())) (Closed (863 :+ ())),Interval (Open (150 :+ ())) (Closed (161 :+ ())),Interval (Closed (35 :+ ())) (Closed (77 :+ ()))]+-- test'' = fromIntervals' . NonEmpty.fromList $ test+-- test = [Interval (Closed (238 :+ ())) (Open (309 :+ ())), Interval (Closed (175 :+ ())) (Closed (269 :+ ())),Interval (Closed (255 :+ ())) (Open (867 :+ ())),Interval (Open (236 :+ ())) (Closed (863 :+ ())),Interval (Open (150 :+ ())) (Closed (161 :+ ())),Interval (Closed (35 :+ ())) (Closed (77 :+ ()))]  --- q =        [78]+-- -- q =        [78] --- test = fromIntervals' . NonEmpty.fromList $ [ closedInterval 0 10---                                             , closedInterval 5 15---                                             , closedInterval 1 4---                                             , closedInterval 3 9---                                             ]-tst = fromIntervals' . NonEmpty.fromList $ [ closedInterval 1 6-                                           , closedInterval 2 6-                                           -- , Interval (Closed $ ext 0) (Open $ ext 1)-                                           ]+-- -- test = fromIntervals' . NonEmpty.fromList $ [ closedInterval 0 10+-- --                                             , closedInterval 5 15+-- --                                             , closedInterval 1 4+-- --                                             , closedInterval 3 9+-- --                                             ]+-- tst = fromIntervals' . NonEmpty.fromList $ [ closedInterval 1 6+--                                            , closedInterval 2 6+--                                            -- , Interval (Closed $ ext 0) (Open $ ext 1)+--                                            ]   -closedInterval a b = ClosedInterval (ext a) (ext b)+-- closedInterval a b = ClosedInterval (ext a) (ext b) -showT :: (Show r, Show v) => SegmentTree v r -> String-showT = drawTree . _unSegmentTree+-- showT :: (Show r, Show v) => SegmentTree v r -> String+-- showT = drawTree . _unSegmentTree  -test' :: (Show r, Num r, Ord r, Enum r) => SegmentTree [I (Interval () r)] r-test' = insert (I $ closedInterval 6 14) $ createTree (NonEmpty.fromList [2,4..20]) []+-- test' :: (Show r, Num r, Ord r, Enum r) => SegmentTree [I (Interval () r)] r+-- test' = insert (I $ closedInterval 6 14) $ createTree (NonEmpty.fromList [2,4..20]) []

@@ -121,11 +121,11 @@   -type instance IntersectionOf (SubLine 2 p r) (Slab o a r) =-  [NoIntersection, SubLine 2 () r]+type instance IntersectionOf (SubLine 2 p s r) (Slab o a r) =+  [NoIntersection, SubLine 2 () s r]  instance (Fractional r, Ord r, HasBoundingLines o) =>-         SubLine 2 a r `IsIntersectableWith` (Slab o a r) where+         SubLine 2 a r r `IsIntersectableWith` (Slab o a r) where    nonEmptyIntersection = defaultNonEmptyIntersection @@ -139,8 +139,7 @@                                  :& RNil)     :& RNil     where-      dropExtra sub = sub&subRange %~ bimap (const ()) id-      singleton p = let x = ext $ toOffset p l in SubLine l (ClosedInterval x x)+      singleton p = let x = ext $ toOffset' p l in SubLine l (ClosedInterval x x)   type instance IntersectionOf (LineSegment 2 p r) (Slab o a r) =@@ -158,5 +157,5 @@   -test :: SubLine 2 () Double-test = (ClosedLineSegment (ext origin) (ext origin))^._SubLine+-- test :: SubLine 2 () Double Double+-- test = (ClosedLineSegment (ext origin) (ext origin))^._SubLine

@@ -3,6 +3,7 @@ module Data.Geometry.SubLine where  import           Control.Lens+import           Data.Bifunctor import           Data.Ext import qualified Data.Foldable as F import           Data.Geometry.Interval@@ -15,39 +16,31 @@ import           Data.Vinyl import           Data.Vinyl.CoRec ++import           Data.Ratio+ --------------------------------------------------------------------------------  -- | Part of a line. The interval is ranged based on the vector of the -- line l, and s.t.t zero is the anchorPoint of l.-data SubLine d p r = SubLine { _line     :: Line d r-                             , _subRange :: Interval p r-                             }--+data SubLine d p s r = SubLine { _line     :: Line d r+                               , _subRange :: Interval p s+                               } makeLenses ''SubLine -type instance Dimension (SubLine d p r) = d-type instance NumType   (SubLine d p r) = r--deriving instance (Show r, Show p, Arity d) => Show (SubLine d p r)-deriving instance (Eq r, Eq p, Arity d)     => Eq (SubLine d p r)-deriving instance Arity d                   => Functor (SubLine d p)-deriving instance Arity d                   => F.Foldable (SubLine d p)-deriving instance Arity d                   => T.Traversable (SubLine d p)--instance Arity d => Bifunctor (SubLine d) where-  bimap f g (SubLine l r) = SubLine (g <$> l) (bimap f g r)+type instance Dimension (SubLine d p s r) = d  --- | Get the point at the given position along line, where 0 corresponds to the--- anchorPoint of the line, and 1 to the point anchorPoint .+^ directionVector-pointAt              :: (Num r, Arity d) => r -> Line d r -> Point d r-pointAt a (Line p v) = p .+^ (a *^ v)-+deriving instance (Show r, Show s, Show p, Arity d) => Show (SubLine d p s r)+-- deriving instance (Read r, Read p, Arity d) => Read (SubLine d p r)+deriving instance (Eq r, Eq s, Fractional r, Eq p, Arity d)     => Eq (SubLine d p s r)+deriving instance Arity d                   => Functor (SubLine d p s)+deriving instance Arity d                   => F.Foldable (SubLine d p s)+deriving instance Arity d                   => T.Traversable (SubLine d p s)   -- | Annotate the subRange with the actual ending points-fixEndPoints    :: (Num r, Arity d) => SubLine d p r -> SubLine d (Point d r :+ p) r+fixEndPoints    :: (Num r, Arity d) => SubLine d p r r -> SubLine d (Point d r :+ p) r r fixEndPoints sl = sl&subRange %~ f   where     ptAt              = flip pointAt (sl^.line)@@ -55,29 +48,40 @@     f ~(Interval l u) = Interval (l&unEndPoint %~ label)                                  (u&unEndPoint %~ label) +-- | forget the extra information stored at the endpoints of the subline.+dropExtra :: SubLine d p s r -> SubLine d () s r+dropExtra = over subRange (first (const ())) --- | given point p on line (Line q v), Get the scalar lambda s.t.--- p = q + lambda v-toOffset              :: (Eq r, Fractional r, Arity d) => Point d r -> Line d r -> r-toOffset p (Line q v) = fromJust' $ scalarMultiple (p .-. q) v-  where-    fromJust' (Just x) = x-    fromJust' _        = error "toOffset: Nothing"+_unBounded :: Prism' (SubLine d p (UnBounded r) r) (SubLine d p r r)+_unBounded = prism' toUnbounded fromUnbounded -type instance IntersectionOf (SubLine 2 p r) (SubLine 2 q r) = [ NoIntersection-                                                               , Point 2 r-                                                               , SubLine 2 p r-                                                               ]+-- | Transform into an subline with a potentially unbounded interval+toUnbounded :: SubLine d p r r -> SubLine d p (UnBounded r) r+toUnbounded = over subRange (fmap Val) +-- | Try to make a potentially unbounded subline into a bounded one.+fromUnbounded               :: SubLine d p (UnBounded r) r -> Maybe (SubLine d p r r)+fromUnbounded (SubLine l i) = SubLine l <$> mapM unBoundedToMaybe i+ -- | given point p, and a Subline l r such that p lies on line l, test if it -- lies on the subline, i.e. in the interval r onSubLine                 :: (Ord r, Fractional r, Arity d)-                          => Point d r -> SubLine d p r -> Bool-onSubLine p (SubLine l r) = toOffset p l `inInterval` r+                          => Point d r -> SubLine d p r r -> Bool+onSubLine p (SubLine l r) = case toOffset p l of+                              Nothing -> False+                              Just x  -> x `inInterval` r  -- | given point p, and a Subline l r such that p lies on line l, test if it -- lies on the subline, i.e. in the interval r-onSubLine2        :: (Ord r, Num r) => Point 2 r -> SubLine 2 p r -> Bool+onSubLineUB                   :: (Ord r, Fractional r)+                              => Point 2 r -> SubLine 2 p (UnBounded r) r -> Bool+p `onSubLineUB` (SubLine l r) = case toOffset p l of+                                  Nothing -> False+                                  Just x  -> Val x `inInterval` r++-- | given point p, and a Subline l r such that p lies on line l, test if it+-- lies on the subline, i.e. in the interval r+onSubLine2        :: (Ord r, Num r) => Point 2 r -> SubLine 2 p r r -> Bool p `onSubLine2` sl = d `inInterval` r   where     -- get the endpoints (a,b) of the subline@@ -87,11 +91,22 @@     d = (p .-. a) `dot` (b .-. a)     -- map to an interval corresponding to the length of the segment     r = Interval (s&unEndPoint.core .~ 0) (e&unEndPoint.core .~ squaredEuclideanDist b a)-      -- note that we take the dist between b and a, so if these are infinity-      -- we get maxInfinity as well ++-- | given point p, and a Subline l r such that p lies on line l, test if it+-- lies on the subline, i.e. in the interval r+onSubLine2UB        :: (Ord r, Fractional r)+                    => Point 2 r -> SubLine 2 p (UnBounded r) r -> Bool+p `onSubLine2UB` sl = p `onSubLineUB` sl+++type instance IntersectionOf (SubLine 2 p s r) (SubLine 2 q s r) = [ NoIntersection+                                                                   , Point 2 r+                                                                   , SubLine 2 p s r+                                                                   ]+ instance (Ord r, Fractional r) =>-         (SubLine 2 p r) `IsIntersectableWith` (SubLine 2 p r) where+         (SubLine 2 p r r) `IsIntersectableWith` (SubLine 2 p r r) where    nonEmptyIntersection = defaultNonEmptyIntersection @@ -107,16 +122,43 @@            )       :& RNil     where-      -- s' = shiftLeft' (toOffset (m^.anchorPoint) l) $ s       s'  = (fixEndPoints sm)^.subRange       s'' = bimap (^.extra) id-          $ s'&start.core .~ toOffset (s'^.start.extra.core) l-              &end.core   .~ toOffset (s'^.end.extra.core)   l+          $ s'&start.core .~ toOffset' (s'^.start.extra.core) l+              &end.core   .~ toOffset' (s'^.end.extra.core)   l -fromLine   :: Arity d => Line d r -> SubLine d () (UnBounded r)-fromLine l = SubLine (fmap Val l) (ClosedInterval (ext MinInfinity) (ext MaxInfinity))+instance (Ord r, Fractional r) =>+         (SubLine 2 p (UnBounded r) r) `IsIntersectableWith` (SubLine 2 p (UnBounded r) r) where+  nonEmptyIntersection = defaultNonEmptyIntersection +  sl@(SubLine l r) `intersect` sm@(SubLine m _) = match (l `intersect` m) $+         (H $ \NoIntersection -> coRec NoIntersection)+      :& (H $ \p@(Point _)    -> if onSubLine2UB p sl && onSubLine2UB p sm+                                 then coRec p+                                 else coRec NoIntersection)+      :& (H $ \_             -> match (r `intersect` s'') $+                                      (H $ \NoIntersection -> coRec NoIntersection)+                                   :& (H $ \i              -> coRec $ SubLine l i)+                                   :& RNil+           )+      :& RNil+    where+      -- convert to points, then convert back to 'r' values (but now w.r.t. l)+      s'  = getEndPointsUnBounded sm+      s'' = second (fmap f) s'+      f = flip toOffset' l +-- | Get the endpoints of an unbounded interval+getEndPointsUnBounded    :: (Num r, Arity d) => SubLine d p (UnBounded r) r+                         -> Interval p (UnBounded (Point d r))+getEndPointsUnBounded sl = second (fmap f) $ sl^.subRange+  where+    f = flip pointAt (sl^.line)++fromLine   :: Arity d => Line d r -> SubLine d () (UnBounded r) r+fromLine l = SubLine l (ClosedInterval (ext MinInfinity) (ext MaxInfinity))++ -- testL :: SubLine 2 () (UnBounded Rational) -- testL = SubLine (horizontalLine 0) (Interval (Closed (only 0)) (Open $ only 10)) @@ -125,3 +167,11 @@   -- test = (testL^.subRange) `intersect` (horL^.subRange)++-- toOffset (Point2 minInfinity minInfinity) (horizontalLine 0)+-- testzz = let f  = bimap (fmap Val) (const ())+--          in++testz :: SubLine 2 () Rational Rational+testz = SubLine (Line (Point2 0 0) (Vector2 10 0))+                (Interval (Closed (0 % 1 :+ ())) (Closed (1 % 1 :+ ())))

@@ -8,6 +8,7 @@ module Data.Geometry.Transformation where  import           Control.Lens (lens,Lens',set)+import           Unsafe.Coerce(unsafeCoerce) import           Data.Geometry.Point import           Data.Geometry.Properties import           Data.Geometry.Vector@@ -16,6 +17,7 @@ import qualified Data.Vector.Fixed as FV import           GHC.TypeLits import           Linear.Matrix ((!*),(!*!))+import qualified Linear.Matrix as Lin  -------------------------------------------------------------------------------- -- * Matrices@@ -34,6 +36,23 @@ mult :: (Arity m, Arity n, Num r) => Matrix n m r -> Vector m r -> Vector n r (Matrix m) `mult` v = m !* v ++class Invertible n r where+  inverse' :: Matrix n n r -> Matrix n n r++instance Fractional r => Invertible 2 r where+  -- >>> inverse' $ Matrix $ Vector2 (Vector2 1 2) (Vector2 3 4.0)+  -- Matrix Vector2 [Vector2 [-2.0,1.0],Vector2 [1.5,-0.5]]+  inverse' (Matrix m) = Matrix . unsafeCoerce . Lin.inv22 . unsafeCoerce $ m++instance Fractional r => Invertible 3 r where+  -- >>> inverse' $ Matrix $ Vector3 (Vector3 1 2 4) (Vector3 4 2 2) (Vector3 1 1 1.0)+  -- Matrix Vector3 [Vector3 [0.0,0.5,-1.0],Vector3 [-0.5,-0.75,3.5],Vector3 [0.5,0.25,-1.5]]+  inverse' (Matrix m) = Matrix . unsafeCoerce . Lin.inv33 . unsafeCoerce $ m++instance Fractional r => Invertible 4 r where+  inverse' (Matrix m) = Matrix . unsafeCoerce . Lin.inv44 . unsafeCoerce $ m+ -------------------------------------------------------------------------------- -- * Transformations @@ -55,6 +74,17 @@ (|.|) :: (Num r, Arity (d + 1)) => Transformation d r -> Transformation d r -> Transformation d r (Transformation f) |.| (Transformation g) = Transformation $ f `multM` g ++-- if it exists?++-- | Compute the inverse transformation+--+-- >>> inverseOf $ translation (Vector2 (10.0) (5.0))+-- Transformation {_transformationMatrix = Matrix Vector3 [Vector3 [1.0,0.0,-10.0],Vector3 [0.0,1.0,-5.0],Vector3 [0.0,0.0,1.0]]}+inverseOf :: (Fractional r, Invertible (d + 1) r)+          => Transformation d r -> Transformation d r+inverseOf = Transformation . inverse' . _transformationMatrix+ -------------------------------------------------------------------------------- -- * Transformable geometry objects @@ -136,3 +166,16 @@ transRow     :: forall n r. (Arity n, Arity (n + 1), Num r)              => Int -> r -> Vector (n + 1) r transRow i x = set (V.element (Proxy :: Proxy n)) x $ mkRow i 1++--------------------------------------------------------------------------------+-- * 3D Rotations++-- | Given three new unit-length basis vectors (u,v,w) that map to (x,y,z),+-- construct the appropriate rotation that does this.+--+--+rotateTo                 :: Num r => Vector 3 (Vector 3 r) -> Transformation 3 r+rotateTo (Vector3 u v w) = Transformation . Matrix $ Vector4 (snoc u        0)+                                                             (snoc v        0)+                                                             (snoc w        0)+                                                             (Vector4 0 0 0 1)

@@ -1,22 +1,38 @@+{-# LANGUAGE ScopedTypeVariables  #-} {-# LANGUAGE DeriveFunctor #-} {-# LANGUAGE UndecidableInstances #-} module Data.Geometry.Triangle where -import Data.Bifunctor-import Control.Lens-import Data.Ext-import Data.Geometry.Point-import Data.Geometry.Vector-import Data.Geometry.Ball-import Data.Geometry.Properties-import Data.Geometry.Transformation-import GHC.TypeLits+import           Control.Lens+import           Data.Bifunctor+import           Data.Either (partitionEithers)+import           Data.Ext+import           Data.Geometry.Ball (Disk, disk)+import           Data.Geometry.Boundary+import           Data.Geometry.HyperPlane+import           Data.Geometry.Line+import           Data.Geometry.LineSegment+import           Data.Geometry.Point+import           Data.Geometry.Properties+import           Data.Geometry.Transformation+import           Data.Geometry.Vector+import qualified Data.Geometry.Vector as V+import qualified Data.List as List+import           Data.Maybe (mapMaybe)+import           Data.Vinyl+import           Data.Vinyl.CoRec+import           GHC.TypeLits ++--------------------------------------------------------------------------------+ data Triangle d p r = Triangle (Point d r :+ p)                                (Point d r :+ p)                                (Point d r :+ p)  deriving instance (Arity d, Show r, Show p) => Show (Triangle d p r)+deriving instance (Arity d, Read r, Read p) => Read (Triangle d p r)+deriving instance (Arity d, Eq r, Eq p)     => Eq (Triangle d p r)  instance Arity d => Functor (Triangle d p) where   fmap f (Triangle p q r) = let f' = first (fmap f) in Triangle (f' p) (f' q) (f' r)@@ -32,6 +48,15 @@   transformBy = transformPointFunctor  +-- | convenience function to construct a triangle without associated data.+triangle'       :: Point d r -> Point d r -> Point d r -> Triangle d () r+triangle' p q r = Triangle (ext p) (ext q) (ext r)+++sideSegments                  :: Triangle d p r -> [LineSegment d p r]+sideSegments (Triangle p q r) =+  [ClosedLineSegment p q, ClosedLineSegment q r, ClosedLineSegment r p]+ -- | Compute the area of a triangle area   :: Fractional r => Triangle 2 p r -> r area t = doubleArea t / 2@@ -47,9 +72,144 @@     Point2 bx by = b^.core     Point2 cx cy = c^.core +-- | Checks if the triangle is degenerate, i.e. has zero area.+isDegenerateTriangle :: (Num r, Eq r) => Triangle 2 p r -> Bool+isDegenerateTriangle = (== 0) . doubleArea  -- | get the inscribed disk. Returns Nothing if the triangle is degenerate, -- i.e. if the points are colinear. inscribedDisk                  :: (Eq r, Fractional r)                                => Triangle 2 p r -> Maybe (Disk () r) inscribedDisk (Triangle p q r) = disk (p^.core) (q^.core) (r^.core)+++instance Num r => HasSupportingPlane (Triangle 3 p r) where+  supportingPlane (Triangle p q r) = from3Points (p^.core) (q^.core) (r^.core)+++-- | Given a point q and a triangle, q inside the triangle, get the baricentric+-- cordinates of q+toBarricentric                                 :: Fractional r+                                               => Point 2 r -> Triangle 2 p r+                                               -> Vector 3 r+toBarricentric (Point2 qx qy) (Triangle a b c) = Vector3 alpha beta gamma+  where+    Point2 ax ay = a^.core+    Point2 bx by = b^.core+    Point2 cx cy = c^.core++    dett  = (by - cy)*(ax - cx) + (cx - bx)*(ay - cy)++    alpha = ((by - cy)*(qx - cx) + (cx - bx)*(qy - cy)) / dett+    beta  = ((cy - ay)*(qx - cx) + (ax - cx)*(qy - cy)) / dett+    gamma = 1 - alpha - beta+    -- see https://en.wikipedia.org/wiki/Barycentric_coordinate_system#Conversion_between_barycentric_and_Cartesian_coordinates++-- | Given a vector of barricentric coordinates and a triangle, get the+-- corresponding point in the same coordinate sytsem as the vertices of the+-- triangle.+fromBarricentric                                  :: (Arity d, Num r)+                                                  => Vector 3 r -> Triangle d p r+                                                  -> Point d r+fromBarricentric (Vector3 a b c) (Triangle p q r) = let f = view (core.vector) in+    Point $ a *^ f p ^+^ b *^ f q ^+^ c *^ f r+++-- | Tests if a point lies inside a triangle, on its boundary, or outside the triangle+inTriangle     :: (Ord r, Fractional r)+                 => Point 2 r -> Triangle 2 p r -> PointLocationResult+inTriangle q t+    | all (`inRange` (OpenRange   0 1)) [a,b,c] = Inside+    | all (`inRange` (ClosedRange 0 1)) [a,b,c] = OnBoundary+    | otherwise                                 = Outside+  where+    Vector3 a b c = toBarricentric q t++-- | Test if a point lies inside or on the boundary of a triangle+onTriangle       :: (Ord r, Fractional r)+                 => Point 2 r -> Triangle 2 p r -> Bool+q `onTriangle` t = let Vector3 a b c = toBarricentric q t+                   in all (`inRange` (ClosedRange 0 1)) [a,b,c]+++-- myQ :: Point 2 Rational+-- myQ = read "Point2 [(-5985) % 16,(-14625) % 1]"+-- myTri :: Triangle 2 () Rational+-- myTri = read "Triangle (Point2 [(-15) % 1,0 % 1] :+ ()) (Point2 [225 % 2,0 % 1] :+ ()) (Point2 [135 % 1,0 % 1] :+ ())"++type instance IntersectionOf (Line 2 r) (Triangle 2 p r) =+  [ NoIntersection, Point 2 r, LineSegment 2 () r ]++instance (Fractional r, Ord r) => (Line 2 r) `IsIntersectableWith` (Triangle 2 p r) where+   nonEmptyIntersection = defaultNonEmptyIntersection++   l `intersect` (Triangle p q r) =+     case first List.nub . partitionEithers . mapMaybe collect $ sides of+       ([],[])   -> coRec NoIntersection+       (_, [s])  -> coRec $ first (const ()) s+       ([a],_)   -> coRec a+       ([a,b],_) -> coRec $ ClosedLineSegment (ext a) (ext b)+       (_,_)     -> error "intersecting a line with a triangle. Triangle is degenerate"+     where+       sides = [ClosedLineSegment p q, ClosedLineSegment q r, ClosedLineSegment r p]++       collect   :: LineSegment 2 p r -> Maybe (Either (Point 2 r) (LineSegment 2 p r))+       collect s = match (s `intersect` l) $+                        (H $ \NoIntersection           -> Nothing)+                     :& (H $ \(a :: Point 2 r)         -> Just $ Left a)+                     :& (H $ \(e :: LineSegment 2 p r) -> Just $ Right e)+                     :& RNil++++type instance IntersectionOf (Line 3 r) (Triangle 3 p r) =+  [ NoIntersection, Point 3 r, LineSegment 3 () r ]++instance (Fractional r, Ord r) => (Line 3 r) `IsIntersectableWith` (Triangle 3 p r) where+   nonEmptyIntersection = defaultNonEmptyIntersection++   l@(Line a v) `intersect` t@(Triangle (p :+ _) (q :+ _) (r :+ _)) =+       match (l `intersect` h) $+            (H $ \NoIntersection   -> coRec NoIntersection)+         :& (H $ \i@(Point3 _ _ _) -> if onTriangle' i then coRec i else coRec NoIntersection)+         :& (H $ \_                -> intersect2d)+         :& RNil+     where+       h@(Plane _ n) = supportingPlane t++       -- 2d triangle and the line in terms of 2d-coordinates wr.t. of a+       -- coordinate system in the supporting plane of t. The origin of this+       -- coordinate system corresponds to the second vertex of t (q)+       t' = Triangle (ext $ project p) (ext origin) (ext $ project r)+       l' = Line (project a) (project' v)++       -- test if the point in terms of its 2d coords lies in side the projected triangle+       onTriangle'                :: Point 3 r -> Bool+       onTriangle' i = (project i) `onTriangle` t'++       -- FIXME! these vectors may not be unit vectors. How do we deal with+       -- that? (and does that really matter here?)+       transf :: Transformation 3 r+       transf = let u = p .-. q+                in rotateTo (Vector3 u (n `cross` u) n) |.| translation ((-1) *^ (toVec q))+       -- inverse of the transformation above.+       invTrans :: Transformation 3 r+       invTrans = inverseOf transf+++       project :: Point 3 r -> Point 2 r+       project = projectPoint . transformBy transf+       project' :: Vector 3 r -> Vector 2 r+       project' = toVec . project . Point++       lift :: Point 2 r -> Point 3 r+       lift = Point . transformBy invTrans . flip V.snoc 0 . toVec+         -- lift a 2d point back into plane coordinates++       intersect2d :: Intersection (Line 3 r) (Triangle 3 p r)+       intersect2d = match (l' `intersect` t') $+            (H $ \NoIntersection    -> coRec NoIntersection)+         :& (H $ \i@(Point2 _ _)    -> coRec $ lift i)+         :& (H $ \(LineSegment s e) -> coRec $ LineSegment (s&unEndPoint.core %~ lift)+                                                           (e&unEndPoint.core %~ lift))+         :& RNil

@@ -1,3 +1,13 @@+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.Vector+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- \(d\)-dimensional vectors.+--+-------------------------------------------------------------------------------- module Data.Geometry.Vector( module Data.Geometry.Vector.VectorFamily                            , module LV                            , C(..)@@ -8,17 +18,19 @@                            , scalarMultiple                            -- reexports                            , FV.replicate-                           , FV.imap,+                           , FV.imap+                           , xComponent, yComponent, zComponent                            ) where  import           Control.Applicative (liftA2)+import           Control.Lens(Lens') import qualified Data.Foldable as F import           Data.Geometry.Properties import           Data.Geometry.Vector.VectorFamily-import           Data.Geometry.Vector.VectorFixed(C(..))+import           Data.Geometry.Vector.VectorFixed (C(..)) import           Data.Maybe-import           Data.Semigroup import qualified Data.Vector.Fixed as FV+import           GHC.TypeLits import           Linear.Affine (Affine(..), qdA, distanceA) import           Linear.Metric (dot,norm,signorm) import           Linear.Vector as LV@@ -101,3 +113,19 @@     g No      = Nothing     g Maybe   = error "scalarMultiple': found a Maybe, which means the vectors either have length zero, or one of them is all Zero!"     g (Yes x) = Just x+++--------------------------------------------------------------------------------+-- * Helper functions specific to two and three dimensional vectors++xComponent :: (1 <= d, Arity d) => Lens' (Vector d r) r+xComponent = element (C :: C 0)+{-# INLINABLE xComponent #-}++yComponent :: (2 <= d, Arity d) => Lens' (Vector d r) r+yComponent = element (C :: C 1)+{-# INLINABLE yComponent #-}++zComponent :: (3 <= d, Arity d) => Lens' (Vector d r) r+zComponent = element (C :: C 2)+{-# INLINABLE zComponent #-}

@@ -1,5 +1,17 @@ {-# LANGUAGE ScopedTypeVariables  #-} {-# LANGUAGE UndecidableInstances #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.Vector.VectorFamily+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Implementation of \(d\)-dimensional vectors. The implementation+-- automatically selects an optimized representation for small (up to size 4)+-- vectors.+--+-------------------------------------------------------------------------------- module Data.Geometry.Vector.VectorFamily where  import           Control.DeepSeq@@ -14,7 +26,6 @@                                                         , VectorFamilyF                                                         , ImplicitArity                                                         )-import           Data.Semigroup import qualified Data.Vector.Fixed as V import           Data.Vector.Fixed.Cont (Peano) import           GHC.TypeLits@@ -24,7 +35,10 @@ import qualified Linear.V3 as L3 import qualified Linear.V4 as L4 import           Linear.Vector-+import           Text.ParserCombinators.ReadP (ReadP, string,pfail)+import           Text.ParserCombinators.ReadPrec (lift)+import           Text.Read (Read(..),readListPrecDefault, readPrec_to_P,minPrec)+import           Data.Proxy  -------------------------------------------------------------------------------- -- * d dimensional Vectors@@ -57,7 +71,10 @@  deriving instance Arity d => Additive (Vector d) deriving instance Arity d => Metric (Vector d)-deriving instance Arity d => Affine (Vector d)+instance Arity d => Affine (Vector d) where+  type Diff (Vector d) = Vector d+  u .-. v = u ^-^ v+  p .+^ v = p ^+^ v  instance Arity d => Ixed (Vector d r) where   ix = element'@@ -71,6 +88,18 @@   show v = mconcat [ "Vector", show $ F.length v , " "                    , show $ F.toList v ] +instance (Read r, Arity d) => Read (Vector d r) where+  readPrec     = lift readVec+  readListPrec = readListPrecDefault++readVec :: forall d r. (Arity d, Read r) => ReadP (Vector d r)+readVec = do let d = natVal (Proxy :: Proxy d)+             _  <- string $ "Vector" <> show d <> " "+             rs <- readPrec_to_P readPrec minPrec+             case vectorFromList rs of+               Just v -> pure v+               _      -> pfail+ deriving instance (FromJSON r, Arity d) => FromJSON (Vector d r) instance (ToJSON r, Arity d) => ToJSON (Vector d r) where   toJSON     = toJSON . _unV@@ -83,18 +112,23 @@  pattern Vector   :: VectorFamilyF (Peano d) r -> Vector d r pattern Vector v = MKVector (VectorFamily v)+{-# COMPLETE Vector #-}  pattern Vector1   :: r -> Vector 1 r pattern Vector1 x = (Vector (Identity x))+{-# COMPLETE Vector1 #-}  pattern Vector2     :: r -> r -> Vector 2 r pattern Vector2 x y = (Vector (L2.V2 x y))+{-# COMPLETE Vector2 #-}  pattern Vector3        :: r -> r -> r -> Vector 3 r pattern Vector3 x y z  = (Vector (L3.V3 x y z))+{-# COMPLETE Vector3 #-}  pattern Vector4         :: r -> r -> r -> r -> Vector 4 r pattern Vector4 x y z w = (Vector (L4.V4 x y z w))+{-# COMPLETE Vector4 #-}  -------------------------------------------------------------------------------- 

@@ -9,13 +9,9 @@ -- 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 qualified Data.Vector.Fixed.Cont as Cont-import           Data.Vector.Fixed.Cont (Peano(..), PeanoNum(..), Fun(..))+import           Data.Vector.Fixed.Cont (PeanoNum(..), Fun(..)) import           GHC.TypeLits import           Linear.Affine (Affine(..)) import           Linear.Metric

@@ -10,7 +10,6 @@ import qualified Data.Vector.Fixed as V import           Data.Vector.Fixed (Arity) import           Data.Vector.Fixed.Boxed-import           Data.Vector.Fixed.Cont (Peano, PeanoNum(..)) import           GHC.Generics (Generic) import           GHC.TypeLits import           Linear.Affine (Affine(..))

@@ -0,0 +1,310 @@+{-# LANGUAGE ScopedTypeVariables #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.LSeq+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+-- Description :  Wrapper around Data.Sequence with type level length annotation.+--+--------------------------------------------------------------------------------+module Data.LSeq( LSeq+                , toSeq+                , empty+                , fromList+                , fromNonEmpty+                , fromSeq+                , toNonEmpty++                , (<|), (|>)+                , (><)+                , eval++                , index+                , adjust+                , partition+                , mapWithIndex+                , take+                , drop+                , unstableSort, unstableSortBy+                , head, last+                , append++                , ViewL(..)+                , viewl+                , pattern (:<|)++                , pattern (:<<)+                , pattern EmptyL++                , ViewR(..)+                , viewr+                , pattern (:|>)+++                , promise+                , forceLSeq+                ) where++import           Control.DeepSeq+import           Control.Lens ((%~), (&), (<&>), (^?), bimap)+import           Control.Lens.At (Ixed(..), Index, IxValue)+import qualified Data.Foldable as F+import qualified Data.List.NonEmpty as NonEmpty+import           Data.Maybe (fromJust)+import qualified Data.Sequence as S+import           GHC.Generics (Generic)+import qualified Data.Traversable as Tr+import           GHC.TypeLits+import           Prelude hiding (drop,take,head,last)++--------------------------------------------------------------------------------++-- $setup+-- >>> :{+-- import Data.Proxy+-- :}++++-- | LSeq n a certifies that the sequence has *at least* n items+newtype LSeq (n :: Nat) a = LSeq (S.Seq a)+                          deriving (Show,Read,Eq,Ord,Foldable,Functor,Traversable+                                   ,Generic,NFData)++toSeq          :: LSeq n a -> S.Seq a+toSeq (LSeq s) = s++instance Semigroup (LSeq n a) where+  (LSeq s) <> (LSeq s') = LSeq (s <> s')++instance Monoid (LSeq 0 a) where+  mempty = empty+  mappend = (<>)++++type instance Index   (LSeq n a) = Int+type instance IxValue (LSeq n a) = a+instance Ixed (LSeq n a) where+  ix i f s@(LSeq xs)+    | 0 <= i && i < S.length xs = f (S.index xs i) <&> \x -> LSeq $ S.update i x xs+    | otherwise                 = pure s++empty :: LSeq 0 a+empty = LSeq S.empty++(<|) :: a -> LSeq n a -> LSeq (1 + n) a+x <| xs = LSeq (x S.<| toSeq xs)++(|>)    :: LSeq n a -> a -> LSeq (1 + n) a+xs |> x = LSeq (toSeq xs S.|> x)++infixr 5 <|+infixl 5 |>++(><) :: LSeq n a -> LSeq m a -> LSeq (n + m) a+xs >< ys = LSeq (toSeq xs <> toSeq ys)++infix 5 ><+++eval :: forall proxy n m a. KnownNat n => proxy n -> LSeq m a -> Maybe (LSeq n a)+eval n (LSeq xs)+  | toInteger (S.length xs) >= natVal n = Just $ LSeq xs+  | otherwise                           = Nothing++++++-- | Promises that the length of this LSeq is actually n. This is not+-- checked.+--+-- This function should be a noop+promise :: LSeq m a -> LSeq n a+promise = LSeq . toSeq+++-- | Forces the first n elements of the LSeq+forceLSeq   :: KnownNat n => proxy n -> LSeq m a -> LSeq n a+forceLSeq n = promise . go (fromInteger $ natVal n)+  where+    -- forces the Lseq for n' positions+    go      :: Int -> LSeq m a -> LSeq m a+    go n' s | n' <= l    = s+            | otherwise  = error msg+      where+        l   = S.length . S.take n' . toSeq $ s+        msg = "forceLSeq: too few elements. expected " <> show n' <> " but found " <> show l++++toNonEmpty :: LSeq (1 + n) a -> NonEmpty.NonEmpty a+toNonEmpty = NonEmpty.fromList . F.toList++-- | appends two sequences.+--+append         :: LSeq n a -> LSeq m a -> LSeq (n + m) a+sa `append` sb = LSeq $ (toSeq sa) <> toSeq sb++--------------------------------------------------------------------------------++-- | get the element with index i, counting from the left and starting at 0.+-- O(log(min(i,n-i)))+index     :: LSeq n a -> Int -> a+index s i = fromJust $ s^?ix i++adjust       :: (a -> a) -> Int -> LSeq n a -> LSeq n a+adjust f i s = s&ix i %~ f+++partition   :: (a -> Bool) -> LSeq n a -> (LSeq 0 a, LSeq 0 a)+partition p = bimap LSeq LSeq . S.partition p . toSeq++mapWithIndex   :: (Int -> a -> b) -> LSeq n a -> LSeq n b+mapWithIndex f = wrapUnsafe (S.mapWithIndex f)++take   :: Int -> LSeq n a -> LSeq 0 a+take i = wrapUnsafe (S.take i)++drop   :: Int -> LSeq n a -> LSeq 0 a+drop i = wrapUnsafe (S.drop i)+++unstableSortBy   :: (a -> a -> Ordering) -> LSeq n a -> LSeq n a+unstableSortBy f = wrapUnsafe (S.unstableSortBy f)++unstableSort :: Ord a => LSeq n a -> LSeq n a+unstableSort = wrapUnsafe (S.unstableSort)+++wrapUnsafe :: (S.Seq a -> S.Seq b) -> LSeq n a -> LSeq m b+wrapUnsafe f = LSeq . f . toSeq++--------------------------------------------------------------------------------++fromSeq :: S.Seq a -> LSeq 0 a+fromSeq = LSeq++fromList :: Foldable f => f a -> LSeq 0 a+fromList = LSeq . S.fromList . F.toList++fromNonEmpty :: NonEmpty.NonEmpty a -> LSeq 1 a+fromNonEmpty = LSeq . S.fromList . F.toList+++--------------------------------------------------------------------------------++data ViewL n a where+  (:<) :: a -> LSeq n a -> ViewL (1 + n) a++infixr 5 :<++instance Semigroup (ViewL n a) where+  (x :< xs) <> (y :< ys) = x :< LSeq (toSeq xs <> (y S.<| toSeq ys))++deriving instance Show a => Show (ViewL n a)+instance Functor (ViewL n) where+  fmap = Tr.fmapDefault+instance Foldable (ViewL n) where+  foldMap = Tr.foldMapDefault+instance Traversable (ViewL n) where+  traverse f (x :< xs) = (:<) <$> f x <*> traverse f xs+instance Eq a => Eq (ViewL n a) where+  s == s' = F.toList s == F.toList s'+instance Ord a => Ord (ViewL n a) where+  s `compare` s' = F.toList s `compare` F.toList s'+++viewl :: LSeq (1 + n) a -> ViewL (1 + n) a+viewl xs = let ~(x S.:< ys) = S.viewl $ toSeq xs in x :< LSeq ys++viewl'    :: LSeq (1 + n) a -> (a, LSeq n a)+viewl' xs = let ~(x S.:< ys) = S.viewl $ toSeq xs in (x,LSeq ys)++infixr 5 :<|++pattern (:<|)    :: a -> LSeq n a -> LSeq (1 + n) a+pattern x :<| xs <- (viewl' -> (x,xs)) -- we need the coerce unfortunately+  where+    x :<| xs = x <| xs+{-# COMPLETE (:<|) #-}++++infixr 5 :<<++pattern (:<<)    :: a -> LSeq 0 a -> LSeq n a+pattern x :<< xs <- (viewLSeq -> Just (x,xs))++pattern EmptyL   :: LSeq n a+pattern EmptyL   <- (viewLSeq -> Nothing)++viewLSeq          :: LSeq n a -> Maybe (a,LSeq 0 a)+viewLSeq (LSeq s) = case S.viewl s of+                      S.EmptyL    -> Nothing+                      (x S.:< xs) -> Just (x,LSeq xs)+++--------------------------------------------------------------------------------++data ViewR n a where+  (:>) :: LSeq n a -> a -> ViewR (1 + n) a++infixl 5 :>++instance Semigroup (ViewR n a) where+  (xs :> x) <> (ys :> y) = LSeq ((toSeq xs S.|> x) <> toSeq ys) :> y++deriving instance Show a => Show (ViewR n a)+instance Functor (ViewR n) where+  fmap = Tr.fmapDefault+instance Foldable (ViewR n) where+  foldMap = Tr.foldMapDefault+instance Traversable (ViewR n) where+  traverse f (xs :> x) = (:>) <$> traverse f xs <*> f x+instance Eq a => Eq (ViewR n a) where+  s == s' = F.toList s == F.toList s'+instance Ord a => Ord (ViewR n a) where+  s `compare` s' = F.toList s `compare` F.toList s'++viewr    :: LSeq (1 + n) a -> ViewR (1 + n) a+viewr xs = let ~(ys S.:> x) = S.viewr $ toSeq xs in LSeq ys :> x++viewr'    :: LSeq (1 + n) a -> (LSeq n a, a)+viewr' xs = let ~(ys S.:> x) = S.viewr $ toSeq xs in (LSeq ys, x)++infixl 5 :|>++pattern (:|>)    :: forall n a. LSeq n a -> a -> LSeq (1 + n) a+pattern xs :|> x <- (viewr' -> (xs,x))+  where+    xs :|> x = xs |> x+{-# COMPLETE (:|>) #-}++--------------------------------------------------------------------------------++-- | Gets the first element of the LSeq+--+-- >>> head $ forceLSeq (Proxy :: Proxy 3) $ fromList [1,2,3]+-- 1+head           :: LSeq (1 + n) a -> a+head (x :<| _) = x++-- s = let (x :< _) = viewl s in x++-- | Get the last element of the LSeq+--+-- >>> last $ forceLSeq (Proxy :: Proxy 3) $ fromList [1,2,3]+-- 3+last           :: LSeq (1 + n) a -> a+last (_ :|> x) = x++-- testL = (eval (Proxy :: Proxy 2) $ fromList [1..5])++-- testL' :: LSeq 2 Integer+-- testL' = fromJust testL++-- test            :: Show a => LSeq (1 + n) a -> String+-- test (x :<| xs) = show x ++ show xs