hgeometry 0.10.0.0 → 0.11.0.0

raw patch · 74 files changed

+3332/−677 lines, 74 filesdep +hashabledep +randomdep ~hgeometry-combinatorialPVP ok

version bump matches the API change (PVP)

Dependencies added: hashable, random

Dependency ranges changed: hgeometry-combinatorial

API changes (from Hackage documentation)

Files

• changelog.org +34/−0
• doctests.hs +5/−0
• hgeometry.cabal +48/−4
• src/Algorithms/Geometry/ConvexHull/DivideAndConquer.hs +0/−1
• src/Algorithms/Geometry/ConvexHull/GrahamScan.hs +72/−3
• src/Algorithms/Geometry/ConvexHull/JarvisMarch.hs +145/−0
• src/Algorithms/Geometry/ConvexHull/Naive.hs +93/−0
• src/Algorithms/Geometry/DelaunayTriangulation/DivideAndConquer.hs +3/−2
• src/Algorithms/Geometry/LineSegmentIntersection/BentleyOttmann.hs +0/−1
• src/Algorithms/Geometry/LineSegmentIntersection/Types.hs +2/−0
• src/Algorithms/Geometry/LinearProgramming/LP2DRIC.hs +2/−2
• src/Algorithms/Geometry/LinearProgramming/Types.hs +5/−5
• src/Algorithms/Geometry/PolyLineSimplification/DouglasPeucker.hs +6/−6
• src/Algorithms/Geometry/SmallestEnclosingBall/Naive.hs +4/−4
• src/Algorithms/Geometry/SmallestEnclosingBall/RIC.hs +27/−10
• src/Algorithms/Geometry/SoS.hs +241/−0
• src/Algorithms/Geometry/SoS/AsPoint.hs +29/−0
• src/Algorithms/Geometry/SoS/Determinant.hs +13/−0
• src/Algorithms/Geometry/SoS/Expr.hs +77/−0
• src/Algorithms/Geometry/SoS/Internal.hs +28/−0
• src/Algorithms/Geometry/SoS/Orientation.hs +83/−0
• src/Algorithms/Geometry/SoS/Sign.hs +30/−0
• src/Algorithms/Geometry/SoS/Symbolic.hs +352/−0
• src/Algorithms/Geometry/WellSeparatedPairDecomposition/Types.hs +1/−0
• src/Algorithms/Geometry/WellSeparatedPairDecomposition/WSPD.hs +2/−2
• src/Data/Geometry/Arrangement/Internal.hs +5/−8
• src/Data/Geometry/Ball.hs +9/−3
• src/Data/Geometry/BezierSpline.hs +169/−0
• src/Data/Geometry/Boundary.hs +7/−2
• src/Data/Geometry/Box.hs +4/−46
• src/Data/Geometry/Box/Corners.hs +70/−0
• src/Data/Geometry/Box/Internal.hs +19/−24
• src/Data/Geometry/Box/Sides.hs +91/−0
• src/Data/Geometry/Directions.hs +49/−0
• src/Data/Geometry/Ellipse.hs +56/−0
• src/Data/Geometry/HalfLine.hs +5/−1
• src/Data/Geometry/HalfSpace.hs +9/−2
• src/Data/Geometry/HyperPlane.hs +41/−1
• src/Data/Geometry/Interval.hs +41/−13
• src/Data/Geometry/IntervalTree.hs +6/−8
• src/Data/Geometry/Line.hs +1/−1
• src/Data/Geometry/Line/Internal.hs +27/−20
• src/Data/Geometry/LineSegment.hs +26/−5
• src/Data/Geometry/Matrix.hs +79/−0
• src/Data/Geometry/Matrix/Internal.hs +14/−0
• src/Data/Geometry/PlanarSubdivision.hs +0/−1
• src/Data/Geometry/PlanarSubdivision/Basic.hs +15/−8
• src/Data/Geometry/Point.hs +9/−392
• src/Data/Geometry/Point/Class.hs +66/−0
• src/Data/Geometry/Point/Internal.hs +252/−0
• src/Data/Geometry/Point/Orientation.hs +31/−0
• src/Data/Geometry/Point/Orientation/Degenerate.hs +150/−0
• src/Data/Geometry/Point/Quadrants.hs +69/−0
• src/Data/Geometry/PolyLine.hs +45/−7
• src/Data/Geometry/Polygon/Core.hs +11/−8
• src/Data/Geometry/PrioritySearchTree.hs +2/−0
• src/Data/Geometry/QuadTree.hs +203/−0
• src/Data/Geometry/QuadTree/Cell.hs +141/−0
• src/Data/Geometry/QuadTree/Quadrants.hs +16/−0
• src/Data/Geometry/QuadTree/Split.hs +32/−0
• src/Data/Geometry/QuadTree/Tree.hs +116/−0
• src/Data/Geometry/RangeTree.hs +1/−1
• src/Data/Geometry/RangeTree/Generic.hs +2/−0
• src/Data/Geometry/RangeTree/Measure.hs +1/−1
• src/Data/Geometry/SegmentTree/Generic.hs +7/−5
• src/Data/Geometry/SubLine.hs +1/−1
• src/Data/Geometry/Transformation.hs +22/−51
• src/Data/Geometry/Triangle.hs +32/−9
• src/Data/Geometry/Vector.hs +25/−8
• src/Data/Geometry/Vector/VectorFamily.hs +19/−1
• src/Data/Geometry/Vector/VectorFamilyPeano.hs +16/−2
• src/Data/PlaneGraph/Core.hs +15/−6
• src/Graphics/Camera.hs +2/−1
• test/Data/Geometry/arrangement.ipe.out.ipe +1/−1

@@ -2,6 +2,40 @@  * Changelog +** 0.11++- Removed Functor instance from Triangle and replaced it with Bifunctor/Bifoldable/Bitraversable+- Testing if a point lies above/below a line is now in a typeclass,+  moreover there now is also an instance of this typeclass for+  planes. Hence, we can test if a point in R^3 lies above or below a+  plane.+- Bugfixes in the incomingEdges and outgoingEdges functions in+  Planar/Plane graphs and Planar subdivisions+- Added separate data types for Sides and Corners of Rectangles.+- More functionality for working with Halfspaces+- Fixed a bug in computing the intersection of overlapping+  linesegments+- PolyLine.fromPoints now returns a Maybe PolyLine rather than a+  Polyine. Use fromPointsUnsafe for the old behavior.+- Interval now no longer exports its constructor. Use the provided+  patterns instead.+- Added an OpenLineSegment pattern/constructor+- The corners and sides functions in Box now return specific types+  representing those rather than four tuples.+- Added a BezierSpline module and data type (Thanks to Maarten).+- Added a QuadTree implementation. It can be built from a set of+  points, and to represent the zeroset of some function.+- Added a Naive implementation of Convex hull in R^3. Note however+  that it works only for points in general position. In particular, no+  four points should be coplanar.+- Added a Data.Geometry.Directions module that defines cardinal and+  InterCardinal directions.+- Added an Ellipse type (mostly so that hgeometry-ipe can read+  ellipses)+- Added FunctorWithIndex, FoldableWithIndex, and TraversableWithIndex+  instances for Vector, and removed specifically exporting imap; we+  can now just use those functions from the Lens package.+ ** 0.10  - renamed the smallest enclosing ball to RIC

@@ -63,6 +63,11 @@   , "Data.Geometry.Polygon"   , "Data.Geometry.Ball"   , "Data.Geometry.Box"+  , "Data.Geometry.HyperPlane"    -- , "Algorithms.Geometry.HiddenSurfaceRemoval.HiddenSurfaceRemoval"+  , "Algorithms.Geometry.ConvexHull.Naive"+  , "Algorithms.Geometry.ConvexHull.JarvisMarch"++  , "Algorithms.Geometry.SoS.Orientation"   ]

@@ -1,5 +1,5 @@ name:                hgeometry-version:             0.10.0.0+version:             0.11.0.0 synopsis:            Geometric Algorithms, Data structures, and Data types. description:   HGeometry provides some basic geometry types, and geometric algorithms and@@ -56,12 +56,17 @@   ghc-options: -O2 -Wall -fno-warn-unticked-promoted-constructors -fno-warn-type-defaults    exposed-modules:+                    -- * Primitives; Simulating General Position+                    Algorithms.Geometry.SoS+                    Algorithms.Geometry.SoS.Symbolic+                     -- * Generic Geometry                     Data.Geometry                     Data.Geometry.Properties                     Data.Geometry.Transformation                     Data.Geometry.Boundary                     Data.Geometry.Duality+                    Data.Geometry.Directions                      -- * Basic Geometry Types                     Data.Geometry.Vector@@ -69,10 +74,13 @@                     Data.Geometry.Vector.VectorFamily                     Data.Geometry.Vector.VectorFamilyPeano +                    Data.Geometry.Matrix+                     -- Data.Geometry.Vector.Vinyl                     Data.Geometry.Interval                     Data.Geometry.Interval.Util                     Data.Geometry.Point+                     Data.Geometry.Line                     Data.Geometry.Line.Internal                     Data.Geometry.LineSegment@@ -86,10 +94,17 @@                     Data.Geometry.Slab                     Data.Geometry.Box                     Data.Geometry.Box.Internal+                    Data.Geometry.Box.Sides+                    Data.Geometry.Box.Corners+                     Data.Geometry.Ball+                    Data.Geometry.Ellipse+                     Data.Geometry.Polygon                     Data.Geometry.Polygon.Convex +                    Data.Geometry.BezierSpline+                     -- * Geometric Data Structures                     Data.Geometry.IntervalTree                     Data.Geometry.SegmentTree@@ -111,13 +126,21 @@                      Data.Geometry.PrioritySearchTree +                    Data.Geometry.QuadTree+                    Data.Geometry.QuadTree.Cell+                    Data.Geometry.QuadTree.Quadrants+                    Data.Geometry.QuadTree.Split+                    Data.Geometry.QuadTree.Tree++                     -- * Algorithms                      -- * Geometric Algorithms                     Algorithms.Geometry.ConvexHull.GrahamScan                     Algorithms.Geometry.ConvexHull.DivideAndConquer                     Algorithms.Geometry.ConvexHull.QuickHull-                    -- Algorithms.Geometry.ConvexHull.JarvisMarch+                    Algorithms.Geometry.ConvexHull.JarvisMarch+                    Algorithms.Geometry.ConvexHull.Naive                      Algorithms.Geometry.LowerEnvelope.DualCH @@ -174,14 +197,32 @@                     Graphics.Render    other-modules:+                    Data.Geometry.Matrix.Internal+                     -- * Implementation Internals of Polygons                     Data.Geometry.Polygon.Core                     Data.Geometry.Polygon.Extremes ++                    Data.Geometry.Point.Internal+                    Data.Geometry.Point.Orientation+                    Data.Geometry.Point.Quadrants+                    Data.Geometry.Point.Orientation.Degenerate+                    Data.Geometry.Point.Class++                    Algorithms.Geometry.SoS.Expr+                    Algorithms.Geometry.SoS.AsPoint+                    Algorithms.Geometry.SoS.Internal+                    Algorithms.Geometry.SoS.Orientation+                    Algorithms.Geometry.SoS.Determinant+                    Algorithms.Geometry.SoS.Sign+++   -- other-extensions:   build-depends:                 base                    >= 4.11      &&     < 5-              , hgeometry-combinatorial >= 0.10.0.0+              , hgeometry-combinatorial >= 0.11.0.0                , bifunctors              >= 4.1               , bytestring              >= 0.10@@ -198,10 +239,13 @@               , deepseq                 >= 1.1               , fingertree              >= 0.1               , MonadRandom             >= 0.5+              , random                  >= 1.1               , QuickCheck              >= 2.5               , quickcheck-instances    >= 0.3               , reflection              >= 2.1               , primitive               >= 0.6.3.0+              , hashable                >= 1.2+               -- , singleton-typelits      >= 0.1.0.0                -- , ghc-typelits-natnormalise >= 0.6@@ -220,7 +264,7 @@               , hspec, QuickCheck, quickcheck-instances  -  hs-source-dirs: src test+  hs-source-dirs: src    default-language:    Haskell2010 

@@ -16,7 +16,6 @@  import           Algorithms.DivideAndConquer import           Control.Arrow ((&&&))-import           Control.Lens ((^.), to) import           Data.Ext import           Data.Geometry.Point import           Data.Geometry.Polygon

@@ -1,6 +1,8 @@ module Algorithms.Geometry.ConvexHull.GrahamScan( convexHull-                                                , upperHull-                                                , lowerHull+                                                , upperHull, upperHull'+                                                , lowerHull, lowerHull'++                                                , upperHullFromSorted, upperHullFromSorted'                                                 ) where  import           Control.Lens ((^.))@@ -23,13 +25,54 @@                        in ConvexPolygon . fromPoints . reverse $ lh ++ uh  -- | Computes the upper hull. The upper hull is given from left to right.+--+-- Specifically. A pair of points defines an edge of the upper hull+-- iff all other points are strictly to the right of its supporting+-- line.+--+-- Note that this definition implies that the segment may be+-- vertical. Use 'upperHull'' if such an edge should not be reported.+--+-- running time: \(O(n\log n)\) upperHull  :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p) upperHull = NonEmpty.reverse . hull id --- | Computes the upper hull. The upper hull is given from left to right+-- | Computes the upper hull, making sure that there are no vertical segments.+--+-- The upper hull is given from left to right+--+upperHull'  :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)+upperHull' = NonEmpty.reverse . dropVertical . hull id++-- | Helper function to remove vertical segments from the hull.+--+-- Tests if the first two points are on a vertical line, if so removes+-- the first point.+dropVertical :: Eq r => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)+dropVertical = \case+  h@(_ :| [])                                            -> h+  h@(p :| (q : rest)) | p^.core.xCoord == q^.core.xCoord -> q :| rest+                      | otherwise                        -> h+++-- | Computes the upper hull. The upper hull is given from left to right.+--+-- Specifically. A pair of points defines an edge of the lower hull+-- iff all other points are strictly to the left of its supporting+-- line.+--+-- Note that this definition implies that the segment may be+-- vertical. Use 'lowerHull'' if such an edge should not be reported.+--+-- running time: \(O(n\log n)\) lowerHull :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p) lowerHull = hull reverse +-- | Computes the lower hull, making sure there are no vertical+-- segments. (Note that the only such segment could be the first+-- segment).+lowerHull' :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)+lowerHull' = dropVertical . hull reverse  -- | Helper function so that that can compute both the upper or the lower hull, depending -- on the function f@@ -43,6 +86,32 @@ incXdecY  :: Ord r => (Point 2 r) :+ p -> (Point 2 r) :+ q -> Ordering incXdecY (Point2 px py :+ _) (Point2 qx qy :+ _) =   compare px qx <> compare qy py+++-- | Given a sequence of points that is sorted on increasing+-- x-coordinate and decreasing y-coordinate, computes the upper+-- hull, in *right to left order*.+--+-- Specifically. A pair of points defines an edge of the upper hull+-- iff all other points are strictly to the right of its supporting+-- line.+--+--+-- Note that In constrast to the 'upperHull' function, the result is+-- returned *from right to left* !!!+--+-- running time: \(O(n)\).+upperHullFromSorted :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)+upperHullFromSorted = \case+  h@(_ :| [])  -> h+  pts          -> hull' $ NonEmpty.toList pts++-- | Computes the upper hull from a sorted input. Removes the last vertical segment.+--+--+-- running time: \(O(n)\).+upperHullFromSorted' :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)+upperHullFromSorted' = dropVertical . upperHullFromSorted   -- | Precondition: The list of input points is sorted

@@ -0,0 +1,145 @@+module Algorithms.Geometry.ConvexHull.JarvisMarch(+    convexHull++  , upperHull, upperHull'+  , lowerHull, lowerHull'+  , steepestCcwFrom, steepestCwFrom+  ) where++import           Control.Lens ((^.))+import           Data.Bifunctor+import           Data.Either (either)+import           Data.Ext+import           Data.Foldable+import           Data.Geometry.Point+import           Data.Geometry.Polygon+import           Data.Geometry.Polygon.Convex (ConvexPolygon(..))+import           Data.Geometry.Vector+import qualified Data.List as List+import           Data.List.NonEmpty (NonEmpty(..), (<|))+import qualified Data.List.NonEmpty as NonEmpty+import           Data.Ord (comparing, Down(..))+import           Data.Semigroup.Foldable++--------------------------------------------------------------------------------++-- | Compute the convexhull using JarvisMarch. The resulting polygon+-- is given in clockwise order.+--+-- running time: \(O(nh)\), where \(n\) is the number of input points+-- and \(h\) is the complexity of the hull.+convexHull            :: (Ord r, Num r)+                      => NonEmpty (Point 2 r :+ p) -> ConvexPolygon p r+convexHull (p :| []) = ConvexPolygon . fromPoints $ [p]+convexHull pts       = ConvexPolygon . fromPoints $ uh <> reverse lh+  where+    lh = case NonEmpty.nonEmpty (NonEmpty.init $ lowerHull pts) of+           Nothing       -> []+           Just (_:|lh') -> lh'+    uh = toList $ upperHull pts++                       -- note that fromList is afe since ps contains at least two elements+  -- where+  --   SP p@(c :+ _) pts = minViewBy incXdecY ps+  --   takeWhile' pf (x :| xs) = x : takeWhile pf xs++upperHull     ::  (Num r, Ord r) =>  NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)+upperHull pts = repeatedly cmp steepestCwFrom s rest+  where+    (s:_ :+ rest) = extractMinimaBy cmp (NonEmpty.toList pts)+    cmp           = comparing (\(Point2 x y :+ _) -> (x, Down y))+                    -- start from the topmost point that has minimum x-coord+                    -- also use cmp as the comparator, so that we also select the last+                    -- vertical segment.++-- | Upepr hull from left to right, without any vertical segments.+upperHull'     ::  (Num r, Ord r) =>  NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)+upperHull' pts = pruneVertical $ repeatedly cmp steepestCwFrom s rest+  where+    (s:_ :+ rest) = extractMinimaBy cmp0 (NonEmpty.toList pts)+    cmp0          = comparing (\(Point2 x y :+ _) -> (x, Down y))+                    -- start from the topmost point that has minimum x-coord+    cmp           = comparing (^.core)+                    -- for the rest select them in normal+                    -- lexicographic order, this causes the last+                    -- vertical segment to be ignored.++-- | Computes the lower hull, from left to right. Includes vertical+-- segments at the start.+--+-- running time: \(O(nh)\), where \(h\) is the complexity of the hull.+lowerHull     ::  (Num r, Ord r) =>  NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)+lowerHull pts = pruneVertical $ repeatedly cmp steepestCcwFrom s rest+  where+    (s:_ :+ rest) = extractMinimaBy cmp0 (NonEmpty.toList pts)+    cmp0          = comparing (\(Point2 x y :+ _) -> (x, Down y))+                    -- start from the topmost point that has minimum x-coord+    cmp           = comparing (^.core)+                    -- for the rest of the comparions use the normal+                    -- lexicographic comparing order.++-- | Jarvis March to compute the lower hull, without any vertical segments.+--+--+-- running time: \(O(nh)\), where \(h\) is the complexity of the hull.+lowerHull'     :: (Num r, Ord r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)+lowerHull' pts = pruneVertical $ repeatedly cmp steepestCcwFrom s rest+  where+    (s:_ :+ rest) = extractMinimaBy cmp (NonEmpty.toList pts)+    cmp           = comparing (^.core)+++-- | Find the next point in counter clockwise order, i.e. the point+-- with minimum slope w.r.t. the given point.+steepestCcwFrom   :: (Ord r, Num r)+               => (Point 2 r :+ a) -> NonEmpty (Point 2 r :+ b)  -> Point 2 r :+ b+steepestCcwFrom p = List.minimumBy (ccwCmpAroundWith (Vector2 0 (-1)) p)++-- | Find the next point in clockwise order, i.e. the point+-- with maximum slope w.r.t. the given point.+steepestCwFrom   :: (Ord r, Num r)+               => (Point 2 r :+ a) -> NonEmpty (Point 2 r :+ b)  -> Point 2 r :+ b+steepestCwFrom p = List.minimumBy (cwCmpAroundWith (Vector2 0 1) p)++repeatedly       :: (a -> a -> Ordering) -> (a -> NonEmpty a -> a) -> a -> [a] -> NonEmpty a+repeatedly cmp f = go+  where+    go m xs' = case NonEmpty.nonEmpty xs' of+      Nothing -> m :| []+      Just xs -> let p = f m xs+                 in m <| go p (NonEmpty.filter (\x -> p `cmp` x == LT) xs)+++-- | Removes the topmost vertical points, if they exist+pruneVertical :: Eq r => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)+pruneVertical = either id id . foldr1With f (\q -> Left $ q:|[])+  where+    f p = \case+      Left (q:|qs) | p^.core.xCoord == q^.core.xCoord -> Left  (p :| qs)+                   | otherwise                        -> Right (p :| q:qs)+      Right pts                                       -> Right (p <| pts)++-- | Foldr, but start by applying some function on the rightmost+-- element to get the starting value.+foldr1With     :: Foldable1 f => (a -> b -> b) -> (a -> b) -> f a -> b+foldr1With f b = go . toNonEmpty+  where+    go (x :| xs) = case NonEmpty.nonEmpty xs of+                     Nothing  -> b x+                     Just xs' -> x `f` (go xs')++-- | extracts all minima from the list. The result consists of the+-- list of minima, and all remaining points. Both lists are returned+-- in the order in which they occur in the input.+--+-- >>> extractMinimaBy compare [1,2,3,0,1,2,3,0,1,2,0,2]+-- [0,0,0] :+ [2,3,1,2,3,1,2,1,2]+extractMinimaBy     :: (a -> a -> Ordering) -> [a] -> [a] :+ [a]+extractMinimaBy cmp = \case+  []     -> [] :+ []+  (x:xs) -> first NonEmpty.toList $ foldr (\y (mins@(m:|_) :+ rest) ->+                                             case m `cmp` y of+                                               LT -> mins :+ y:rest+                                               EQ -> (y NonEmpty.<| mins) :+ rest+                                               GT -> (y:|[]) :+ NonEmpty.toList mins <> rest+                                          ) ((x:|[]) :+ []) xs

@@ -0,0 +1,93 @@+module Algorithms.Geometry.ConvexHull.Naive( ConvexHull+                                           , lowerHull', lowerHullAll++                                           , isValidTriangle, upperHalfSpaceOf+                                           ) where++import           Control.Lens+import           Data.Ext+import           Data.Foldable (toList)+import           Data.Geometry.HalfSpace+import           Data.Geometry.HyperPlane+import           Data.Geometry.Line+import           Data.Geometry.Point+import           Data.Geometry.Triangle+import           Data.Geometry.Vector+import           Data.Intersection(intersects)+import qualified Data.List as List+import           Data.List.NonEmpty (NonEmpty(..))+import qualified Data.List.NonEmpty as NonEmpty+import           Data.Maybe (listToMaybe, isNothing)+import           Data.Util+--------------------------------------------------------------------------------++type ConvexHull d p r = [Triangle 3 p r]++-- | Computes the lower hull without its vertical triangles.+--+-- pre: The points are in general position. In particular, no four+-- points should be coplanar.+--+-- running time: \(O(n^4)\)+lowerHull' :: forall r p. (Ord r, Fractional r, Show r)+           => NonEmpty (Point 3 r :+ p) -> ConvexHull 3 p r+lowerHull' = filter (not . isVertical) . lowerHullAll+  where+    isVertical (Triangle p q r) =+      ccw' (p&core %~ projectPoint) (q&core %~ projectPoint) (r&core %~ projectPoint) == CoLinear++-- | Generates a set of triangles to be used to construct a complete+-- convex hull. In particular, it may contain vertical triangles.+--+-- pre: The points are in general position. In particular, no four+-- points should be coplanar.+--+-- running time: \(O(n^4)\)+lowerHullAll                 :: forall r p. (Ord r, Fractional r, Show r)+                             => NonEmpty (Point 3 r :+ p) -> ConvexHull 3 p r+lowerHullAll (toList -> pts) = let mkT (Three p q r) = Triangle p q r in+    [ t | t <- mkT <$> uniqueTriplets pts, isNothing (isValidTriangle t pts) ]++++killOverlapping :: (Ord r, Fractional r) => [Triangle 3 p r] -> [Triangle 3 p r]+killOverlapping = foldr keepIfNotOverlaps []+  where+    keepIfNotOverlaps t ts | any (t `overlaps`) ts = ts+                           | otherwise             = t:ts+++t1@(Triangle p q r) `overlaps` t2@(Triangle a b c) = upperHalfSpaceOf t1 == upperHalfSpaceOf t2+                                                  && False++++-- | Tests if this is a valid triangle for the lower envelope. That+-- is, if all point lie above the plane through these points. Returns+-- a Maybe; if the result is a Nothing the triangle is valid, if not+-- it returns a counter example.+--+-- >>> let t = (Triangle (ext origin) (ext $ Point3 1 0 0) (ext $ Point3 0 1 0))+-- >>> isValidTriangle t [ext $ Point3 5 5 0]+-- Nothing+-- >>> let t = (Triangle (ext origin) (ext $ Point3 1 0 0) (ext $ Point3 0 1 0))+-- >>> isValidTriangle t [ext $ Point3 5 5 (-10)]+-- Just (Point3 [5,5,-10] :+ ())+isValidTriangle   :: (Num r, Ord r)+                  => Triangle 3 p r -> [Point 3 r :+ q] -> Maybe (Point 3 r :+ q)+isValidTriangle t = listToMaybe . filter (\a -> not $ (a^.core) `intersects` h)+  where+    h = upperHalfSpaceOf t+++-- | Computes the halfspace above the triangle.+--+-- >>> upperHalfSpaceOf (Triangle (ext $ origin) (ext $ Point3 10 0 0) (ext $ Point3 0 10 0))+-- HalfSpace {_boundingPlane = HyperPlane {_inPlane = Point3 [0,0,0], _normalVec = Vector3 [0,0,100]}}+upperHalfSpaceOf                  :: (Ord r, Num r) => Triangle 3 p r -> HalfSpace 3 r+upperHalfSpaceOf (Triangle p q r) = HalfSpace h+  where+    h' = from3Points (p^.core) (q^.core) (r^.core)+    c  = p&core.zCoord -~ 1+    h  = if (c^.core) `liesBelow` h' then h' else h'&normalVec %~ ((-1) *^)+    a `liesBelow` plane = (a `onSideUpDown` plane) == Below

@@ -8,21 +8,22 @@ 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 qualified Data.CircularSeq as CS 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.Geometry.Polygon.Convex as Convex 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           Data.Measured.Size import qualified Data.Vector as V  -------------------------------------------------------------------------------

@@ -26,7 +26,6 @@ import qualified Data.Map as M import           Data.Maybe import           Data.Ord (Down(..), comparing)-import           Data.OrdSeq (Compare) import qualified Data.Set as SS -- status struct import qualified Data.Set.Util as SS -- status struct import qualified Data.Set as EQ -- event queue

@@ -16,6 +16,8 @@  -------------------------------------------------------------------------------- +type Compare a = a -> a -> Ordering+ -- get the endpoints of a line segment endPoints'   :: (HasEnd s, HasStart s) => s -> (StartCore s, EndCore s) endPoints' s = (s^.start.core,s^.end.core)

@@ -86,8 +86,8 @@                            , _current :: !(Point d r)                            } -deriving instance (Arity d, Show r)   => Show    (LPState d r)-deriving instance (Arity d, Eq r)     => Eq      (LPState d r)+deriving instance (Arity d, Show r)             => Show    (LPState d r)+deriving instance (Arity d, Eq r, Fractional r) => Eq      (LPState d r)  obj     :: Lens' (LPState d r) (Vector d r) obj     = lens _obj     (\(LPState _ s p) o -> LPState o s p)

@@ -26,14 +26,14 @@                     | UnBounded (HalfLine d r) makePrisms ''LPSolution -deriving instance (Arity d, Show r)   => Show    (LPSolution d r)-deriving instance (Arity d, Eq r)     => Eq      (LPSolution d r)+deriving instance (Arity d, Show r)             => Show    (LPSolution d r)+deriving instance (Arity d, Eq r, Fractional r) => Eq      (LPSolution d r)  data LinearProgram d r = LinearProgram { _objective   :: !(Vector d r)                                        , _constraints :: [HalfSpace d r]                                        } makeLenses ''LinearProgram -deriving instance Arity d             => Functor (LinearProgram d)-deriving instance (Arity d, Show r)   => Show    (LinearProgram d r)-deriving instance (Arity d, Eq r)     => Eq      (LinearProgram d r)+deriving instance Arity d                       => Functor (LinearProgram d)+deriving instance (Arity d, Show r)             => Show    (LinearProgram d r)+deriving instance (Arity d, Fractional r, Eq r) => Eq      (LinearProgram d r)

@@ -22,15 +22,15 @@ douglasPeucker         :: (Ord r, Fractional r, Arity d)                        => r -> PolyLine d p r -> PolyLine d p r douglasPeucker eps pl-    | dst <= (eps*eps) = fromPoints [a,b]+    | dst <= (eps*eps) = fromPointsUnsafe [a,b] -- at least two points, so we are fine.     | otherwise        = douglasPeucker eps pref `merge` douglasPeucker eps subf   where-    pts     = pl^.points-    a       = LSeq.head pts-    b       = LSeq.last pts-    (i,dst)             = maxDist pts (ClosedLineSegment a b)+    pts         = pl^.points+    a           = LSeq.head pts+    b           = LSeq.last pts+    (i,dst)     = maxDist pts (ClosedLineSegment a b) -    (pref,subf)         = split i pl+    (pref,subf) = split i pl  -------------------------------------------------------------------------------- -- * Internal functions

@@ -22,8 +22,8 @@ import Data.List (minimumBy) import Data.Function (on) import Data.Maybe (fromMaybe)-import Data.Util(STR(..),SP(..), uniquePairs, uniqueTriplets)-+import Data.Util(uniquePairs, uniqueTriplets)+import qualified Data.Util as Util --------------------------------------------------------------------------------  -- | Horrible O(n^4) implementation that simply tries all disks, checks if they@@ -38,11 +38,11 @@  pairs     :: Fractional r => [Point 2 r :+ p] -> [DiskResult p r] pairs pts = [ DiskResult (fromDiameter (a^.core) (b^.core)) (Two a b)-            | SP a b <- uniquePairs pts]+            | Util.Two 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)-               | STR a b c <- uniqueTriplets pts]+               | Util.Three a b c <- uniqueTriplets pts]  disk'       :: (Ord r, Fractional r)             => Point 2 r :+ p -> Point 2 r :+ p -> Point 2 r :+ p -> Disk () r

@@ -31,6 +31,7 @@ import           Data.Ord (comparing) import           System.Random.Shuffle (shuffle) +import Data.RealNumber.Rational import Debug.Trace  --------------------------------------------------------------------------------@@ -149,16 +150,32 @@  -------------------------------------------------------------------------------- -test :: Maybe (DiskResult () Rational)-test = smallestEnclosingDiskWithPoints p q myPts-  where-    p = ext $ Point2 0 (-6)-    q = ext $ Point2 0 6+-- test :: Maybe (DiskResult () Rational)+-- test = smallestEnclosingDiskWithPoints p q myPts+--   where+--     p = ext $ Point2 0 (-6)+--     q = ext $ Point2 0 6  -myPts = map ext [Point2 5 1, Point2 3 3, Point2 (-2) 2, Point2 (-4) 5]+-- myPts = map ext [Point2 5 1, Point2 3 3, Point2 (-2) 2, Point2 (-4) 5] -disk'' r = (r:+) <$> disk (p^.core) (q^.core) (r^.core)-  where-    p = ext $ Point2 0 (-6)-    q = ext $ Point2 0 6+-- disk'' r = (r:+) <$> disk (p^.core) (q^.core) (r^.core)+--   where+--     p = ext $ Point2 0 (-6)+--     q = ext $ Point2 0 6+++-- maartenBug :: DiskResult () Double+-- maartenBug = let (p:q:rest) = maartenBug'+--              in smallestEnclosingDisk' p q rest++-- maartenBug' :: [Point 2 Double :+ ()]+-- maartenBug' = [ Point2 (7.2784424e-3) (249.23) :+ ()+--               , Point2 (-5.188493   ) (249.23) :+ ()+--               , Point2 (-10.382694  ) (249.23) :+ ()+--               , Point2 (-15.575621  ) (249.23) :+ ()+--               , Point2 (0.0         ) (249.23) :+ ()+--               , Point2 (0.0         ) (239.9031) :+ ()+--               , Point2 (0.0         ) (230.37791) :+ ()+--               , Point2 (0.0         ) (220.67882) :+ ()+--               ]

@@ -0,0 +1,241 @@+--------------------------------------------------------------------------------+-- |+-- Module      :  Algorithms.Geometry.SoS+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Implementation of+-- Simulation of Simplicity: A Technique to Cope with Degenerate Cases in Geometric Algorithms+--+-- By+-- Herbert Edelsbrunner and Ernst Peter Mucke+--+--------------------------------------------------------------------------------+module Algorithms.Geometry.SoS+  ( module Algorithms.Geometry.SoS.Sign+  , module Algorithms.Geometry.SoS.Orientation+  , module Algorithms.Geometry.SoS.Determinant+  ) where++-- import Algorithms.Geometry.SoS.Internal+import Algorithms.Geometry.SoS.Orientation+import Algorithms.Geometry.SoS.Determinant+import Algorithms.Geometry.SoS.Sign+import Control.CanAquire+import Control.Lens+import Data.Ext+import Data.Geometry.Point.Internal+import Data.Geometry.Properties+import Data.Geometry.Vector++--------------------------------------------------------------------------------++--------------------------------------------------------------------------------++++-- sideTest'             :: ( SoS p, Dimension p ~ 2, r ~ NumType p+--                          , Eq r, Num r+--                          ) => [p] -> Sign+-- sideTest' (q:p1:p2:_) = sideTest q (Vector2 p1 p2)+++++++--------------------------------------------------------------------------------+++----------------------------------------+--------------------------------------------------------------------------------++--------------------------------------------------------------------------------++++++-- instance (i `CanAquire` Point d r, Arity d) => P i d r `CanAquire` Point d (R i) where+--   aquire (P i) = Point $ pure ()+++++--------------------------------------------------------------------------------++++--------------------------------------------------------------------------------+++--------------------------------------------------------------------------------+++-- -- TODO: Remove this one+-- instance HasIndex (Point d r :+ Int) where+--   indexOf = view extra+++-- test1 :: Sign+-- test1 = sideTest (Point1 1 :+ 0 :: Point 1 Int :+ Int) (Vector1 $ Point1 5 :+ 1)++-- test2 :: Sign+-- test2 = sideTest (Point1 5 :+ 0 :: Point 1 Int :+ Int) (Vector1 $ Point1 5 :+ 1)+++-- test3 :: Sign+-- test3 = sideTest (Point2 (-1) 5 :+ 0 :: Point 2 Int :+ Int) (Vector2 (Point2 0 0  :+ 1)+--                                                                      (Point2 0 10 :+ 2)+--                                                             )+++-- pattern Point1 x = Point (Vector1 x)+++-- testV :: Sign+-- testV = simulateSimplicity sideTest' [ Point2 (-1) 5+--                                      , Point2 0 0+--                                      , Point2 0 10+--                                      ]++++++--------------------------------------------------------------------------------+++++++++-- cmpSignificance                   :: Ord k => Bag k -> Bag k -> Ordering+-- cmpSignificance (Bag e1) (Bag e2) = e1e2 `compare` e2e1+--   where+--     e1e2 = fmap fst . Map.lookupMax $ e1 `Map.difference` e2+--     e2e1 = fmap fst . Map.lookupMax $ e2 `Map.difference` e1++++-- -- | Represents a Sum of terms, i.e. a value that has the form:+-- --+-- -- \[+-- --   \sum c \Pi_{(i,j)} \varepsilon(i,j)+-- -- \]+-- newtype Symbolic i j r = Symbolic [Term i j r] deriving (Show,Eq,Functor)++-- instance (Ord i, Ord j, Num r) => Num (Symbolic i j r) where+--   (Symbolic ts) + (Symbolic ts') = Symbolic (ts `addTerms` ts')+--   negate = fmap negate+--   (Symbolic ts) * (Symbolic ts') = Symbolic $ multiplyTerms ts ts'+--   fromInteger x = constant (fromInteger x)+--   -- abs x | signum x == -1 = (-1)*x+--   --       | oterwise       = x++--   -- signum = undefined+++++++++++-- -- | Adds two lists of terms+-- addTerms        :: forall i j r. (Ord i, Ord j, Num r)+--                 => [Term i j r] -> [Term i j r] -> [Term i j r]+-- addTerms ts ts' = (\(eps,c) -> Term c eps) <$> Map.toList m+--   where+--     m :: Map.Map (EpsFold i j) r+--     m = Map.fromListWith (+) [ (eps,c) | (Term c eps) <- ts <> ts' ]++-- multiplyTerms        :: forall i j r. (Ord i, Ord j, Num r)+--                      => [Term i j r] -> [Term i j r] -> [Term i j r]+-- multiplyTerms ts ts' = (\(eps,c) -> Term c eps) <$> Map.toList m+--   where+--     m :: Map.Map (EpsFold i j) r+--     m = Map.fromListWith (+) [ (es <> es',c*d) | (Term c es) <- ts, (Term d es') <- ts' ]+++++-- orderedTerms               :: (Ord i, Ord j) => Symbolic i j r -> [Term i j r]+-- orderedTerms (Symbolic ts) = List.sortBy (\(Term _ e1) (Term _ e2) -> cmpSignificance e1 e2) ts++++++++++++++++++  -- zipWith (\j x -> Term x $ singleton (i,j)) [0..] . toList+++++++-- orderTerms               :: (Ord i, Ord j) => Symbolic i j r -> Symbolic i j r+-- orderTerms (Symbolic ts) = Symbolic $ List.sortBy cmpSignificance ts++++-- fromPoint'   :: Foldable f => i -> f r -> Symbolic i Int r+-- fromPoint' i = Symbolic . zipWith (\j x -> Term x [(i,j)]) [0..] . toList++++-- testZ :: Symbolic Int Int Int+-- testZ = (5 + 6) *++++++  --   case sign i of+  --                   (-1) -> Negative $ fromInteger i+  --                   0    -> Zero+  --                   _    -> Positive $ fromInteger i+  -- negate        = \case+  --   Negative c -> Positive c+  --   Positive c -> Negative c+++-- newtype N = N String deriving (Show,Eq)+++-- instance Num N where+--   (N x) + (N y) = N $ x <> "+" <> y+--   (N x) * (N y) = N $ x <> y+--   negate  (N x) = N ("negate(" <> x <> ")")+--   fromInteger = N . show+++-- n       :: (Ord i, Ord j) => String -> i -> j -> Symbolic i j N+-- n x i j = Symbolic [Term (N x) mempty, Term 1 (singleton (i,j))]++++++-- testM3 = det33 $ V3 (fromPoint' [N "px", N "py"] <> 1)+--                     (fromPoint' [N "px", N "py"] <> 1)+--    (fromPoint' [N "px", N "py"] <> 1)+-- -- (V3 (N "qx") (N "qy") 1)+-- --                     (V3 (N "rx") (N "ry") 1)

@@ -0,0 +1,29 @@+module Algorithms.Geometry.SoS.AsPoint where++import           Control.CanAquire+import           Control.Lens+import           Data.Ext+import           Data.Geometry.Point.Class+import           Data.Geometry.Point.Internal+import           Data.Geometry.Properties+import           Data.Geometry.Vector++--------------------------------------------------------------------------------+-- | a P is a 'read only' point in d dimensions+newtype P i d r = P i deriving (Eq, Show)++-- | Indxec type that can disambiguate points+newtype SoSIndex i = SoSIndex i deriving (Show,Eq,Ord)++instance HasIndex (P i d r) i where+  indexOf (P i) = i++instance Int `CanAquire` (Point d r) => (P Int d r) `CanAquire` (Point d r) where+  aquire (P i) = aquire i++type instance NumType   (P i d r) = r+type instance Dimension (P i d r) = d++asPointWithIndex       :: (Arity d, i `CanAquire` Point d r)+                       => P i d r -> Point d r :+ SoSIndex i+asPointWithIndex (P i) = aquire i :+ (SoSIndex i)

@@ -0,0 +1,13 @@+module Algorithms.Geometry.SoS.Determinant where++import           Algorithms.Geometry.SoS.Sign+import           Algorithms.Geometry.SoS.Symbolic+import           Data.Geometry.Matrix+++-- | pre: computes the sign of the determinant+signDet   :: (HasDeterminant d, Ord i, Num r, Ord r) => Matrix d d (Symbolic i r) -> Sign+signDet m = case det m `compare` 0 of+              LT -> Negative+              GT -> Positive+              EQ -> error "signDet: determinant is zero! this should not happen!"

@@ -0,0 +1,77 @@+{-# LANGUAGE TemplateHaskell #-}+module Algorithms.Geometry.SoS.Expr where++import           Control.Lens+import qualified Data.List as List+import           Data.List.NonEmpty (NonEmpty(..),nonEmpty)++--------------------------------------------------------------------------------++data Expr v r = Constant r+              | Negate (Expr v r)+              | Sum  [Expr v r]+              | Prod [Expr v r]+              | Var v+              deriving (Show,Eq)+makePrisms ''Expr+++foldExpr :: (r -> b) -> (b -> b) -> ([b] -> b) -> ([b] -> b) -> (v -> b) -> Expr v r -> b+foldExpr con' neg' sum' prod' var' = go+  where+    go = \case+      Constant c -> con' c+      Negate e   -> neg'  $ go e+      Sum es     -> sum'  $ map go es+      Prod es    -> prod' $ map go es+      Var v      -> var' v++-- | Test if the expression has any variables.+hasVariables :: Expr v r -> Bool+hasVariables = foldExpr (const False)+                        id+                        or+                        or+                        (const True)++instance (Num r) => Num (Expr i r) where+  fromInteger = Constant . fromInteger+  negate      = \case+    Negate e -> e+    e        -> Negate e++  (Sum es) + (Sum es') = Sum $ es <> es'+  (Sum es) + e         = Sum (e:es)+  e        + (Sum es)  = Sum (e:es)+  e        + e'        = Sum [e,e']++  (Prod es) * (Prod es') = Prod $ es <> es'+  (Prod es) * e          = Prod (e:es)+  e         * (Prod es)   = Prod (e:es)+  e         * e'          = Prod [e,e']+++simplify :: (Num r, Eq r) => Expr v r -> Expr v r+simplify = \case+  Prod es  -> case filter (isn't $ _Constant.only 1) es of+                []  -> Constant 1+                es' -> Prod $ map simplify es'+  Sum  es  -> case filter (isn't $ _Constant.only 0) es of+                []  -> Constant 0+                es' -> Sum $ map simplify es'+  Negate e -> Negate $ simplify e+  e        -> e++prettyP :: (Show r, Show v) => Expr v r -> String+prettyP = \case+  Constant c  -> show c+  Negate e    -> "(-1)*(" <> prettyP e <> ")"+  Prod es     -> mconcat [ "("+                            , List.intercalate ")*(" (prettyP <$> es)+                            , ")"+                            ]+  Sum es     -> mconcat [ "("+                        , List.intercalate ") + (" (prettyP <$> es)+                        , ")"+                        ]+  Var v -> show v

@@ -0,0 +1,28 @@+module Algorithms.Geometry.SoS.Internal where++import           Algorithms.Geometry.SoS.AsPoint+import           Algorithms.Geometry.SoS.Orientation+import           Control.CanAquire+import           Data.Geometry.Point.Internal++--------------------------------------------------------------------------------++-- simulateSimplicity :: forall t d r b. (Traversable t, SoSD d)+--                    => (forall p. ( AsPoint p, HasIndex p+--                                  , d ~ Dimension p, r ~ NumType p+--                                  ) => t p -> b)+--                    -> t (Point d r) -> b+-- simulateSimplicity = simulateSimplicity'+++-- | The actual implementation of SoS+simulateSimplicity'     :: forall t d r b. (Traversable t, SoS d)+                        => (forall i. ( CanAquire i (Point d r)+                                      , SoS d+                                      ) => t (P i d r) -> b)+                        -> t (Point d r) -> b+simulateSimplicity' alg = runAcquire alg'+  where+    alg' :: forall i. CanAquire i (Point d r) => t i -> b+    alg' = alg . fmap (P @i @d @r)+      -- ideally the fmap would just be a coerce, but GHC does not want to do that.

@@ -0,0 +1,83 @@+module Algorithms.Geometry.SoS.Orientation( SoS++                                          , sideTest+                                          , sideTest'++                                          , toSymbolic+                                          ) where++import Algorithms.Geometry.SoS.Determinant+import Algorithms.Geometry.SoS.Sign+import Algorithms.Geometry.SoS.Symbolic+import Control.Lens hiding (snoc,cons)+import Data.Ext+import Data.Geometry.Matrix+import Data.Geometry.Point+import Data.Geometry.Vector+import GHC.TypeNats++--------------------------------------------------------------------------------++++-- | A dimension d has support for SoS when we can: compute a+-- dterminant of a d+1 by d+1 dimensional matrix.+type SoS d = (Arity d, HasDeterminant (d+1))++-- | Given a query point q, and a vector of d points defining a+-- hyperplane test if q lies above or below the hyperplane. Each point+-- is assumed to have an unique index of type i that can be used to+-- disambiguate it in case of degeneracies.+--+-- some 1D examples:+--+-- >>> sideTest (Point1 0 :+ 0) (Vector1 $ Point1 2 :+ 1)+-- Negative+-- >>> sideTest (Point1 10 :+ 0) (Vector1 $ Point1 2 :+ 1)+-- Positive+-- >>> sideTest (Point1 2 :+ 0) (Vector1 $ Point1 2 :+ 1)+-- Positive+-- >>> sideTest (Point1 2 :+ 3) (Vector1 $ Point1 2 :+ 1)+-- Negative+--+-- some 2D examples:+--+-- >>> sideTest (Point2 1 2 :+ 0) $ Vector2 (Point2 0 0 :+ 1) (Point2 2 2 :+ 3)+-- Positive+-- >>> sideTest (Point2 1 (-2) :+ 0) $ Vector2 (Point2 0 0 :+ 1) (Point2 2 2 :+ 3)+-- Negative+-- >>> sideTest (Point2 1 1 :+ 0) $ Vector2 (Point2 0 0 :+ 1) (Point2 2 2 :+ 3)+-- Positive+-- >>> sideTest (Point2 1 1 :+ 10) $ Vector2 (Point2 0 0 :+ 1) (Point2 2 2 :+ 3)+-- Negative+-- >>> sideTest (Point2 1 1 :+ 10) $ Vector2 (Point2 0 0 :+ 3) (Point2 2 2 :+ 1)+-- Negative+sideTest      :: (SoS d, Num r, Ord r, Ord i)+              => Point d r :+ i -> Vector d (Point d r :+ i) -> Sign+sideTest q ps = sideTest'' . fmap toSymbolic $ cons q ps++-- | Given an input point, transform its number type to include+-- symbolic $\varepsilon$ expressions so that we can use SoS.+toSymbolic          :: (Ord i, Arity d) => Point d r :+ i -> Point d (Symbolic (i,Int) r)+toSymbolic (p :+ i) = p&vector' %~ imap (\j x -> symbolic x (i,j))++-- | Given a point q and a vector of d points defining a hyperplane,+-- test on which side of the hyperplane q lies.+--+-- TODO: Specify what the sign means+sideTest'      :: (Num r, Ord r, Ord i, HasDeterminant (d+1), Arity d, Arity (d+1))+               => Point d (Symbolic i r) -> Vector d (Point d (Symbolic i r)) -> Sign+sideTest' q ps = sideTest'' $ cons q ps++-- | Given a vector of points, tests if the point encoded in the first+-- row is above/below the hyperplane defined by the remaining points+-- (rows).+sideTest'' :: (Num r, Ord r, Ord i, HasDeterminant (d+1), Arity d, Arity (d+1))+           => Vector (d+1) (Point d (Symbolic i r)) -> Sign+sideTest'' = signDet . Matrix . fmap mkLambdaRow++-- | Given a point produces the vector/row corresponding to this point+-- in a homogeneous matrix represetnation. I.e. we add a 1 as an+-- additonal column at the end.+mkLambdaRow :: (Num r, Arity d, Arity (d+1)) => Point d r -> Vector (d+1) r+mkLambdaRow = flip snoc 1 . view vector'

@@ -0,0 +1,30 @@+module Algorithms.Geometry.SoS.Sign where++import qualified Data.List as List+import           Data.Maybe++--------------------------------------------------------------------------------++-- | The sign of an expression+data Sign = Negative | Positive deriving (Show,Eq,Ord,Enum,Bounded)++flipSign :: Sign -> Sign+flipSign = \case+  Negative -> Positive+  Positive -> Negative++--------------------------------------------------------------------------------++-- | Given the terms, in decreasing order of significance, computes the sign+--+-- i.e. expects a list of terms, we base the sign on the sign of the first non-zero term.+--+-- pre: the list contains at least one such a term.+signFromTerms :: (Num r, Eq r) => [r] -> Sign+signFromTerms = List.head . mapMaybe signum'+  where+    signum' x = case signum x of+                  -1    -> Just Negative+                  0     -> Nothing+                  1     -> Just Positive+                  _     -> error "signum': absurd"

@@ -0,0 +1,352 @@+module Algorithms.Geometry.SoS.Symbolic(+    EpsFold+  , eps, mkEpsFold+  , hasNoPertubation+  , factors+  , suitableBase++  , Term(..), term, constantFactor++  , Symbolic+  , constant, symbolic, perturb++  , toTerms+  , signOf+  ) where++import           Algorithms.Geometry.SoS.Sign (Sign(..))+import           Control.Lens+import           Data.Foldable (toList)+import qualified Data.List as List+import qualified Data.Map as Map+import qualified Data.Map.Merge.Strict as Map+import           Data.Maybe (isNothing)+import           Data.Word+import           Test.QuickCheck (Arbitrary(..), listOf, suchThat)+import           Test.QuickCheck.Instances ()++--------------------------------------------------------------------------------+-- * EpsFolds++{-+Let \(\mathcal{I}\) be a bag with indices, let \(c\) be an upper+bound on the number of times a single item may occur in+\(\mathcal{I}\), and let \(\varepsilon\) be a function mapping indices+to real numbers that satisfies:++1. \(0 < \varepsilon(j) < 1\), for all \(1 \leq j\),+2. \(\prod_{0 \leq i \leq j} \varepsilon(i)^c > \varepsilon(k)\), for all \(1 \leq j < k\)++Note that such a function exists:++\begin{lemma}+  \label{lem:condition_2}+  Let \(\delta \in (0,1)\) and \(d \geq c+1\). The function+  \(\varepsilon(i) = \delta^{d^i}\) satisfies condition 2.+\end{lemma}++\begin{proof}+  By transitivity it suffices to argue this for \(k=j+1\):++  \begin{align*}+           &\qquad \prod_{0 \leq i \leq j} \varepsilon(i)^c > \varepsilon(j+1) \\+    \equiv &\qquad \prod_{0 \leq i \leq j} (\delta^{d^i})^c > \delta^{d^{j+1}}\\+    \equiv &\qquad \prod_{0 \leq i \leq j} \delta^{cd^i}    > \delta^{d^{j+1}} \\+    \equiv &\qquad \delta^{\sum_{0 \leq i \leq j} cd^i} > \delta^{d^{j+1}} &+                                                                    \text{using+                                                                    }+                                                                    \delta \in (0,1)\\+    \equiv &\qquad \sum_{0 \leq i \leq j} cd^i < d^{j+1} \\+    \equiv &\qquad c\sum_{0 \leq i \leq j} d^i < d^{j+1} \\+  \end{align*}++  We prove this by induction.++  For the base case \(j=0\): we have \(0 < 1\), which is trivially true.++  For the step case we have the induction hypothesis+  \(c\sum_{0 \leq i \leq j} d^i < d^{j+1}\), and we have to prove that+  \(c\sum_{0 \leq i \leq j+1} d^i < d^{j+2}\):++  \begin{align*}+    c\sum_{0 \leq i \leq j+1} d^i+    &= cd^{j+1} + c\sum_{0 \leq i \leq j} d^i \\+    &< cd^{j+1} + d^{j+1}   & \text{using IH}  \\+    &= (c+1)d^{j+1}        & \text{using that } c+1 \leq d \\+    &\leq dd^{j+1}  \\+    &=d^{j+2}+  \end{align*}+  This completes the proof.+\end{proof}+++++++An EpsFold now represents the term++\[ \prod_{i \in \mathcal{I}} \varepsilon(i) \]++for some bag \(\mathcal{I}\).+++Let \(\mathcal{J}\) be some sub-bag of \(\mathcal{I}\). Note that+condition 2 implies that:++\(\prod_{i \in \mathcal{J}} \varepsilon(i) > \varepsilon(k)\), for all \(1 \leq j < k\)++This means that when comparing two EpsFolds, say \(e_1\) and \(e_2\),+representing bags \(\mathcal{I}_1\) and \(\mathcal{I}_2\),+respectively. It suffices to compare the largest index+\(j \in \mathcal{I}_1\setminus\mathcal{I}_2\) with the largest index+\(k \in \mathcal{I}_2\setminus\mathcal{I}_1\). We have that++\(e_1 > e_2\) if and only if \(j < k\).+-}+newtype EpsFold i = Pi (Bag i) deriving (Semigroup,Monoid)++-- | Gets the factors+factors         :: EpsFold i -> Bag i+factors (Pi is) = is++-- | Creates the term \(\varepsilon(i)\)+eps :: i -> EpsFold i+eps = Pi . singleton++mkEpsFold :: Ord i => [i] -> EpsFold i+mkEpsFold = Pi . foldMap singleton++++-- | computes a base 'd' that can be used as:+--+-- \( \varepsilon(i) = \varepsilon^{d^i} \)+suitableBase :: EpsFold i -> Int+suitableBase = max 2 . (1+) . maxMultiplicity . factors++instance Show i => Show (EpsFold i) where+  showsPrec d (Pi b) = showParen (d > app_prec) $+                         showString "Pi " . showsPrec d (toList b)+    where+      app_prec = 10+++instance Ord i => Eq (EpsFold i) where+  e1 == e2 = (e1 `compare` e2) == EQ++instance Ord i => Ord (EpsFold i) where+  (Pi e1) `compare` (Pi e2) = k `compare` j -- note that k and j are flipped here+    where+      j = maximum' $ e1 `difference` e2+      k = maximum' $ e2 `difference` e1+    -- note: If the terms are all the same, the difference of the bags is empty+    -- and thus both e1e2 and e2e1 are Nothing and thus equal.++    -- otherwise, let j be the largest term that is in e1 but not in e2.+    -- If e2 does not have any terms at all (Nothing) it will be bigger than e1+    --+    -- if e2 does have a term, let k be the largest one, then the+    -- biggest of those terms is the pair whose indices comes first.++instance (Arbitrary i, Ord i) => Arbitrary (EpsFold i) where+  arbitrary = (mkEpsFold . take 4) <$> listOf arbitrary+++-- | Test if the epsfold has no pertubation at all (i.e. if it is \(\Pi_{\emptyset}\)+hasNoPertubation        :: EpsFold i -> Bool+hasNoPertubation (Pi b) = null b+++--------------------------------------------------------------------------------+-- * Terms++-- | A term 'Term c es' represents a term:+--+-- \[ c \Pi_{i \in es} \varepsilon(i)+-- \]+--+-- for a constant c and an arbitrarily small value \(\varepsilon\),+-- parameterized by i.+data Term i r = Term r (EpsFold i) deriving (Eq,Functor)++-- | Lens to access the constant 'c' in the term.+constantFactor :: Lens' (Term i r) r+constantFactor = lens (\(Term c _) -> c) (\(Term _ es) c -> Term c es)+++instance (Show i, Show r) => Show (Term i r) where+  showsPrec d (Term c es) = showParen (d > up_prec) $+                               showsPrec (up_prec + 1) c+                             . showString " * "+                             . showsPrec (up_prec + 1) es+    where+      up_prec = 5+++-- | Creates a singleton term+term     :: r -> i -> Term i r+term r i = Term r $ eps i++instance (Ord i, Ord r, Num r) => Ord (Term i r) where+  (Term c e1) `compare` (Term d e2) = case (hasNoPertubation e1, hasNoPertubation e2) of+                                        (True,True) -> c    `compare` d+                                        _           -> case (signum c, signum d) of+                                                         (-1,-1) -> e2 `compare` e1+                                                         (0,0)   -> e1 `compare` e2+                                                         (1,1)   -> e1 `compare` e2+                                                         (-1,_)  -> LT+                                                         (_,-1)  -> GT+                                                         _       -> error "SoS: Term.ord absurd"+  -- If both the eps folds are zero, and thus we just have constants+  -- then we should compare the individual terms.++  -- if *one* of the two has an eps term, then we can choose eps to be+  -- arbitrarily small, i.e. small enough so that that terms is+  -- actually smaller than the other term.  this is reflected since+  -- findMax will then return a Noting, which is smaller than anything+  -- else++  -- if both terms have epsilon terms, we first look at the sign. If+  -- they have non-negative signs we compare the eps-folds as in the+  -- paper. (Lemma 3.3). If both are negative, that reverses the+  -- ordering. If the signs are different then we can base the+  -- ordering on that.++instance (Arbitrary r, Arbitrary (EpsFold i), Ord i) => Arbitrary (Term i r) where+  arbitrary = Term <$> arbitrary <*> arbitrary++--------------------------------------------------------------------------------+-- * Symbolic++-- | Represents a Sum of terms, i.e. a value that has the form:+--+-- \[+--   \sum c \Pi_i \varepsilon(i)+-- \]+--+-- The terms are represented in order of decreasing significance.+--+-- The main idea in this type is that, if symbolic values contains+-- \(\varepsilon(i)\) terms we can always order them. That is, two+-- Symbolic terms will be equal only if:+--+-- - they contain *only* a constant term (that is equal)+-- - they contain the exact same \(\varepsilon\)-fold.+--+newtype Symbolic i r = Sum (Map.Map (EpsFold i) r) deriving (Functor)++-- | Produces a list of terms, in decreasing order of significance+toTerms         :: Symbolic i r -> [Term i r]+toTerms (Sum m) = map (\(i,c) -> Term c i) . Map.toDescList $ m++-- | Computing the Sign of an expression. (Nothing represents zero)+signOf   :: (Num r, Eq r) => Symbolic i r -> Maybe Sign+signOf e = case List.dropWhile (== 0) . map (\(Term c _) -> signum c) $ toTerms e of+             []     -> Nothing+             (-1:_) -> Just Negative+             _      -> Just Positive++instance (Ord i, Eq r, Num r) => Eq (Symbolic i r) where+  e1 == e2 = isNothing $ signOf (e1 - e2)++instance (Ord i, Ord r, Num r) => Ord (Symbolic i r) where+  e1 `compare` e2 = case signOf (e1 - e2) of+                      Nothing       -> EQ+                      Just Negative -> LT+                      Just Positive -> GT++instance (Ord i, Num r, Eq r) => Num (Symbolic i r) where+  (Sum e1) + (Sum e2) = Sum $ Map.merge Map.preserveMissing -- insert things only in e1+                                        Map.preserveMissing -- insert things only in e2+                                        combine+                                        e1 e2+    where+      -- if things are in both e1 and e2, we add the constant terms. If they are non-zero+      -- we use this value in the map. Otherwise we drop it.+      combine = Map.zipWithMaybeMatched+                (\_ c d -> let x = c + d in if x /= 0 then Just x else Nothing)+    -- Symbolic $ Map.unionWith (+) ts ts'++  negate = fmap negate++  (Sum ts) * (Sum ts') = Sum $ Map.fromListWith (+) [ (es <> es',c*d)+                                                    | (es, c) <- Map.toList ts+                                                    , (es',d) <- Map.toList ts'+                                                    , c*d /= 0+                                                    ]++  fromInteger x = constant (fromInteger x)++  signum s = case signOf s of+               Nothing       -> 0+               Just Negative -> (-1)+               Just Positive -> 1++  abs x | signum x == -1 = (-1)*x+        | otherwise      = x+++instance (Show i, Show r) => Show (Symbolic i r) where+  showsPrec d s = showParen (d > app_prec) $+                    showString "Sum " . showsPrec d (toTerms s)+    where+      app_prec = 10++instance (Arbitrary r, Ord i, Arbitrary (EpsFold i)) => Arbitrary (Symbolic i r) where+  arbitrary = Sum <$> arbitrary++----------------------------------------++-- | Creates a constant symbolic value+constant   :: Ord i => r -> Symbolic i r+constant c = Sum $ Map.singleton mempty c++-- | Creates a symbolic vlaue with a single indexed term. If you just need a constant (i.e. non-indexed), use 'constant'+symbolic     :: Ord i => r -> i -> Symbolic i r+symbolic r i = Sum $ Map.singleton (eps i) r++-- | given the value c and the index i, creates the perturbed value+-- \(c + \varepsilon(i)\)+perturb      :: (Num r, Ord i) => r -> i -> Symbolic i r+perturb c i = Sum $ Map.fromAscList [ (mempty,c) , (eps i,1) ]+++--------------------------------------------------------------------------------++-- | The word specifiies how many *duplicates* there are. I.e. If the+-- Bag maps k to i, then k has multiplicity i+1.+newtype Bag a = Bag (Map.Map a Int) deriving (Show,Eq,Ord,Arbitrary)++singleton   :: k -> Bag k+singleton x = Bag $ Map.singleton x 0+++instance Foldable Bag where+  -- ^ Takes multiplicity into account.+  foldMap f (Bag m) =+    Map.foldMapWithKey (\k d -> foldMap f (List.replicate (fromIntegral d+1) k)) m+  null (Bag m) = Map.null m++instance Ord k => Semigroup (Bag k) where+  (Bag m) <> (Bag m') = Bag $ Map.unionWith (\d d' -> d + d' + 1) m m'++instance Ord k => Monoid (Bag k) where+  mempty = Bag $ Map.empty++-- | Computes the difference of the two maps+difference                   :: Ord a => Bag a -> Bag a -> Bag a+difference (Bag m1) (Bag m2) = Bag $ Map.differenceWith updateCount m1 m2+  where+    updateCount i j = let d = i - j -- note that we should actually compare (i+1) and (j+1)+                      in if d <= 0 then Nothing -- we have no copies left+                                   else Just $ d - 1+++maximum'         :: Bag b -> Maybe b+maximum' (Bag m) = fmap fst . Map.lookupMax $ m+++-- | maximum multiplicity of an element in the bag+maxMultiplicity         :: Bag a -> Int+maxMultiplicity (Bag m) = maximum . (0:) . map (1+) . Map.elems $ m

@@ -19,6 +19,7 @@ import           Data.Geometry.Point import           Data.Geometry.Vector import qualified Data.LSeq as LSeq+import           Data.Measured.Class import qualified Data.Sequence as S import qualified Data.Traversable as Tr 

@@ -50,7 +50,7 @@                   => NonEmpty.NonEmpty (Point d r :+ p) -> SplitTree d p r () fairSplitTree pts = foldUp node' Leaf $ fairSplitTree' n pts'   where-    pts' = GV.imap sortOn . pure . g $ pts+    pts' = imap sortOn . pure . g $ pts     n    = length $ pts'^.GV.element (C :: C 0)      sortOn' i = NonEmpty.sortWith (^.core.unsafeCoord i)@@ -354,7 +354,7 @@ -- -- 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 ->+extends = imap (\i pts ->                      ClosedRange ((LSeq.head pts)^.core.unsafeCoord (i + 1))                                  ((LSeq.last pts)^.core.unsafeCoord (i + 1))) 

@@ -11,6 +11,7 @@ -------------------------------------------------------------------------------- module Data.Geometry.Arrangement.Internal where +import           Algorithms.BinarySearch import           Control.Lens import qualified Data.CircularSeq as CSeq import           Data.Ext@@ -25,7 +26,6 @@ import qualified Data.List as List import           Data.Maybe import           Data.Ord (Down(..))-import           Data.Sequence.Util import qualified Data.Vector as V import           Data.Vinyl.CoRec @@ -192,13 +192,10 @@ 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+    Sides t r b l       = fmap sideIntersections'      $ sides   rect+    Corners tl tr br bl = fmap ((,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)@@ -270,8 +267,8 @@     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+    Sides t r b l = sides'' $ arr^.boundedArea+    sides''       = fmap (\(ClosedLineSegment a c) -> LineSegment (Closed a) (Open c)) . sides      findSide q = fmap fst . List.find (onSegment q . snd) $ zip [1..] [t,r,b,l] 

@@ -123,21 +123,27 @@ pattern Sphere c r = Boundary (Ball c r) {-# COMPLETE Sphere #-} --+-- |+_BallSphere :: Iso (Disk p r) (Disk p s) (Circle p r) (Circle p s)+_BallSphere = _Boundary  -------------------------------------------------------------------------------- -- * Disks and Circles, aka 2-dimensional Balls and Spheres  type Disk p r = Ball 2 p r +-- | Given the center and the squared radius, constructs a disk pattern Disk     :: Point 2 r :+ p -> r -> Disk p r pattern Disk c r = Ball c r {-# COMPLETE Disk #-} - type Circle p r = Sphere 2 p r +-- | Iso for converting between Disks and Circles, i.e. forgetting the boundary+_DiskCircle  :: Iso (Disk p r) (Disk p s) (Circle p r) (Circle p s)+_DiskCircle = _BallSphere++-- | Given the center and the squared radius, constructs a circle pattern Circle     :: Point 2 r :+ p ->  r -> Circle p r pattern Circle c r = Sphere c r {-# COMPLETE Circle #-}

@@ -0,0 +1,169 @@+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE TemplateHaskell #-}+module Data.Geometry.BezierSpline(+    BezierSpline (BezierSpline)+  , controlPoints+  , fromPointSeq++  , evaluate+  , split+  , subBezier+  , tangent+  , approximate+  , parameterOf+  , snap++  , pattern Bezier2, pattern Bezier3+  ) where++import           Control.Lens hiding (Empty)+import qualified Data.Foldable as F+import           Data.Geometry.Point+import           Data.Geometry.Properties+import           Data.Geometry.Transformation+import           Data.Geometry.Vector+import           Data.LSeq (LSeq)+import qualified Data.LSeq as LSeq+import           Data.Sequence (Seq(..))+import qualified Data.Sequence as Seq+import           Data.Traversable (fmapDefault,foldMapDefault)+import           GHC.TypeNats+import qualified Test.QuickCheck as QC++--------------------------------------------------------------------------------++-- | Datatype representing a Bezier curve of degree \(n\) in \(d\)-dimensional space.+newtype BezierSpline n d r = BezierSpline { _controlPoints :: LSeq (1+n) (Point d r) }+makeLenses ''BezierSpline++-- | Quadratic Bezier Spline+pattern Bezier2      :: Point d r -> Point d r -> Point d r -> BezierSpline 2 d r+pattern Bezier2 p q r <- ((F.toList . LSeq.take 3 . _controlPoints) -> [p,q,r])+  where+    Bezier2 p q r = fromPointSeq . Seq.fromList $ [p,q,r]+{-# COMPLETE Bezier2 #-}++-- | Cubic Bezier Spline+pattern Bezier3         :: Point d r -> Point d r -> Point d r -> Point d r -> BezierSpline 3 d r+pattern Bezier3 p q r s <- ((F.toList . LSeq.take 4 . _controlPoints) -> [p,q,r,s])+  where+    Bezier3 p q r s = fromPointSeq . Seq.fromList $ [p,q,r,s]+{-# COMPLETE Bezier3 #-}++deriving instance (Arity d, Eq r) => Eq (BezierSpline n d r)++type instance Dimension (BezierSpline n d r) = d+type instance NumType   (BezierSpline n d r) = r+++instance (Arity n, Arity d, QC.Arbitrary r) => QC.Arbitrary (BezierSpline n d r) where+  arbitrary = fromPointSeq . Seq.fromList <$> QC.vector (fromIntegral . (1+) . natVal $ C @n)++-- | Constructs the Bezier Spline from a given sequence of points.+fromPointSeq :: Seq (Point d r) -> BezierSpline n d r+fromPointSeq = BezierSpline . LSeq.promise . LSeq.fromSeq+++instance (Arity d, Show r) => Show (BezierSpline n d r) where+  show (BezierSpline ps) =+    mconcat [ "BezierSpline", show $ length ps - 1, " ", show (F.toList ps) ]++instance Arity d => Functor (BezierSpline n d) where+  fmap = fmapDefault++instance Arity d => Foldable (BezierSpline n d) where+  foldMap = foldMapDefault++instance Arity d => Traversable (BezierSpline n d) where+  traverse f (BezierSpline ps) = BezierSpline <$> traverse (traverse f) ps++instance (Fractional r, Arity d, Arity (d + 1), Arity n)+          => IsTransformable (BezierSpline n d r) where+  transformBy = transformPointFunctor++instance PointFunctor (BezierSpline n d) where+  pmap f = over controlPoints (fmap f)++-- | Evaluate a BezierSpline curve at time t in [0, 1]+--+-- pre: \(t \in [0,1]\)+evaluate    :: (Arity d, Ord r, Num r) => BezierSpline n d r -> r -> Point d r+evaluate b t = evaluate' (b^.controlPoints.to LSeq.toSeq)+  where+    evaluate' = \case+      (p :<| Empty)  -> p+      pts@(_ :<| tl) -> let (ini :|> _) = pts in evaluate' $ Seq.zipWith blend ini tl+      _              -> error "evaluate: absurd"++    blend p q = p .+^ t *^ (q .-. p)+++tangent   :: (Arity d, Num r, 1 <= n) => BezierSpline n d r -> Vector d r+tangent b = b^?!controlPoints.ix 1  .-. b^?!controlPoints.ix 0++-- | Restrict a Bezier curve to th,e piece between parameters t < u in [0, 1].+subBezier     :: (KnownNat n, Arity d, Ord r, Num r)+              => r -> r -> BezierSpline n d r -> BezierSpline n d r+subBezier t u = fst . split u . snd . split t++-- | Split a Bezier curve at time t in [0, 1] into two pieces.+split :: forall n d r. (KnownNat n, Arity d, Ord r, Num r)+      => r -> BezierSpline n d r -> (BezierSpline n d r, BezierSpline n d r)+split t b | t < 0 || t > 1 = error "Split parameter out of bounds."+          | otherwise      = let n  = fromIntegral $ natVal (C @n)+                                 ps = collect t $ b^.controlPoints+                             in ( fromPointSeq . Seq.take (n + 1) $ ps+                                , fromPointSeq . Seq.drop n       $ ps+                                )++collect   :: (Arity d, Ord r, Num r) => r -> LSeq n (Point d r) -> Seq (Point d r)+collect t = go . LSeq.toSeq+  where+    go = \case+      ps@(_ :<| Empty) -> ps+      ps@(p :<| tl)    -> let (ini :|> q) = ps in (p :<| go (Seq.zipWith blend ini tl)) :|> q+      _                -> error "collect: absurd"++    blend p q = p .+^ t *^ (q .-. p)++-- {-++-- -- | Merge to Bezier pieces. Assumes they can be merged into a single piece of the same degree+-- --   (as would e.g. be the case for the result of a 'split' operation).+-- --   Does not test whether this is the case!+-- merge :: (Arity d, Ord r, Num r) => (Bezier d r, Bezier d r) -> Bezier d r++-- -}++-- | Approximate Bezier curve by Polyline with given resolution.+approximate :: forall n d r. (KnownNat n, Arity d, Ord r, Fractional r)+            => r -> BezierSpline n d r -> [Point d r]+approximate eps b+    | squaredEuclideanDist p q < eps^2 = [p,q]+    | otherwise                        = let (b1, b2) = split 0.5 b+                                         in approximate eps b1 ++ tail (approximate eps b2)+  where+    p = b^.controlPoints.to LSeq.head+    q = b^.controlPoints.to LSeq.last++-- | Given a point on (or close to) a Bezier curve, return the corresponding parameter value.+--   (For points far away from the curve, the function will return the parameter value of+--   an approximate locally closest point to the input point.)+parameterOf      :: (Arity d, Ord r, Fractional r) => BezierSpline n d r -> Point d r -> r+parameterOf b p = binarySearch (qdA p . evaluate b) treshold (1 - treshold)+  where treshold = 0.0001++binarySearch                                    :: (Ord r, Fractional r) => (r -> r) -> r -> r -> r+binarySearch f l r | abs (f l - f r) < treshold = m+                   | derivative f m  > 0        = binarySearch f l m+                   | otherwise                  = binarySearch f m r+  where m = (l + r) / 2+        treshold = 0.0001++derivative     :: Fractional r => (r -> r) -> r -> r+derivative f x = (f (x + delta) - f x) / delta+  where delta = 0.00001++-- | Snap a point close to a Bezier curve to the curve.+snap   :: (Arity d, Ord r, Fractional r) => BezierSpline n d r -> Point d r -> Point d r+snap b = evaluate b . parameterOf b

@@ -1,7 +1,8 @@ module Data.Geometry.Boundary where -import           Data.Geometry.Properties-import           Data.Geometry.Transformation+import Control.Lens (iso,Iso)+import Data.Geometry.Properties+import Data.Geometry.Transformation  -------------------------------------------------------------------------------- @@ -12,6 +13,10 @@  type instance NumType (Boundary g)   = NumType g type instance Dimension (Boundary g) = Dimension g++-- | Iso for converting between things with a boundary and without its boundary+_Boundary :: Iso g h (Boundary g) (Boundary h)+_Boundary = iso Boundary (\(Boundary b) -> b)   -- | Result of a query that asks if something is Inside a g, *on* the boundary

@@ -1,5 +1,3 @@-{-# LANGUAGE TemplateHaskell  #-}-{-# LANGUAGE ScopedTypeVariables  #-} {-# LANGUAGE UndecidableInstances  #-} {-# LANGUAGE DeriveAnyClass  #-} {-# OPTIONS_GHC -fno-warn-orphans #-}@@ -14,56 +12,16 @@ -- -------------------------------------------------------------------------------- module Data.Geometry.Box( module Data.Geometry.Box.Internal-                        , topSide, leftSide, bottomSide, rightSide-                        , sides, sides'+                        , module Data.Geometry.Box.Corners+                        , module Data.Geometry.Box.Sides                         ) where  import Control.DeepSeq+import Data.Geometry.Box.Corners import Data.Geometry.Box.Internal-import Data.Geometry.LineSegment+import Data.Geometry.Box.Sides import Data.Geometry.Vector  --------------------------------------------------------------------------------  deriving instance (NFData p, NFData r, Arity d) => NFData (Box d p r)---topSide :: Num r => Rectangle p r -> LineSegment 2 p r-topSide = (\(l,r,_,_) -> ClosedLineSegment l r) . corners---- | Oriented from *left to right*-bottomSide :: Num r => Rectangle p r -> LineSegment 2 p r-bottomSide = (\(_,_,r,l) -> ClosedLineSegment l r) . corners-----leftSide  :: Num r => Rectangle p r -> LineSegment 2 p r-leftSide = (\(t,_,_,b) -> ClosedLineSegment b t) . corners---- | The right side, oriented from *bottom* to top-rightSide :: Num r => Rectangle p r -> LineSegment 2 p r-rightSide = (\(_,t,b,_) -> ClosedLineSegment b t) . corners----- | The sides of the rectangle, in order (Top, Right, Bottom, Left). The sides--- themselves are also oriented in clockwise order. If, you want them in the--- same order as the functions `topSide`, `bottomSide`, `leftSide`, and--- `rightSide`, use `sides'` instead.-sides :: Num r => Rectangle p r -> ( LineSegment 2 p r-                                   , LineSegment 2 p r-                                   , LineSegment 2 p r-                                   , LineSegment 2 p r-                                   )-sides = (\(t,r,b,l) -> (t,flipSegment r,flipSegment b,l)) . sides'----- | The sides of the rectangle. The order of the segments is (Top, Right,--- Bottom, Left).  Note that the segments themselves, are oriented as described--- by the functions topSide, bottomSide, leftSide, rightSide (basically: from--- left to right, and from bottom to top). If you want the segments oriented--- along the boundary of the rectangle, use the `sides` function instead.-sides'   :: Num r => Rectangle p r -> ( LineSegment 2 p r-                                      , LineSegment 2 p r-                                      , LineSegment 2 p r-                                      , LineSegment 2 p r-                                      )-sides' r = (topSide r, rightSide r, bottomSide r, leftSide r)

@@ -0,0 +1,70 @@+{-# LANGUAGE TemplateHaskell  #-}+module Data.Geometry.Box.Corners( Corners(Corners), northWest, northEast, southEast, southWest+                                , corners, cornersInDirection+                                ) where++import Control.Lens (makeLenses,Ixed(..),Index, IxValue,(%~),(&),(^?!))+import Data.Ext+import Data.Functor.Apply+import Data.Geometry.Box.Internal+import Data.Geometry.Directions+import Data.Geometry.Point+import Data.Semigroup.Foldable.Class+import Data.Semigroup.Traversable.Class+import Data.Util+import GHC.Generics (Generic)++--------------------------------------------------------------------------------++-- | A Quadrant data type+data Corners a = Corners { _northWest  :: !a+                         , _northEast  :: !a+                         , _southEast  :: !a+                         , _southWest  :: !a+                         } deriving (Show,Eq,Ord,Generic,Functor,Foldable,Traversable)+makeLenses ''Corners+++type instance Index   (Corners a) = InterCardinalDirection+type instance IxValue (Corners a) = a++instance Ixed (Corners a) where+  ix = \case+    NorthWest -> northWest+    NorthEast -> northEast+    SouthEast -> southEast+    SouthWest -> southWest++instance Foldable1 Corners+instance Traversable1 Corners where+  traverse1 f (Corners a b c d) = Corners <$> f a <.> f b <.> f c <.> f d++instance Applicative Corners where+  pure x = Corners x x x x+  (Corners f g h i) <*> (Corners a b c d) = Corners (f a) (g b) (h c) (i d)++instance Semigroup a => Semigroup (Corners a) where+  s <> s' = (<>) <$> s <*> s'+instance Monoid a => Monoid (Corners a) where+  mempty = pure mempty+++--------------------------------------------------------------------------------++-- | Get the corners of a rectangle, the order is:+-- (TopLeft, TopRight, BottomRight, BottomLeft).+-- The extra values in the Top points are taken from the Top point,+-- the extra values in the Bottom points are taken from the Bottom point+corners :: Num r => Rectangle p r -> Corners (Point 2 r :+ p)+corners r     = let w = width r+                    p = (_maxP r)&core %~ _cwMax+                    q = (_minP r)&core %~ _cwMin+                in Corners (p&core.xCoord %~ (subtract w)) p+                           (q&core.xCoord %~ (+ w))        q+++--------------------------------------------------------------------------------++-- | Gets the corners in a particular direction+cornersInDirection     :: CardinalDirection -> Corners p -> Two p+cornersInDirection d c = (\icd -> c^?!ix icd) <$> interCardinalsOf d

@@ -15,8 +15,11 @@  import           Control.DeepSeq import           Control.Lens+import           Data.Bifoldable import           Data.Bifunctor+import           Data.Bitraversable import           Data.Ext+import qualified Data.Foldable as F import           Data.Geometry.Point import           Data.Geometry.Properties import           Data.Geometry.Transformation@@ -25,12 +28,11 @@ import qualified Data.List.NonEmpty as NE import qualified Data.Range as R import qualified Data.Semigroup.Foldable as F-import qualified Data.Foldable as F import qualified Data.Vector.Fixed as FV import           Data.Vinyl.CoRec (asA) import           GHC.Generics (Generic) import           GHC.TypeLits-import           Test.QuickCheck(Arbitrary(..))+import           Test.QuickCheck (Arbitrary(..))  -------------------------------------------------------------------------------- @@ -60,6 +62,10 @@                      } deriving Generic makeLenses ''Box ++++ -- | Given the point with the lowest coordinates and the point with highest -- coordinates, create a box. box          :: Point d r :+ p -> Point d r :+ p -> Box d p r@@ -108,12 +114,14 @@       r `intersect'` s = asA @(R.Range r) $ r `intersect` s  instance Arity d => Bifunctor (Box d) where-  bimap :: forall p q r s. (p -> q) -> (r -> s) -> Box d p r -> Box d q s-  bimap f g (Box mi ma) = Box (bimap g' f mi) (bimap g' f ma)+  bimap = bimapDefault+instance Arity d => Bifoldable (Box d) where+  bifoldMap = bifoldMapDefault+instance Arity d => Bitraversable (Box d) where+  bitraverse f g (Box mi ma) = Box <$> bitraverse (tr g) f mi <*> bitraverse (tr g) f ma     where-      g' :: Functor g => g (Point d r) -> g (Point d s)-      g' = fmap (fmap g)-+      tr    :: (Traversable t, Applicative f) => (r -> f s) -> t (Point d r) -> f (t (Point d s))+      tr g' = traverse $ traverse g'  -- -- In principle this should also just work for Boxes in higher dimensions. It is just -- -- that we need a better way to compute their corners@@ -240,23 +248,7 @@ height = widthIn (C :: C 2)  --- | Get the corners of a rectangle, the order is:--- (TopLeft, TopRight, BottomRight, BottomLeft).--- The extra values in the Top points are taken from the Top point,--- the extra values in the Bottom points are taken from the Bottom point-corners :: Num r => Rectangle p r -> ( Point 2 r :+ p-                                     , Point 2 r :+ p-                                     , Point 2 r :+ p-                                     , Point 2 r :+ p-                                     )-corners r     = let w = width r-                    p = (_maxP r)&core %~ _cwMax-                    q = (_minP r)&core %~ _cwMin-                in ( p&core.xCoord %~ (subtract w)-                   , p-                   , q&core.xCoord %~ (+ w)-                   , q-                   )+--------------------------------------------------------------------------------  -------------------------------------------------------------------------------- -- * Constructing bounding boxes@@ -282,3 +274,6 @@  instance IsBoxable (Box d p r) where   boundingBox (Box m m') = Box (m&extra .~ ()) (m'&extra .~ ())++instance IsBoxable c => IsBoxable (c :+ e) where+  boundingBox = boundingBox . view core

@@ -0,0 +1,91 @@+{-# LANGUAGE TemplateHaskell  #-}+module Data.Geometry.Box.Sides( Sides(Sides), north, east, south, west+                              , topSide, bottomSide, leftSide, rightSide+                              , sides, sides'++                              , sideDirections+                              ) where++import Data.Geometry.Directions+import Data.Geometry.Box.Internal+import Data.Geometry.Box.Corners+import Data.Geometry.LineSegment+import Data.Functor.Apply+import Data.Semigroup.Foldable.Class+import Data.Semigroup.Traversable.Class+import GHC.Generics (Generic)+import Control.Lens(makeLenses, Ixed(..), Index, IxValue)++--------------------------------------------------------------------------------++-- | The four sides of a rectangle+data Sides a = Sides { _north :: !a+                     , _east  :: !a+                     , _south :: !a+                     , _west  :: !a+                     } deriving (Show,Read,Eq,Generic,Ord,Foldable,Functor,Traversable)+makeLenses ''Sides++instance Applicative Sides where+  pure x = Sides x x x x+  (Sides f g h i) <*> (Sides a b c d) = Sides (f a) (g b) (h c) (i d)++instance Foldable1 Sides+instance Traversable1 Sides where+  traverse1 f (Sides a b c d) = Sides <$> f a <.> f b <.> f c <.> f d++instance Semigroup a => Semigroup (Sides a) where+  s <> s' = (<>) <$> s <*> s'+instance Monoid a => Monoid (Sides a) where+  mempty = pure mempty+++type instance Index   (Sides a) = CardinalDirection+type instance IxValue (Sides a) = a++instance Ixed (Sides a) where+  ix = \case+    North -> north+    East  -> east+    South -> south+    West  -> west++-- | Constructs a Sides value that indicates the appropriate+-- direction.+sideDirections :: Sides CardinalDirection+sideDirections = Sides North East South West++--------------------------------------------------------------------------------++topSide :: Num r => Rectangle p r -> LineSegment 2 p r+topSide = (\(Corners l r _ _) -> ClosedLineSegment l r) . corners++-- | Oriented from *left to right*+bottomSide :: Num r => Rectangle p r -> LineSegment 2 p r+bottomSide = (\(Corners _ _ r l) -> ClosedLineSegment l r) . corners++--+leftSide  :: Num r => Rectangle p r -> LineSegment 2 p r+leftSide = (\(Corners t _ _ b) -> ClosedLineSegment b t) . corners++-- | The right side, oriented from *bottom* to top+rightSide :: Num r => Rectangle p r -> LineSegment 2 p r+rightSide = (\(Corners _ t b _) -> ClosedLineSegment b t) . corners+++-- | The sides of the rectangle, in order (Top, Right, Bottom, Left). The sides+-- themselves are also oriented in clockwise order. If, you want them in the+-- same order as the functions `topSide`, `bottomSide`, `leftSide`, and+-- `rightSide`, use `sides'` instead.+sides   :: Num r => Rectangle p r -> Sides (LineSegment 2 p r)+sides r = let Corners nw ne se sw = corners r+          in Sides (ClosedLineSegment nw ne) (ClosedLineSegment ne se)+                   (ClosedLineSegment se sw) (ClosedLineSegment sw nw)++-- | The sides of the rectangle. The order of the segments is (Top, Right,+-- Bottom, Left).  Note that the segments themselves, are oriented as described+-- by the functions topSide, bottomSide, leftSide, rightSide (basically: from+-- left to right, and from bottom to top). If you want the segments oriented+-- along the boundary of the rectangle, use the `sides` function instead.+sides'   :: Num r => Rectangle p r -> Sides (LineSegment 2 p r)+sides' r = Sides (topSide r) (rightSide r) (bottomSide r) (leftSide r)

@@ -0,0 +1,49 @@+{-# LANGUAGE TemplateHaskell  #-}+module Data.Geometry.Directions( CardinalDirection(..)+                               , _North, _East, _South, _West+                               , oppositeDirection++                                , InterCardinalDirection(..)+                                , _NorthWest, _NorthEast, _SouthEast, _SouthWest++                                , interCardinalsOf+                                ) where++import Control.Lens (makePrisms)+import Data.Util+import GHC.Generics (Generic)++--------------------------------------------------------------------------------++data CardinalDirection = North | East | South | West deriving (Show,Read,Eq,Ord,Enum,Bounded)+makePrisms ''CardinalDirection++--------------------------------------------------------------------------------+-- * Functions on Cardinal Directions++-- | Computes the direction opposite to the given one.+oppositeDirection :: CardinalDirection -> CardinalDirection+oppositeDirection = \case+  North -> South+  East  -> West+  South -> North+  West  -> East++--------------------------------------------------------------------------------++-- | Intercardinal directions+data InterCardinalDirection = NorthWest | NorthEast | SouthEast | SouthWest+  deriving (Show,Read,Eq,Ord,Enum,Generic)+makePrisms ''InterCardinalDirection++--------------------------------------------------------------------------------+-- * Functions on InterCardinal Directions++-- | Get the two intercardinal directions, in increasing order,+-- corresponding to the cardinal direction.+interCardinalsOf :: CardinalDirection -> Two InterCardinalDirection+interCardinalsOf = \case+  North -> Two NorthWest NorthEast+  East  -> Two NorthEast SouthEast+  South -> Two SouthEast SouthWest+  West  -> Two SouthWest NorthWest

@@ -0,0 +1,56 @@+{-# LANGUAGE TemplateHaskell #-}+module Data.Geometry.Ellipse(+    Ellipse(Ellipse)+  , affineTransformation+  , ellipseMatrix+  , unitEllipse+  , circleToEllipse, ellipseToCircle, _EllipseCircle+  ) where++import Control.Lens+import Data.Ext+import Data.Geometry.Ball+import Data.Geometry.Matrix+import Data.Geometry.Transformation+import Data.Geometry.Point+import Data.Geometry.Properties+import Data.Geometry.Vector++--------------------------------------------------------------------------------++-- | A typre representing planar ellipses+newtype Ellipse r = Ellipse { _affineTransformation :: Transformation 2 r }+                   deriving (Show,Eq,Functor,Foldable,Traversable)+makeLenses ''Ellipse++type instance Dimension (Ellipse r) = 2+type instance NumType   (Ellipse r) = r++instance Num r => IsTransformable (Ellipse r) where+  transformBy t (Ellipse t') = Ellipse $ t |.| t'+++ellipseMatrix :: Iso (Ellipse r) (Ellipse s) (Matrix 3 3 r) (Matrix 3 3 s)+ellipseMatrix = affineTransformation.transformationMatrix++-- | Ellipse representing the unit circle+unitEllipse :: Num r => Ellipse r+unitEllipse = Ellipse $ Transformation identityMatrix++--------------------------------------------------------------------------------+-- | Converting between ellipses and circles++_EllipseCircle :: (Floating r, Eq r) => Prism' (Ellipse r) (Circle () r)+_EllipseCircle = prism' circleToEllipse ellipseToCircle++ellipseToCircle   :: (Num r, Eq r) => Ellipse r -> Maybe (Circle () r)+ellipseToCircle e = case e^.ellipseMatrix of+      Matrix (Vector3 (Vector3 sx 0 x)+                      (Vector3 0 sy y)+                      (Vector3 0 0  1)+             )+           | sx == sy -> Just $ Circle (ext $ Point2 x y) (sx*sx)+      _               -> Nothing++circleToEllipse                            :: Floating r => Circle p r -> Ellipse r+circleToEllipse (Circle (Point v :+ _) rr) = Ellipse $ translation v |.| uniformScaling (sqrt rr)

@@ -31,7 +31,6 @@ makeLenses ''HalfLine  deriving instance (Show r, Arity d)   => Show    (HalfLine d r)-deriving instance (Eq r, Arity d)     => Eq      (HalfLine d r) deriving instance (NFData r, Arity d) => NFData  (HalfLine d r)  deriving instance Arity d           => Functor (HalfLine d)@@ -40,6 +39,11 @@  type instance Dimension (HalfLine d r) = d type instance NumType   (HalfLine d r) = r+++instance (Eq r, Fractional r, Arity d) => Eq (HalfLine d r) where+  (HalfLine p u) == (HalfLine q v) = let lam = scalarMultiple u v+                                     in p == q && (signum <$> lam) == Just 1  instance HasStart (HalfLine d r) where   type StartCore  (HalfLine d r) = Point d r

@@ -38,13 +38,14 @@  -------------------------------------------------------------------------------- --- | A Halfspace in \(d\) dimensions.+-- | A Halfspace in \(d\) dimensions. Note that the intended side of+-- the halfspace is already indicated by the normal vector of the+-- bounding plane. newtype HalfSpace d r = HalfSpace { _boundingPlane :: HyperPlane d  r }                        deriving Generic makeLenses ''HalfSpace  deriving instance (Arity d, Show r)   => Show    (HalfSpace d r)-deriving instance (Arity d, Eq r)     => Eq      (HalfSpace d r) -- deriving instance (NFData r, Arity d) => NFData  (HalfSpace d r) deriving instance Arity d => Functor     (HalfSpace d) deriving instance Arity d => Foldable    (HalfSpace d)@@ -54,6 +55,12 @@ type instance Dimension (HalfSpace d r) = d  deriving instance (Arity d, Arity (d + 1), Fractional r) => IsTransformable (HalfSpace d r)++instance (Arity d, Eq r, Fractional r) => Eq (HalfSpace d r) where+  (HalfSpace h) == (HalfSpace h') = let u = h^.normalVec+                                        v = h'^.normalVec+                                        d = quadrance (u ^+^ v) - (quadrance u)+                                    in h == h' && signum d == 1  -------------------------------------------------------------------------------- 

@@ -25,12 +25,15 @@ 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, Eq r, Fractional r) => Eq (HyperPlane d r) where+  (HyperPlane p u) == h@(HyperPlane _ v) = p `intersects` h && u `isScalarMultipleOf` v++ 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) @@ -68,12 +71,29 @@  pattern Plane     :: Point 3 r -> Vector 3 r -> Plane r pattern Plane p n = HyperPlane p n+{-# COMPLETE Plane #-} +-- | Produces a plane. If r lies counter clockwise of q w.r.t. p then+-- the normal vector of the resulting plane is pointing "upwards".+--+-- >>> from3Points origin (Point3 1 0 0) (Point3 0 1 0)+-- HyperPlane {_inPlane = Point3 [0,0,0], _normalVec = Vector3 [0,0,1]} 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) +instance OnSideUpDownTest (Plane r) where+  -- >>> (Point3 5 5 5) `onSideUpDown` from3Points origin (Point3 1 0 0) (Point3 0 1 0)+  -- Above+  -- >>> (Point3 5 5 (-5)) `onSideUpDown` from3Points origin (Point3 1 0 0) (Point3 0 1 0)+  -- Below+  -- >>> (Point3 5 5 0) `onSideUpDown` from3Points origin (Point3 1 0 0) (Point3 0 1 0)+  -- On+  q `onSideUpDown` (Plane p n) = let v = q .-. p in case (n `dot` v) `compare` 0 of+                                   LT -> Below+                                   EQ -> On+                                   GT -> Above  type instance IntersectionOf (Line 3 r) (Plane r) = [NoIntersection, Point 3 r, Line 3 r] @@ -107,3 +127,23 @@  instance HasSupportingPlane (HyperPlane d r) where   supportingPlane = id+++-- | Given+-- * a plane,+-- * a unit vector in the plane that will represent the y-axis (i.e. the "view up" vector), and+-- * a point in the plane,+--+-- computes the plane coordinates of the given point, using the+-- inPlane point as the origin, the normal vector of the plane as the+-- unit vector in the "z-direction" and the view up vector as the+-- y-axis.+--+-- >>> planeCoordinatesWith (Plane origin (Vector3 0 0 1)) (Vector3 0 1 0) (Point3 10 10 0)+-- Point2 [10.0,10.0]+planeCoordinatesWith       :: Fractional r => Plane r -> Vector 3 r -> Point 3 r -> Point 2 r+planeCoordinatesWith h vup = projectPoint . transformBy (planeCoordinatesTransform h vup)++planeCoordinatesTransform                    :: Num r => Plane r -> Vector 3 r -> Transformation 3 r+planeCoordinatesTransform (HyperPlane o n) v =   rotateTo (Vector3 (v `cross` n) v n)+                                             |.| translation ((-1) *^ toVec o)

@@ -1,7 +1,10 @@ {-# LANGUAGE TemplateHaskell  #-} module Data.Geometry.Interval(                              -- * 1 dimensional Intervals-                               Interval(..)+                               Interval+                             , fromRange, toRange+                             , _Range+                              , pattern OpenInterval                              , pattern ClosedInterval                              , pattern Interval@@ -12,12 +15,13 @@                              , inInterval                              , shiftLeft' +                             , asProperInterval, flipInterval+                              , module Data.Range-                             )-       where+                             ) where  import           Control.DeepSeq-import           Control.Lens (makeLenses, (^.),(%~),(&), Lens')+import           Control.Lens (lens, (^.),(%~),(&), Lens') import           Data.Bifunctor import           Data.Bitraversable import           Data.Ext@@ -34,20 +38,31 @@ --------------------------------------------------------------------------------  -- | An Interval is essentially a 'Data.Range' but with possible payload-newtype Interval a r = GInterval { _unInterval :: Range (r :+ a) }+--+-- We can think of an interval being defined as:+--+-- >>> data Interval a r = Interval (EndPoint (r :+ a)) (EndPoint (r :+ a))+newtype Interval a r = GInterval { toRange :: Range (r :+ a) }                      deriving (Eq,Generic,Arbitrary)-makeLenses ''Interval +_Range :: Lens' (Interval a r) (Range (r :+ a))+_Range = lens toRange (const GInterval)+{-# INLINE _Range #-}++-- | Constrct an interval from a Range+fromRange :: Range (r :+ a) -> Interval a r+fromRange = GInterval+ deriving instance (NFData a, NFData r) => NFData (Interval a r)  instance (Show a, Show r) => Show (Interval a r) where   show ~(Interval l u) = concat [ "Interval (", show l, ") (", show u,")"]  instance Functor (Interval a) where-  fmap = T.fmapDefault+  fmap f (GInterval r) = GInterval $ fmap (first f) r  instance F.Foldable (Interval a) where-  foldMap = T.foldMapDefault+  foldMap f (GInterval r) = foldMap (f . (^.core)) r  instance T.Traversable (Interval a) where   traverse f (GInterval r) = GInterval <$> T.traverse f' r@@ -63,7 +78,7 @@ --  inInterval and inRange is that the extra value is *not* used in the --  comparison with inInterval, whereas it is in inRange. inInterval       :: Ord r => r -> Interval a r -> Bool-x `inInterval` r = x `inRange` (fmap (^.core) $ r^.unInterval )+x `inInterval` r = x `inRange` (fmap (^.core) $ r^._Range )   pattern OpenInterval       :: (r :+ a) -> (r :+ a) -> Interval a r@@ -87,7 +102,8 @@ instance HasStart (Interval a r) where   type StartCore (Interval a r) = r   type StartExtra (Interval a r) = a-  start = unInterval.lower.unEndPoint+  start = _Range.lower.unEndPoint+  {-# INLINE start #-}  class HasEnd t where   type EndCore t@@ -97,7 +113,8 @@ instance HasEnd (Interval a r) where   type EndCore (Interval a r) = r   type EndExtra (Interval a r) = a-  end = unInterval.upper.unEndPoint+  end = _Range.upper.unEndPoint+  {-# INLINE end #-}  type instance Dimension (Interval a r) = 1 type instance NumType   (Interval a r) = r@@ -121,6 +138,17 @@        g (Arg _ x) = x +-- | Shifts the interval to the left by delta+shiftLeft'       :: Num r => r -> Interval a r -> Interval a r+shiftLeft' delta = fmap (subtract delta) -shiftLeft'   :: Num r => r -> Interval a r -> Interval a r-shiftLeft' x = fmap (subtract x)++-- | Makes sure the start and endpoint are oriented such that the+-- starting value is smaller than the ending value.+asProperInterval                                     :: Ord r => Interval p r -> Interval p r+asProperInterval i | (i^.start.core) > (i^.end.core) = flipInterval i+                   | otherwise                       = i++-- | Flips the start and endpoint of the interval.+flipInterval :: Interval a r -> Interval a r+flipInterval = _Range %~ \(Range s t) -> Range t s

@@ -55,7 +55,7 @@                  => [i] -> IntervalTree i r fromIntervals is = foldr insert (createTree pts) is   where-    endPoints (toRange -> Range' a b) = [a,b]+    endPoints (asRange -> Range' a b) = [a,b]     pts = List.sort . concatMap endPoints $ is  -- | Lists the intervals. We don't guarantee anything about the order@@ -100,7 +100,7 @@                           => i -> IntervalTree i r -> IntervalTree i r insert i (IntervalTree t) = IntervalTree $ insert' t   where-    ri@(Range a b) = toRange i+    ri@(Range a b) = asRange i      insert' Nil = Nil     insert' (Internal l nd@(_splitPoint -> m) r)@@ -119,7 +119,7 @@           => i -> IntervalTree i r -> IntervalTree i r delete i (IntervalTree t) = IntervalTree $ delete' t   where-    ri@(Range a b) = toRange i+    ri@(Range a b) = asRange i      delete' Nil = Nil     delete' (Internal l nd@(_splitPoint -> m) r)@@ -137,15 +137,13 @@  -- | Anything that looks like an interval class IntervalLike i where-  toRange :: i -> Range (NumType i)+  asRange :: i -> Range (NumType i)  instance IntervalLike (Range r) where-  toRange = id+  asRange = id  instance IntervalLike (Interval p r) where-  toRange = fmap (^.core) . _unInterval--+  asRange = fmap (^.core) . toRange  -------------------------------------------------------------------------------- 

@@ -80,7 +80,7 @@         _     -> coRec NoIntersection       _     -> error "intersect; line x boundary rect; absurd"     where-      (t,r,b,l) = sides' rect+      Sides t r b l = sides' rect       ints = map (\s -> sl `intersect` toSL s) [t,r,b,l]        nub' = map L.head . L.group . L.sort

@@ -53,6 +53,8 @@ instance (Arity d, Eq r, Fractional r) => Eq (Line d r) where   l@(Line p _) == m = l `isParallelTo` m && p `onLine` m ++ instance (Arbitrary r, Arity d, Num r, Eq r) => Arbitrary (Line d r) where   arbitrary = do p <- arbitrary                  q <- suchThat arbitrary (/= p)@@ -215,29 +217,34 @@     :& (H $ \_              -> Nothing)    -- l is a vertical line (through x=0)     :& RNil + -- | Result of a side test data SideTestUpDown = Below | On | Above deriving (Show,Read,Eq,Ord) --- | Given a point q and a line l, compute to which side of l q lies. For--- vertical lines the left side of the line is interpeted as below.------ >>> Point2 10 10 `onSideUpDown` (lineThrough origin $ Point2 10 5)--- Above--- >>> Point2 10 10 `onSideUpDown` (lineThrough origin $ Point2 (-10) 5)--- Above--- >>> Point2 5 5 `onSideUpDown` (verticalLine 10)--- Below--- >>> Point2 5 5 `onSideUpDown` (lineThrough origin $ Point2 (-3) (-3))--- On-onSideUpDown                :: (Ord r, Num r) => Point 2 r -> Line 2 r -> SideTestUpDown-q `onSideUpDown` (Line p v) = let r    =  p .+^ v-                                  f z         = (z^.xCoord, -z^.yCoord)-                                  minBy g a b = F.minimumBy (comparing g) [a,b]-                                  maxBy g a b = F.maximumBy (comparing g) [a,b]-                              in case ccw (minBy f p r) (maxBy f p r) q of-                                   CCW      -> Above-                                   CW       -> Below-                                   CoLinear -> On+class OnSideUpDownTest t where+  onSideUpDown :: (d ~ Dimension t, r ~ NumType t, Ord r, Num r)+               => Point d r -> t -> SideTestUpDown++instance OnSideUpDownTest (Line 2 r) where+  -- | Given a point q and a line l, compute to which side of l q lies. For+  -- vertical lines the left side of the line is interpeted as below.+  --+  -- >>> Point2 10 10 `onSideUpDown` (lineThrough origin $ Point2 10 5)+  -- Above+  -- >>> Point2 10 10 `onSideUpDown` (lineThrough origin $ Point2 (-10) 5)+  -- Above+  -- >>> Point2 5 5 `onSideUpDown` (verticalLine 10)+  -- Below+  -- >>> Point2 5 5 `onSideUpDown` (lineThrough origin $ Point2 (-3) (-3))+  -- On+  q `onSideUpDown` (Line p v) = let r    =  p .+^ v+                                    f z         = (z^.xCoord, -z^.yCoord)+                                    minBy g a b = F.minimumBy (comparing g) [a,b]+                                    maxBy g a b = F.maximumBy (comparing g) [a,b]+                                in case ccw (minBy f p r) (maxBy f p r) q of+                                     CCW      -> Above+                                     CW       -> Below+                                     CoLinear -> On  -- | Result of a side test data SideTest = LeftSide | OnLine | RightSide deriving (Show,Read,Eq,Ord)

@@ -14,6 +14,7 @@                                 , pattern LineSegment                                 , pattern LineSegment'                                 , pattern ClosedLineSegment+                                , pattern OpenLineSegment                                 , endPoints                                  , _SubLine@@ -26,6 +27,8 @@                                 , segmentLength                                 , sqDistanceToSeg, sqDistanceToSegArg                                 , flipSegment++                                , interpolate                                 ) where  import           Control.Arrow ((&&&))@@ -45,7 +48,7 @@ import           Data.Vinyl import           Data.Vinyl.CoRec import           GHC.TypeLits-import           Test.QuickCheck+import           Test.QuickCheck(Arbitrary(..))  -------------------------------------------------------------------------------- -- * d-dimensional LineSegments@@ -77,12 +80,16 @@ pattern LineSegment'      s t <- ((^.start) &&& (^.end) -> (s,t)) {-# COMPLETE LineSegment' #-} -pattern ClosedLineSegment     :: Point d r :+ p-                              -> Point d r :+ p-                              -> LineSegment d p r+pattern ClosedLineSegment     :: Point d r :+ p -> Point d r :+ p -> LineSegment d p r pattern ClosedLineSegment s t = GLineSegment (ClosedInterval s t) {-# COMPLETE ClosedLineSegment #-} +pattern OpenLineSegment     :: Point d r :+ p -> Point d r :+ p -> LineSegment d p r+pattern OpenLineSegment s t = GLineSegment (OpenInterval s t)+{-# COMPLETE OpenLineSegment #-}+++ type instance Dimension (LineSegment d p r) = d type instance NumType   (LineSegment d p r) = r @@ -187,7 +194,6 @@       :& (H $ coRec . subLineToSegment)       :& RNil - instance (Ord r, Fractional r) =>          (LineSegment 2 p r) `IsIntersectableWith` (Line 2 r) where   nonEmptyIntersection = defaultNonEmptyIntersection@@ -292,3 +298,18 @@ --                   (q&unEndPoint %~ ff)  -- ss'' = ss'^._SubLine++-- | Linearly interpolate the two endpoints with a value in the range [0,1]+--+-- >>> interpolate 0.5 $ ClosedLineSegment (ext $ origin) (ext $ Point2 10.0 10.0)+-- Point2 [5.0,5.0]+-- >>> interpolate 0.1 $ ClosedLineSegment (ext $ origin) (ext $ Point2 10.0 10.0)+-- Point2 [1.0,1.0]+-- >>> interpolate 0 $ ClosedLineSegment (ext $ origin) (ext $ Point2 10.0 10.0)+-- Point2 [0.0,0.0]+-- >>> interpolate 1 $ ClosedLineSegment (ext $ origin) (ext $ Point2 10.0 10.0)+-- Point2 [10.0,10.0]+interpolate                      :: (Fractional r, Arity d) => r -> LineSegment d p r -> Point d r+interpolate t (LineSegment' p q) = Point $ (asV p ^* (1-t)) ^+^ (asV q ^* t)+  where+    asV = (^.core.vector)

@@ -0,0 +1,79 @@+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.Matrix+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- type-indexed matrices.+--+--------------------------------------------------------------------------------+module Data.Geometry.Matrix(+    Matrix(Matrix)+  , identityMatrix++  , multM+  , mult++  , Invertible(..)+  , HasDeterminant(..)+  ) where++import           Control.Lens (imap)+import           Data.Geometry.Matrix.Internal (mkRow)+import           Data.Geometry.Vector+import           Linear.Matrix ((!*),(!*!))+import qualified Linear.Matrix as Lin+import           Unsafe.Coerce (unsafeCoerce)++--------------------------------------------------------------------------------+-- * Matrices++-- | a matrix of n rows, each of m columns, storing values of type r+newtype Matrix n m r = Matrix (Vector n (Vector m r))++deriving instance (Show r, Arity n, Arity m) => Show (Matrix n m r)+deriving instance (Eq r, Arity n, Arity m)   => Eq (Matrix n m r)+deriving instance (Ord r, Arity n, Arity m)  => Ord (Matrix n m r)+deriving instance (Arity n, Arity m)         => Functor (Matrix n m)+deriving instance (Arity n, Arity m)         => Foldable (Matrix n m)+deriving instance (Arity n, Arity m)         => Traversable (Matrix n m)++multM :: (Arity r, Arity c, Arity c', Num a) => Matrix r c a -> Matrix c c' a -> Matrix r c' a+(Matrix a) `multM` (Matrix b) = Matrix $ a !*! b++mult :: (Arity m, Arity n, Num r) => Matrix n m r -> Vector m r -> Vector n r+(Matrix m) `mult` v = m !* v++-- | Produces the Identity Matrix+identityMatrix :: (Arity d, Num r) => Matrix d d r+identityMatrix = Matrix $ imap mkRow (pure 1)++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+++class Arity d => HasDeterminant d where+  det :: Num r => Matrix d d r -> r++instance HasDeterminant 1 where+  det (Matrix (Vector1 (Vector1 x))) = x+instance HasDeterminant 2 where+  det = Lin.det22 . unsafeCoerce+instance HasDeterminant 3 where+  det = Lin.det33 . unsafeCoerce+instance HasDeterminant 4 where+  det = Lin.det44 . unsafeCoerce

@@ -0,0 +1,14 @@+{-# LANGUAGE Unsafe #-}+module Data.Geometry.Matrix.Internal where++import           Control.Lens (set)+import           Data.Geometry.Vector+import qualified Data.Vector.Fixed as FV++--------------------------------------------------------------------------------+-- * Helper functions to easily create matrices++-- | Creates a row with zeroes everywhere, except at position i, where the+-- value is the supplied value.+mkRow     :: forall d r. (Arity d, Num r) => Int -> r -> Vector d r+mkRow i x = set (FV.element i) x zero

@@ -24,7 +24,6 @@ import           Data.Geometry.PlanarSubdivision.Basic import           Data.Geometry.PlanarSubdivision.Merge import           Data.Geometry.Polygon-import qualified Data.PlaneGraph as PG import           Data.Proxy  

@@ -423,13 +423,23 @@                             d''      = PG.nextIncidentEdge d' g                         in g^.dataOf d'' --- | All incoming edges incident to vertex v, in counterclockwise order around v.+-- | All edges incident to vertex v in incoming direction+-- (i.e. pointing into v) in counterclockwise order around v.+--+-- running time: \(O(k)\), where \(k) is the total number of incident edges of v incomingEdges      :: VertexId' s -> PlanarSubdivision s v e f r -> V.Vector (Dart s)-incomingEdges v ps = V.filter (not . isPositive) $ incidentEdges v ps+incomingEdges v ps = orient <$> incidentEdges v ps+  where+    orient d = if headOf d ps == v then d else twin d --- | All outgoing edges incident to vertex v, in counterclockwise order around v.+-- | All edges incident to vertex v in outgoing direction+-- (i.e. pointing away from v) in counterclockwise order around v.+--+-- running time: \(O(k)\), where \(k) is the total number of incident edges of v outgoingEdges      :: VertexId' s -> PlanarSubdivision s v e f r  -> V.Vector (Dart s)-outgoingEdges v ps = V.filter isPositive $ incidentEdges v ps+outgoingEdges v ps = orient <$> incidentEdges v ps+  where+    orient d = if tailOf d ps == v then d else twin d   -- | Gets the neighbours of a particular vertex, in counterclockwise order@@ -437,10 +447,7 @@ -- -- running time: \(O(k)\), where \(k\) is the output size neighboursOf      :: VertexId' s -> PlanarSubdivision s v e f r -> V.Vector (VertexId' s)-neighboursOf v ps = otherVtx <$> incidentEdges v ps-  where-    otherVtx d = let u = tailOf d ps in if u == v then headOf d ps else u-+neighboursOf v ps = flip tailOf ps <$> incomingEdges v ps  -- | The face to the left of the dart --

@@ -13,18 +13,17 @@ module Data.Geometry.Point( Point(..)                           , origin, vector                           , pointFromList--                          , coord , unsafeCoord-                           , projectPoint +                          , pattern Point1                           , pattern Point2                           , pattern Point3                           , xCoord, yCoord, zCoord                            , PointFunctor(..) -                          , CCW(..), ccw, ccw'+                          , CCW, ccw, ccw'+                          , pattern CCW, pattern CW, pattern CoLinear                            , ccwCmpAround, cwCmpAround, ccwCmpAroundWith, cwCmpAroundWith                           , sortAround, insertIntoCyclicOrder@@ -34,394 +33,12 @@                           , cmpByDistanceTo                            , squaredEuclideanDist, euclideanDist-                          ) where -import           Control.DeepSeq-import           Control.Lens-import           Data.Aeson-import qualified Data.CircularList as C-import qualified Data.CircularList.Util as CU-import           Data.Ext-import qualified Data.Foldable as F-import           Data.Geometry.Properties-import           Data.Geometry.Vector-import qualified Data.Geometry.Vector as Vec-import qualified Data.List as L-import           Data.Ord (comparing)-import           Data.Proxy-import           GHC.Generics (Generic)-import           GHC.TypeLits-import           Test.QuickCheck (Arbitrary)-import           Text.ParserCombinators.ReadP (ReadP, string,pfail)-import           Text.ParserCombinators.ReadPrec (lift)-import           Text.Read (Read(..),readListPrecDefault, readPrec_to_P,minPrec) ------------------------------------------------------------------------------------- $setup--- >>> :{--- let myVector :: Vector 3 Int---     myVector = Vector3 1 2 3---     myPoint = Point myVector--- :}-------------------------------------------------------------------------------------- * A d-dimensional Point---- | A d-dimensional point.-newtype Point d r = Point { toVec :: Vector d r } deriving (Generic)--instance (Show r, Arity d) => Show (Point d r) where-  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)-deriving instance Arity d                => Functor (Point d)-deriving instance Arity d                => Foldable (Point d)-deriving instance Arity d                => Traversable (Point d)-deriving instance (Arity d, NFData r)    => NFData (Point d r)-deriving instance (Arity d, Arbitrary r) => Arbitrary (Point d r)--type instance NumType (Point d r) = r-type instance Dimension (Point d r) = d--instance Arity d =>  Affine (Point d) where-  type Diff (Point d) = Vector d--  p .-. q = toVec p ^-^ toVec q-  p .+^ v = Point $ toVec p ^+^ v--instance (FromJSON r, Arity d, KnownNat d) => FromJSON (Point d r) where-  parseJSON = fmap Point . parseJSON--instance (ToJSON r, Arity d) => ToJSON (Point d r) where-  toJSON     = toJSON     . toVec-  toEncoding = toEncoding . toVec---- | Point representing the origin in d dimensions------ >>> origin :: Point 4 Int--- Point4 [0,0,0,0]-origin :: (Arity d, Num r) => Point d r-origin = Point $ pure 0----- ** Accessing points---- | Lens to access the vector corresponding to this point.------ >>> (Point3 1 2 3) ^. vector--- Vector3 [1,2,3]--- >>> origin & vector .~ Vector3 1 2 3--- Point3 [1,2,3]-vector :: Lens' (Point d r) (Vector d r)-vector = lens toVec (const Point)----- | Get the coordinate in a given dimension. This operation is unsafe in the--- sense that no bounds are checked. Consider using `coord` instead.--------- >>> Point3 1 2 3 ^. unsafeCoord 2--- 2-unsafeCoord   :: Arity d => Int -> Lens' (Point d r) r-unsafeCoord i = vector . singular (ix (i-1))-                -- Points are 1 indexed, vectors are 0 indexed---- | Get the coordinate in a given dimension------ >>> Point3 1 2 3 ^. coord (C :: C 2)--- 2--- >>> Point3 1 2 3 & coord (C :: C 1) .~ 10--- Point3 [10,2,3]--- >>> Point3 1 2 3 & coord (C :: C 3) %~ (+1)--- Point3 [1,2,4]-coord   :: forall proxy i d r. (1 <= i, i <= d, ((i - 1) + 1) ~ i-                               , Arity (i - 1), Arity d-                               ) => proxy i -> Lens' (Point d r) r-coord _ = vector . Vec.element (Proxy :: Proxy (i-1))-{-# INLINABLE coord #-}----- somehow these rules don't fire--- {-# SPECIALIZE coord :: C 1 -> Lens' (Point 2 r) r#-}--- {-# SPECIALIZE coord :: C 2 -> Lens' (Point 2 r) r#-}----- | Constructs a point from a list of coordinates------ >>> pointFromList [1,2,3] :: Maybe (Point 3 Int)--- Just Point3 [1,2,3]-pointFromList :: Arity d => [r] -> Maybe (Point d r)-pointFromList = fmap Point . Vec.vectorFromList----- | Project a point down into a lower dimension.-projectPoint :: (Arity i, Arity d, i <= d) => Point d r -> Point i r-projectPoint = Point . prefix . toVec------------------------------------------------------------------------------------- * Convenience functions to construct 2 and 3 dimensional points----- | We provide pattern synonyms Point2 and Point3 for 2 and 3 dimensional points. i.e.--- we can write:------ >>> :{---   let---     f              :: Point 2 r -> r---     f (Point2 x y) = x---   in f (Point2 1 2)--- :}--- 1------ if we want.-pattern Point2       :: r -> r -> Point 2 r-pattern Point2 x y = Point (Vector2 x y)-{-# COMPLETE Point2 #-}---- | Similarly, we can write:------ >>> :{---   let---     g                :: Point 3 r -> r---     g (Point3 x y z) = z---   in g myPoint--- :}--- 3-pattern Point3       :: r -> r -> r -> Point 3 r-pattern Point3 x y z = (Point (Vector3 x y z))-{-# COMPLETE Point3 #-}---- | Shorthand to access the first coordinate C 1------ >>> Point3 1 2 3 ^. xCoord--- 1--- >>> Point2 1 2 & xCoord .~ 10--- Point2 [10,2]-xCoord :: (1 <= d, Arity d) => Lens' (Point d r) r-xCoord = coord (C :: C 1)-{-# INLINABLE xCoord #-}---- | Shorthand to access the second coordinate C 2------ >>> Point2 1 2 ^. yCoord--- 2--- >>> Point3 1 2 3 & yCoord %~ (+1)--- Point3 [1,3,3]-yCoord :: (2 <= d, Arity d) => Lens' (Point d r) r-yCoord = coord (C :: C 2)-{-# INLINABLE yCoord #-}---- | Shorthand to access the third coordinate C 3------ >>> Point3 1 2 3 ^. zCoord--- 3--- >>> Point3 1 2 3 & zCoord %~ (+1)--- Point3 [1,2,4]-zCoord :: (3 <= d, Arity d) => Lens' (Point d r) r-zCoord = coord (C :: C 3)-{-# INLINABLE zCoord #-}-------------------------------------------------------------------------------------- * Point Functors---- | Types that we can transform by mapping a function on each point in the structure-class PointFunctor g where-  pmap :: (Point (Dimension (g r)) r -> Point (Dimension (g s)) s) -> g r -> g s--  -- pemap :: (d ~ Dimension (g r)) => (Point d r :+ p -> Point d s :+ p) -> g r -> g s-  -- pemap =--instance PointFunctor (Point d) where-  pmap f = f-------------------------------------------------------------------------------------- * Functions specific to Two Dimensional points--data CCW = CCW | CoLinear | CW-         deriving (Show,Eq)---- | Given three points p q and r determine the orientation when going from p to r via q.-ccw :: (Ord r, Num r) => Point 2 r -> Point 2 r -> Point 2 r -> CCW-ccw p q r = case z `compare` 0 of-              LT -> CW-              GT -> CCW-              EQ -> CoLinear-     where-       Vector2 ux uy = q .-. p-       Vector2 vx vy = r .-. p-       z             = ux * vy - uy * vx---- | Given three points p q and r determine the orientation when going from p to r via q.-ccw' :: (Ord r, Num r) => Point 2 r :+ a -> Point 2 r :+ b -> Point 2 r :+ c -> CCW-ccw' p q r = ccw (p^.core) (q^.core) (r^.core)---- | Sort the points arround the given point p in counter clockwise order with--- respect to the rightward horizontal ray starting from p.  If two points q--- and r are colinear with p, the closest one to p is reported first.--- running time: O(n log n)-sortAround   :: (Ord r, Num r)-             => Point 2 r :+ q -> [Point 2 r :+ p] -> [Point 2 r :+ p]-sortAround c = L.sortBy (ccwCmpAround c <> cmpByDistanceTo c)----- | Quadrants of two dimensional points. in CCW order-data Quadrant = TopRight | TopLeft | BottomLeft | BottomRight-              deriving (Show,Read,Eq,Ord,Enum,Bounded)---- | Quadrants around point c; quadrants are closed on their "previous"--- boundary (i..e the boundary with the previous quadrant in the CCW order),--- open on next boundary. The origin itself is assigned the topRight quadrant-quadrantWith                   :: (Ord r, 1 <= d, 2 <= d, Arity d)-                               => Point d r :+ q -> Point d r :+ p -> Quadrant-quadrantWith (c :+ _) (p :+ _) = case ( (c^.xCoord) `compare` (p^.xCoord)-                                      , (c^.yCoord) `compare` (p^.yCoord) ) of-                                   (EQ, EQ) -> TopRight-                                   (LT, EQ) -> TopRight-                                   (LT, LT) -> TopRight-                                   (EQ, LT) -> TopLeft-                                   (GT, LT) -> TopLeft-                                   (GT, EQ) -> BottomLeft-                                   (GT, GT) -> BottomLeft-                                   (EQ, GT) -> BottomRight-                                   (LT, GT) -> BottomRight---- | Quadrants with respect to the origin-quadrant :: (Ord r, Num r, 1 <= d, 2 <= d, Arity d) => Point d r :+ p -> Quadrant-quadrant = quadrantWith (ext origin)---- | Given a center point c, and a set of points, partition the points into--- quadrants around c (based on their x and y coordinates). The quadrants are--- reported in the order topLeft, topRight, bottomLeft, bottomRight. The points--- are in the same order as they were in the original input lists.--- Points with the same x-or y coordinate as p, are "rounded" to above.-partitionIntoQuadrants       :: (Ord r, 1 <= d, 2 <= d, Arity d)-                             => Point d r :+ q-                             -> [Point d r :+ p]-                             -> ( [Point d r :+ p], [Point d r :+ p]-                                , [Point d r :+ p], [Point d r :+ p]-                                )-partitionIntoQuadrants c pts = (topL, topR, bottomL, bottomR)-  where-    (below',above')   = L.partition (on yCoord) pts-    (bottomL,bottomR) = L.partition (on xCoord) below'-    (topL,topR)       = L.partition (on xCoord) above'--    on l q       = q^.core.l < c^.core.l------ | Given a zero vector z, a center c, and two points p and q,--- compute the ccw ordering of p and q around c with this vector as zero--- direction.------ pre: the points p,q /= c-ccwCmpAroundWith                              :: (Ord r, Num r)-                                              => Vector 2 r-                                              -> Point 2 r :+ c-                                              -> Point 2 r :+ a -> Point 2 r :+ b-                                              -> Ordering-ccwCmpAroundWith z@(Vector2 zx zy) (c :+ _) (q :+ _) (r :+ _) =-    case (ccw c a q, ccw c a r) of-      (CCW,CCW)      -> cmp-      (CCW,CW)       -> LT-      (CCW,CoLinear) | onZero r  -> GT-                     | otherwise -> LT--      (CW, CCW)      -> GT-      (CW, CW)       -> cmp-      (CW, CoLinear) -> GT--      (CoLinear, CCW) | onZero q  -> LT-                      | otherwise -> GT--      (CoLinear, CW)      -> LT-      (CoLinear,CoLinear) -> case (onZero q, onZero r) of-                               (True, True)   -> EQ-                               (False, False) -> EQ-                               (True, False)  -> LT-                               (False, True)  -> GT-  where-    a = c .+^ z-    b = c .+^ Vector2 (-zy) zx-    -- b is on a perpendicular vector to z--    -- test if the point lies on the ray defined by z, starting in c-    onZero d = case ccw c b d of-                 CCW      -> False-                 CW       -> True-                 CoLinear -> True -- this shouldh appen only when you ask for c itself--    cmp = case ccw c q r of-            CCW      -> LT-            CW       -> GT-            CoLinear -> EQ---- | Given a zero vector z, a center c, and two points p and q,--- compute the cw ordering of p and q around c with this vector as zero--- direction.------ pre: the points p,q /= c-cwCmpAroundWith     :: (Ord r, Num r)-                    => Vector 2 r-                    -> Point 2 r :+ a-                    -> Point 2 r :+ b -> Point 2 r :+ c-                    -> Ordering-cwCmpAroundWith z c = flip (ccwCmpAroundWith z c)------ | Compare by distance to the first argument-cmpByDistanceTo              :: (Ord r, Num r, Arity d)-                             => Point d r :+ c -> Point d r :+ p -> Point d r :+ q -> Ordering-cmpByDistanceTo (c :+ _) p q = comparing (squaredEuclideanDist c) (p^.core) (q^.core)----- | Counter clockwise ordering of the points around c. Points are ordered with--- respect to the positive x-axis.-ccwCmpAround :: (Num r, Ord r)-             => Point 2 r :+ qc -> Point 2 r :+ p -> Point 2 r :+ q -> Ordering-ccwCmpAround = ccwCmpAroundWith (Vector2 1 0)---- | Clockwise ordering of the points around c. Points are ordered with--- respect to the positive x-axis.-cwCmpAround :: (Num r, Ord r)-            => Point 2 r :+ qc -> Point 2 r :+ p -> Point 2 r :+ q -> Ordering-cwCmpAround = cwCmpAroundWith (Vector2 1 0)----- | Given a center c, a new point p, and a list of points ps, sorted in--- counter clockwise order around c. Insert p into the cyclic order. The focus--- of the returned cyclic list is the new point p.------ running time: O(n)-insertIntoCyclicOrder   :: (Ord r, Num r)-                        => Point 2 r :+ q -> Point 2 r :+ p-                        -> C.CList (Point 2 r :+ p) -> C.CList (Point 2 r :+ p)-insertIntoCyclicOrder c = CU.insertOrdBy (ccwCmpAround c <> cmpByDistanceTo c)----- | Squared Euclidean distance between two points-squaredEuclideanDist :: (Num r, Arity d) => Point d r -> Point d r -> r-squaredEuclideanDist = qdA+                          , AsAPoint(..), coord, unsafeCoord, vector'+                          ) where --- | Euclidean distance between two points-euclideanDist :: (Floating r, Arity d) => Point d r -> Point d r -> r-euclideanDist = distanceA+import Data.Geometry.Point.Class+import Data.Geometry.Point.Internal hiding (coord, unsafeCoord)+import Data.Geometry.Point.Orientation.Degenerate+import Data.Geometry.Point.Quadrants

@@ -0,0 +1,66 @@+module Data.Geometry.Point.Class where++import           Control.Lens+import           Data.Geometry.Point.Internal (Point)+import qualified Data.Geometry.Point.Internal as Internal+import           Data.Geometry.Vector+import           GHC.TypeNats++--------------------------------------------------------------------------------++-- $setup+-- >>> import Data.Geometry.Point.Internal (pattern Point2, pattern Point3)++class ToAPoint point d r where+  toPoint   :: Prism' (point d r) (Point d r)++class AsAPoint p where+  asAPoint :: Lens (p d r) (p d' r') (Point d r) (Point d' r')++vector' :: AsAPoint p => Lens (p d r) (p d r') (Vector d r) (Vector d r')+vector' = asAPoint . lens Internal.toVec (const Internal.Point)++coord   :: (1 <= i, i <= d, KnownNat i, Arity d, AsAPoint p) => proxy i -> Lens' (p d r) r+coord i = asAPoint.Internal.coord i++unsafeCoord   :: (Arity d, AsAPoint p) => Int -> Lens' (p d r) r+unsafeCoord i = asAPoint.Internal.unsafeCoord i++instance ToAPoint Point d r where+  toPoint = prism' id Just++instance AsAPoint Point where+  asAPoint = id+++++-- | Shorthand to access the first coordinate C 1+--+-- >>> Point3 1 2 3 ^. xCoord+-- 1+-- >>> Point2 1 2 & xCoord .~ 10+-- Point2 [10,2]+xCoord :: (1 <= d, Arity d, AsAPoint point) => Lens' (point d r) r+xCoord = coord (C :: C 1)+{-# INLINABLE xCoord #-}++-- | Shorthand to access the second coordinate C 2+--+-- >>> Point2 1 2 ^. yCoord+-- 2+-- >>> Point3 1 2 3 & yCoord %~ (+1)+-- Point3 [1,3,3]+yCoord :: (2 <= d, Arity d, AsAPoint point) => Lens' (point d r) r+yCoord = coord (C :: C 2)+{-# INLINABLE yCoord #-}++-- | Shorthand to access the third coordinate C 3+--+-- >>> Point3 1 2 3 ^. zCoord+-- 3+-- >>> Point3 1 2 3 & zCoord %~ (+1)+-- Point3 [1,2,4]+zCoord :: (3 <= d, Arity d,AsAPoint point) => Lens' (point d r) r+zCoord = coord (C :: C 3)+{-# INLINABLE zCoord #-}

@@ -0,0 +1,252 @@+{-# LANGUAGE ScopedTypeVariables  #-}+{-# LANGUAGE UndecidableInstances #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.Point+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- \(d\)-dimensional points.+--+--------------------------------------------------------------------------------+module Data.Geometry.Point.Internal+  ( Point(..)+  , origin, vector+  , pointFromList++  , coord , unsafeCoord++  , projectPoint++  , pattern Point1+  , pattern Point2+  , pattern Point3+  , PointFunctor(..)++  , cmpByDistanceTo+  , squaredEuclideanDist, euclideanDist+  ) where++import           Control.DeepSeq+import           Control.Lens+import           Data.Aeson+import           Data.Ext+import qualified Data.Foldable as F+import           Data.Geometry.Properties+import           Data.Geometry.Vector+import qualified Data.Geometry.Vector as Vec+import           Data.Hashable+import           Data.Ord (comparing)+import           Data.Proxy+import           GHC.Generics (Generic)+import           GHC.TypeLits+import           System.Random (Random(..))+import           Test.QuickCheck (Arbitrary)+import           Text.ParserCombinators.ReadP (ReadP, string,pfail)+import           Text.ParserCombinators.ReadPrec (lift)+import           Text.Read (Read(..),readListPrecDefault, readPrec_to_P,minPrec)+++--------------------------------------------------------------------------------+-- $setup+-- >>> :{+-- let myVector :: Vector 3 Int+--     myVector = Vector3 1 2 3+--     myPoint = Point myVector+-- :}+++--------------------------------------------------------------------------------+-- * A d-dimensional Point++-- | A d-dimensional point.+newtype Point d r = Point { toVec :: Vector d r } deriving (Generic)++instance (Show r, Arity d) => Show (Point d r) where+  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)+deriving instance Arity d                => Functor (Point d)+deriving instance Arity d                => Foldable (Point d)+deriving instance Arity d                => Traversable (Point d)+deriving instance (Arity d, NFData r)    => NFData (Point d r)+deriving instance (Arity d, Arbitrary r) => Arbitrary (Point d r)+deriving instance (Arity d, Hashable r)  => Hashable (Point d r)+deriving instance (Arity d, Random r)    => Random (Point d r)+++type instance NumType (Point d r) = r+type instance Dimension (Point d r) = d++instance Arity d =>  Affine (Point d) where+  type Diff (Point d) = Vector d++  p .-. q = toVec p ^-^ toVec q+  p .+^ v = Point $ toVec p ^+^ v++instance (FromJSON r, Arity d, KnownNat d) => FromJSON (Point d r) where+  parseJSON = fmap Point . parseJSON++instance (ToJSON r, Arity d) => ToJSON (Point d r) where+  toJSON     = toJSON     . toVec+  toEncoding = toEncoding . toVec++-- | Point representing the origin in d dimensions+--+-- >>> origin :: Point 4 Int+-- Point4 [0,0,0,0]+origin :: (Arity d, Num r) => Point d r+origin = Point $ pure 0+++-- ** Accessing points++-- | Lens to access the vector corresponding to this point.+--+-- >>> (Point3 1 2 3) ^. vector+-- Vector3 [1,2,3]+-- >>> origin & vector .~ Vector3 1 2 3+-- Point3 [1,2,3]+vector :: Lens' (Point d r) (Vector d r)+vector = lens toVec (const Point)+{-# INLINABLE vector #-}++-- | Get the coordinate in a given dimension. This operation is unsafe in the+-- sense that no bounds are checked. Consider using `coord` instead.+--+--+-- >>> Point3 1 2 3 ^. unsafeCoord 2+-- 2+unsafeCoord   :: Arity d => Int -> Lens' (Point d r) r+unsafeCoord i = vector . singular (ix (i-1))+                -- Points are 1 indexed, vectors are 0 indexed+{-# INLINABLE unsafeCoord #-}++-- | Get the coordinate in a given dimension+--+-- >>> Point3 1 2 3 ^. coord (C :: C 2)+-- 2+-- >>> Point3 1 2 3 & coord (C :: C 1) .~ 10+-- Point3 [10,2,3]+-- >>> Point3 1 2 3 & coord (C :: C 3) %~ (+1)+-- Point3 [1,2,4]+coord   :: forall proxy i d r. (1 <= i, i <= d, Arity d, KnownNat i)+        => proxy i -> Lens' (Point d r) r+coord _ = unsafeCoord $ fromIntegral (natVal $ C @i)+{-# INLINABLE coord #-}++ -- somehow these rules don't fire+-- {-# SPECIALIZE coord :: C 1 -> Lens' (Point 2 r) r#-}+-- {-# SPECIALIZE coord :: C 2 -> Lens' (Point 2 r) r#-}+-- {-# SPECIALIZE coord :: C 3 -> Lens' (Point 3 r) r#-}+++-- | Constructs a point from a list of coordinates. The length of the+-- list has to match the dimension exactly.+--+-- >>> pointFromList [1,2,3] :: Maybe (Point 3 Int)+-- Just Point3 [1,2,3]+-- >>> pointFromList [1] :: Maybe (Point 3 Int)+-- Nothing+-- >>> pointFromList [1,2,3,4] :: Maybe (Point 3 Int)+-- Nothing+pointFromList :: Arity d => [r] -> Maybe (Point d r)+pointFromList = fmap Point . Vec.vectorFromList+++-- | Project a point down into a lower dimension.+projectPoint :: (Arity i, Arity d, i <= d) => Point d r -> Point i r+projectPoint = Point . prefix . toVec++--------------------------------------------------------------------------------+-- * Convenience functions to construct 1, 2 and 3 dimensional points++-- | We provide pattern synonyms for 1, 2 and 3 dimensional points. i.e.+-- we can write:+--+--+-- >>> :{+--   let+--     f            :: Num r => Point 1 r -> r+--     f (Point1 x) = x + 1+--   in f (Point1 1)+-- :}+-- 2+pattern Point1   :: r -> Point 1 r+pattern Point1 x = Point (Vector1 x)+{-# COMPLETE Point1 #-}+++-- | Pattern synonym for 2 dimensional points+--+-- >>> :{+--   let+--     f              :: Point 2 r -> r+--     f (Point2 x y) = x+--   in f (Point2 1 2)+-- :}+-- 1+pattern Point2       :: r -> r -> Point 2 r+pattern Point2 x y = Point (Vector2 x y)+{-# COMPLETE Point2 #-}++-- | Similarly, we can write:+--+-- >>> :{+--   let+--     g                :: Point 3 r -> r+--     g (Point3 x y z) = z+--   in g myPoint+-- :}+-- 3+pattern Point3       :: r -> r -> r -> Point 3 r+pattern Point3 x y z = (Point (Vector3 x y z))+{-# COMPLETE Point3 #-}++--------------------------------------------------------------------------------+-- * Point Functors++-- | Types that we can transform by mapping a function on each point in the structure+class PointFunctor g where+  pmap :: (Point (Dimension (g r)) r -> Point (Dimension (g s)) s) -> g r -> g s++  -- pemap :: (d ~ Dimension (g r)) => (Point d r :+ p -> Point d s :+ p) -> g r -> g s+  -- pemap =++instance PointFunctor (Point d) where+  pmap f = f+++--------------------------------------------------------------------------------+-- * Functions specific to Two Dimensional points++-- | Compare by distance to the first argument+cmpByDistanceTo              :: (Ord r, Num r, Arity d)+                             => Point d r :+ c -> Point d r :+ p -> Point d r :+ q -> Ordering+cmpByDistanceTo (c :+ _) p q = comparing (squaredEuclideanDist c) (p^.core) (q^.core)+++++-- | Squared Euclidean distance between two points+squaredEuclideanDist :: (Num r, Arity d) => Point d r -> Point d r -> r+squaredEuclideanDist = qdA++-- | Euclidean distance between two points+euclideanDist :: (Floating r, Arity d) => Point d r -> Point d r -> r+euclideanDist = distanceA

@@ -0,0 +1,31 @@+module Data.Geometry.Point.Orientation where++import Algorithms.Geometry.SoS.Orientation+import Algorithms.Geometry.SoS.Sign+import Data.Ext+import Data.Geometry.Point.Internal+import Data.Geometry.Vector++--------------------------------------------------------------------------------++--------------------------------------------------------------------------------++newtype StrictCCW = SCCW Sign deriving Eq++pattern CCW :: StrictCCW+pattern CCW = SCCW Negative++pattern CW  :: StrictCCW+pattern CW  = SCCW Positive+{-# COMPLETE CCW, CW #-}++instance Show StrictCCW where+  show = \case+    CCW -> "CCW"+    CW  -> "CW"+++-- | Given three points p q and r determine the orientation when going from p to r via q.+ccw       :: (Ord r, Num r, Ord i)+          => Point 2 r :+ i -> Point 2 r :+ i -> Point 2 r :+ i -> StrictCCW+ccw p q r = SCCW $ sideTest r (Vector2 p q)

@@ -0,0 +1,150 @@+module Data.Geometry.Point.Orientation.Degenerate(+    CCW(..)+  , pattern CCW, pattern CW, pattern CoLinear++  , ccw, ccw'++  , sortAround++  , ccwCmpAroundWith, cwCmpAroundWith+  , ccwCmpAround, cwCmpAround++  , insertIntoCyclicOrder+  ) where++import           Control.Lens+import qualified Data.CircularList as C+import qualified Data.CircularList.Util as CU+import           Data.Ext+import           Data.Geometry.Point.Internal+import           Data.Geometry.Vector+import qualified Data.List as L++--------------------------------------------------------------------------------++-- | Data type for expressing the orientation of three points, with+-- the option of allowing Colinearities.+newtype CCW = CCWWrap Ordering deriving Eq++pattern CCW      :: CCW+pattern CCW      = CCWWrap GT++pattern CW       :: CCW+pattern CW       = CCWWrap LT++pattern CoLinear :: CCW+pattern CoLinear = CCWWrap EQ+{-# COMPLETE CCW, CW, CoLinear #-}++instance Show CCW where+  show = \case+    CCW      -> "CCW"+    CW       -> "CW"+    CoLinear -> "CoLinear"+++-- | Given three points p q and r determine the orientation when going from p to r via q.+ccw :: (Ord r, Num r) => Point 2 r -> Point 2 r -> Point 2 r -> CCW+ccw p q r = CCWWrap $ z `compare` 0+            -- case z `compare` 0 of+            --   LT -> CW+            --   GT -> CCW+            --   EQ -> CoLinear+     where+       Vector2 ux uy = q .-. p+       Vector2 vx vy = r .-. p+       z             = ux * vy - uy * vx++-- | Given three points p q and r determine the orientation when going from p to r via q.+ccw' :: (Ord r, Num r) => Point 2 r :+ a -> Point 2 r :+ b -> Point 2 r :+ c -> CCW+ccw' p q r = ccw (p^.core) (q^.core) (r^.core)++-- | Sort the points arround the given point p in counter clockwise order with+-- respect to the rightward horizontal ray starting from p.  If two points q+-- and r are colinear with p, the closest one to p is reported first.+-- running time: O(n log n)+sortAround   :: (Ord r, Num r)+             => Point 2 r :+ q -> [Point 2 r :+ p] -> [Point 2 r :+ p]+sortAround c = L.sortBy (ccwCmpAround c <> cmpByDistanceTo c)+++-- | Given a zero vector z, a center c, and two points p and q,+-- compute the ccw ordering of p and q around c with this vector as zero+-- direction.+--+-- pre: the points p,q /= c+ccwCmpAroundWith                              :: (Ord r, Num r)+                                              => Vector 2 r+                                              -> Point 2 r :+ c+                                              -> Point 2 r :+ a -> Point 2 r :+ b+                                              -> Ordering+ccwCmpAroundWith z@(Vector2 zx zy) (c :+ _) (q :+ _) (r :+ _) =+    case (ccw c a q, ccw c a r) of+      (CCW,CCW)      -> cmp+      (CCW,CW)       -> LT+      (CCW,CoLinear) | onZero r  -> GT+                     | otherwise -> LT++      (CW, CCW)      -> GT+      (CW, CW)       -> cmp+      (CW, CoLinear) -> GT++      (CoLinear, CCW) | onZero q  -> LT+                      | otherwise -> GT++      (CoLinear, CW)      -> LT+      (CoLinear,CoLinear) -> case (onZero q, onZero r) of+                               (True, True)   -> EQ+                               (False, False) -> EQ+                               (True, False)  -> LT+                               (False, True)  -> GT+  where+    a = c .+^ z+    b = c .+^ Vector2 (-zy) zx+    -- b is on a perpendicular vector to z++    -- test if the point lies on the ray defined by z, starting in c+    onZero d = case ccw c b d of+                 CCW      -> False+                 CW       -> True+                 CoLinear -> True -- this shouldh appen only when you ask for c itself++    cmp = case ccw c q r of+            CCW      -> LT+            CW       -> GT+            CoLinear -> EQ++-- | Given a zero vector z, a center c, and two points p and q,+-- compute the cw ordering of p and q around c with this vector as zero+-- direction.+--+-- pre: the points p,q /= c+cwCmpAroundWith     :: (Ord r, Num r)+                    => Vector 2 r+                    -> Point 2 r :+ a+                    -> Point 2 r :+ b -> Point 2 r :+ c+                    -> Ordering+cwCmpAroundWith z c = flip (ccwCmpAroundWith z c)++-- | Counter clockwise ordering of the points around c. Points are ordered with+-- respect to the positive x-axis.+ccwCmpAround :: (Num r, Ord r)+             => Point 2 r :+ qc -> Point 2 r :+ p -> Point 2 r :+ q -> Ordering+ccwCmpAround = ccwCmpAroundWith (Vector2 1 0)++-- | Clockwise ordering of the points around c. Points are ordered with+-- respect to the positive x-axis.+cwCmpAround :: (Num r, Ord r)+            => Point 2 r :+ qc -> Point 2 r :+ p -> Point 2 r :+ q -> Ordering+cwCmpAround = cwCmpAroundWith (Vector2 1 0)+++-- | Given a center c, a new point p, and a list of points ps, sorted in+-- counter clockwise order around c. Insert p into the cyclic order. The focus+-- of the returned cyclic list is the new point p.+--+-- running time: O(n)+insertIntoCyclicOrder   :: (Ord r, Num r)+                        => Point 2 r :+ q -> Point 2 r :+ p+                        -> C.CList (Point 2 r :+ p) -> C.CList (Point 2 r :+ p)+insertIntoCyclicOrder c = CU.insertOrdBy (ccwCmpAround c <> cmpByDistanceTo c)

@@ -0,0 +1,69 @@+module Data.Geometry.Point.Quadrants where++import           Control.DeepSeq+import           Control.Lens+import           Data.Aeson+import           Data.Ext+import qualified Data.Foldable as F+import           Data.Geometry.Point.Class+import           Data.Geometry.Point.Internal+import           Data.Geometry.Properties+import           Data.Geometry.Vector+import qualified Data.Geometry.Vector as Vec+import           Data.Hashable+import qualified Data.List as L+import           Data.Ord (comparing)+import           Data.Proxy+import           GHC.Generics (Generic)+import           GHC.TypeLits+import           System.Random (Random(..))+import           Test.QuickCheck (Arbitrary)+import           Text.ParserCombinators.ReadP (ReadP, string,pfail)+import           Text.ParserCombinators.ReadPrec (lift)+import           Text.Read (Read(..),readListPrecDefault, readPrec_to_P,minPrec)++--------------------------------------------------------------------------------++-- | Quadrants of two dimensional points. in CCW order+data Quadrant = TopRight | TopLeft | BottomLeft | BottomRight+              deriving (Show,Read,Eq,Ord,Enum,Bounded)++-- | Quadrants around point c; quadrants are closed on their "previous"+-- boundary (i..e the boundary with the previous quadrant in the CCW order),+-- open on next boundary. The origin itself is assigned the topRight quadrant+quadrantWith                   :: (Ord r, 1 <= d, 2 <= d, Arity d)+                               => Point d r :+ q -> Point d r :+ p -> Quadrant+quadrantWith (c :+ _) (p :+ _) = case ( (c^.xCoord) `compare` (p^.xCoord)+                                      , (c^.yCoord) `compare` (p^.yCoord) ) of+                                   (EQ, EQ) -> TopRight+                                   (LT, EQ) -> TopRight+                                   (LT, LT) -> TopRight+                                   (EQ, LT) -> TopLeft+                                   (GT, LT) -> TopLeft+                                   (GT, EQ) -> BottomLeft+                                   (GT, GT) -> BottomLeft+                                   (EQ, GT) -> BottomRight+                                   (LT, GT) -> BottomRight++-- | Quadrants with respect to the origin+quadrant :: (Ord r, Num r, 1 <= d, 2 <= d, Arity d) => Point d r :+ p -> Quadrant+quadrant = quadrantWith (ext origin)++-- | Given a center point c, and a set of points, partition the points into+-- quadrants around c (based on their x and y coordinates). The quadrants are+-- reported in the order topLeft, topRight, bottomLeft, bottomRight. The points+-- are in the same order as they were in the original input lists.+-- Points with the same x-or y coordinate as p, are "rounded" to above.+partitionIntoQuadrants       :: (Ord r, 1 <= d, 2 <= d, Arity d)+                             => Point d r :+ q+                             -> [Point d r :+ p]+                             -> ( [Point d r :+ p], [Point d r :+ p]+                                , [Point d r :+ p], [Point d r :+ p]+                                )+partitionIntoQuadrants c pts = (topL, topR, bottomL, bottomR)+  where+    (below',above')   = L.partition (on yCoord) pts+    (bottomL,bottomR) = L.partition (on xCoord) below'+    (topL,topR)       = L.partition (on xCoord) above'++    on l q       = q^.core.l < c^.core.l

@@ -5,7 +5,9 @@  import           Control.Lens import           Data.Aeson+import           Data.Bifoldable import           Data.Bifunctor+import           Data.Bitraversable import           Data.Ext import qualified Data.Foldable as F import           Data.Geometry.Box@@ -17,10 +19,18 @@ import           Data.LSeq (LSeq, pattern (:<|)) import qualified Data.LSeq as LSeq import qualified Data.List.NonEmpty as NE-import           GHC.Generics(Generic)+import           GHC.Generics (Generic) import           GHC.TypeLits  --------------------------------------------------------------------------------++--------------------------------------------------------------------------------+-- $setup+-- >>> :{+-- let myPolyLine = fromPointsUnsafe $ map ext [origin, Point2 10.0 10.0, Point2 10.0 20.0]+-- :}++-------------------------------------------------------------------------------- -- * d-dimensional Polygonal Lines (PolyLines)  -- | A Poly line in R^d has at least 2 vertices@@ -50,25 +60,33 @@   pmap f = over points (fmap (first f))  instance Arity d => Bifunctor (PolyLine d) where-  bimap f g (PolyLine pts) = PolyLine $ fmap (bimap (fmap g) f) pts+  bimap = bimapDefault+instance Arity d => Bifoldable (PolyLine d) where+  bifoldMap = bifoldMapDefault+instance Arity d => Bitraversable (PolyLine d) where+  bitraverse f g (PolyLine pts) = PolyLine <$> traverse (bitraverse (traverse g) f) pts  instance (ToJSON p, ToJSON r, Arity d) => ToJSON (PolyLine d p r) where     toEncoding = genericToEncoding defaultOptions instance (FromJSON p, FromJSON r, Arity d, KnownNat d) => FromJSON (PolyLine d p r) +-- | Builds a Polyline from a list of points, if there are sufficiently many points+fromPoints :: [Point d r :+ p] -> Maybe (PolyLine d p r)+fromPoints = fmap PolyLine . LSeq.eval (C @ 2) . LSeq.fromList+ -- | pre: The input list contains at least two points-fromPoints :: [Point d r :+ p] -> PolyLine d p r-fromPoints = PolyLine . LSeq.forceLSeq (C  :: C 2) . LSeq.fromList+fromPointsUnsafe :: [Point d r :+ p] -> PolyLine d p r+fromPointsUnsafe = PolyLine . LSeq.forceLSeq (C @ 2) . LSeq.fromList  -- | pre: The input list contains at least two points. All extra vields are -- initialized with mempty.-fromPoints' :: (Monoid p) => [Point d r] -> PolyLine d p r-fromPoints' = fromPoints . map (\p -> p :+ mempty)+fromPointsUnsafe' :: (Monoid p) => [Point d r] -> PolyLine d p r+fromPointsUnsafe' = fromPointsUnsafe . map (\p -> p :+ mempty)   -- | We consider the line-segment as closed. fromLineSegment                     :: LineSegment d p r -> PolyLine d p r-fromLineSegment ~(LineSegment' p q) = fromPoints [p,q]+fromLineSegment ~(LineSegment' p q) = fromPointsUnsafe [p,q]  -- | Convert to a closed line segment by taking the first two points. asLineSegment                            :: PolyLine d p r -> LineSegment d p r@@ -80,3 +98,23 @@ asLineSegment' (PolyLine pts) = case F.toList pts of                                   [p,q] -> Just $ ClosedLineSegment p q                                   _     -> Nothing++-- | Computes the edges, as linesegments, of an LSeq+edgeSegments    :: Arity d => PolyLine d p r -> LSeq 1 (LineSegment d p r)+edgeSegments pl = let vs = pl^.points+                  in LSeq.zipWith ClosedLineSegment (LSeq.init vs) (LSeq.tail vs)+++-- | Linearly interpolate the polyline with a value in the range+-- \([0,n-1]\), where \(n\) is the number of vertices of the polyline.+--+-- running time: \(O(\log n)\)+--+-- >>> interpolatePoly 0.5 myPolyLine+-- Point2 [5.0,5.0]+-- >>> interpolatePoly 1.5 myPolyLine+-- Point2 [10.0,15.0]+interpolatePoly      :: (RealFrac r, Arity d) => r -> PolyLine d p r -> Point d r+interpolatePoly t pl = let i = floor t in case edgeSegments pl^?ix i of+                         Nothing -> pl^.points.to LSeq.last.core+                         Just e  -> interpolate (t-fromIntegral i) e

@@ -73,12 +73,15 @@ import           Data.Util import           Data.Vinyl.CoRec (asA) +-- import Data.RealNumber.Rational+ --------------------------------------------------------------------------------  {- $setup+>>> import Data.RealNumber.Rational >>> :{ -- import qualified Data.CircularSeq as C-let simplePoly :: SimplePolygon () Rational+let simplePoly :: SimplePolygon () (RealNumber 10)     simplePoly = SimplePolygon . C.fromList . map ext $ [ Point2 0 0                                                         , Point2 10 0                                                         , Point2 10 10@@ -268,17 +271,17 @@ -- -- -- >>> mapM_ print . polygonVertices $ withIncidentEdges simplePoly--- Point2 [0 % 1,0 % 1] :+ SP LineSegment (Closed (Point2 [1 % 1,11 % 1] :+ ())) (Closed (Point2 [0 % 1,0 % 1] :+ ())) LineSegment (Closed (Point2 [0 % 1,0 % 1] :+ ())) (Closed (Point2 [10 % 1,0 % 1] :+ ()))--- Point2 [10 % 1,0 % 1] :+ SP LineSegment (Closed (Point2 [0 % 1,0 % 1] :+ ())) (Closed (Point2 [10 % 1,0 % 1] :+ ())) LineSegment (Closed (Point2 [10 % 1,0 % 1] :+ ())) (Closed (Point2 [10 % 1,10 % 1] :+ ()))--- Point2 [10 % 1,10 % 1] :+ SP LineSegment (Closed (Point2 [10 % 1,0 % 1] :+ ())) (Closed (Point2 [10 % 1,10 % 1] :+ ())) LineSegment (Closed (Point2 [10 % 1,10 % 1] :+ ())) (Closed (Point2 [5 % 1,15 % 1] :+ ()))--- Point2 [5 % 1,15 % 1] :+ SP LineSegment (Closed (Point2 [10 % 1,10 % 1] :+ ())) (Closed (Point2 [5 % 1,15 % 1] :+ ())) LineSegment (Closed (Point2 [5 % 1,15 % 1] :+ ())) (Closed (Point2 [1 % 1,11 % 1] :+ ()))--- Point2 [1 % 1,11 % 1] :+ SP LineSegment (Closed (Point2 [5 % 1,15 % 1] :+ ())) (Closed (Point2 [1 % 1,11 % 1] :+ ())) LineSegment (Closed (Point2 [1 % 1,11 % 1] :+ ())) (Closed (Point2 [0 % 1,0 % 1] :+ ()))+-- Point2 [0,0] :+ V2 LineSegment (Closed (Point2 [1,11] :+ ())) (Closed (Point2 [0,0] :+ ())) LineSegment (Closed (Point2 [0,0] :+ ())) (Closed (Point2 [10,0] :+ ()))+-- Point2 [10,0] :+ V2 LineSegment (Closed (Point2 [0,0] :+ ())) (Closed (Point2 [10,0] :+ ())) LineSegment (Closed (Point2 [10,0] :+ ())) (Closed (Point2 [10,10] :+ ()))+-- Point2 [10,10] :+ V2 LineSegment (Closed (Point2 [10,0] :+ ())) (Closed (Point2 [10,10] :+ ())) LineSegment (Closed (Point2 [10,10] :+ ())) (Closed (Point2 [5,15] :+ ()))+-- Point2 [5,15] :+ V2 LineSegment (Closed (Point2 [10,10] :+ ())) (Closed (Point2 [5,15] :+ ())) LineSegment (Closed (Point2 [5,15] :+ ())) (Closed (Point2 [1,11] :+ ()))+-- Point2 [1,11] :+ V2 LineSegment (Closed (Point2 [5,15] :+ ())) (Closed (Point2 [1,11] :+ ())) LineSegment (Closed (Point2 [1,11] :+ ())) (Closed (Point2 [0,0] :+ ())) withIncidentEdges                    :: Polygon t p r                                      -> Polygon t (Two (LineSegment 2 p r)) r withIncidentEdges (SimplePolygon vs) =       SimplePolygon $ C.zip3LWith f (C.rotateL vs) vs (C.rotateR vs)   where-    f p c n = c&extra .~ SP (ClosedLineSegment p c) (ClosedLineSegment c n)+    f p c n = c&extra .~ Two (ClosedLineSegment p c) (ClosedLineSegment c n) withIncidentEdges (MultiPolygon vs hs) = MultiPolygon vs' hs'   where     (SimplePolygon vs') = withIncidentEdges $ SimplePolygon vs@@ -565,7 +568,7 @@ -- will be numbered last, in the same order. -- -- >>> numberVertices simplePoly--- SimplePolygon (CSeq [Point2 [0 % 1,0 % 1] :+ SP 0 (),Point2 [10 % 1,0 % 1] :+ SP 1 (),Point2 [10 % 1,10 % 1] :+ SP 2 (),Point2 [5 % 1,15 % 1] :+ SP 3 (),Point2 [1 % 1,11 % 1] :+ SP 4 ()])+-- SimplePolygon (CSeq [Point2 [0,0] :+ SP 0 (),Point2 [10,0] :+ SP 1 (),Point2 [10,10] :+ SP 2 (),Point2 [5,15] :+ SP 3 (),Point2 [1,11] :+ SP 4 ()]) numberVertices :: Polygon t p r -> Polygon t (SP Int p) r numberVertices = snd . bimapAccumL (\a p -> (a+1,SP a p)) (\a r -> (a,r)) 0   -- TODO: Make sure that this does not have the same issues as foldl vs foldl'

@@ -26,6 +26,8 @@ import           Data.Geometry.Point import           Data.List.NonEmpty (NonEmpty(..)) import qualified Data.List.NonEmpty as NonEmpty+import           Data.Measured.Class ()+import           Data.Measured.Size import           Data.Ord (comparing, Down(..)) import           Data.Range import qualified Data.Set as Set

@@ -0,0 +1,203 @@+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE TypeApplications #-}+module Data.Geometry.QuadTree-- ( module Data.Geometry.QuadTree.Cell+                             -- , module Data.Geometry.QuadTree.Quadrants+                             -- , module Data.Geometry.QuadTree.Split+                             -- , QuadTree(..)+                             -- , leaves+                             -- , withCells+                             -- )+                             where+++import           Control.Lens (makeLenses, (^.), (.~), (&), (^?!), ix, view)+import           Data.Ext+import qualified Data.Foldable as F+import           Data.Geometry.Box+import           Data.Geometry.Point+import           Data.Geometry.QuadTree.Cell+import           Data.Geometry.QuadTree.Quadrants+import           Data.Geometry.QuadTree.Split+import           Data.Geometry.QuadTree.Tree (Tree(..))+import qualified Data.Geometry.QuadTree.Tree as Tree+import           Data.Geometry.Vector+import           Data.Intersection+import           Data.List.NonEmpty (NonEmpty(..))+import           Data.Tree.Util (TreeNode(..), levels)+import           GHC.Generics (Generic)+--------------------------------------------------------------------------------++-- | QuadTree on the starting cell+data QuadTree v p r = QuadTree { _startingCell  :: !(Cell r)+                               , _tree          :: !(Tree v p)+                               }+                    deriving (Show,Eq,Generic,Functor,Foldable,Traversable)+makeLenses ''QuadTree++--------------------------------------------------------------------------------+-- * Functions operating on the QuadTree (in terms of the 'Tree' type)++withCells    :: (Fractional r, Ord r) => QuadTree v p r -> QuadTree (v :+ Cell r) (p :+ Cell r) r+withCells qt = qt&tree .~ withCellsTree qt++withCellsTree                :: (Fractional r, Ord r)+                             => QuadTree v p r -> Tree (v :+ Cell r) (p :+ Cell r)+withCellsTree (QuadTree c t) = Tree.withCells c t++leaves :: QuadTree v p r -> NonEmpty p+leaves = Tree.leaves . view tree++perLevel :: QuadTree v p r -> NonEmpty (NonEmpty (TreeNode v p))+perLevel = levels . Tree.toRoseTree . view tree+++--------------------------------------------------------------------------------++-- | Given a starting cell, a Tree builder, and some input required by+-- the builder, constructs a quadTree.+buildOn            :: Cell r -> (Cell r -> i -> Tree v p) -> i -> QuadTree v p r+buildOn c0 builder = QuadTree c0 . builder c0++-- | The Equivalent of Tree.build for constructing a QuadTree+build     :: (Fractional r, Ord r) => (Cell r -> i -> Split i v p) -> Cell r -> i -> QuadTree v p r+build f c = buildOn c (Tree.build f)++-- | Build a QuadtTree from a set of points.+--+-- pre: the points lie inside the initial given cell.+--+-- running time: \(O(nh)\), where \(n\) is the number of points and+-- \(h\) is the height of the resulting quadTree.+fromPointsBox   :: (Fractional r, Ord r)+                 => Cell r -> [Point 2 r :+ p] -> QuadTree () (Maybe (Point 2 r :+ p)) r+fromPointsBox c = buildOn c Tree.fromPoints++fromPoints     :: (RealFrac r, Ord r)+               => NonEmpty (Point 2 r :+ p) -> QuadTree () (Maybe (Point 2 r :+ p)) r+fromPoints pts = buildOn c Tree.fromPoints (F.toList pts)+  where+    c = fitsRectangle $ boundingBoxList (view core <$> pts)++-- | Locates the cell containing the given point, if it exists.+--+-- running time: \(O(h)\), where \(h\) is the height of the quadTree+findLeaf                                       :: (Fractional r, Ord r)+                                               => Point 2 r -> QuadTree v p r -> Maybe (p :+ Cell r)+findLeaf q (QuadTree c0 t) | q `intersects` c0  = Just $ findLeaf' c0 t+                           | otherwise          = Nothing+  where+    -- |+    -- pre: p intersects c+    findLeaf' c = \case+      Leaf p    -> p :+ c+      Node _ qs -> let quad = quadrantOf q c+                   in findLeaf' ((splitCell c)^?!ix quad) (qs^?!ix quad)++--------------------------------------------------------------------------------+++fromZeros :: (Fractional r, Ord r, Num a, Eq a, v ~ Quadrants Sign)+          => Cell r -> (Point 2 r -> a) -> QuadTree v (Either v Sign) r+fromZeros = fromZerosWith (limitWidthTo (-1))+++fromZerosWith            ::  (Fractional r, Ord r, Eq a, Num a)+                         => Limiter r (Corners Sign) (Corners Sign) Sign+                         -> Cell r+                         -> (Point 2 r -> a)+                         -> QuadTree (Quadrants Sign) (Signs Sign) r+fromZerosWith limit c0 f = fromZerosWith' limit c0 (fromSignum f)+++type Signs sign = Either (Corners sign) sign+++fromZerosWith'           :: (Eq sign, Fractional r, Ord r)+                         => Limiter r (Corners sign) (Corners sign) sign+                         -> Cell r+                         -> (Point 2 r -> sign)+                         -> QuadTree (Quadrants sign) (Signs sign) r+fromZerosWith' limit c0 f = build (limit $ shouldSplitZeros f) c0 (f <$> cellCorners c0)++++-- type Sign = Ordering++-- pattern Negative :: Sign+-- pattern Negative = LT+-- pattern Zero :: Sign+-- pattern Zero     = EQ+-- pattern Positive :: Sign+-- pattern Positive = GT+-- {-# COMPLETE Negative, Zero, Positive #-}++-- fromOrdering :: Ordering -> Sign+-- fromOrdering = id+++data Sign = Negative | Zero | Positive deriving (Show,Eq,Ord)++++-- | Interpret an ordering result as a Sign+fromOrdering :: Ordering -> Sign+fromOrdering = \case+    LT -> Negative+    EQ -> Zero+    GT -> Positive++fromSignum   :: (Num a, Eq a) => (b -> a) -> b -> Sign+fromSignum f = \x -> case signum (f x) of+                       -1 -> Negative+                       0  -> Zero+                       1  -> Positive+                       _  -> error "absurd: fromSignum"++-- | Splitter that determines if we should split a cell based on the+-- sign of the corners.+shouldSplitZeros :: forall r sign. (Fractional r, Eq sign)+                 => (Point 2 r -> sign) -- ^ The function we are evaluating+                 -> Splitter r+                             (Quadrants sign) -- the input are the signs of the corners+                             (Quadrants sign) -- at internal nodes we store signs of corners+                             sign+shouldSplitZeros f (Cell w' p) qs@(Quadrants nw ne se sw) | all sameSign qs = No ne+                                                          | otherwise       = Yes qs qs'+  where+    m = fAt rr rr+    n = fAt rr ww+    e = fAt ww rr+    s = fAt rr 0+    w = fAt 0  rr++    sameSign = (== ne)++    -- signs at the new corners+    qs' = Quadrants (Quadrants nw n m w)+                    (Quadrants n ne e m)+                    (Quadrants m e se s)+                    (Quadrants w m s sw)++    r     = w' - 1+    rr    = pow r+    ww    = pow w'++    fAt x y = f $ p .+^ Vector2 x y+++isZeroCell   :: (Eq sign) => sign -- ^ the zero value+             -> Either v sign -> Bool+isZeroCell z = \case+    Left _  -> True -- if we kept splitting then we must have a sign transition+    Right s -> s == z++--------------------------------------------------------------------------------++++-- | Constructs an empty/complete tree from the starting width+completeTree    :: (Fractional r, Ord r) => Cell r -> QuadTree () () r+completeTree c0 =+    build (\_ w -> if w == 0 then No () else Yes () (pure $ w - 1)) c0 (c0^.cellWidthIndex)++--------------------------------------------------------------------------------

@@ -0,0 +1,141 @@+{-# LANGUAGE TemplateHaskell #-}+module Data.Geometry.QuadTree.Cell where++import Control.Lens (makeLenses, (^.),(&),(%~),ix, to)+import Data.Ext+import Data.Geometry.Box+import Data.Geometry.Directions+import Data.Geometry.LineSegment+import Data.Geometry.Point+import Data.Geometry.Properties+import Data.Geometry.QuadTree.Quadrants+import Data.Geometry.Vector++--------------------------------------------------------------------------------++-- | side lengths will be 2^i for some integer i+type WidthIndex = Int++-- | A Cell corresponding to a node in the QuadTree+data Cell r = Cell { _cellWidthIndex :: {-# UNPACK #-} !WidthIndex+                   , _lowerLeft      ::                !(Point 2 r)+                   } deriving (Show,Eq,Functor,Foldable,Traversable)+makeLenses ''Cell++-- | Computes a cell that contains the given rectangle+fitsRectangle   :: (RealFrac r, Ord r) => Rectangle p r -> Cell r+fitsRectangle b = Cell w ((b^.to minPoint.core) .-^ Vector2 1 1)+  where+    w = lg' . ceiling . (1+) . maximum . size $ b++    -- "approximate log" that over approximates by a factor of at most two.+    lg'   :: Integer -> WidthIndex+    lg' n = go 1+      where+        go i | floor (pow i) <= n = go (i+1) -- note that the floor does not really do anything+                                             -- since i is integral and >= 1.+             | otherwise  = i++type instance Dimension (Cell r) = 2+type instance NumType   (Cell r) = r++type instance IntersectionOf (Point 2 r) (Cell r) = '[ NoIntersection, Point 2 r]++instance (Ord r, Fractional r) => (Point 2 r) `IsIntersectableWith` (Cell r) where+  nonEmptyIntersection = defaultNonEmptyIntersection+  p `intersect` c = p `intersect` toBox c++pow   :: Fractional r => WidthIndex -> r+pow i = case i `compare` 0 of+          LT -> 1 / (2 ^ (-1*i))+          EQ -> 1+          GT -> 2 ^ i++cellWidth            :: Fractional r => Cell r -> r+cellWidth (Cell w _) = pow w++toBox            :: Fractional r => Cell r -> Box 2 () r+toBox (Cell w p) = box (ext $ p) (ext $ p .+^ Vector2 (pow w) (pow w))++inCell            :: (Fractional r, Ord r) => Point 2 r :+ p -> Cell r -> Bool+inCell (p :+ _) c = p `inBox` toBox c++cellCorners :: Fractional r => Cell r -> Quadrants (Point 2 r)+cellCorners = fmap (^.core) . corners . toBox++-- | Sides are open+cellSides :: Fractional r => Cell r -> Sides (LineSegment 2 () r)+cellSides = fmap (\(ClosedLineSegment p q) -> OpenLineSegment p q) . sides . toBox++splitCell            :: (Num r, Fractional r) => Cell r -> Quadrants (Cell r)+splitCell (Cell w p) = Quadrants (Cell r $ f 0 rr)+                                 (Cell r $ f rr rr)+                                 (Cell r $ f rr 0)+                                 (Cell r p)+  where+    r     = w - 1+    rr    = pow r+    f x y = p .+^ Vector2 x y+++midPoint            :: Fractional r => Cell r -> Point 2 r+midPoint (Cell w p) = let rr = pow (w - 1) in p .+^ Vector2 rr rr+++--------------------------------------------------------------------------------++-- | Partitions the points into quadrants. See 'quadrantOf' for the+-- precise rules.+partitionPoints   :: (Fractional r, Ord r)+                  => Cell r -> [Point 2 r :+ p] -> Quadrants [Point 2 r :+ p]+partitionPoints c = foldMap (\p -> let q = quadrantOf (p^.core) c in mempty&ix q %~ (p:))++-- | Computes the quadrant of the cell corresponding to the current+-- point. Note that we decide the quadrant solely based on the+-- midpoint. If the query point lies outside the cell, it is still+-- assigned a quadrant.+--+-- - The northEast quadrants includes its bottom and left side+-- - The southEast quadrant  includes its            left side+-- - The northWest quadrant  includes its bottom          side+-- - The southWest quadrants does not include any of its sides.+--+--+-- >>> quadrantOf (Point2 9 9) (Cell 4 origin)+-- NorthEast+-- >>> quadrantOf (Point2 8 9) (Cell 4 origin)+-- NorthEast+-- >>> quadrantOf (Point2 8 8) (Cell 4 origin)+-- NorthEast+-- >>> quadrantOf (Point2 8 7) (Cell 4 origin)+-- SouthEast+-- >>> quadrantOf (Point2 4 7) (Cell 4 origin)+-- SouthWest+-- >>> quadrantOf (Point2 4 10) (Cell 4 origin)+-- NorthWest+-- >>> quadrantOf (Point2 4 40) (Cell 4 origin)+-- NorthEast+-- >>> quadrantOf (Point2 4 40) (Cell 4 origin)+-- NorthWest+quadrantOf     :: forall r. (Fractional r, Ord r)+               => Point 2 r -> Cell r -> InterCardinalDirection+quadrantOf q c = let m = midPoint c+                 in case (q^.xCoord < m^.xCoord, q^.yCoord < m^.yCoord) of+                      (False,False) -> NorthEast+                      (False,True)  -> SouthEast+                      (True,False)  -> NorthWest+                      (True,True)   -> SouthWest++++-- | Given two cells c and me, compute on which side of `me` the cell+-- `c` is.+--+-- pre: c and me are non-overlapping+relationTo        :: (Fractional r, Ord r)+                  => (p :+ Cell r) -> Cell r -> Sides (Maybe (p :+ Cell r))+c `relationTo` me = f <$> Sides b l t r <*> cellSides me+  where+    Sides t r b l = cellSides (c^.extra)+    f e e' | e `intersects` e' = Just c+           | otherwise         = Nothing

@@ -0,0 +1,16 @@+module Data.Geometry.QuadTree.Quadrants( pattern Quadrants+                                       , Quadrants+                                       , module Data.Geometry.Box.Corners+                                       ) where++import           Data.Geometry.Box.Corners++--------------------------------------------------------------------------------++type Quadrants = Corners++pattern Quadrants         :: a -> a -> a -> a -> Corners a+pattern Quadrants a b c d = Corners a b c d+{-# COMPLETE Quadrants #-}++--------------------------------------------------------------------------------

@@ -0,0 +1,32 @@+{-# LANGUAGE TemplateHaskell #-}+module Data.Geometry.QuadTree.Split where++import Control.Lens (makePrisms,(^.))+import Data.Geometry.QuadTree.Cell+import Data.Geometry.QuadTree.Quadrants++--------------------------------------------------------------------------------++-- | Data Type to Decide if we should continue splitting the current cell+data Split i v p = No !p | Yes !v (Quadrants i) deriving (Show,Eq,Ord)+makePrisms ''Split++-- | A splitter is a function that determines weather or not we should the given cell+-- corresponding to the given input (i).+type Splitter r i v p = Cell r -> i -> Split i v p++-- | Transformer that limits the depth of a splitter+type Limiter r i v p = Splitter r i v p+                    -> Splitter r i v (Either i p)++-- | Split only when the Cell-width is at least wMin+limitWidthTo        :: WidthIndex -- ^ smallest allowed width of a cell (i.e. width of a leaf)+                    -> Limiter r i v p+limitWidthTo wMin f = \c pts -> case f c pts of+                                  No p                                -> No (Right p)+                                  Yes v qs | wMin < c^.cellWidthIndex -> Yes v qs+                                           | otherwise                -> No (Left pts)+  -- note that it is important that we still evaluate the function so+  -- that we can distinguish at the last level i.e. between a regular+  -- " we are done splitting (No (Right p))" and a "we are no longer+  -- allowed to split further (No (Left p))"

@@ -0,0 +1,116 @@+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE TypeApplications #-}+module Data.Geometry.QuadTree.Tree where+++import           Control.Lens (makePrisms)+import           Data.Bifoldable+import           Data.Bifunctor+import           Data.Bitraversable+import           Data.Ext+import qualified Data.Foldable as F+import           Data.Functor.Apply+import           Data.Geometry.Point+import           Data.Geometry.QuadTree.Cell+import           Data.Geometry.QuadTree.Quadrants+import           Data.Geometry.QuadTree.Split+import           Data.List.NonEmpty (NonEmpty(..))+import qualified Data.List.NonEmpty as NonEmpty+import           Data.Semigroup.Foldable.Class+import           Data.Semigroup.Traversable.Class+import qualified Data.Tree as RoseTree+import           Data.Tree.Util (TreeNode(..))++--------------------------------------------------------------------------------++-- | Our cells use Rational numbers as their numeric type+-- type CellR = Cell (RealNumber 10)++-- | The Actual Tree type representing a quadTree+data Tree v p = Leaf !p+              | Node !v (Quadrants (Tree v p)) -- quadrants are stored lazily on purpose+              deriving (Show,Eq)+makePrisms ''Tree++instance Bifunctor Tree where+  bimap = bimapDefault++instance Bifoldable Tree where+  bifoldMap = bifoldMapDefault++instance Bitraversable Tree where+  bitraverse f g = \case+    Leaf p    -> Leaf <$> g p+    Node v qs -> Node <$> f v <*> traverse (bitraverse f g) qs++instance Bifoldable1 Tree+instance Bitraversable1 Tree where+  bitraverse1 f g = \case+    Leaf p    -> Leaf <$> g p+    Node v qs -> Node <$> f v <.> traverse1 (bitraverse1 f g) qs++-- | Fold on the Tree type.+foldTree     :: (p -> b) -> (v -> Quadrants b -> b) -> Tree v p -> b+foldTree f g = go+  where+    go = \case+      Leaf p    -> f p+      Node v qs -> g v (go <$> qs)++-- | Produce a list of all leaves of a quad tree+leaves :: Tree v p -> NonEmpty p+leaves = NonEmpty.fromList . bifoldMap (const []) (:[])++-- | Converts into a RoseTree+toRoseTree :: Tree v p -> RoseTree.Tree (TreeNode v p)+toRoseTree = foldTree (\p    -> RoseTree.Node (LeafNode p)     [])+                      (\v qs -> RoseTree.Node (InternalNode v) (F.toList qs))++-- | Computes the height of the quadtree+height :: Tree v p -> Integer+height = foldTree (const 1) (\_ -> (1 +) . maximum)+++--------------------------------------------------------------------------------++--------------------------------------------------------------------------------+-- * Functions operating on the QuadTree (in temrs of the 'Tree' type)++-- | Builds a QuadTree+build             :: Fractional r+                  => Splitter r pts v p -> Cell r -> pts -> Tree v p+build shouldSplit = build'+  where+    build' cc pts = case shouldSplit cc pts of+                      No p     -> Leaf p+                      Yes v qs -> Node v $ build' <$> splitCell cc <*> qs++-- | Annotate the tree with its corresponing cells+withCells :: Fractional r => Cell r -> Tree v p -> Tree (v :+ Cell r) (p :+ Cell r)+withCells c0 = \case+  Leaf p    -> Leaf (p :+ c0)+  Node v qs -> Node (v :+ c0) (withCells <$> splitCell c0 <*> qs)+++--------------------------------------------------------------------------------+++-- | Build a QuadtTree from a set of points.+--+-- pre: the points lie inside the initial given cell.+--+-- running time: \(O(nh)\), where \(n\) is the number of points and+-- \(h\) is the height of the resulting quadTree.+fromPoints :: (Fractional r, Ord r)+           => Cell r -> [Point 2 r :+ p]+           -> Tree () (Maybe (Point 2 r :+ p))+fromPoints = build fromPointsF++-- | The function that can be used to build a quadTree 'fromPoints'+fromPointsF   :: (Fractional r, Ord r)+              => Splitter r [Point 2 r :+ p] () (Maybe (Point 2 r :+ p))+fromPointsF c = \case+      []   -> No Nothing+      [p]  -> No (Just p)+      pts  -> Yes () $ partitionPoints c pts+        -- (\cell -> filter (`inCell` cell) pts) <$> splitCell c

@@ -2,7 +2,6 @@ module Data.Geometry.RangeTree where  import           Control.Lens hiding (element)-import           Data.BinaryTree (Measured(..)) import           Data.Ext import qualified Data.Foldable as F import           Data.Geometry.Point@@ -11,6 +10,7 @@ import           Data.Geometry.Vector import           Data.List.NonEmpty (NonEmpty(..)) import qualified Data.List.NonEmpty as NonEmpty+import           Data.Measured.Class import           Data.Proxy import           Data.Range import           Data.Semigroup.Foldable

@@ -9,6 +9,8 @@ import           Data.List.NonEmpty (NonEmpty(..)) import qualified Data.List.NonEmpty as NonEmpty import           Data.Range+import           Data.Measured.Class+import           Data.Measured.Size import           Data.Semigroup import           Data.Semigroup.Foldable import qualified Data.Set as Set

@@ -1,6 +1,6 @@ module Data.Geometry.RangeTree.Measure where -import           Data.BinaryTree(Measured(..))+import Data.Measured.Class import Data.Functor.Product(Product(..)) import Data.Functor.Classes 

@@ -34,6 +34,8 @@ import qualified Data.List as List import           Data.List.NonEmpty (NonEmpty) import qualified Data.List.NonEmpty as NonEmpty+import           Data.Measured.Class+import           Data.Measured.Size import           GHC.Generics (Generic)  --------------------------------------------------------------------------------@@ -73,7 +75,7 @@ --                         AtomicRange -> OpenRange   MinInfinity MaxInfinity  -data BuildLeaf a = LeafSingleton a | LeafRange a a deriving (Show,Eq)+data BuildLeaf a = LeafSingleton !a | LeafRange !a !a deriving (Show,Eq)  -- | Given a sorted list of endpoints, without duplicates, construct a segment tree --@@ -118,7 +120,7 @@                    -> NonEmpty (Interval p r) -> SegmentTree v r fromIntervals f is = foldr (insert . f) (createTree pts mempty) is   where-    endPoints (toRange -> Range' a b) = [a,b]+    endPoints (asRange -> Range' a b) = [a,b]     pts = nub' . NonEmpty.sort . NonEmpty.fromList . concatMap endPoints $ is     nub' = fmap NonEmpty.head . NonEmpty.group1 @@ -186,7 +188,7 @@                          => i -> SegmentTree v r -> SegmentTree v r insert i (SegmentTree t) = SegmentTree $ insertRoot t   where-    ri@(Range a b) = toRange i+    ri@(Range a b) = asRange i     insertRoot t' = maybe t' (flip insert' t') $ getRange t'      insert' inR         lf@(Leaf nd@(LeafData rr _))@@ -209,7 +211,7 @@           => i -> SegmentTree v r -> SegmentTree v r delete i (SegmentTree t) = SegmentTree $ delete' t   where-    (Range _ b) = toRange i+    (Range _ b) = asRange i      delete' (Leaf ld) = Leaf $ ld&leafAssoc %~ deleteAssoc i     delete' (Node l nd@(_splitPoint -> m) r)@@ -254,7 +256,7 @@   instance IntervalLike a => IntervalLike (I a) where-  toRange = toRange . _unI+  asRange = asRange . _unI   fromIntervals' :: (Eq p, Ord r)

@@ -136,7 +136,7 @@       :& RNil     where       s'  = (fixEndPoints sm)^.subRange-      s'' = bimap (^.extra) id+      s'' = asProperInterval . first (^.extra)           $ s'&start.core .~ toOffset' (s'^.start.extra.core) l               &end.core   .~ toOffset' (s'^.end.extra.core)   l 

@@ -1,69 +1,37 @@+{-# LANGUAGE Unsafe #-} {-# LANGUAGE UndecidableInstances #-} module Data.Geometry.Transformation where -import           Control.Lens (lens,Lens',set)-import           Unsafe.Coerce(unsafeCoerce)+import           Control.Lens (iso,set,Iso,imap)+import           Data.Geometry.Matrix+import           Data.Geometry.Matrix.Internal (mkRow) import           Data.Geometry.Point import           Data.Geometry.Properties import           Data.Geometry.Vector import qualified Data.Geometry.Vector as V import           Data.Proxy-import qualified Data.Vector.Fixed as FV import           GHC.TypeLits-import           Linear.Matrix ((!*),(!*!))-import qualified Linear.Matrix as Lin ------------------------------------------------------------------------------------ * Matrices --- | a matrix of n rows, each of m columns, storing values of type r-newtype Matrix n m r = Matrix (Vector n (Vector m r))--deriving instance (Show r, Arity n, Arity m) => Show (Matrix n m r)-deriving instance (Eq r, Arity n, Arity m)   => Eq (Matrix n m r)-deriving instance (Ord r, Arity n, Arity m)  => Ord (Matrix n m r)-deriving instance (Arity n, Arity m)         => Functor (Matrix n m)--multM :: (Arity r, Arity c, Arity c', Num a) => Matrix r c a -> Matrix c c' a -> Matrix r c' a-(Matrix a) `multM` (Matrix b) = Matrix $ a !*! b--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  -- | A type representing a Transformation for d dimensional objects newtype Transformation d r = Transformation { _transformationMatrix :: Matrix (d + 1) (d + 1) r } -transformationMatrix :: Lens' (Transformation d r) (Matrix (d + 1) (d + 1) r)-transformationMatrix = lens _transformationMatrix (const Transformation)+transformationMatrix :: Iso (Transformation d r)       (Transformation d       s)+                            (Matrix (d + 1) (d + 1) r) (Matrix (d + 1) (d + 1) s)+transformationMatrix = iso _transformationMatrix Transformation  deriving instance (Show r, Arity (d + 1)) => Show (Transformation d r) deriving instance (Eq r, Arity (d + 1))   => Eq (Transformation d r) deriving instance (Ord r, Arity (d + 1))  => Ord (Transformation d r) deriving instance Arity (d + 1)           => Functor (Transformation d)+deriving instance Arity (d + 1)           => Foldable (Transformation d)+deriving instance Arity (d + 1)           => Traversable (Transformation d)  type instance NumType (Transformation d r) = r - -- | Compose transformations (right to left) (|.|) :: (Num r, Arity (d + 1)) => Transformation d r -> Transformation d r -> Transformation d r (Transformation f) |.| (Transformation g) = Transformation $ f `multM` g@@ -112,16 +80,17 @@  translation   :: (Num r, Arity d, Arity (d + 1))               => Vector d r -> Transformation d r-translation v = Transformation . Matrix $ V.imap transRow (snoc v 1)+translation v = Transformation . Matrix $ imap transRow (snoc v 1)   scaling   :: (Num r, Arity d, Arity (d + 1))           => Vector d r -> Transformation d r-scaling v = Transformation . Matrix $ V.imap mkRow (snoc v 1)+scaling v = Transformation . Matrix $ imap mkRow (snoc v 1)  uniformScaling :: (Num r, Arity d, Arity (d + 1)) => r -> Transformation d r uniformScaling = scaling . pure + -------------------------------------------------------------------------------- -- * Functions that execute transformations @@ -142,14 +111,6 @@ scaleUniformlyBy = transformBy  . uniformScaling  ------------------------------------------------------------------------------------ * Helper functions to easily create matrices---- | Creates a row with zeroes everywhere, except at position i, where the--- value is the supplied value.-mkRow     :: forall d r. (Arity d, Num r) => Int -> r -> Vector d r-mkRow i x = set (FV.element i) x zero- -- | Row in a translation matrix -- transRow     :: forall n r. ( Arity n, Arity (n- 1), ((n - 1) + 1) ~ n --                             , Num r) => Int -> r -> Vector n r@@ -171,3 +132,13 @@                                                              (snoc v        0)                                                              (snoc w        0)                                                              (Vector4 0 0 0 1)++--------------------------------------------------------------------------------+-- * 2D Transformations++-- | Skew transformation that keeps the y-coordinates fixed and shifts+-- the x coordinates.+skewX        :: Num r => r -> Transformation 2 r+skewX lambda = Transformation . Matrix $ Vector3 (Vector3 1 lambda 0)+                                                 (Vector3 0 1      0)+                                                 (Vector3 0 0      1)

@@ -3,8 +3,11 @@ {-# LANGUAGE UndecidableInstances #-} module Data.Geometry.Triangle where +import           Control.DeepSeq import           Control.Lens+import           Data.Bifoldable import           Data.Bifunctor+import           Data.Bitraversable import           Data.Either (partitionEithers) import           Data.Ext import           Data.Geometry.Ball (Disk, disk)@@ -19,28 +22,48 @@ import qualified Data.Geometry.Vector as V import qualified Data.List as List import           Data.Maybe (mapMaybe)+import           Data.Util import           Data.Vinyl import           Data.Vinyl.CoRec+import           GHC.Generics (Generic) import           GHC.TypeLits - --------------------------------------------------------------------------------  -- | Triangles in \(d\)-dimensional space.-data Triangle d p r = Triangle (Point d r :+ p)-                               (Point d r :+ p)-                               (Point d r :+ p)+data Triangle d p r = Triangle !(Point d r :+ p)+                               !(Point d r :+ p)+                               !(Point d r :+ p)+                      deriving (Generic) -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)+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)+instance (Arity d, NFData r, NFData p) => NFData (Triangle d p r) +instance Arity d => Bifunctor  (Triangle d) where bimap = bimapDefault+instance Arity d => Bifoldable (Triangle d) where bifoldMap = bifoldMapDefault +instance Arity d => Bitraversable (Triangle d) where+  bitraverse f g (Triangle p q r) = let tr = bitraverse (traverse g) f in+    Triangle <$> tr p <*> tr q <*> tr 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)++instance Field1 (Triangle d p r) (Triangle d p r) (Point d r :+ p) (Point d r :+ p) where+  _1 = lens (\(Triangle p _ _) -> p) (\(Triangle _ q r) p -> Triangle p q r)+instance Field2 (Triangle d p r) (Triangle d p r) (Point d r :+ p) (Point d r :+ p) where+  _2 = lens (\(Triangle _ q _) -> q) (\(Triangle p _ r) q -> Triangle p q r)+instance Field3 (Triangle d p r) (Triangle d p r) (Point d r :+ p) (Point d r :+ p) where+  _3 = lens (\(Triangle _ _ r) -> r) (\(Triangle p q _) r -> Triangle p q r)+ type instance NumType   (Triangle d p r) = r type instance Dimension (Triangle d p r) = d++_TriangleThreePoints :: Iso' (Triangle d p r) (Three (Point d r :+ p))+_TriangleThreePoints = iso (\(Triangle p q r) -> Three p q r) (\(Three p q r) -> Triangle p q r)  instance PointFunctor (Triangle d p) where   pmap f (Triangle p q r) = Triangle (p&core %~ f) (q&core %~ f) (r&core %~ f)

@@ -14,29 +14,29 @@                            , module LV                            , C(..)                            , Affine(..)-                           , qdA, distanceA+                           , quadrance, qdA, distanceA                            , dot, norm, signorm                            , isScalarMultipleOf                            , scalarMultiple                            -- reexports                            , FV.replicate-                           , FV.imap                            , xComponent, yComponent, zComponent                            ) where  import           Control.Applicative (liftA2)-import           Control.Lens(Lens')+import           Control.Lens (Lens')+import           Control.Monad.State import qualified Data.Foldable as F import           Data.Geometry.Properties import           Data.Geometry.Vector.VectorFamily import           Data.Geometry.Vector.VectorFixed (C(..))-import           Data.Maybe 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-import           Test.QuickCheck+import           Linear.Metric (dot,norm,signorm,quadrance)+import           Linear.Vector as LV hiding (E(..))+import           System.Random (Random(..))+import           Test.QuickCheck (Arbitrary(..),infiniteList)  -------------------------------------------------------------------------------- @@ -46,13 +46,21 @@ instance (Arbitrary r, Arity d) => Arbitrary (Vector d r) where   arbitrary = vectorFromListUnsafe <$> infiniteList +instance (Random r, Arity d) => Random (Vector d r) where+  randomR (lows,highs) g0 = flip runState g0 $+                            FV.zipWithM (\l h -> state $ randomR (l,h)) lows highs+  random g0 = flip runState g0 $ FV.replicateM (state random)  -- | 'isScalarmultipleof u v' test if v is a scalar multiple of u. -- -- >>> Vector2 1 1 `isScalarMultipleOf` Vector2 10 10 -- True+-- >>> Vector3 1 1 2 `isScalarMultipleOf` Vector3 10 10 20+-- True -- >>> Vector2 1 1 `isScalarMultipleOf` Vector2 10 1 -- False+-- >>> Vector2 1 1 `isScalarMultipleOf` Vector2 (-1) (-1)+-- True -- >>> Vector2 1 1 `isScalarMultipleOf` Vector2 11.1 11.1 -- True -- >>> Vector2 1 1 `isScalarMultipleOf` Vector2 11.1 11.2@@ -63,11 +71,20 @@ -- True -- >>> Vector2 2 1 `isScalarMultipleOf` Vector2 4 0 -- False+-- >>> Vector3 2 1 0 `isScalarMultipleOf` Vector3 4 0 5+-- False+-- >>> Vector3 0 0 0 `isScalarMultipleOf` Vector3 4 0 5+-- True isScalarMultipleOf       :: (Eq r, Fractional r, Arity d)                          => Vector d r -> Vector d r -> Bool-u `isScalarMultipleOf` v = isJust $ scalarMultiple u v+u `isScalarMultipleOf` v = let d = u `dot` v+                               num = quadrance u * quadrance v+                           in num == 0 || 1 == d*d / num+-- u `isScalarMultipleOf` v = isJust $ scalarMultiple u v {-# SPECIALIZE     isScalarMultipleOf :: (Eq r, Fractional r) => Vector 2 r -> Vector 2 r -> Bool  #-}+{-# SPECIALIZE+    isScalarMultipleOf :: (Eq r, Fractional r) => Vector 3 r -> Vector 3 r -> Bool  #-}  -- | scalarMultiple u v computes the scalar labmda s.t. v = lambda * u (if it exists) scalarMultiple     :: (Eq r, Fractional r, Arity d)

@@ -39,6 +39,7 @@ import           Text.ParserCombinators.ReadPrec (lift) import           Text.Read (Read(..),readListPrecDefault, readPrec_to_P,minPrec) import           Data.Proxy+import           Data.Hashable  -------------------------------------------------------------------------------- -- * d dimensional Vectors@@ -59,8 +60,11 @@ unV = lens _unV (const MKVector) {-# INLINE unV #-} -type Arity d = (ImplicitArity (Peano d), KnownNat d)+-- type Arity d = (ImplicitArity (Peano d), KnownNat d)+class (ImplicitArity (Peano d), KnownNat d) => Arity d+instance (ImplicitArity (Peano d), KnownNat d) => Arity d + deriving instance (Eq r,  Arity d) => Eq  (Vector d r) deriving instance (Ord r, Arity d) => Ord (Vector d r) @@ -69,6 +73,15 @@ deriving instance Arity d => Traversable (Vector d) deriving instance Arity d => Applicative (Vector d) +++instance Arity d => FunctorWithIndex     Int (Vector d) where+  imap = V.imap+instance Arity d => FoldableWithIndex    Int (Vector d)+instance Arity d => TraversableWithIndex Int (Vector d) where+  itraverse = V.imapM++ deriving instance Arity d => Additive (Vector d) deriving instance Arity d => Metric (Vector d) instance Arity d => Affine (Vector d) where@@ -76,6 +89,8 @@   u .-. v = u ^-^ v   p .+^ v = p ^+^ v +deriving instance (Arity d, Hashable r) => Hashable (Vector d r)+ instance Arity d => Ixed (Vector d r) where   ix = element' @@ -167,6 +182,9 @@  -------------------------------------------------------------------------------- -- * Snoccing and consindg++cons   :: (Arity d, Arity (d+1)) => r -> Vector d r -> Vector (d + 1) r+cons x = vectorFromListUnsafe . (x:) . F.toList  -- | Add an element at the back of the vector snoc     :: (Arity (d + 1), Arity d) => Vector d r -> r -> Vector (d + 1) r

@@ -5,7 +5,7 @@ import           Control.Applicative (liftA2) import           Control.DeepSeq import           Control.Lens hiding (element)-import           Data.Aeson(FromJSON(..),ToJSON(..))+import           Data.Aeson (FromJSON(..),ToJSON(..)) -- import           Data.Aeson (ToJSON(..),FromJSON(..)) import qualified Data.Foldable as F import qualified Data.Geometry.Vector.VectorFixed as FV@@ -19,6 +19,7 @@ import qualified Linear.V3 as L3 import qualified Linear.V4 as L4 import           Linear.Vector+import           Data.Hashable  -------------------------------------------------------------------------------- -- * Natural number stuff@@ -77,7 +78,9 @@ unVF = lens _unVF (const VectorFamily) {-# INLINE unVF #-} -type ImplicitArity d = (ImplicitPeano d, V.Arity (FromPeano d))+-- type ImplicitArity d = (ImplicitPeano d, V.Arity (FromPeano d))+class (ImplicitPeano d, V.Arity (FromPeano d)) => ImplicitArity d+instance (ImplicitPeano d, V.Arity (FromPeano d)) => ImplicitArity d  instance (Eq r, ImplicitArity d) => Eq (VectorFamily d r) where   (VectorFamily u) == (VectorFamily v) = case (implicitPeano :: SingPeano d) of@@ -187,6 +190,17 @@                            (SS (SS (SS (SS SZ))))     -> rnf v                            (SS (SS (SS (SS (SS _))))) -> rnf v   {-# INLINE rnf #-}+++instance (ImplicitPeano d, Hashable r) => Hashable (VectorFamily d r) where+  hashWithSalt = case (implicitPeano :: SingPeano d) of+                   SZ                         -> hashWithSalt+                   (SS SZ)                    -> hashWithSalt+                   (SS (SS SZ))               -> hashWithSalt+                   (SS (SS (SS SZ)))          -> hashWithSalt+                   (SS (SS (SS (SS SZ))))     -> hashWithSalt+                   (SS (SS (SS (SS (SS _))))) -> hashWithSalt+  instance ImplicitArity d => Ixed (VectorFamily d r) where   ix = element'

@@ -13,7 +13,8 @@ -- embedding. -- ---------------------------------------------------------------------------------module Data.PlaneGraph.Core( PlaneGraph(PlaneGraph), graph+module Data.PlaneGraph.Core( -- $setup+                             PlaneGraph(PlaneGraph), graph                            , PlanarGraph                            , VertexData(VertexData), vData, location, vtxDataToExt                            , fromSimplePolygon, fromConnectedSegments@@ -87,7 +88,6 @@  -------------------------------------------------------------------------------- - -- $setup -- >>> import Data.Proxy -- >>> import Data.PlaneGraph.AdjRep(Gr(Gr),Face(Face),Vtx(Vtx))@@ -390,17 +390,26 @@ incidentEdges   :: VertexId' s -> PlaneGraph s v e f r -> V.Vector (Dart s) incidentEdges v = PG.incidentEdges v . _graph --- | All incoming edges incident to vertex v, in counterclockwise order around v.++-- | All edges incident to vertex v in incoming direction+-- (i.e. pointing into v) in counterclockwise order around v. --+-- running time: \(O(k)\), where \(k) is the total number of incident edges of v+-- -- >>> incomingEdges (VertexId 1) smallG--- [Dart (Arc 1) -1]+-- [Dart (Arc 1) +1,Dart (Arc 4) -1,Dart (Arc 3) -1] incomingEdges   :: VertexId' s -> PlaneGraph s v e f r -> V.Vector (Dart s) incomingEdges v = PG.incomingEdges v . _graph --- | All outgoing edges incident to vertex v, in counterclockwise order around v.+++-- | All edges incident to vertex v in incoming direction+-- (i.e. pointing into v) in counterclockwise order around v. --+-- running time: \(O(k)\), where \(k) is the total number of incident edges of v+-- -- >>> outgoingEdges (VertexId 1) smallG--- [Dart (Arc 4) +1,Dart (Arc 3) +1]+-- [Dart (Arc 1) -1,Dart (Arc 4) +1,Dart (Arc 3) +1] outgoingEdges   :: VertexId' s -> PlaneGraph s v e f r  -> V.Vector (Dart s) outgoingEdges v = PG.outgoingEdges v . _graph 

@@ -22,9 +22,10 @@                       ) where  import Control.Lens+import Data.Geometry.Matrix import Data.Geometry.Point-import Data.Geometry.Vector import Data.Geometry.Transformation+import Data.Geometry.Vector  -------------------------------------------------------------------------------- 

@@ -133,7 +133,7 @@ <textstyle name="center" begin="\begin{center}" end="\end{center}"/> <textstyle name="itemize" begin="\begin{itemize}" end="\end{itemize}"/> <textstyle name="item" begin="\begin{itemize}\item{}" end="\end{itemize}"/>-</ipestyle><page><group><use pos="172.658536585365 637.048780487804" name="mark/disk(sx)"/><use pos="172.658536585365 637.215447154471" name="mark/disk(sx)"/><use pos="172.658536585365 637.898780487804" name="mark/disk(sx)"/><use pos="172.658536585365 746.000060937506" name="mark/disk(sx)"/><use pos="173.658536585365 638.048780487804" name="mark/disk(sx)"/><use pos="195.199902499989 746.000060937506" name="mark/disk(sx)"/><use pos="252.176396921200 637.048780487804" name="mark/disk(sx)"/><use pos="256.914285714285 707.428571428571" name="mark/disk(sx)"/><use pos="258.512437254286 650.776865588142" name="mark/disk(sx)"/><use pos="278.447793072843 693.970129329472" name="mark/disk(sx)"/><use pos="302.000073125007 745.000060937506" name="mark/disk(sx)"/><use pos="302.461611663469 746.000060937506" name="mark/disk(sx)"/><use pos="303.200073125007 746.000060937506" name="mark/disk(sx)"/><use pos="330.322580645161 661.548387096774" name="mark/disk(sx)"/><use pos="331.322580645161 637.048780487804" name="mark/disk(sx)"/><use pos="331.322580645161 660.923387096774" name="mark/disk(sx)"/><use pos="331.322580645161 661.698387096774" name="mark/disk(sx)"/><use pos="331.322580645161 746.000060937506" name="mark/disk(sx)"/><path>172.658536585365 637.215447154471 m+</ipestyle><page><layer name="alpha"/><view layers="alpha" active="alpha"/><group><use pos="172.658536585365 637.048780487804" name="mark/disk(sx)"/><use pos="172.658536585365 637.215447154471" name="mark/disk(sx)"/><use pos="172.658536585365 637.898780487804" name="mark/disk(sx)"/><use pos="172.658536585365 746.000060937506" name="mark/disk(sx)"/><use pos="173.658536585365 638.048780487804" name="mark/disk(sx)"/><use pos="195.199902499989 746.000060937506" name="mark/disk(sx)"/><use pos="252.176396921200 637.048780487804" name="mark/disk(sx)"/><use pos="256.914285714285 707.428571428571" name="mark/disk(sx)"/><use pos="258.512437254286 650.776865588142" name="mark/disk(sx)"/><use pos="278.447793072843 693.970129329472" name="mark/disk(sx)"/><use pos="302.000073125007 745.000060937506" name="mark/disk(sx)"/><use pos="302.461611663469 746.000060937506" name="mark/disk(sx)"/><use pos="303.200073125007 746.000060937506" name="mark/disk(sx)"/><use pos="330.322580645161 661.548387096774" name="mark/disk(sx)"/><use pos="331.322580645161 637.048780487804" name="mark/disk(sx)"/><use pos="331.322580645161 660.923387096774" name="mark/disk(sx)"/><use pos="331.322580645161 661.698387096774" name="mark/disk(sx)"/><use pos="331.322580645161 746.000060937506" name="mark/disk(sx)"/><path>172.658536585365 637.215447154471 m 173.658536585365 638.048780487804 l </path><path>173.658536585365 638.048780487804 m 256.914285714285 707.428571428571 l