hgeometry 0.9.0.0 → 0.10.0.0

raw patch · 45 files changed

+1516/−997 lines, 45 filesdep +hspecdep +primitivedep −singletonsdep ~QuickCheckdep ~hgeometry-combinatorialdep ~mtlbinary-addedPVP ok

version bump matches the API change (PVP)

Dependencies added: hspec, primitive

Dependencies removed: singletons

Dependency ranges changed: QuickCheck, hgeometry-combinatorial, mtl, quickcheck-instances

API changes (from Hackage documentation)

Files

• changelog.org +18/−1
• docs/Data/PlaneGraph/small.png binary
• hgeometry.cabal +21/−10
• src/Algorithms/Geometry/ClosestPair/DivideAndConquer.hs +5/−26
• src/Algorithms/Geometry/ClosestPair/Naive.hs +29/−18
• src/Algorithms/Geometry/ConvexHull/DivideAndConquer.hs +59/−15
• src/Algorithms/Geometry/ConvexHull/GrahamScan.hs +3/−2
• src/Algorithms/Geometry/ConvexHull/QuickHull.hs +109/−0
• src/Algorithms/Geometry/DelaunayTriangulation/DivideAndConquer.hs +7/−5
• src/Algorithms/Geometry/DelaunayTriangulation/Naive.hs +1/−1
• src/Algorithms/Geometry/Diameter.hs +0/−35
• src/Algorithms/Geometry/Diameter/Naive.hs +32/−0
• src/Algorithms/Geometry/LowerEnvelope/DualCH.hs +2/−2
• src/Algorithms/Geometry/PolygonTriangulation/MakeMonotone.hs +43/−21
• src/Algorithms/Geometry/PolygonTriangulation/Triangulate.hs +7/−6
• src/Algorithms/Geometry/PolygonTriangulation/TriangulateMonotone.hs +5/−5
• src/Algorithms/Geometry/PolygonTriangulation/Types.hs +3/−1
• src/Algorithms/Geometry/SmallestEnclosingBall/Naive.hs +3/−2
• src/Algorithms/Geometry/SmallestEnclosingBall/RIC.hs +164/−0
• src/Algorithms/Geometry/SmallestEnclosingBall/RandomizedIncrementalConstruction.hs +0/−84
• src/Algorithms/Geometry/SmallestEnclosingBall/Types.hs +1/−2
• src/Algorithms/Geometry/WellSeparatedPairDecomposition/WSPD.hs +2/−5
• src/Data/Geometry/Ball.hs +7/−7
• src/Data/Geometry/Box/Internal.hs +16/−15
• src/Data/Geometry/Duality.hs +1/−1
• src/Data/Geometry/Line/Internal.hs +45/−23
• src/Data/Geometry/LineSegment.hs +10/−10
• src/Data/Geometry/PlanarSubdivision/Basic.hs +1/−1
• src/Data/Geometry/PlanarSubdivision/Merge.hs +7/−11
• src/Data/Geometry/Point.hs +120/−83
• src/Data/Geometry/PolyLine.hs +6/−1
• src/Data/Geometry/Polygon.hs +37/−492
• src/Data/Geometry/Polygon/Convex.hs +85/−57
• src/Data/Geometry/Polygon/Core.hs +571/−0
• src/Data/Geometry/Polygon/Extremes.hs +39/−0
• src/Data/Geometry/PrioritySearchTree.hs +1/−1
• src/Data/Geometry/Slab.hs +2/−2
• src/Data/Geometry/Transformation.hs +0/−8
• src/Data/Geometry/Triangle.hs +5/−3
• src/Data/Geometry/Vector.hs +3/−3
• src/Data/Geometry/Vector/VectorFamilyPeano.hs +3/−3
• src/Data/PlaneGraph/Core.hs +2/−1
• src/Data/PlaneGraph/IO.hs +10/−3
• test/Data/PlaneGraph/myPlaneGraph.yaml +18/−18
• test/Data/PlaneGraph/small.yaml +13/−13

@@ -2,7 +2,24 @@  * Changelog -** 0.9 (unreleased)+** 0.10++- renamed the smallest enclosing ball to RIC+- improved tangency finding on convex hulls/chains+- changes to how we order points in ccwCmpAround and cwCmpAround;+  these will report EQ if points appear at the same angle from the+  center point.+- new functions ccwCmpAroundWith and cwCmpAroundWith that allow you to+  specify the direction corresponding to "zero".+- bugfixes, in particular triangulating a polygon with holes now works properly.+- removed some unused dependencies+- we are no longer depending on ghc-plugins; as a result hgeometry+  now also compiles with ghcjs+- more ToJSON/FromJSON instances.+- removed the 'point2' and 'point3' functions in favor of the pattern+  synonyms Point2 and Point3.++** 0.9  - Implemented 2D Linear Programming using randomized incremental   construction (in \(O(n)\) expected time). This allows us to solve

binary file changed (absent → 24478 bytes)

@@ -1,5 +1,5 @@ name:                hgeometry-version:             0.9.0.0+version:             0.10.0.0 synopsis:            Geometric Algorithms, Data structures, and Data types. description:   HGeometry provides some basic geometry types, and geometric algorithms and@@ -44,7 +44,7 @@ extra-source-files:  README.md                      changelog.org -Extra-doc-files:+Extra-doc-files:     docs/Data/PlaneGraph/small.png                      -- docs/**/*.png  cabal-version:       2.0@@ -116,11 +116,13 @@                     -- * Geometric Algorithms                     Algorithms.Geometry.ConvexHull.GrahamScan                     Algorithms.Geometry.ConvexHull.DivideAndConquer+                    Algorithms.Geometry.ConvexHull.QuickHull+                    -- Algorithms.Geometry.ConvexHull.JarvisMarch                      Algorithms.Geometry.LowerEnvelope.DualCH                      Algorithms.Geometry.SmallestEnclosingBall.Types-                    Algorithms.Geometry.SmallestEnclosingBall.RandomizedIncrementalConstruction+                    Algorithms.Geometry.SmallestEnclosingBall.RIC                     Algorithms.Geometry.SmallestEnclosingBall.Naive                      Algorithms.Geometry.DelaunayTriangulation.Types@@ -134,7 +136,7 @@                     Algorithms.Geometry.WellSeparatedPairDecomposition.WSPD                     Algorithms.Geometry.WellSeparatedPairDecomposition.Types -                    Algorithms.Geometry.Diameter+                    Algorithms.Geometry.Diameter.Naive                      -- Algorithms.Geometry.Sweep @@ -160,6 +162,7 @@                      Algorithms.Geometry.FrechetDistance.Discrete +                     -- * Embedded Planar Graphs                     Data.PlaneGraph                     Data.PlaneGraph.Core@@ -171,12 +174,14 @@                     Graphics.Render    other-modules:-+                    -- * Implementation Internals of Polygons+                    Data.Geometry.Polygon.Core+                    Data.Geometry.Polygon.Extremes    -- other-extensions:   build-depends:                 base                    >= 4.11      &&     < 5-              , hgeometry-combinatorial >= 0.9.0.0+              , hgeometry-combinatorial >= 0.10.0.0                , bifunctors              >= 4.1               , bytestring              >= 0.10@@ -185,7 +190,7 @@               , lens                    >= 4.2               , semigroupoids           >= 5               , semigroups              >= 0.18-              , singletons              >= 2.0+              -- , singletons              >= 2.0               , linear                  >= 1.10               , fixed-vector            >= 1.0               , vector-builder          >= 0.3.7@@ -196,7 +201,12 @@               , QuickCheck              >= 2.5               , quickcheck-instances    >= 0.3               , reflection              >= 2.1+              , primitive               >= 0.6.3.0+              -- , singleton-typelits      >= 0.1.0.0 +              -- , ghc-typelits-natnormalise >= 0.6+              -- , ghc-typelits-knownnat     >= 0.6+               , vector                  >= 0.11               , data-clist              >= 0.1.2.3               , text                    >= 1.1.1.0@@ -204,12 +214,14 @@               , aeson                   >= 1.0               , yaml                    >= 0.8 -              , mtl+              , mtl                     >= 2.2               , template-haskell +              , hspec, QuickCheck, quickcheck-instances -  hs-source-dirs: src +  hs-source-dirs: src test+   default-language:    Haskell2010    default-extensions: TypeFamilies@@ -234,7 +246,6 @@                     , DeriveFoldable                     , DeriveTraversable                     , DeriveGeneric-                    , AutoDeriveTypeable                       , FlexibleInstances

@@ -14,15 +14,12 @@                                                        , CP                                                        , CCP(..)                                                        , mergePairs-                                                       , mergeSortedBy-                                                       , mergeSortedListsBy                                                        ) where +import           Algorithms.DivideAndConquer import           Control.Lens-import           Data.BinaryTree import           Data.Ext-import qualified Data.Foldable as F import           Data.Geometry.Point import           Data.LSeq (LSeq) import qualified Data.LSeq as LSeq@@ -30,7 +27,7 @@ import           Data.List.NonEmpty (NonEmpty(..)) import qualified Data.List.NonEmpty as NonEmpty import           Data.Ord (comparing)-import           Data.Semigroup.Foldable(foldMap1, toNonEmpty)+import           Data.Semigroup.Foldable (toNonEmpty) import           Data.UnBounded import           Data.Util @@ -40,11 +37,11 @@ -- \(n\) points. -- -- running time: \(O(n)\)-closestPair :: ( Ord r, Num r) => LSeq 2 (Point 2 r :+ p) -> Two (Point 2 r :+ p)-closestPair = f . foldMap1 mkCCP . asBalancedBinLeafTree . toNonEmpty+closestPair :: (Ord r, Num r) => LSeq 2 (Point 2 r :+ p) -> Two (Point 2 r :+ p)+closestPair = f . divideAndConquer1 mkCCP . toNonEmpty             . LSeq.unstableSortBy (comparing (^.core))   where-    mkCCP (Elem p) = CCP (p :| []) Top+    mkCCP p = CCP (p :| []) Top     f = \case           CCP _ (ValT (SP cp _)) -> cp           CCP _ Top              -> error "closestPair: absurd."@@ -68,24 +65,6 @@       cmp p q = comparing (^.core.yCoord) p q <> comparing (^.core.xCoord) p q  --- | Given an ordering and two nonempty sequences ordered according to that--- ordering, merge them-mergeSortedBy           :: (a -> a -> Ordering) -> NonEmpty a -> NonEmpty a -> NonEmpty a-mergeSortedBy cmp ls rs = NonEmpty.fromList-                        $ mergeSortedListsBy cmp (F.toList ls) (F.toList rs)----- | Given an ordering and two nonempty sequences ordered according to that--- ordering, merge them-mergeSortedListsBy     :: (a -> a -> Ordering) -> [a] -> [a] -> [a]-mergeSortedListsBy cmp = go-  where-    go []         ys     = ys-    go xs         []     = xs-    go xs@(x:xs') ys@(y:ys') = case x `cmp` y of-                                 LT -> x : go xs' ys-                                 EQ -> x : go xs' ys-                                 GT -> y : go xs  ys'  -- | Function that does the actual merging work mergePairs            :: forall p r. (Ord r, Num r)

@@ -9,12 +9,16 @@ -- \(n\) points in \(\mathbb{R}^d\). -- ---------------------------------------------------------------------------------module Algorithms.Geometry.ClosestPair.Naive where+module Algorithms.Geometry.ClosestPair.Naive( closestPair+                                            , closestPairWith+                                            , DistanceFunction+                                            ) where +import           Control.Lens ((^.),_1) import           Data.Ext import qualified Data.Foldable as F-import           Data.Geometry (qdA) import           Data.Geometry.Point+import           Data.Geometry.Properties (NumType) import           Data.Geometry.Vector (Arity) import           Data.LSeq (LSeq) import qualified Data.List.NonEmpty as NonEmpty@@ -23,24 +27,31 @@  -------------------------------------------------------------------------------- --- | Naive algorithm to compute the closest pair in \(d\) dimensions. Runs in--- \(O(n^2)\) time (for any constant \(d\)). Note that we need at least two elements--- for there to be a closest pair.-closestPair :: ( Ord r, Arity d, Num r-               ) => LSeq 2 (Point d r :+ p) -> Two (Point d r :+ p)-closestPair = getVal . getMin . sconcat . fmap (uncurry' mkPair) . pairs-  where-    uncurry' f (Two a b) = f a b-    getVal (Arg _ x) = x+-- | Naive algorithm to compute the closest pair according to the+-- (squared) Euclidean distance in \(d\) dimensions. Note that we need+-- at least two elements for there to be a closest pair.+--+-- running time: \(O(dn^2)\) time.+closestPair :: ( Ord r, Arity d, Num r)+            => LSeq 2 (Point d r :+ p) -> Two (Point d r :+ p)+closestPair = (^._1) . closestPairWith (\p q -> squaredEuclideanDist (p^.core) (q^.core)) --- | A pair of points-type PP d p r = ArgMin r (Two (Point d r :+ p)) --- | Create a pair of points-mkPair                         :: (Arity d, Num r)-                               => Point d r :+ p -> Point d r :+ p -> PP d p r-mkPair pp@(p :+ _) qq@(q :+ _) = let dst = qdA p q-                                 in Min (Arg dst (Two pp qq))+type DistanceFunction g = g -> g -> NumType g++-- | Naive algorithm to compute the closest pair of points (and the+-- distance realized by those points) given a distance function.  Note+-- that we need at least two elements for there to be a closest pair.+--+-- running time: \(O(T(d)n^2)\), where \(T(d)\) is the time required+-- to evaluate the distance between two points in \(\mathbb{R}^d\).+closestPairWith   :: Ord r+                  => DistanceFunction (Point d r :+ p)+                  -> LSeq 2 (Point d r :+ p) -> SP (Two (Point d r :+ p)) r+closestPairWith d = getVal . getMin . sconcat . fmap mkPair . pairs+  where+    getVal (Arg dist x) = SP x dist+    mkPair (Two p q)    = Min (Arg (d p q) (Two p q))  -- | Produce all lists from a vec of elements. Since the Vec contains at least two -- elements, the resulting list is non-empty

@@ -9,30 +9,74 @@ -- of a set of \(n\) points in \(\mathbb{R}^2\). -- ---------------------------------------------------------------------------------module Algorithms.Geometry.ConvexHull.DivideAndConquer(convexHull) where+module Algorithms.Geometry.ConvexHull.DivideAndConquer( convexHull+                                                      , upperHull+                                                      , lowerHull+                                                      ) where -import           Control.Lens ((^.))-import           Data.BinaryTree+import           Algorithms.DivideAndConquer+import           Control.Arrow ((&&&))+import           Control.Lens ((^.), to) import           Data.Ext-import           Data.Function (on) import           Data.Geometry.Point import           Data.Geometry.Polygon import           Data.Geometry.Polygon.Convex+import           Data.List.NonEmpty (NonEmpty(..)) import qualified Data.List.NonEmpty as NonEmpty-import           Data.Semigroup.Foldable-+import           Data.Util --------------------------------------------------------------------------------  -- | \(O(n \log n)\) time ConvexHull using divide and conquer. The resulting polygon is -- given in clockwise order.-convexHull :: (Ord r, Num r)-           => NonEmpty.NonEmpty (Point 2 r :+ p) -> ConvexPolygon p r-convexHull = unMerge-           . foldMap1 (Merge . ConvexPolygon . fromPoints . (:[]) . _unElem)-           . asBalancedBinLeafTree-           . NonEmpty.sortBy (compare `on` (^.core))+convexHull           :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> ConvexPolygon p r+convexHull (p :| []) = ConvexPolygon . fromPoints $ [p]+convexHull pts       = combine . (upperHull' &&& lowerHull') . NonEmpty.sortBy incXdecY $ pts+  where+    combine (l:|uh,_:|lh) = ConvexPolygon . fromPoints $ l : uh <> reverse (init lh) -newtype Merge r p = Merge { unMerge :: ConvexPolygon p r }+----------------------------------------+-- * Computing a lower hull -instance (Num r, Ord r) => Semigroup (Merge r p) where-  (Merge lp) <> (Merge rp) = let (ch,_,_) = merge lp rp in Merge ch+-- | \(O(n \log n)\) time LowerHull using divide and conquer. The resulting Hull is+-- given from left to right, i.e. in counter clockwise order.+lowerHull :: (Ord r, Num r)+          => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)+lowerHull = lowerHull' . NonEmpty.sortBy incXdecY++lowerHull' :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)+lowerHull' = unLH . divideAndConquer1 (LH . (:|[]))++newtype LH r p = LH { unLH :: NonEmpty (Point 2 r :+ p) } deriving (Eq,Show)++instance (Num r, Ord r) => Semigroup (LH r p) where+  (LH lh) <> (LH rh) = LH $ hull lowerTangent' lh rh++----------------------------------------+-- * Computing an upper hull++-- | \(O(n \log n)\) time UpperHull using divide and conquer. The resulting Hull is+-- given from left to right, i.e. in clockwise order.+upperHull :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)+upperHull = upperHull' . NonEmpty.sortBy incXdecY++upperHull' :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p)+upperHull' = unUH . divideAndConquer1 (UH . (:|[]))++newtype UH r p = UH { unUH :: NonEmpty (Point 2 r :+ p) }++instance (Num r, Ord r) => Semigroup (UH r p) where+  (UH lh) <> (UH rh) = UH $ hull upperTangent' lh rh++----------------------------------------++-- | The function that does the actual merging part+hull               :: (NonEmpty p -> NonEmpty p -> Two (p :+ [p]))+                   -> NonEmpty p -> NonEmpty p -> NonEmpty p+hull tangent lh rh = let Two (l :+ lh') (r :+ rh') = tangent (NonEmpty.reverse lh) rh+                     in NonEmpty.fromList $ (reverse lh') <> [l,r] <> rh'++--------------------------------------------------------------------------------++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

@@ -22,10 +22,11 @@                            lh  = NonEmpty.tail . hull' $ reverse ps'                        in ConvexPolygon . fromPoints . reverse $ lh ++ uh +-- | Computes the upper hull. 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 = hull id-+upperHull = NonEmpty.reverse . hull id +-- | Computes the upper hull. The upper hull is given from left to right lowerHull :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p) lowerHull = hull reverse 

@@ -0,0 +1,109 @@+module Algorithms.Geometry.ConvexHull.QuickHull( convexHull ) where++import           Control.Lens ((^.),(&),(.~))+import           Data.Ext+import qualified Data.Foldable as F+import           Data.Geometry.Line+import           Data.Geometry.Point+import           Data.Geometry.Polygon+import           Data.Geometry.Polygon.Convex (ConvexPolygon(..))+import           Data.Geometry.Triangle+import qualified Data.List as List+import           Data.List.NonEmpty (NonEmpty(..))+import           Data.Ord (comparing)+import           Data.Util+++-- import Data.Ratio+-- import qualified Data.List.NonEmpty as NonEmpty+-- import qualified Algorithms.Geometry.ConvexHull.GrahamScan as GC+-- import Debug.Trace+--------------------------------------------------------------------------------++-- | ConvexHull using Quickhull. The resulting polygon is given in+-- clockwise order.+--+-- running time: \(O(n^2)\)+convexHull            :: (Ord r, Fractional r, Show r, Show p)+                      => NonEmpty (Point 2 r :+ p) -> ConvexPolygon p r+convexHull (p :| []) = ConvexPolygon . fromPoints $ [p]+convexHull ps        = ConvexPolygon . fromPoints+                     $ [l] <> hull l r above <> [r] <> (reverse $ hull l r below)+  where+    STR l r mids  = findExtremes ps+    m             = lineThrough (l^.core) (r^.core)+    (above,below) = List.partition (\(p :+ _) -> p `liesAbove` m) mids++-- | Finds the leftmost and rightmost point in the list+findExtremes            :: Ord r+                        => NonEmpty (Point 2 r :+ q)+                        -> STR (Point 2 r :+ q) (Point 2 r :+ q) [Point 2 r :+ q]+findExtremes (p :| pts ) = foldr f (STR p p []) pts+  where+    f q (STR l r ms) = case (incXdecY q l, incXdecY q r) of+                         (LT,_)  -> STR q r (addIfNot r l ms)+                         (EQ,_)  -> STR l r ms -- ditch q; it is the same as l+                         (GT,GT) -> STR l q (addIfNot l r ms)+                         (GT,EQ) -> STR l r ms -- ditch q; it is the same as r+                         (GT,LT) -> STR l r (q:ms)++    addIfNot y x xs | (x^.core) /= (y^.core) = x:xs+                    | otherwise              = xs++-- findExtremesBy         :: (a -> a -> Ordering)+--                        -> NonEmpty a+--                        -> STR a a [a]+-- findExtremesBy cmp pts = let l = F.minimumBy cmp pts+--                              r = F.maximumBy cmp pts+--                              a /=. b = a `cmp` b /= EQ+--                          in STR l r [p | p <- F.toList pts, p /=. l, p /=. r]+++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++-- | include neigher left or right+--+hull         :: (Fractional r, Ord r)+             => Point 2 r :+ p -> Point 2 r :+ p -> [Point 2 r :+ p] -> [Point 2 r :+ p]+hull _ _ []  = []+hull l r pts = hull l mid ls <> [mid] <> hull mid r rs+  where+    m       = lineThrough (l^.core) (r^.core)+    mid     = F.maximumBy (comparing dist) pts++    dist (p :+ _) = p `sqDistanceTo` m+    t       = Triangle l mid r+    -- line through l and mid, which splits the remaining points in a left half and a right half.+    splitL   = lineThrough (l^.core) (mid^.core)+    rightSide = (r^.core) `onSide` splitL -- define the side containing r the right side++    (ls,rs) = List.partition (\(p :+ _) -> p `onSide` splitL /= rightSide)+            . filter (\(p :+ _) -> not $ p `onTriangle` t) $ pts+++-- mPoint2 [x,y] = Point2 x y++-- -- testPoints = NonEmpty.fromList+-- --   [ mPoint2 [22536303956634 % 7570647828779,(-5816376064439) % 1228319866920] :+ 1+-- --   , mPoint2 [(-3136920648983) % 824638230353,(-14583744643665) % 9604445576558] :+ 2+-- --   , mPoint2 [(-11653462784667) % 6525086575987,(-598434515815) % 1364557986096] :+ 3+-- --   , mPoint2 [(-7841595901661) % 3282967141364,(-207167076115) % 482378191549] :+ 4+-- --   ]+++-- testPoints :: NonEmpty (Point 2 Rational :+ Int)+-- testPoints = read "(Point2 [(-11199966464450) % 1365514034959,4065659138075 % 2296468530516] :+ 1) :| [Point2 [86686001553073 % 2736621704548,(-63774454571048) % 1880665081093] :+ 2,Point2 [(-77322324895231) % 8260610289790,(-41165682123514) % 2291705829063] :+ 3,Point2 [2292642905947 % 2659329735076,87045289726355 % 2752214350419] :+ 4]"+++-- toDouble          :: Point 2 Rational :+ a -> Point 2 Double :+ a+-- toDouble (p :+ x) = (realToFrac <$> p) :+ x+++--        -- Falsified (after 36 tests and 4 shrinks):++++-- testPoints2 :: NonEmpty (Point 2 Rational :+ Int)+-- testPoints2 = read "(Point2 [(-9876468593885) % 9254762894818,(-34982972348716) % 7450362538495] :+ 1) :| [Point2 [(-11974177119403) % 7705693443554,(-37634868551543) % 9311528788922] :+ 2,Point2 [(-32383659855458) % 9565531378857,20253950785876 % 8268868939819] :+ 3,Point2 [42425655100996 % 8786996213535,(-7972873491283) % 1604043452399] :+ 4]"

@@ -162,6 +162,7 @@                    True  -> (,False) <$> rotateR' l r r1 (pred' r1)                    False -> pure (r1,True) + -- | The code that does the actual rotating rotateR'     :: (Ord r, Fractional r)              => VertexID -> VertexID -> Vertex -> Vertex -> Merge p r Vertex@@ -231,7 +232,8 @@                       . IM.adjustWithKey (insert'' u) v   where     -- inserts b into the adjacency list of a-    insert'' bi ai = CU.insertOrdBy (cwCmpAround (ptMap V.! ai) `on` (ptMap V.!)) bi+    insert'' bi ai = CU.insertOrdBy (cwCmpAround' (ptMap V.! ai) `on` (ptMap V.!)) bi+    cwCmpAround' c p q = cwCmpAround c p q <> cmpByDistanceTo c p q   -- | Deletes an edge@@ -241,21 +243,21 @@     delete' x = CL.filterL (/= x) -- should we rotate left or right if it is the focus?  -- -- | Lifted version of Convex.IsLeftOf isLeftOf           :: (Ord r, Num r)                    => VertexID -> (VertexID, VertexID) -> Merge p r Bool p `isLeftOf` (l,r) = withPtMap . snd . fst <$> ask   where-    withPtMap ptMap = (ptMap V.! p) `Convex.isLeftOf` (ptMap V.! l, ptMap V.! r)+    withPtMap ptMap = (ptMap V.! p) `isLeftOf'` (ptMap V.! l, ptMap V.! r)+    a `isLeftOf'` (b,c) = ccw' b c a == CCW  -- | Lifted version of Convex.IsRightOf isRightOf           :: (Ord r, Num r)                     => VertexID -> (VertexID, VertexID) -> Merge p r Bool p `isRightOf` (l,r) = withPtMap . snd . fst <$> ask   where-    withPtMap ptMap = (ptMap V.! p) `Convex.isRightOf` (ptMap V.! l, ptMap V.! r)+    withPtMap ptMap = (ptMap V.! p) `isRightOf'` (ptMap V.! l, ptMap V.! r)+    a `isRightOf'` (b,c) = ccw' b c a == CW  -------------------------------------------------------------------------------- -- * Some Helper functions

@@ -59,7 +59,7 @@ -- running time: O(m log m), where m=|vs| is the number of vertices to sort. sortAround'               :: (Num r, Ord r)                           => Mapping p r -> VertexID -> [VertexID] -> [VertexID]-sortAround' (_,ptsV) u vs = reverse . map (^.extra) $ sortArround (f u) (map f vs)+sortAround' (_,ptsV) u vs = reverse . map (^.extra) $ sortAround (f u) (map f vs)   where     f v = (ptsV V.! v)&extra .~ v 

@@ -1,35 +0,0 @@-module Algorithms.Geometry.Diameter where--import Control.Lens-import Data.Ext-import Data.Geometry-import Data.List(maximumBy)-------------------------------------------------------------------------------------diameterNaive :: (Ord r, Floating r, Arity d) => [Point d r :+ p] -> r-diameterNaive = maybe 0 (\(p,q) -> euclideanDist (p^.core) (q^.core))-              . diametralPairNaive---- | Computes the Euclidean diametral pair by naively trying all pairs.------ running time: \(O(n^2)\)-diametralPairNaive :: (Ord r, Num r, Arity d)-                   => [Point d r :+ p] -> Maybe (Point d r :+ p, Point d r :+ p)-diametralPairNaive = diametralPairWithNaive squaredEuclideanDist----- | Given a distance function and a list of points pts, computes the diametral--- pair by naively trying all pairs.------ running time: \(O(n^2)\)-diametralPairWithNaive               :: Ord r-                                     => (Point d r -> Point d r -> r)-                                     -> [Point d r :+ p]-                                     -> Maybe (Point d r :+ p, Point d r :+ p)-diametralPairWithNaive f pts@(_:_:_) = Just $ maximumBy cmp [ (p,q) | p <- pts, q <- pts ]-  where-    f' (p,q) = f (p^.core) (q^.core)-    tp `cmp` tq = f' tp `compare` f' tq-diametralPairWithNaive _ _           = Nothing

@@ -0,0 +1,32 @@+module Algorithms.Geometry.Diameter.Naive where++import Control.Lens+import Data.Ext+import Data.Geometry+import Data.List(maximumBy)++--------------------------------------------------------------------------------++diameter :: (Ord r, Floating r, Arity d) => [Point d r :+ p] -> r+diameter = maybe 0 (\(p,q) -> euclideanDist (p^.core) (q^.core)) . diametralPair++-- | Computes the Euclidean diametral pair by naively trying all pairs.+--+-- running time: \(O(n^2)\)+diametralPair :: (Ord r, Num r, Arity d)+                   => [Point d r :+ p] -> Maybe (Point d r :+ p, Point d r :+ p)+diametralPair = diametralPairWith squaredEuclideanDist++-- | Given a distance function and a list of points pts, computes the diametral+-- pair by naively trying all pairs.+--+-- running time: \(O(n^2)\)+diametralPairWith               :: Ord r+                                     => (Point d r -> Point d r -> r)+                                     -> [Point d r :+ p]+                                     -> Maybe (Point d r :+ p, Point d r :+ p)+diametralPairWith f pts@(_:_:_) = Just $ maximumBy cmp [ (p,q) | p <- pts, q <- pts ]+  where+    f' (p,q) = f (p^.core) (q^.core)+    tp `cmp` tq = f' tp `compare` f' tq+diametralPairWith _ _           = Nothing

@@ -16,11 +16,11 @@ type Envelope a r = NonEmpty (Line 2 r :+ a)  -- | Given a list of non-vertical lines, computes the lower envelope using--- duality.+-- duality. The lines are given in left to right order. -- -- \(O(n\log n)\) lowerEnvelope :: (Ord r, Fractional r) => NonEmpty (Line 2 r :+ a) -> Envelope a r-lowerEnvelope = lowerEnvelopeWith upperHull+lowerEnvelope = NonEmpty.reverse . lowerEnvelopeWith upperHull   type UpperHullAlgorithm a r = NonEmpty (Point 2 r :+ a) -> NonEmpty (Point 2 r :+ a)

@@ -1,7 +1,13 @@ {-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE ScopedTypeVariables #-}-module Algorithms.Geometry.PolygonTriangulation.MakeMonotone where+module Algorithms.Geometry.PolygonTriangulation.MakeMonotone( makeMonotone+                                                            , computeDiagonals ++                                                            , VertexType(..)+                                                            , classifyVertices+                                                            ) where+ import           Algorithms.Geometry.LineSegmentIntersection.BentleyOttmann ( xCoordAt                                                                             , ordAt) import           Algorithms.Geometry.PolygonTriangulation.Types@@ -29,16 +35,15 @@ import qualified Data.Vector.Mutable as MV  +-- import Debug.Trace+-- import qualified          Data.CircularSeq as CC ----------------------------------------------------------------------------------  data VertexType = Start | Merge | Split | End | Regular deriving (Show,Read,Eq) ---- How about the hole vertices?- -- | assigns a vertex type to each vertex ----- pre: the polygon is given in CCW order+-- pre: Both the outer boundary and the inner boundary of the polygon are given in CCW order. -- -- running time: \(O(n)\). classifyVertices                     :: (Num r, Ord r)@@ -119,7 +124,7 @@     f = first (\i -> vertexInfo^.ix' i._2)      pg :: Polygon t (SP Int (p :+ VertexType)) r-    pg = numberVertices . classifyVertices . toCounterClockWiseOrder $ p'+    pg = numberVertices . holesToCW . classifyVertices . toCCW $ p'     vertexInfo :: V.Vector (STR (Point 2 r) p VertexType)     vertexInfo = let vs = polygonVertices pg                      n  = F.length vs@@ -137,9 +142,18 @@     sweep'' :: NonEmpty.NonEmpty (Event r) -> Sweep p r ()     sweep'' = mapM_ handle +    -- make everything counterclockwise+    toCCW p = (toCounterClockWiseOrder' p)&polygonHoles'.traverse %~ toCounterClockWiseOrder'+    -- make the holes clockwise:+    holesToCW p = p&polygonHoles'.traverse %~ toClockwiseOrder'+++ -- | Computes a set of diagionals that decompose the polygon into y-monotone -- pieces. --+-- pre: the polygon boundary is given in counterClockwise order.+-- -- running time: \(O(n\log n)\) makeMonotone      :: (Fractional r, Ord r)                   => proxy s -> Polygon t p r@@ -277,21 +291,21 @@   -- testPolygon :: SimplePolygon Int Rational--- testPolygon = fromPoints [ point2 20 20 :+ 1---                          , point2 18 19 :+ 2---                          , point2 16 25 :+ 3---                          , point2 13 23 :+ 4---                          , point2 10 24 :+ 5---                          , point2 6  22 :+ 6---                          , point2 8  21 :+ 7---                          , point2 7  18 :+ 8---                          , point2 2  19 :+ 9---                          , point2 1  10 :+ 10---                          , point2 3  5  :+ 11---                          , point2 11 7  :+ 12---                          , point2 15 1  :+ 13---                          , point2 12 15 :+ 14---                          , point2 15 12 :+ 15+-- testPolygon = fromPoints [ Point2 20 20 :+ 1+--                          , Point2 18 19 :+ 2+--                          , Point2 16 25 :+ 3+--                          , Point2 13 23 :+ 4+--                          , Point2 10 24 :+ 5+--                          , Point2 6  22 :+ 6+--                          , Point2 8  21 :+ 7+--                          , Point2 7  18 :+ 8+--                          , Point2 2  19 :+ 9+--                          , Point2 1  10 :+ 10+--                          , Point2 3  5  :+ 11+--                          , Point2 11 7  :+ 12+--                          , Point2 15 1  :+ 13+--                          , Point2 12 15 :+ 14+--                          , Point2 15 12 :+ 15 --                          ]  -- vertexTypes = [Start,Merge,Start,Merge,Start,Regular,Regular,Merge,Start,Regular,End,Split,End,Split,End]@@ -307,3 +321,11 @@ --                      ] --                outFile = "/Users/frank/tmp/out.ipe" --            writeIpeFile outFile . singlePageFromContent $ out+++-- myPoly :: Polygon Multi () Rational+-- myPoly = MultiPolygon (CC.fromList $ read "[Point2 [16 % 1,80 % 1] :+ (),Point2 [16 % 1,16 % 1] :+ (),Point2 [144 % 1,16 % 1] :+ (),Point2 [144 % 1,80 % 1] :+ ()]"+--                       )+--   [ fromPoints $ read "[Point2 [88 % 1,48 % 1] :+ (),Point2 [112 % 1,40 % 1] :+ (),Point2 [112 % 1,48 % 1] :+ (),Point2 [80 % 1,56 % 1] :+ ()]"+--   , fromPoints $ read "[Point2 [32 % 1,64 % 1] :+ (),Point2 [32 % 1,32 % 1] :+ (),Point2 [64 % 1,32 % 1] :+ (),Point2 [64 % 1,64 % 1] :+ ()]"+--   ]

@@ -31,8 +31,8 @@ -- -- running time: \(O(n \log n)\) triangulate'        :: (Ord r, Fractional r)-                   => proxy s -> Polygon t p r-                   -> PlaneGraph s p PolygonEdgeType PolygonFaceData r+                    => proxy s -> Polygon t p r+                    -> PlaneGraph s p PolygonEdgeType PolygonFaceData r triangulate' px pg' = constructGraph px e es diags   where     (pg, diags)   = computeDiagonals' pg'@@ -56,16 +56,17 @@ computeDiagonals' pg' = (pg, monotoneDiags <> extraDiags)   where     pg            = toCounterClockWiseOrder pg'-    monotoneP     = MM.makeMonotone (Identity pg') pg -- use some arbitrary proxy type+    monotoneP     = MM.makeMonotone (Identity pg) pg -- use some arbitrary proxy type     -- outerFaceId'  = outerFaceId monotoneP      monotoneDiags = map (^._2.core) . filter (\e' -> e'^._2.extra == Diagonal)                   . F.toList . edgeSegments $ monotoneP     extraDiags    = concatMap (TM.computeDiagonals . toCounterClockWiseOrder')                   . lefts . map (^._2.core)+                  . filter (\mp -> mp^._2.extra == Inside) -- triangulate only the insides                   -- . filter (\f -> f^._1 /= outerFaceId')                   . F.toList . rawFacePolygons $ monotoneP -    -- we alredy know we get the polgyons in *clockwise* order, so skip the-    -- check if it is counter clockwise-    toCounterClockWiseOrder' = reverseOuterBoundary+    -- -- we alredy know we get the polgyons in *clockwise* order, so skip the+    -- -- check if it is counter clockwise+    -- toCounterClockWiseOrder'' = reverseOuterBoundary

@@ -145,11 +145,11 @@  -------------------------------------------------------------------------------- --- testPolygon = fromPoints . map ext $ [ point2 10 10---                                      , point2 5 20---                                      , point2 3 14---                                      , point2 1 1---                                      , point2 8 8 ]+-- testPolygon = fromPoints . map ext $ [ Point2 10 10+--                                      , Point2 5 20+--                                      , Point2 3 14+--                                      , Point2 1 1+--                                      , Point2 8 8 ]   

@@ -12,6 +12,7 @@ import qualified Data.PlaneGraph as PG import qualified Data.Vector as V import qualified Data.Vector.Mutable as MV+ --------------------------------------------------------------------------------  data PolygonEdgeType = Original | Diagonal@@ -83,7 +84,7 @@ constructGraph px e origs diags =     subdiv & PG.vertexData.traverse  %~ NonEmpty.head            & PG.faceData             .~ faceData'-           & PG.rawDartData.traverse  %~ snd+           & PG.rawDartData.traverse %~ snd   where     subdiv :: PG.PlaneGraph s (NonEmpty p) (Bool,PolygonEdgeType) () r     subdiv = PG.fromConnectedSegments px $ e' : origs' <> diags'@@ -101,6 +102,7 @@      -- the interior faces     intFaces = flip PG.leftFace subdiv <$> queryDarts+     faceData' :: V.Vector PolygonFaceData     faceData' = V.create $ do                   v' <- MV.replicate (PG.numFaces subdiv) Outside

@@ -9,7 +9,9 @@ -- points in \(\mathbb{R}^2\) -- ---------------------------------------------------------------------------------module Algorithms.Geometry.SmallestEnclosingBall.Naive where+module Algorithms.Geometry.SmallestEnclosingBall.Naive( smallestEnclosingDisk+                                                      , enclosesAll+                                                      ) where  -- just for the types import Control.Lens@@ -55,7 +57,6 @@                            => [Point 2 r :+ p] -> [DiskResult p r] -> DiskResult p r smallestEnclosingDisk' pts = minimumBy (compare `on` (^.enclosingDisk.squaredRadius))                            . filter (flip enclosesAll pts)-  -- | check if a disk encloses all points enclosesAll   :: (Num r, Ord r) => DiskResult p r -> [Point 2 r :+ q] -> Bool

@@ -0,0 +1,164 @@+{-# LANGUAGE DeriveFunctor  #-}+{-# LANGUAGE TemplateHaskell  #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Algorithms.Geometry.SmallestEnclosingBall.RIC+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- An randomized algorithm to compute the smallest enclosing disk of a set of+-- \(n\) points in \(\mathbb{R}^2\). The expected running time is \(O(n)\).+--+--------------------------------------------------------------------------------+module Algorithms.Geometry.SmallestEnclosingBall.RIC(+    smallestEnclosingDisk'+  , smallestEnclosingDisk+  , smallestEnclosingDiskWithPoint+  , smallestEnclosingDiskWithPoints+  ) where++import           Algorithms.Geometry.SmallestEnclosingBall.Types+import           Control.Lens+import           Control.Monad.Random.Class+import           Data.Ext+import qualified Data.Foldable as F+import           Data.Geometry+import           Data.Geometry.Ball+import qualified Data.List as List+import           Data.List.NonEmpty(NonEmpty(..))+import           Data.Maybe (fromMaybe, mapMaybe, catMaybes)+import           Data.Ord (comparing)+import           System.Random.Shuffle (shuffle)++import Debug.Trace++--------------------------------------------------------------------------------++-- | Compute the smallest enclosing disk of a set of points,+-- implemented using randomized incremental construction.+--+-- pre: the input has at least two points.+--+-- running time: expected \(O(n)\) time, where \(n\) is the number of input points.+smallestEnclosingDisk           :: (Ord r, Fractional r, MonadRandom m+                                               -- , Show r, Show p+                                   )+                                => [Point 2 r :+ p]+                                -> m (DiskResult p r)++smallestEnclosingDisk pts@(_:_:_) = ((\(p:q:pts') -> smallestEnclosingDisk' p q pts')+                                    . F.toList) <$> shuffle pts+smallestEnclosingDisk _           = error "smallestEnclosingDisk: Too few points"++-- | Smallest enclosing disk.+smallestEnclosingDisk'     :: (Ord r, Fractional r+                                               -- , Show r, Show p++                              )+                           => Point 2 r :+ p -> Point 2 r :+ p -> [Point 2 r :+ p]+                           -> DiskResult p r+smallestEnclosingDisk' a b = foldr addPoint (initial a b) . List.tails+  where+    -- The empty case occurs only initially+    addPoint []      br            = br+    addPoint (p:pts) br@(DiskResult d _)+      | (p^.core) `inClosedBall` d = br+      | otherwise                  = fromJust' $ smallestEnclosingDiskWithPoint p (a :| (b : pts))+    fromJust' = fromMaybe (error "smallestEncosingDisk' : fromJust, absurd")++-- | Smallest enclosing disk, given that p should be on it.+smallestEnclosingDiskWithPoint              :: (Ord r, Fractional r+                                               -- , Show r, Show p+                                               )+                                            => Point 2 r :+ p -> NonEmpty (Point 2 r :+ p)+                                            -> Maybe (DiskResult p r)+smallestEnclosingDiskWithPoint p (a :| pts) = foldr addPoint (Just $ initial p a) $ List.tails pts+  where+    addPoint []       br   = br+    addPoint (q:pts') br@(Just (DiskResult d _))+      | (q^.core) `inClosedBall` d = br+      | otherwise                  = smallestEnclosingDiskWithPoints p q (a:pts')+    addPoint _        br           = br+++-- | Smallest enclosing disk, given that p and q should be on it+--+-- running time: \(O(n)\)+smallestEnclosingDiskWithPoints        :: (Ord r, Fractional r+                                          -- , Show r, Show p+                                          )+                                       => Point 2 r :+ p -> Point 2 r :+ p -> [Point 2 r :+ p]+                                       -> Maybe (DiskResult p r)+smallestEnclosingDiskWithPoints p q ps = minimumOn (^.enclosingDisk.squaredRadius)+                                       $ catMaybes [mkEnclosingDisk dl, mkEnclosingDisk dr, mdc]+  where+    centers = mapMaybe disk' ps+    -- generate a disk with p q and r+    disk' r = (r:+) <$> disk (p^.core) (q^.core) (r^.core)++    -- partition the points in to those on the left and those on the+    -- right.  Note that centers still contains only those points (and+    -- disks) for which the three points are not colinear. So the+    -- points are either on the left or on the right.+    (leftCenters,rightCenters) = List.partition (\(r :+ _) -> ccw' p q r == CCW) centers+    -- note that we consider 'leftmost' with respect to going from p+    -- to q. This does not really have a global meaning.++    -- we need to find the leftmost and rightmost center on the+    -- bisector. In case there are left-centers, this means that among+    -- the left centers we want to find the point that is furthest way+    -- from p (or q). If there are no left-centers, we with to find+    -- the closest one among the right-centers.+    leftDist z = let c = z^.extra.center+                     s = if ccw' p q c == CCW then 1 else -1+                 in s * squaredEuclideanDist (p^.core) (c^.core)++    dl = maximumOn leftDist leftCenters  -- disk that has the "leftmost" center+    dr = minimumOn leftDist rightCenters -- disk that has the "rightmost" center++    -- diameteral disk+    dd = fromDiameter (p^.core) (q^.core)+    mdc | isEnclosingDisk dd ps = Just $ DiskResult dd (Two p q)+        | otherwise             = Nothing++    -- test if d is an enclosing disk.+    mkEnclosingDisk  md = md >>= mkEnclosingDisk'+    mkEnclosingDisk' (r :+ d) | isEnclosingDisk d ps = Just (DiskResult d (Three p q r))+                              | otherwise            = Nothing+++isEnclosingDisk   :: (Foldable t, Ord r, Num r)+                  => Disk p r -> t (Point 2 r :+ extra) -> Bool+isEnclosingDisk d = all (\s -> (s^.core) `inClosedBall` d)++-- | Constructs the initial 'DiskResult' from two points+initial     :: Fractional r => Point 2 r :+ p -> Point 2 r :+ p -> DiskResult p r+initial p q = DiskResult (fromDiameter (p^.core) (q^.core)) (Two p q)++maximumOn   :: Ord b => (a -> b) -> [a] -> Maybe a+maximumOn f = \case+    [] -> Nothing+    xs -> Just $ List.maximumBy (comparing f) xs++minimumOn   :: Ord b => (a -> b) -> [a] -> Maybe a+minimumOn f = \case+    [] -> Nothing+    xs -> Just $ List.minimumBy (comparing f) xs+++--------------------------------------------------------------------------------++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]++disk'' r = (r:+) <$> disk (p^.core) (q^.core) (r^.core)+  where+    p = ext $ Point2 0 (-6)+    q = ext $ Point2 0 6

@@ -1,84 +0,0 @@-{-# LANGUAGE DeriveFunctor  #-}-{-# LANGUAGE TemplateHaskell  #-}------------------------------------------------------------------------------------ |--- Module      :  Algorithms.Geometry.SmallestEnclosingBall.RandomizedIncrementalConstruction--- Copyright   :  (C) Frank Staals--- License     :  see the LICENSE file--- Maintainer  :  Frank Staals------ An randomized algorithm to compute the smallest enclosing disk of a set of--- \(n\) points in \(\mathbb{R}^2\). The expected running time is \(O(n)\).-------------------------------------------------------------------------------------module Algorithms.Geometry.SmallestEnclosingBall.RandomizedIncrementalConstruction where--import           Algorithms.Geometry.SmallestEnclosingBall.Types-import           Control.Lens-import           Control.Monad.Random.Class-import           Data.Ext-import qualified Data.Foldable as F-import           Data.Geometry-import           Data.Geometry.Ball-import qualified Data.List as L-import           Data.List.NonEmpty-import           Data.Maybe (fromMaybe)-import           System.Random.Shuffle (shuffle)-------------------------------------------------------------------------------------- | O(n) expected time algorithm to compute the smallest enclosing disk of a--- set of points. we need at least two points.--- implemented using randomized incremental construction-smallestEnclosingDisk           :: (Ord r, Fractional r, MonadRandom m)-                                => [Point 2 r :+ p]-                                -> m (DiskResult p r)--smallestEnclosingDisk pts@(_:_:_) = ((\(p:q:pts') -> smallestEnclosingDisk' p q pts')-                                    . F.toList) <$> shuffle pts-smallestEnclosingDisk _           = error "smallestEnclosingDisk: Too few points"---- | Smallest enclosing disk.-smallestEnclosingDisk'     :: (Ord r, Fractional r)-                           => Point 2 r :+ p -> Point 2 r :+ p -> [Point 2 r :+ p]-                           -> DiskResult p r-smallestEnclosingDisk' a b = foldr addPoint (initial a b) . L.tails-  where-    -- The epty case occurs only initially-    addPoint []      br   = br-    addPoint (p:pts) br@(DiskResult d _)-      | (p^.core) `inClosedBall` d = br-      | otherwise                  = smallestEnclosingDiskWithPoint p (a :| (b : pts))----- | Smallest enclosing disk, given that p should be on it.-smallestEnclosingDiskWithPoint              :: (Ord r, Fractional r)-                                            => Point 2 r :+ p -> NonEmpty (Point 2 r :+ p)-                                            -> DiskResult p r-smallestEnclosingDiskWithPoint p (a :| pts) = foldr addPoint (initial p a) $ L.tails pts-  where-    addPoint []       br   = br-    addPoint (q:pts') br@(DiskResult d _)-      | (q^.core) `inClosedBall` d = br-      | otherwise                  = smallestEnclosingDiskWithPoints p q (a:pts')------ | Smallest enclosing disk, given that p and q should be on it-smallestEnclosingDiskWithPoints     :: (Ord r, Fractional r)-                                    => Point 2 r :+ p -> Point 2 r :+ p -> [Point 2 r :+ p]-                                    -> DiskResult p r-smallestEnclosingDiskWithPoints p q = foldr addPoint (initial p q)-  where-    addPoint r br@(DiskResult d _)-      | (r^.core) `inClosedBall` d = br-      | otherwise                  = DiskResult (circle' r) (Three p q r)--    circle' r = fromMaybe degen $ disk (p^.core) (q^.core) (r^.core)-    degen = error "smallestEnclosingDisk: Unhandled degeneracy, three points on a line"-    -- TODO: handle degenerate case----- | Constructs the initial 'DiskResult' from two points-initial     :: Fractional r => Point 2 r :+ p -> Point 2 r :+ p -> DiskResult p r-initial p q = DiskResult (fromDiameter (p^.core) (q^.core)) (Two p q)

@@ -1,4 +1,3 @@-{-# LANGUAGE DeriveFunctor  #-} {-# LANGUAGE TemplateHaskell  #-} -------------------------------------------------------------------------------- -- |@@ -41,5 +40,5 @@ --    and the points defining it data DiskResult p r = DiskResult { _enclosingDisk  :: Disk () r                                  , _definingPoints :: TwoOrThree (Point 2 r :+ p)-                                 }+                                 } deriving (Show,Eq) makeLenses ''DiskResult

@@ -209,7 +209,7 @@ -- | Assigns the points to a level. Returns the list of levels used. The first -- level in the list is the level assigned to the rest of the nodes. Their -- level is actually still set to Nothing in the underlying array.-assignLevels                  :: (Fractional r, Ord r, Arity d, KnownNat d+assignLevels                  :: (Fractional r, Ord r, Arity d                                  , Show r, Show p                                  )                               => Int -- ^ Number of items we need to collect@@ -362,9 +362,6 @@ -------------------------------------------------------------------------------- -- * Finding Well Separated Pairs --- type AlwaysTrueWSPD d = ( Arity d, KnownNat d---                         , AlwaysTruePFT d, AlwaysTrueTransformation d)- findPairs                     :: (Floating r, Ord r, Arity d, Arity (d + 1))                               => r -> SplitTree d p r a -> SplitTree d p r a                               -> [WSP d p r a]@@ -436,7 +433,7 @@   -- | Computes the maximum width of a splitTree-maxWidth                             :: (Arity d, KnownNat d, Num r)+maxWidth                             :: (Arity d, Num r)                                      => SplitTree d p r a -> r maxWidth (Leaf _)                    = 0 maxWidth (Node _ (NodeData i b _) _) = fromJust $ widthIn' i b

@@ -84,11 +84,11 @@  -- | Test if a point lies strictly inside a ball ----- >>> (point2 0.5 0.0) `insideBall` unitBall+-- >>> (Point2 0.5 0.0) `insideBall` unitBall -- True--- >>> (point2 1 0) `insideBall` unitBall+-- >>> (Point2 1 0) `insideBall` unitBall -- False--- >>> (point2 2 0) `insideBall` unitBall+-- >>> (Point2 2 0) `insideBall` unitBall -- False insideBall       :: (Arity d, Ord r, Num r)                  => Point d r -> Ball d p r -> Bool@@ -104,9 +104,9 @@  -- | Test if a point lies on the boundary of a ball. ----- >>> (point2 1 0) `onBall` unitBall+-- >>> (Point2 1 0) `onBall` unitBall -- True--- >>> (point3 1 1 0) `onBall` unitBall+-- >>> (Point3 1 1 0) `onBall` unitBall -- False onBall       :: (Arity d, Ord r, Num r)              => Point d r -> Ball d p r -> Bool@@ -145,7 +145,7 @@ -- | Given three points, get the disk through the three points. If the three -- input points are colinear we return Nothing ----- >>> disk (point2 0 10) (point2 10 0) (point2 (-10) 0)+-- >>> disk (Point2 0 10) (Point2 10 0) (Point2 (-10) 0) -- Just (Ball {_center = Point2 [0.0,0.0] :+ (), _squaredRadius = 100.0}) disk       :: (Eq r, Fractional r)            => Point 2 r -> Point 2 r -> Point 2 r -> Maybe (Disk () r)@@ -177,7 +177,7 @@     ynom   = det33 $ V3 (fy px py) (fy qx qy) (fy sx sy)      denom  = (2 *) . det33 $ V3 (V3 px py 1) (V3 qx qy 1) (V3 sx sy 1)-    c      = point2 (xnom / denom) (ynom / denom)+    c      = Point2 (xnom / denom) (ynom / denom)   

@@ -25,6 +25,7 @@ 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)@@ -171,9 +172,9 @@  -- | Check if a point lies a box ----- >>> origin `inBox` (boundingBoxList' [point3 1 2 3, point3 10 20 30] :: Box 3 () Int)+-- >>> origin `inBox` (boundingBoxList' [Point3 1 2 3, Point3 10 20 30] :: Box 3 () Int) -- False--- >>> origin `inBox` (boundingBoxList' [point3 (-1) (-2) (-3), point3 10 20 30] :: Box 3 () Int)+-- >>> origin `inBox` (boundingBoxList' [Point3 (-1) (-2) (-3), Point3 10 20 30] :: Box 3 () Int) -- True inBox :: (Arity d, Ord r) => Point d r -> Box d p r -> Bool p `inBox` b = FV.and . FV.zipWith R.inRange (toVec p) . extent $ b@@ -182,7 +183,7 @@ -- resulting vector is 0 indexed whereas one would normally count dimensions -- starting at zero. ----- >>> extent (boundingBoxList' [point3 1 2 3, point3 10 20 30] :: Box 3 () Int)+-- >>> extent (boundingBoxList' [Point3 1 2 3, Point3 10 20 30] :: Box 3 () Int) -- Vector3 [Range (Closed 1) (Closed 10),Range (Closed 2) (Closed 20),Range (Closed 3) (Closed 30)] extent                                 :: Arity d                                        => Box d p r -> Vector d (R.Range r)@@ -191,16 +192,16 @@ -- | Get the size of the box (in all dimensions). Note that the resulting vector is 0 indexed -- whereas one would normally count dimensions starting at zero. ----- >>> size (boundingBoxList' [origin, point3 1 2 3] :: Box 3 () Int)+-- >>> size (boundingBoxList' [origin, Point3 1 2 3] :: Box 3 () Int) -- Vector3 [1,2,3] size :: (Arity d, Num r) => Box d p r -> Vector d r size = fmap R.width . extent  -- | Given a dimension, get the width of the box in that dimension. Dimensions are 1 indexed. ----- >>> widthIn (C :: C 1) (boundingBoxList' [origin, point3 1 2 3] :: Box 3 () Int)+-- >>> widthIn (C :: C 1) (boundingBoxList' [origin, Point3 1 2 3] :: Box 3 () Int) -- 1--- >>> widthIn (C :: C 3) (boundingBoxList' [origin, point3 1 2 3] :: Box 3 () Int)+-- >>> widthIn (C :: C 3) (boundingBoxList' [origin, Point3 1 2 3] :: Box 3 () Int) -- 3 widthIn   :: forall proxy p i d r. (Arity d, Arity (i - 1), Num r, ((i-1)+1) <= d)           => proxy i -> Box d p r -> r@@ -209,13 +210,13 @@  -- | Same as 'widthIn' but with a runtime int instead of a static dimension. ----- >>> widthIn' 1 (boundingBoxList' [origin, point3 1 2 3] :: Box 3 () Int)+-- >>> widthIn' 1 (boundingBoxList' [origin, Point3 1 2 3] :: Box 3 () Int) -- Just 1--- >>> widthIn' 3 (boundingBoxList' [origin, point3 1 2 3] :: Box 3 () Int)+-- >>> widthIn' 3 (boundingBoxList' [origin, Point3 1 2 3] :: Box 3 () Int) -- Just 3--- >>> widthIn' 10 (boundingBoxList' [origin, point3 1 2 3] :: Box 3 () Int)+-- >>> widthIn' 10 (boundingBoxList' [origin, Point3 1 2 3] :: Box 3 () Int) -- Nothing-widthIn'   :: (Arity d, KnownNat d, Num r) => Int -> Box d p r -> Maybe r+widthIn'   :: (Arity d, Num r) => Int -> Box d p r -> Maybe r widthIn' i = preview (V.element' (i-1)) . size  @@ -224,14 +225,14 @@  type Rectangle = Box 2 --- >>> width (boundingBoxList' [origin, point2 1 2] :: Rectangle () Int)+-- >>> width (boundingBoxList' [origin, Point2 1 2] :: Rectangle () Int) -- 1 -- >>> width (boundingBoxList' [origin] :: Rectangle () Int) -- 0 width :: Num r => Rectangle p r -> r width = widthIn (C :: C 1) --- >>> height (boundingBoxList' [origin, point2 1 2] :: Rectangle () Int)+-- >>> height (boundingBoxList' [origin, Point2 1 2] :: Rectangle () Int) -- 2 -- >>> height (boundingBoxList' [origin] :: Rectangle () Int) -- 0@@ -270,9 +271,9 @@   -- | Unsafe version of boundingBoxList, that does not check if the list is non-empty-boundingBoxList' :: (IsBoxable g, Ord (NumType g), Arity (Dimension g))-                 => [g] -> Box (Dimension g) () (NumType g)-boundingBoxList' = boundingBoxList . NE.fromList+boundingBoxList' :: (IsBoxable g, Foldable c, Ord (NumType g), Arity (Dimension g))+                 => c g -> Box (Dimension g) () (NumType g)+boundingBoxList' = boundingBoxList . NE.fromList . F.toList  ---------------------------------------- 

@@ -13,7 +13,7 @@ -- | Returns Nothing if the input line is vertical -- Maps a line l: y = ax + b to a point (a,-b) dualPoint   :: (Fractional r, Eq r) => Line 2 r -> Maybe (Point 2 r)-dualPoint l = (\(a,b) -> point2 a (-b)) <$> toLinearFunction l+dualPoint l = (\(a,b) -> Point2 a (-b)) <$> toLinearFunction l  -- | Pre: the input line is not vertical dualPoint' :: (Fractional r, Eq r) => Line 2 r -> Point 2 r

@@ -68,10 +68,10 @@ lineThrough p q = Line p (q .-. p)  verticalLine   :: Num r => r -> Line 2 r-verticalLine x = Line (point2 x 0) (Vector2 0 1)+verticalLine x = Line (Point2 x 0) (Vector2 0 1)  horizontalLine   :: Num r => r -> Line 2 r-horizontalLine y = Line (point2 0 y) (Vector2 1 0)+horizontalLine y = Line (Point2 0 y) (Vector2 1 0)  -- | Given a line l with anchor point p and vector v, get the line -- perpendicular to l that also goes through p. The resulting line m is@@ -94,9 +94,9 @@  -- | Test if the two lines are parallel. ----- >>> lineThrough origin (point2 1 0) `isParallelTo` lineThrough (point2 1 1) (point2 2 1)+-- >>> lineThrough origin (Point2 1 0) `isParallelTo` lineThrough (Point2 1 1) (Point2 2 1) -- True--- >>> lineThrough origin (point2 1 0) `isParallelTo` lineThrough (point2 1 1) (point2 2 2)+-- >>> lineThrough origin (Point2 1 0) `isParallelTo` lineThrough (Point2 1 1) (Point2 2 2) -- False isParallelTo                         :: (Eq r, Fractional r, Arity d)                                      => Line d r -> Line d r -> Bool@@ -106,11 +106,11 @@  -- | Test if point p lies on line l ----- >>> origin `onLine` lineThrough origin (point2 1 0)+-- >>> origin `onLine` lineThrough origin (Point2 1 0) -- True--- >>> point2 10 10 `onLine` lineThrough origin (point2 2 2)+-- >>> Point2 10 10 `onLine` lineThrough origin (Point2 2 2) -- True--- >>> point2 10 5 `onLine` lineThrough origin (point2 2 2)+-- >>> Point2 10 5 `onLine` lineThrough origin (Point2 2 2) -- False onLine                :: (Eq r, Fractional r, Arity d) => Point d r -> Line d r -> Bool p `onLine` (Line q v) = p == q || (p .-. q) `isScalarMultipleOf` v@@ -203,7 +203,7 @@  -- | Create a line from the linear function ax + b fromLinearFunction     :: Num r => r -> r -> Line 2 r-fromLinearFunction a b = Line (point2 0 b) (Vector2 1 a)+fromLinearFunction a b = Line (Point2 0 b) (Vector2 1 a)  -- | get values a,b s.t. the input line is described by y = ax + b. -- returns Nothing if the line is vertical@@ -216,34 +216,56 @@     :& RNil  -- | Result of a side test-data SideTest = Below | On | Above deriving (Show,Read,Eq,Ord)+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 `onSide` (lineThrough origin $ point2 10 5)+-- >>> Point2 10 10 `onSideUpDown` (lineThrough origin $ Point2 10 5) -- Above--- >>> point2 10 10 `onSide` (lineThrough origin $ point2 (-10) 5)+-- >>> Point2 10 10 `onSideUpDown` (lineThrough origin $ Point2 (-10) 5) -- Above--- >>> point2 5 5 `onSide` (verticalLine 10)+-- >>> Point2 5 5 `onSideUpDown` (verticalLine 10) -- Below--- >>> point2 5 5 `onSide` (lineThrough origin $ point2 (-3) (-3))+-- >>> Point2 5 5 `onSideUpDown` (lineThrough origin $ Point2 (-3) (-3)) -- On-onSide                :: (Ord r, Num r) => Point 2 r -> Line 2 r -> SideTest-q `onSide` (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+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 +-- | Result of a side test+data SideTest = LeftSide | OnLine | RightSide 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 `onSide` (lineThrough origin $ Point2 10 5)+-- LeftSide+-- >>> Point2 10 10 `onSide` (lineThrough origin $ Point2 (-10) 5)+-- RightSide+-- >>> Point2 5 5 `onSide` (verticalLine 10)+-- LeftSide+-- >>> Point2 5 5 `onSide` (lineThrough origin $ Point2 (-3) (-3))+-- OnLine+onSide                :: (Ord r, Num r) => Point 2 r -> Line 2 r -> SideTest+q `onSide` (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 p r q of+                          CCW      -> LeftSide+                          CW       -> RightSide+                          CoLinear -> OnLine  -- | Test if the query point q lies (strictly) above line l liesAbove       :: (Ord r, Num r) => Point 2 r -> Line 2 r -> Bool-q `liesAbove` l = q `onSide` l == Above+q `liesAbove` l = q `onSideUpDown` l == Above   -- | Get the bisector between two points

@@ -203,30 +203,30 @@  -- | Test if a point lies on a line segment. ----- >>> (point2 1 0) `onSegment` (ClosedLineSegment (origin :+ ()) (point2 2 0 :+ ()))+-- >>> (Point2 1 0) `onSegment` (ClosedLineSegment (origin :+ ()) (Point2 2 0 :+ ())) -- True--- >>> (point2 1 1) `onSegment` (ClosedLineSegment (origin :+ ()) (point2 2 0 :+ ()))+-- >>> (Point2 1 1) `onSegment` (ClosedLineSegment (origin :+ ()) (Point2 2 0 :+ ())) -- False--- >>> (point2 5 0) `onSegment` (ClosedLineSegment (origin :+ ()) (point2 2 0 :+ ()))+-- >>> (Point2 5 0) `onSegment` (ClosedLineSegment (origin :+ ()) (Point2 2 0 :+ ())) -- False--- >>> (point2 (-1) 0) `onSegment` (ClosedLineSegment (origin :+ ()) (point2 2 0 :+ ()))+-- >>> (Point2 (-1) 0) `onSegment` (ClosedLineSegment (origin :+ ()) (Point2 2 0 :+ ())) -- False--- >>> (point2 1 1) `onSegment` (ClosedLineSegment (origin :+ ()) (point2 3 3 :+ ()))+-- >>> (Point2 1 1) `onSegment` (ClosedLineSegment (origin :+ ()) (Point2 3 3 :+ ())) -- True -- -- Note that the segments are assumed to be closed. So the end points lie on the segment. ----- >>> (point2 2 0) `onSegment` (ClosedLineSegment (origin :+ ()) (point2 2 0 :+ ()))+-- >>> (Point2 2 0) `onSegment` (ClosedLineSegment (origin :+ ()) (Point2 2 0 :+ ())) -- True--- >>> origin `onSegment` (ClosedLineSegment (origin :+ ()) (point2 2 0 :+ ()))+-- >>> origin `onSegment` (ClosedLineSegment (origin :+ ()) (Point2 2 0 :+ ())) -- True -- -- -- This function works for arbitrary dimensons. ----- >>> (point3 1 1 1) `onSegment` (ClosedLineSegment (origin :+ ()) (point3 3 3 3 :+ ()))+-- >>> (Point3 1 1 1) `onSegment` (ClosedLineSegment (origin :+ ()) (Point3 3 3 3 :+ ())) -- True--- >>> (point3 1 2 1) `onSegment` (ClosedLineSegment (origin :+ ()) (point3 3 3 3 :+ ()))+-- >>> (Point3 1 2 1) `onSegment` (ClosedLineSegment (origin :+ ()) (Point3 3 3 3 :+ ())) -- False onSegment       :: (Ord r, Fractional r, Arity d)                 => Point d r -> LineSegment d p r -> Bool@@ -277,7 +277,7 @@                 in (s&start .~ q)&end .~ p  -- testSeg :: LineSegment 2 () Rational--- testSeg = LineSegment (Open $ ext origin)  (Closed $ ext (point2 10 0))+-- testSeg = LineSegment (Open $ ext origin)  (Closed $ ext (Point2 10 0))  -- horL' :: Line 2 Rational -- horL' = horizontalLine 0

@@ -180,7 +180,7 @@     oF@(FaceId (VertexId of')) = PG.leftFace ofD g      mkFaceIdx i | i == 0    = Nothing-                | otherwise = Just (c,mkFaceId i)+                | otherwise = Just (c,mkFaceId . flipID $ i)      -- at index i we are storing the outerface     mkFaceData                 :: Int -> f -> FaceData (Dart s) f

@@ -15,10 +15,13 @@                                             , embedAsHolesIn                                             ) where +import           Algorithms.DivideAndConquer import           Control.Lens hiding (holes)-import           Data.BinaryTree (asBalancedBinLeafTree, foldUp, Elem(..))+import           Data.Ext import           Data.Geometry.PlanarSubdivision.Basic import           Data.Geometry.PlanarSubdivision.Raw+import           Data.Geometry.Point+import           Data.Geometry.Polygon import           Data.PlanarGraph.Dart import           Data.PlaneGraph ( Dart, VertexId(..), FaceId(..)                                 , VertexId', FaceId'@@ -26,14 +29,7 @@ import qualified Data.PlaneGraph as PG import           Data.Semigroup.Foldable import qualified Data.Vector as V---- import Data.Coerce-import           Unsafe.Coerce(unsafeCoerce)---import           Data.Ext-import           Data.Geometry.Point-import           Data.Geometry.Polygon+import           Unsafe.Coerce (unsafeCoerce)  -------------------------------------------------------------------------------- -- * Embedding one subdivision in another one@@ -118,7 +114,7 @@                => (f -> f -> f)                -> t (PlanarSubdivision s v e f r)                -> PlanarSubdivision s v e f r-mergeAllWith f = foldUp (\l _ r -> mergeWith f l r) _unElem . asBalancedBinLeafTree . toNonEmpty+mergeAllWith f = divideAndConquer1With (mergeWith f) id . toNonEmpty  -- | Merge a pair of *disjoint* planar subdivisions, unifying their -- outer face. For the outerface data it simply takes the data of the@@ -153,7 +149,7 @@         -- we have to shift the number of the *Arcs*. Since every dart         -- consists of two arcs, we have to shift by numDarts / 2         -- Furthermore, we take numFaces - 1 since we want the first-        -- *internal* face of p2 (the one with FaceId 1) to correspond with the first free+        -- /internal/ face of p2 (the one with FaceId 1) to correspond with the first free         -- position (at index numFaces)      cs = p1^.components <> p2'^.components

@@ -10,8 +10,32 @@ -- \(d\)-dimensional points. -- ---------------------------------------------------------------------------------module Data.Geometry.Point where+module Data.Geometry.Point( Point(..)+                          , origin, vector+                          , pointFromList +                          , coord , unsafeCoord++                          , projectPoint++                          , pattern Point2+                          , pattern Point3+                          , xCoord, yCoord, zCoord++                          , PointFunctor(..)++                          , CCW(..), ccw, ccw'++                          , ccwCmpAround, cwCmpAround, ccwCmpAroundWith, cwCmpAroundWith+                          , sortAround, insertIntoCyclicOrder++                          , Quadrant(..), quadrantWith, quadrant, partitionIntoQuadrants++                          , cmpByDistanceTo++                          , squaredEuclideanDist, euclideanDist+                          ) where+ import           Control.DeepSeq import           Control.Lens import           Data.Aeson@@ -23,13 +47,14 @@ 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)-import           Test.QuickCheck(Arbitrary)   --------------------------------------------------------------------------------@@ -99,7 +124,7 @@  -- | Lens to access the vector corresponding to this point. ----- >>> (point3 1 2 3) ^. vector+-- >>> (Point3 1 2 3) ^. vector -- Vector3 [1,2,3] -- >>> origin & vector .~ Vector3 1 2 3 -- Point3 [1,2,3]@@ -111,7 +136,7 @@ -- sense that no bounds are checked. Consider using `coord` instead. -- ----- >>> point3 1 2 3 ^. unsafeCoord 2+-- >>> Point3 1 2 3 ^. unsafeCoord 2 -- 2 unsafeCoord   :: Arity d => Int -> Lens' (Point d r) r unsafeCoord i = vector . singular (ix (i-1))@@ -119,11 +144,11 @@  -- | Get the coordinate in a given dimension ----- >>> point3 1 2 3 ^. coord (C :: C 2)+-- >>> Point3 1 2 3 ^. coord (C :: C 2) -- 2--- >>> point3 1 2 3 & coord (C :: C 1) .~ 10+-- >>> 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 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@@ -160,15 +185,13 @@ --   let --     f              :: Point 2 r -> r --     f (Point2 x y) = x---   in f (point2 1 2)+--   in f (Point2 1 2) -- :} -- 1 -- -- if we want. pattern Point2       :: r -> r -> Point 2 r-pattern Point2 x y   <- (_point2 -> (x,y))-  where-    Point2 x y = point2 x y+pattern Point2 x y = Point (Vector2 x y) {-# COMPLETE Point2 #-}  -- | Similarly, we can write:@@ -181,47 +204,14 @@ -- :} -- 3 pattern Point3       :: r -> r -> r -> Point 3 r-pattern Point3 x y z <- (_point3 -> (x,y,z))-  where-    Point3 x y z = point3 x y z+pattern Point3 x y z = (Point (Vector3 x y z)) {-# COMPLETE Point3 #-} --- | Construct a 2 dimensional point------ >>> point2 1 2--- Point2 [1,2]-point2     :: r -> r -> Point 2 r-point2 x y = Point $ Vector2 x y---- | Destruct a 2 dimensional point------ >>> _point2 $ point2 1 2--- (1,2)-_point2 :: Point 2 r -> (r,r)-_point2 = (\(Vector2 x y) -> (x,y)) . toVec------ | Construct a 3 dimensional point------ >>> point3 1 2 3--- Point3 [1,2,3]-point3       :: r -> r -> r -> Point 3 r-point3 x y z = Point $ Vector3 x y z---- | Destruct a 3 dimensional point------ >>> _point3 $ point3 1 2 3--- (1,2,3)-_point3 :: Point 3 r -> (r,r,r)-_point3 = (\(Vector3 x y z) -> (x,y,z)) . toVec-- -- | Shorthand to access the first coordinate C 1 ----- >>> point3 1 2 3 ^. xCoord+-- >>> Point3 1 2 3 ^. xCoord -- 1--- >>> point2 1 2 & xCoord .~ 10+-- >>> Point2 1 2 & xCoord .~ 10 -- Point2 [10,2] xCoord :: (1 <= d, Arity d) => Lens' (Point d r) r xCoord = coord (C :: C 1)@@ -229,9 +219,9 @@  -- | Shorthand to access the second coordinate C 2 ----- >>> point2 1 2 ^. yCoord+-- >>> Point2 1 2 ^. yCoord -- 2--- >>> point3 1 2 3 & yCoord %~ (+1)+-- >>> 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)@@ -239,9 +229,9 @@  -- | Shorthand to access the third coordinate C 3 ----- >>> point3 1 2 3 ^. zCoord+-- >>> Point3 1 2 3 ^. zCoord -- 3--- >>> point3 1 2 3 & zCoord %~ (+1)+-- >>> 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)@@ -287,9 +277,9 @@ -- 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)-sortArround   :: (Ord r, Num r)-               => Point 2 r :+ q -> [Point 2 r :+ p] -> [Point 2 r :+ p]-sortArround c = L.sortBy (ccwCmpAround c)+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@@ -336,38 +326,85 @@      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.--- Points nearer to the center come before--- points further away.-ccwCmpAround       :: (Num r, Ord r)-                   => Point 2 r :+ qc -> Point 2 r :+ p -> Point 2 r :+ q -> Ordering-ccwCmpAround c q r = case (quadrantWith c q `compare` quadrantWith c r) of-                       EQ -> case ccw (c^.core) (q^.core) (r^.core) of-                         CCW      -> LT-                         CW       -> GT-                         CoLinear -> qdA (c^.core) (q^.core)-                                     `compare`-                                     qdA (c^.core) (r^.core)-                       x -> x -- if the quadrant differs, use the order-                              -- specified by the quadrant.+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. Points nearer to the center come before--- points further away.-cwCmpAround       :: (Num r, Ord r)-                  => Point 2 r :+ qc -> Point 2 r :+ p -> Point 2 r :+ q -> Ordering-cwCmpAround c q r = case (quadrantWith c q `compare` quadrantWith c r) of-                       EQ -> case ccw (c^.core) (q^.core) (r^.core) of-                         CCW      -> GT-                         CW       -> LT-                         CoLinear -> qdA (c^.core) (q^.core)-                                     `compare`-                                     qdA (c^.core) (r^.core)-                       LT -> GT-                       GT -> LT -- if the quadrant differs, use the order-                                -- specified by the quadrant.-+-- 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@@ -378,7 +415,7 @@ 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)+insertIntoCyclicOrder c = CU.insertOrdBy (ccwCmpAround c <> cmpByDistanceTo c)   -- | Squared Euclidean distance between two points

@@ -4,6 +4,7 @@ module Data.Geometry.PolyLine where  import           Control.Lens+import           Data.Aeson import           Data.Bifunctor import           Data.Ext import qualified Data.Foldable as F@@ -16,13 +17,14 @@ import           Data.LSeq (LSeq, pattern (:<|)) import qualified Data.LSeq as LSeq import qualified Data.List.NonEmpty as NE+import           GHC.Generics(Generic) import           GHC.TypeLits  -------------------------------------------------------------------------------- -- * d-dimensional Polygonal Lines (PolyLines)  -- | A Poly line in R^d has at least 2 vertices-newtype PolyLine d p r = PolyLine { _points :: LSeq 2 (Point d r :+ p) }+newtype PolyLine d p r = PolyLine { _points :: LSeq 2 (Point d r :+ p) } deriving (Generic) makeLenses ''PolyLine  deriving instance (Show r, Show p, Arity d) => Show    (PolyLine d p r)@@ -50,6 +52,9 @@ instance Arity d => Bifunctor (PolyLine d) where   bimap f g (PolyLine pts) = PolyLine $ fmap (bimap (fmap 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)  -- | pre: The input list contains at least two points fromPoints :: [Point d r :+ p] -> PolyLine d p r

@@ -1,3 +1,4 @@+{-# LANGUAGE TemplateHaskell #-} -------------------------------------------------------------------------------- -- | -- Module      :  Data.Geometry.Polygon@@ -8,514 +9,58 @@ -- A Polygon data type and some basic functions to interact with them. -- ---------------------------------------------------------------------------------module Data.Geometry.Polygon where--import           Algorithms.Geometry.LinearProgramming.LP2DRIC-import           Algorithms.Geometry.LinearProgramming.Types-import           Control.DeepSeq-import           Control.Lens hiding (Simple)-import           Control.Monad.Random.Class-import           Data.Bifoldable-import           Data.Bifunctor-import           Data.Bitraversable-import qualified Data.CircularSeq as C-import           Data.Ext-import qualified Data.Foldable as F-import           Data.Geometry.Boundary-import           Data.Geometry.Box-import           Data.Geometry.Line-import           Data.Geometry.HalfSpace(rightOf)-import           Data.Geometry.LineSegment-import           Data.Geometry.Point-import           Data.Geometry.Properties-import           Data.Geometry.Transformation-import           Data.Geometry.Triangle (Triangle(..), inTriangle)-import           Data.Geometry.Vector-import qualified Data.List as List-import           Data.List.NonEmpty (NonEmpty(..))-import qualified Data.List.NonEmpty as NonEmpty-import           Data.Maybe (mapMaybe, catMaybes)-import           Data.Ord (comparing)-import           Data.Semigroup (sconcat)-import           Data.Semigroup.Foldable-import qualified Data.Sequence as Seq-import           Data.Util-import           Data.Vinyl.CoRec (asA)-------------------------------------------------------------------------------------- * Polygons--{- $setup->>> :{--- import qualified Data.CircularSeq as C-let simplePoly :: SimplePolygon () Rational-    simplePoly = SimplePolygon . C.fromList . map ext $ [ point2 0 0-                                                        , point2 10 0-                                                        , point2 10 10-                                                        , point2 5 15-                                                        , point2 1 11-                                                        ]-:} -}---- | We distinguish between simple polygons (without holes) and Polygons with holes.-data PolygonType = Simple | Multi---data Polygon (t :: PolygonType) p r where-  SimplePolygon :: C.CSeq (Point 2 r :+ p)                         -> Polygon Simple p r-  MultiPolygon  :: C.CSeq (Point 2 r :+ p) -> [Polygon Simple p r] -> Polygon Multi  p r--instance Bifunctor (Polygon t) where-  bimap = bimapDefault--instance Bifoldable (Polygon t) where-  bifoldMap = bifoldMapDefault--instance Bitraversable (Polygon t) where-  bitraverse f g p = case p of-    SimplePolygon vs   -> SimplePolygon <$> bitraverseVertices f g vs-    MultiPolygon vs hs -> MultiPolygon  <$> bitraverseVertices f g vs-                                        <*> traverse (bitraverse f g) hs--instance (NFData p, NFData r) => NFData (Polygon t p r) where-  rnf (SimplePolygon vs)   = rnf vs-  rnf (MultiPolygon vs hs) = rnf (vs,hs)--bitraverseVertices     :: (Applicative f, Traversable t) => (p -> f q) -> (r -> f s)-                  -> t (Point 2 r :+ p) -> f (t (Point 2 s :+ q))-bitraverseVertices f g = traverse (bitraverse (traverse g) f)--type SimplePolygon = Polygon Simple--type MultiPolygon  = Polygon Multi---- | Either a simple or multipolygon-type SomePolygon p r = Either (Polygon Simple p r) (Polygon Multi p r)--type instance Dimension (SomePolygon p r) = 2-type instance NumType   (SomePolygon p r) = r---- | Polygons are per definition 2 dimensional-type instance Dimension (Polygon t p r) = 2-type instance NumType   (Polygon t p r) = r--instance (Show p, Show r) => Show (Polygon t p r) where-  show (SimplePolygon vs)   = "SimplePolygon " <> show vs-  show (MultiPolygon vs hs) = "MultiPolygon " <> show vs <> " " <> show hs--instance (Eq p, Eq r) => Eq (Polygon t p r) where-  (SimplePolygon vs)   == (SimplePolygon vs')    = vs == vs'-  (MultiPolygon vs hs) == (MultiPolygon vs' hs') = vs == vs' && hs == hs'--instance PointFunctor (Polygon t p) where-  pmap f (SimplePolygon vs)   = SimplePolygon (fmap (first f) vs)-  pmap f (MultiPolygon vs hs) = MultiPolygon  (fmap (first f) vs) (map (pmap f) hs)--instance Fractional r => IsTransformable (Polygon t p r) where-  transformBy = transformPointFunctor--instance IsBoxable (Polygon t p r) where-  boundingBox = boundingBoxList' . toListOf (outerBoundary.traverse.core)--type instance IntersectionOf (Line 2 r) (Boundary (Polygon t p r)) =-  '[Seq.Seq (Either (Point 2 r) (LineSegment 2 () r))]--type instance IntersectionOf (Point 2 r) (Polygon t p r) = [NoIntersection, Point 2 r]--instance (Fractional r, Ord r) => (Point 2 r) `IsIntersectableWith` (Polygon t p r) where-  nonEmptyIntersection = defaultNonEmptyIntersection-  q `intersects` pg = q `inPolygon` pg /= Outside-  q `intersect` pg | q `intersects` pg = coRec q-                   | otherwise         = coRec NoIntersection---- instance IsIntersectableWith (Line 2 r) (Boundary (Polygon t p r)) where---   nonEmptyIntersection _ _ (CoRec xs) = null xs---   l `intersect` (Boundary (SimplePolygon vs)) =---     undefined-  -- l `intersect` (Boundary (MultiPolygon vs hs)) = coRec .-  --    Seq.sortBy f . Seq.fromList-  --     . concatMap (unpack . (l `intersect`) . Boundary)-  --     $ SimplePolygon vs : hs-  --   where-  --     unpack (CoRec x) = x-  --     f = undefined------- * Functions on Polygons--outerBoundary :: forall t p r. Lens' (Polygon t p r) (C.CSeq (Point 2 r :+ p))-outerBoundary = lens g s-  where-    g                     :: Polygon t p r -> C.CSeq (Point 2 r :+ p)-    g (SimplePolygon vs)  = vs-    g (MultiPolygon vs _) = vs--    s                           :: Polygon t p r -> C.CSeq (Point 2 r :+ p)-                                -> Polygon t p r-    s (SimplePolygon _)      vs = SimplePolygon vs-    s (MultiPolygon  _   hs) vs = MultiPolygon vs hs--polygonHoles :: forall p r. Lens' (Polygon Multi p r) [Polygon Simple p r]-polygonHoles = lens g s-  where-    g                     :: Polygon Multi p r -> [Polygon Simple p r]-    g (MultiPolygon _ hs) = hs-    s                     :: Polygon Multi p r -> [Polygon Simple p r]-                          -> Polygon Multi p r-    s (MultiPolygon vs _) = MultiPolygon vs----- | Access the i^th vertex on the outer boundary-outerVertex   :: Int -> Lens' (Polygon t p r) (Point 2 r :+ p)-outerVertex i = outerBoundary.C.item i---- running time: \(O(\log i)\)-outerBoundaryEdge     :: Int -> Polygon t p r -> LineSegment 2 p r-outerBoundaryEdge i p = let u = p^.outerVertex i-                            v = p^.outerVertex (i+1)-                        in LineSegment (Closed u) (Open v)----- | Get all holes in a polygon-holeList                     :: Polygon t p r -> [Polygon Simple p r]-holeList (SimplePolygon _)   = []-holeList (MultiPolygon _ hs) = hs----- | The vertices in the polygon. No guarantees are given on the order in which--- they appear!-polygonVertices                      :: Polygon t p r-                                     -> NonEmpty.NonEmpty (Point 2 r :+ p)-polygonVertices (SimplePolygon vs)   = toNonEmpty vs-polygonVertices (MultiPolygon vs hs) =-  sconcat $ toNonEmpty vs NonEmpty.:| map polygonVertices hs----- | Creates a simple polygon from the given list of vertices.------ pre: the input list constains no repeated vertices.-fromPoints :: [Point 2 r :+ p] -> SimplePolygon p r-fromPoints = SimplePolygon . C.fromList----- | The edges along the outer boundary of the polygon. The edges are half open.------ running time: \(O(n)\)-outerBoundaryEdges :: Polygon t p r -> C.CSeq (LineSegment 2 p r)-outerBoundaryEdges = toEdges . (^.outerBoundary)---- | Lists all edges. The edges on the outer boundary are given before the ones--- on the holes. However, no other guarantees are given on the order.------ running time: \(O(n)\)-listEdges    :: Polygon t p r -> [LineSegment 2 p r]-listEdges pg = let f = F.toList . outerBoundaryEdges-               in  f pg <> concatMap f (holeList pg)---- | Pairs every vertex with its incident edges. The first one is its--- predecessor edge, the second one its successor edge.------ >>> 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] :+ ()))-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)-withIncidentEdges (MultiPolygon vs hs) = MultiPolygon vs' hs'-  where-    (SimplePolygon vs') = withIncidentEdges $ SimplePolygon vs-    hs' = map withIncidentEdges hs---- -- | Gets the i^th edge on the outer boundary of the polygon, that is the edge----- with vertices i and i+1 with respect to the current focus. All indices--- -- modulo n.--- ------ | Given the vertices of the polygon. Produce a list of edges. The edges are--- half-open.-toEdges    :: C.CSeq (Point 2 r :+ p) -> C.CSeq (LineSegment 2 p r)-toEdges vs = C.zipLWith (\p q -> LineSegment (Closed p) (Open q)) vs (C.rotateR vs)-  -- let vs' = F.toList vs in-  -- C.fromList $ zipWith (\p q -> LineSegment (Closed p) (Open q)) vs' (tail vs' ++ vs')----- | Test if q lies on the boundary of the polygon. Running time: O(n)------ >>> point2 1 1 `onBoundary` simplePoly--- False--- >>> point2 0 0 `onBoundary` simplePoly--- True--- >>> point2 10 0 `onBoundary` simplePoly--- True--- >>> point2 5 13 `onBoundary` simplePoly--- False--- >>> point2 5 10 `onBoundary` simplePoly--- False--- >>> point2 10 5 `onBoundary` simplePoly--- True--- >>> point2 20 5 `onBoundary` simplePoly--- False------ TODO: testcases multipolygon-onBoundary        :: (Fractional r, Ord r) => Point 2 r -> Polygon t p r -> Bool-q `onBoundary` pg = any (q `onSegment`) es-  where-    out = SimplePolygon $ pg^.outerBoundary-    es = concatMap (F.toList . outerBoundaryEdges) $ out : holeList pg---- | Check if a point lies inside a polygon, on the boundary, or outside of the polygon.--- Running time: O(n).------ >>> point2 1 1 `inPolygon` simplePoly--- Inside--- >>> point2 0 0 `inPolygon` simplePoly--- OnBoundary--- >>> point2 10 0 `inPolygon` simplePoly--- OnBoundary--- >>> point2 5 13 `inPolygon` simplePoly--- Inside--- >>> point2 5 10 `inPolygon` simplePoly--- Inside--- >>> point2 10 5 `inPolygon` simplePoly--- OnBoundary--- >>> point2 20 5 `inPolygon` simplePoly--- Outside------ TODO: Add some testcases with multiPolygons--- TODO: Add some more onBoundary testcases-inPolygon                                :: forall t p r. (Fractional r, Ord r)-                                         => Point 2 r -> Polygon t p r-                                         -> PointLocationResult-q `inPolygon` pg-    | q `onBoundary` pg                             = OnBoundary-    | odd kl && odd kr && not (any (q `inHole`) hs) = Inside-    | otherwise                                     = Outside-  where-    l = horizontalLine $ q^.yCoord--    -- Given a line segment, compute the intersection point (if a point) with the-    -- line l-    intersectionPoint = asA @(Point 2 r) . (`intersect` l)--    -- Count the number of intersections that the horizontal line through q-    -- maxes with the polygon, that are strictly to the left and strictly to-    -- the right of q. If these numbers are both odd the point lies within the polygon.-    ---    ---    -- note that: - by the asA (Point 2 r) we ignore horizontal segments (as desired)-    --            - by the filtering, we effectively limit l to an open-half line, starting-    --               at the (open) point q.-    --            - by using half-open segments as edges we avoid double counting-    --               intersections that coincide with vertices.-    --            - If the point is outside, and on the same height as the-    --              minimum or maximum coordinate of the polygon. The number of-    --              intersections to the left or right may be one. Thus-    --              incorrectly classifying the point as inside. To avoid this,-    --              we count both the points to the left *and* to the right of-    --              p. Only if both are odd the point is inside.  so that if-    --              the point is outside, and on the same y-coordinate as one-    --              of the extermal vertices (one ofth)-    ---    -- See http://geomalgorithms.com/a03-_inclusion.html for more information.-    SP kl kr = count (\p -> (p^.xCoord) `compare` (q^.xCoord))-             . mapMaybe intersectionPoint . F.toList . outerBoundaryEdges $ pg--    -- For multi polygons we have to test if we do not lie in a hole .-    inHole = insidePolygon-    hs     = holeList pg--    count   :: (a -> Ordering) -> [a] -> SP Int Int-    count f = foldr (\x (SP lts gts) -> case f x of-                             LT -> SP (lts + 1) gts-                             EQ -> SP lts       gts-                             GT -> SP lts       (gts + 1)) (SP 0 0)----- | Test if a point lies strictly inside the polgyon.-insidePolygon        :: (Fractional r, Ord r) => Point 2 r -> Polygon t p r -> Bool-q `insidePolygon` pg = q `inPolygon` pg == Inside----- testQ = map (`inPolygon` testPoly) [ point2 1 1    -- Inside---                                    , point2 0 0    -- OnBoundary---                                    , point2 5 14   -- Inside---                                    , point2 5 10   -- Inside---                                    , point2 10 5   -- OnBoundary---                                    , point2 20 5   -- Outside---                                    ]---- testPoly :: SimplePolygon () Rational--- testPoly = SimplePolygon . C.fromList . map ext $ [ point2 0 0---                                                   , point2 10 0---                                                   , point2 10 10---                                                   , point2 5 15---                                                   , point2 1 11---                                                   ]---- | Compute the area of a polygon-area                        :: Fractional r => Polygon t p r -> r-area poly@(SimplePolygon _) = abs $ signedArea poly-area (MultiPolygon vs hs)   = area (SimplePolygon vs) - sum [area h | h <- hs]----- | Compute the signed area of a simple polygon. The the vertices are in--- clockwise order, the signed area will be negative, if the verices are given--- in counter clockwise order, the area will be positive.-signedArea      :: Fractional r => SimplePolygon p r -> r-signedArea poly = x / 2-  where-    x = sum [ p^.core.xCoord * q^.core.yCoord - q^.core.xCoord * p^.core.yCoord-            | LineSegment' p q <- F.toList $ outerBoundaryEdges poly  ]----- | Compute the centroid of a simple polygon.-centroid      :: Fractional r => SimplePolygon p r -> Point 2 r-centroid poly = Point $ sum' xs ^/ (6 * signedArea poly)-  where-    xs = [ (toVec p ^+^ toVec q) ^* (p^.xCoord * q^.yCoord - q^.xCoord * p^.yCoord)-         | LineSegment' (p :+ _) (q :+ _) <- F.toList $ outerBoundaryEdges poly  ]--    sum' = F.foldl' (^+^) zero----- | Pick a  point that is inside the polygon.------ (note: if the polygon is degenerate; i.e. has <3 vertices, we report a--- vertex of the polygon instead.)------ pre: the polygon is given in CCW order------ running time: \(O(n)\)-pickPoint    :: (Ord r, Fractional r) => Polygon p t r -> Point 2 r-pickPoint pg | isTriangle pg = centroid . SimplePolygon $ pg^.outerBoundary-             | otherwise     = let LineSegment' (p :+ _) (q :+ _) = findDiagonal pg-                               in p .+^ (0.5 *^ (q .-. p))---- | Test if the polygon is a triangle------ running time: \(O(1)\)-isTriangle :: Polygon p t r -> Bool-isTriangle = \case-    SimplePolygon vs   -> go vs-    MultiPolygon vs [] -> go vs-    MultiPolygon _  _  -> False-  where-    go vs = case toNonEmpty vs of-              (_ :| [_,_]) -> True-              _            -> False---- | Find a diagonal of the polygon.------ pre: the polygon is given in CCW order------ running time: \(O(n)\)-findDiagonal    :: (Ord r, Fractional r) => Polygon t p r -> LineSegment 2 p r-findDiagonal pg = List.head . catMaybes . F.toList $ diags-     -- note that a diagonal is guaranteed to exist, so the usage of head is safe.-  where-    vs      = pg^.outerBoundary-    diags   = C.zip3LWith f (C.rotateL vs) vs (C.rotateR vs)-    f u v w = case ccw (u^.core) (v^.core) (w^.core) of-                CCW      -> Just $ findDiag u v w-                            -- v is a convex vertex, so find a diagonal-                            -- (either uw) or from v to a point inside the-                            -- triangle-                CW       -> Nothing -- v is a reflex vertex-                CoLinear -> Nothing -- colinear vertex!?--    -- we test if uw is a diagonal by figuring out if there is a vertex-    -- strictly inside the triangle t. If there is no such vertex then uw must-    -- be a diagonal (i.e. uw intersects the polygon boundary iff there is a-    -- vtx inside t).  If there are vertices inside the triangle, we find the-    -- one z furthest from the line(segment) uw. It then follows that vz is a-    -- diagonal. Indeed this is pretty much the argument used to prove that any-    -- polygon can be triangulated. See BKOS Chapter 3 for details.-    findDiag u v w = let t  = Triangle u v w-                         uw = ClosedLineSegment u w-                     in maybe uw (ClosedLineSegment v)-                      . safeMaximumOn (distTo $ supportingLine uw)-                      . filter (\(z :+ _) -> z `inTriangle` t == Inside)-                      . F.toList . polygonVertices-                      $ pg--    distTo l (z :+ _) = sqDistanceTo z l-+module Data.Geometry.Polygon( PolygonType(..)+                            , Polygon(..)+                            , _SimplePolygon, _MultiPolygon+                            , SimplePolygon, MultiPolygon, SomePolygon -safeMaximumOn   :: Ord b => (a -> b) -> [a] -> Maybe a-safeMaximumOn f = \case-  [] -> Nothing-  xs -> Just $ List.maximumBy (comparing f) xs+                            , fromPoints +                            , polygonVertices, listEdges --- | Test if the outer boundary of the polygon is in clockwise or counter--- clockwise order.------ running time: \(O(n)\)----isCounterClockwise :: (Eq r, Fractional r) => Polygon t p r -> Bool-isCounterClockwise = (\x -> x == abs x) . signedArea-                   . fromPoints . F.toList . (^.outerBoundary)+                            , outerBoundary, outerBoundaryEdges+                            , outerVertex, outerBoundaryEdge +                            , polygonHoles, polygonHoles'+                            , holeList --- | Orient the outer boundary to clockwise order-toClockwiseOrder         :: (Eq r, Fractional r) => Polygon t p r -> Polygon t p r-toClockwiseOrder p-  | isCounterClockwise p = reverseOuterBoundary p-  | otherwise            = p+                            , inPolygon, insidePolygon, onBoundary --- | Orient the outer boundary to counter clockwise order-toCounterClockWiseOrder    :: (Eq r, Fractional r) => Polygon t p r -> Polygon t p r-toCounterClockWiseOrder p-  | not $ isCounterClockwise p = reverseOuterBoundary p-  | otherwise                  = p+                            , area, signedArea -reverseOuterBoundary   :: Polygon t p r -> Polygon t p r-reverseOuterBoundary p = p&outerBoundary %~ C.reverseDirection+                            , centroid+                            , pickPoint +                            , isTriangle, isStarShaped --- | Convert a Polygon to a simple polygon by forgetting about any holes.-asSimplePolygon                        :: Polygon t p r -> SimplePolygon p r-asSimplePolygon poly@(SimplePolygon _) = poly-asSimplePolygon (MultiPolygon vs _)    = SimplePolygon vs+                            , isCounterClockwise+                            , toCounterClockWiseOrder, toCounterClockWiseOrder'+                            , toClockwiseOrder, toClockwiseOrder'+                            , reverseOuterBoundary +                            , findDiagonal --- | Comparison that compares which point is 'larger' in the direction given by--- the vector u.-cmpExtreme       :: (Num r, Ord r)-                 => Vector 2 r -> Point 2 r :+ p -> Point 2 r :+ q -> Ordering-cmpExtreme u p q = u `dot` (p^.core .-. q^.core) `compare` 0+                            , withIncidentEdges, numberVertices +                            , asSimplePolygon+                            , extremesLinear, cmpExtreme+                            ) where --- | Finds the extreme points, minimum and maximum, in a given direction------ running time: \(O(n)\)-extremesLinear     :: (Ord r, Num r) => Vector 2 r -> Polygon t p r-                   -> (Point 2 r :+ p, Point 2 r :+ p)-extremesLinear u p = let vs = p^.outerBoundary-                         f  = cmpExtreme u-                     in (F.minimumBy f vs, F.maximumBy f vs)+import           Algorithms.Geometry.LinearProgramming.LP2DRIC+import           Algorithms.Geometry.LinearProgramming.Types+import           Control.Lens hiding (Simple)+import           Control.Monad.Random.Class+import           Data.Ext+import qualified Data.Foldable as F+import           Data.Geometry.HalfSpace (rightOf)+import           Data.Geometry.Line+import           Data.Geometry.Point+import           Data.Geometry.Polygon.Core+import           Data.Geometry.Polygon.Extremes  --- | assigns unique integer numbers to all vertices. Numbers start from 0, and--- are increasing along the outer boundary. The vertices of holes--- will be numbered last, in the same order.------ >>> numberVertices simplePoly--- SimplePolygon 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 ()]-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'- --------------------------------------------------------------------------------+-- * Polygons  -- | Test if a Simple polygon is star-shaped. Returns a point in the kernel -- (i.e. from which the entire polygon is visible), if it exists.

@@ -12,8 +12,8 @@ -------------------------------------------------------------------------------- module Data.Geometry.Polygon.Convex( ConvexPolygon(..), simplePolygon                                    , merge-                                   , lowerTangent, upperTangent-                                   , isLeftOf, isRightOf+                                   , lowerTangent, lowerTangent'+                                   , upperTangent, upperTangent'                                     , extremes                                    , maxInDirection@@ -26,22 +26,27 @@  import           Control.DeepSeq import           Control.Lens hiding ((:<), (:>))-import           Data.CircularSeq (focus,CSeq)+import           Data.CircularSeq (CSeq) import qualified Data.CircularSeq as C import           Data.Ext import qualified Data.Foldable as F-import           Data.Function (on, )+import           Data.Function (on) import           Data.Geometry.Box (IsBoxable(..)) import           Data.Geometry.LineSegment import           Data.Geometry.Point-import           Data.Geometry.Polygon (fromPoints, SimplePolygon, cmpExtreme, outerBoundary)+import           Data.Geometry.Polygon.Core (fromPoints, SimplePolygon, outerBoundary)+import           Data.Geometry.Polygon.Extremes(cmpExtreme) import           Data.Geometry.Properties import           Data.Geometry.Transformation import           Data.Geometry.Vector+import           Data.List.NonEmpty (NonEmpty(..))+import qualified Data.List.NonEmpty as NonEmpty import           Data.Maybe (fromJust) import           Data.Ord (comparing)+import           Data.Semigroup.Foldable (Foldable1(..)) import           Data.Sequence (viewl,viewr, ViewL(..), ViewR(..)) import qualified Data.Sequence as S+import           Data.Util  -- import           Data.Geometry.Ipe -- import           Debug.Trace@@ -88,7 +93,7 @@ --            -- print $ toPlaneGraph (Proxy :: Proxy DT) dt --            -- writeIpeFile outFile . singlePageFromContent $ out --            -- mapM_ (print . extremesNaive (v2 1 0)) polies---            pure $ map (flip rightTangent (point2 80 528)) polies+--            pure $ map (flip rightTangent (Point2 80 528)) polies   @@ -220,11 +225,12 @@ rotateTo'   :: Eq a => (a :+ b) -> CSeq (a :+ b) -> CSeq (a :+ b) rotateTo' x = fromJust . C.findRotateTo (coreEq x) - coreEq :: Eq a => (a :+ b) -> (a :+ b) -> Bool coreEq = (==) `on` (^.core)  +--------------------------------------------------------------------------------+-- * Computing Tangents  -- | Compute the lower tangent of the two polgyons --@@ -234,33 +240,45 @@ --        - The vertices of the polygons are given in clockwise order -- -- Running time: O(n+m), where n and m are the sizes of the two polygons respectively-lowerTangent                                       :: (Num r, Ord r)-                                                   => ConvexPolygon p r-                                                   -> ConvexPolygon p r-                                                   -> LineSegment 2 p r-lowerTangent (getVertices -> l) (getVertices -> r) = rotate xx yy zz zz''+lowerTangent       :: (Num r, Ord r)+                   => ConvexPolygon p r+                   -> ConvexPolygon p r+                   -> LineSegment 2 p r+lowerTangent lp rp = ClosedLineSegment l r   where-    xx = rightMost l-    yy = leftMost r+    mkH f = NonEmpty.fromList . F.toList . f . getVertices+    lh = mkH (C.rightElements . rightMost) lp+    rh = mkH (C.leftElements  . leftMost)  rp+    (Two (l :+ _) (r :+ _)) = lowerTangent' lh rh -    zz   = pred' yy-    zz'' = succ' xx+-- | Compute the lower tangent of the two convex chains lp and rp+--+--   pre: - the chains lp and rp have at least 1 vertex+--        - lp and rp are disjoint, and there is a vertical line+--          having lp on the left and rp on the right.+--        - The vertices in the left-chain are given in clockwise order, (right to left)+--        - The vertices in the right chain are given in counterclockwise order (left-to-right)+--+-- The result returned is the two endpoints l and r of the tangents,+-- and the remainders lc and rc of the chains (i.e.)  such that the lower hull+-- of both chains is: (reverse lc) ++ [l,h] ++ rc+--+-- Running time: \(O(n+m)\), where n and m are the sizes of the two chains+-- respectively+lowerTangent'       :: (Ord r, Num r, Foldable1 f)+                    => f (Point 2 r :+ p) -> f (Point 2 r :+ p)+                    -> Two ((Point 2 r :+ p) :+ [Point 2 r :+ p])+lowerTangent' l0 r0 = go (toNonEmpty l0) (toNonEmpty r0)+  where+    ne = NonEmpty.fromList+    isRight' []    _ _ = False+    isRight' (x:_) l r = ccw' l r x /= CCW +    go lh@(l:|ls) rh@(r:|rs) | isRight' rs l r = go lh      (ne rs)+                             | isRight' ls l r = go (ne ls) rh+                             | otherwise       = Two (l :+ ls) (r :+ rs) -    rotate x y z z''-      | focus z   `isRightOf` (focus x, focus y) = rotate x   z (pred' z) z''-                                                      -- rotate the right polygon CCW-      | focus z'' `isRightOf` (focus x, focus y) = rotate z'' y z         (succ' z'')-                                                      -- rotate the left polygon CW-      | otherwise                                = ClosedLineSegment (focus x)-                                                                     (focus y) -succ' :: CSeq a -> CSeq a-succ' = C.rotateR--pred' :: CSeq a -> CSeq a-pred' = C.rotateL- -- | Compute the upper tangent of the two polgyons -- --   pre: - polygons lp and rp have at least 1 vertex@@ -269,33 +287,43 @@ --        - The vertices of the polygons are given in clockwise order -- -- Running time: O(n+m), where n and m are the sizes of the two polygons respectively-upperTangent                                       :: (Num r, Ord r)-                                                   => ConvexPolygon p r-                                                   -> ConvexPolygon p r-                                                   -> LineSegment 2 p r-upperTangent (getVertices -> l) (getVertices -> r) = rotate xx yy zz zz'+upperTangent       :: (Num r, Ord r)+                   => ConvexPolygon p r+                   -> ConvexPolygon p r+                   -> LineSegment 2 p r+upperTangent lp rp = ClosedLineSegment l r   where-    xx = rightMost l-    yy = leftMost r--    zz  = succ' yy-    zz' = pred' xx--    rotate x y z z'-      | focus z  `isLeftOf` (focus x, focus y) = rotate x  z (succ' z) z'-                                                    -- rotate the right polygon CW-      | focus z' `isLeftOf` (focus x, focus y) = rotate z' y z        (pred' z')-                                                    -- rotate the left polygon CCW-      | otherwise                              = ClosedLineSegment (focus x)-                                                                   (focus y)+    mkH f = NonEmpty.fromList . F.toList . f . getVertices+    lh = mkH (C.leftElements  . rightMost) lp+    rh = mkH (C.rightElements . leftMost)  rp+    (Two (l :+ _) (r :+ _)) = upperTangent' lh rh -isRightOf           :: (Num r, Ord r)-                    => Point 2 r :+ p -> (Point 2 r :+ p', Point 2 r :+ p'') -> Bool-a `isRightOf` (b,c) = ccw (b^.core) (c^.core) (a^.core) == CW+-- | Compute the upper tangent of the two convex chains lp and rp+--+--   pre: - the chains lp and rp have at least 1 vertex+--        - lp and rp are disjoint, and there is a vertical line+--          having lp on the left and rp on the right.+--        - The vertices in the left-chain are given in clockwise order, (right to left)+--        - The vertices in the right chain are given in counterclockwise order (left-to-right)+--+-- The result returned is the two endpoints l and r of the tangents,+-- and the remainders lc and rc of the chains (i.e.)  such that the upper hull+-- of both chains is: (reverse lc) ++ [l,h] ++ rc+--+-- Running time: \(O(n+m)\), where n and m are the sizes of the two chains+-- respectively+upperTangent'       :: (Ord r, Num r, Foldable1 f)+                    => f (Point 2 r :+ p) -> f (Point 2 r :+ p)+                    -> Two ((Point 2 r :+ p) :+ [Point 2 r :+ p])+upperTangent' l0 r0 = go (toNonEmpty l0) (toNonEmpty r0)+  where+    ne = NonEmpty.fromList+    isLeft' []    _ _ = False+    isLeft' (x:_) l r = ccw' l r x /= CW -isLeftOf            :: (Num r, Ord r)-                    => Point 2 r :+ p -> (Point 2 r :+ p', Point 2 r :+ p'') -> Bool-a `isLeftOf` (b,c) = ccw (b^.core) (c^.core) (a^.core) == CCW+    go lh@(l:|ls) rh@(r:|rs) | isLeft' rs l r = go lh      (ne rs)+                             | isLeft' ls l r = go (ne ls) rh+                             | otherwise      = Two (l :+ ls) (r :+ rs)  -------------------------------------------------------------------------------- @@ -363,13 +391,13 @@ --               in if p cur nxt then go xs' else xs  -- test1 :: Num r => ConvexPolygon () r--- test1 = ConvexPolygon . fromPoints . map ext . reverse $ [origin, point2 1 4, point2 5 6, point2 10 3]+-- test1 = ConvexPolygon . fromPoints . map ext . reverse $ [origin, Point2 1 4, Point2 5 6, Point2 10 3]  -- test2 :: Num r => ConvexPolygon () r--- test2 = ConvexPolygon . fromPoints . map ext . reverse $ [point2 11 6, point2 10 10, point2 15 18, point2 12 5]+-- test2 = ConvexPolygon . fromPoints . map ext . reverse $ [Point2 11 6, Point2 10 10, Point2 15 18, Point2 12 5]  -- testA :: Num r => ConvexPolygon () r--- testA = ConvexPolygon . fromPoints . map ext $ [origin, point2 5 1, point2 2 2]+-- testA = ConvexPolygon . fromPoints . map ext $ [origin, Point2 5 1, Point2 2 2]  -- testB :: Num r => ConvexPolygon () r--- testB = ConvexPolygon . fromPoints . map ext $ [origin, point2 5 3, point2 (-2) 2, point2 (-2) 1]+-- testB = ConvexPolygon . fromPoints . map ext $ [origin, Point2 5 3, Point2 (-2) 2, Point2 (-2) 1]

@@ -0,0 +1,571 @@+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.Polygon.Core+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- A Polygon data type and some basic functions to interact with them.+--+--------------------------------------------------------------------------------+module Data.Geometry.Polygon.Core( PolygonType(..)+                                 , Polygon(..)+                                 , _SimplePolygon, _MultiPolygon+                                 , SimplePolygon, MultiPolygon, SomePolygon+++                                 , fromPoints++                                 , polygonVertices, listEdges++                                 , outerBoundary, outerBoundaryEdges+                                 , outerVertex, outerBoundaryEdge++                                 , polygonHoles, polygonHoles'+                                 , holeList++                                 , inPolygon, insidePolygon, onBoundary++                                 , area, signedArea++                                 , centroid+                                 , pickPoint++                                 , isTriangle++                                 , isCounterClockwise+                                 , toCounterClockWiseOrder, toCounterClockWiseOrder'+                                 , toClockwiseOrder, toClockwiseOrder'+                                 , reverseOuterBoundary++                                 , findDiagonal++                                 , withIncidentEdges, numberVertices++                                 , asSimplePolygon+                                 ) where++import           Control.DeepSeq+import           Control.Lens hiding (Simple)+import           Data.Bifoldable+import           Data.Bifunctor+import           Data.Bitraversable+import qualified Data.CircularSeq as C+import           Data.Ext+import qualified Data.Foldable as F+import           Data.Geometry.Boundary+import           Data.Geometry.Box+import           Data.Geometry.Line+import           Data.Geometry.LineSegment+import           Data.Geometry.Point+import           Data.Geometry.Properties+import           Data.Geometry.Transformation+import           Data.Geometry.Triangle (Triangle(..), inTriangle)+import           Data.Geometry.Vector+import qualified Data.List as List+import           Data.List.NonEmpty (NonEmpty(..))+import qualified Data.List.NonEmpty as NonEmpty+import           Data.Maybe (mapMaybe, catMaybes)+import           Data.Ord (comparing)+import           Data.Semigroup (sconcat)+import           Data.Semigroup.Foldable+import qualified Data.Sequence as Seq+import           Data.Util+import           Data.Vinyl.CoRec (asA)++--------------------------------------------------------------------------------++{- $setup+>>> :{+-- import qualified Data.CircularSeq as C+let simplePoly :: SimplePolygon () Rational+    simplePoly = SimplePolygon . C.fromList . map ext $ [ Point2 0 0+                                                        , Point2 10 0+                                                        , Point2 10 10+                                                        , Point2 5 15+                                                        , Point2 1 11+                                                        ]+:} -}++-- | We distinguish between simple polygons (without holes) and Polygons with holes.+data PolygonType = Simple | Multi+++data Polygon (t :: PolygonType) p r where+  SimplePolygon :: C.CSeq (Point 2 r :+ p)                         -> Polygon Simple p r+  MultiPolygon  :: C.CSeq (Point 2 r :+ p) -> [Polygon Simple p r] -> Polygon Multi  p r++-- | Prism to 'test' if we are a simple polygon+_SimplePolygon :: Prism' (Polygon Simple p r) (C.CSeq (Point 2 r :+ p))+_SimplePolygon = prism' SimplePolygon (\(SimplePolygon vs) -> Just vs)++-- | Prism to 'test' if we are a Multi polygon+_MultiPolygon :: Prism' (Polygon Multi p r) (C.CSeq (Point 2 r :+ p), [Polygon Simple p r])+_MultiPolygon = prism' (uncurry MultiPolygon) (\(MultiPolygon vs hs) -> Just (vs,hs))++instance Bifunctor (Polygon t) where+  bimap = bimapDefault++instance Bifoldable (Polygon t) where+  bifoldMap = bifoldMapDefault++instance Bitraversable (Polygon t) where+  bitraverse f g p = case p of+    SimplePolygon vs   -> SimplePolygon <$> bitraverseVertices f g vs+    MultiPolygon vs hs -> MultiPolygon  <$> bitraverseVertices f g vs+                                        <*> traverse (bitraverse f g) hs++instance (NFData p, NFData r) => NFData (Polygon t p r) where+  rnf (SimplePolygon vs)   = rnf vs+  rnf (MultiPolygon vs hs) = rnf (vs,hs)++bitraverseVertices     :: (Applicative f, Traversable t) => (p -> f q) -> (r -> f s)+                  -> t (Point 2 r :+ p) -> f (t (Point 2 s :+ q))+bitraverseVertices f g = traverse (bitraverse (traverse g) f)++type SimplePolygon = Polygon Simple++type MultiPolygon  = Polygon Multi++-- | Either a simple or multipolygon+type SomePolygon p r = Either (Polygon Simple p r) (Polygon Multi p r)++type instance Dimension (SomePolygon p r) = 2+type instance NumType   (SomePolygon p r) = r++-- | Polygons are per definition 2 dimensional+type instance Dimension (Polygon t p r) = 2+type instance NumType   (Polygon t p r) = r++instance (Show p, Show r) => Show (Polygon t p r) where+  show (SimplePolygon vs)   = "SimplePolygon (" <> show vs <> ")"+  show (MultiPolygon vs hs) = "MultiPolygon (" <> show vs <> ") (" <> show hs <> ")"++-- instance (Read p, Read r) => Show (Polygon t p r) where+--   show (SimplePolygon vs)   = "SimplePolygon (" <> show vs <> ")"+--   show (MultiPolygon vs hs) = "MultiPolygon (" <> show vs <> ") (" <> show hs <> ")"++++instance (Eq p, Eq r) => Eq (Polygon t p r) where+  (SimplePolygon vs)   == (SimplePolygon vs')    = vs == vs'+  (MultiPolygon vs hs) == (MultiPolygon vs' hs') = vs == vs' && hs == hs'++instance PointFunctor (Polygon t p) where+  pmap f (SimplePolygon vs)   = SimplePolygon (fmap (first f) vs)+  pmap f (MultiPolygon vs hs) = MultiPolygon  (fmap (first f) vs) (map (pmap f) hs)++instance Fractional r => IsTransformable (Polygon t p r) where+  transformBy = transformPointFunctor++instance IsBoxable (Polygon t p r) where+  boundingBox = boundingBoxList' . toListOf (outerBoundary.traverse.core)++type instance IntersectionOf (Line 2 r) (Boundary (Polygon t p r)) =+  '[Seq.Seq (Either (Point 2 r) (LineSegment 2 () r))]++type instance IntersectionOf (Point 2 r) (Polygon t p r) = [NoIntersection, Point 2 r]++instance (Fractional r, Ord r) => (Point 2 r) `IsIntersectableWith` (Polygon t p r) where+  nonEmptyIntersection = defaultNonEmptyIntersection+  q `intersects` pg = q `inPolygon` pg /= Outside+  q `intersect` pg | q `intersects` pg = coRec q+                   | otherwise         = coRec NoIntersection++-- instance IsIntersectableWith (Line 2 r) (Boundary (Polygon t p r)) where+--   nonEmptyIntersection _ _ (CoRec xs) = null xs+--   l `intersect` (Boundary (SimplePolygon vs)) =+--     undefined+  -- l `intersect` (Boundary (MultiPolygon vs hs)) = coRec .+  --    Seq.sortBy f . Seq.fromList+  --     . concatMap (unpack . (l `intersect`) . Boundary)+  --     $ SimplePolygon vs : hs+  --   where+  --     unpack (CoRec x) = x+  --     f = undefined+++++-- * Functions on Polygons++outerBoundary :: forall t p r. Lens' (Polygon t p r) (C.CSeq (Point 2 r :+ p))+outerBoundary = lens g s+  where+    g                     :: Polygon t p r -> C.CSeq (Point 2 r :+ p)+    g (SimplePolygon vs)  = vs+    g (MultiPolygon vs _) = vs++    s                           :: Polygon t p r -> C.CSeq (Point 2 r :+ p)+                                -> Polygon t p r+    s (SimplePolygon _)      vs = SimplePolygon vs+    s (MultiPolygon  _   hs) vs = MultiPolygon vs hs++polygonHoles :: forall p r. Lens' (Polygon Multi p r) [Polygon Simple p r]+polygonHoles = lens g s+  where+    g                     :: Polygon Multi p r -> [Polygon Simple p r]+    g (MultiPolygon _ hs) = hs+    s                     :: Polygon Multi p r -> [Polygon Simple p r]+                          -> Polygon Multi p r+    s (MultiPolygon vs _) = MultiPolygon vs++polygonHoles' :: Traversal' (Polygon t p r) [Polygon Simple p r]+polygonHoles' = \f -> \case+  p@(SimplePolygon _)  -> pure p+  (MultiPolygon vs hs) -> MultiPolygon vs <$> f hs++-- | Access the i^th vertex on the outer boundary+outerVertex   :: Int -> Lens' (Polygon t p r) (Point 2 r :+ p)+outerVertex i = outerBoundary.C.item i++-- running time: \(O(\log i)\)+outerBoundaryEdge     :: Int -> Polygon t p r -> LineSegment 2 p r+outerBoundaryEdge i p = let u = p^.outerVertex i+                            v = p^.outerVertex (i+1)+                        in LineSegment (Closed u) (Open v)+++-- | Get all holes in a polygon+holeList                     :: Polygon t p r -> [Polygon Simple p r]+holeList (SimplePolygon _)   = []+holeList (MultiPolygon _ hs) = hs+++-- | The vertices in the polygon. No guarantees are given on the order in which+-- they appear!+polygonVertices                      :: Polygon t p r+                                     -> NonEmpty.NonEmpty (Point 2 r :+ p)+polygonVertices (SimplePolygon vs)   = toNonEmpty vs+polygonVertices (MultiPolygon vs hs) =+  sconcat $ toNonEmpty vs NonEmpty.:| map polygonVertices hs+++-- | Creates a simple polygon from the given list of vertices.+--+-- pre: the input list constains no repeated vertices.+fromPoints :: [Point 2 r :+ p] -> SimplePolygon p r+fromPoints = SimplePolygon . C.fromList+++-- | The edges along the outer boundary of the polygon. The edges are half open.+--+-- running time: \(O(n)\)+outerBoundaryEdges :: Polygon t p r -> C.CSeq (LineSegment 2 p r)+outerBoundaryEdges = toEdges . (^.outerBoundary)++-- | Lists all edges. The edges on the outer boundary are given before the ones+-- on the holes. However, no other guarantees are given on the order.+--+-- running time: \(O(n)\)+listEdges    :: Polygon t p r -> [LineSegment 2 p r]+listEdges pg = let f = F.toList . outerBoundaryEdges+               in  f pg <> concatMap f (holeList pg)++-- | Pairs every vertex with its incident edges. The first one is its+-- predecessor edge, the second one its successor edge (in terms of+-- the ordering along the boundary).+--+--+-- >>> mapM_ print . polygonVertices $ withIncidentEdges simplePoly+-- Point2 [0 % 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] :+ ()))+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)+withIncidentEdges (MultiPolygon vs hs) = MultiPolygon vs' hs'+  where+    (SimplePolygon vs') = withIncidentEdges $ SimplePolygon vs+    hs' = map withIncidentEdges hs++-- -- | Gets the i^th edge on the outer boundary of the polygon, that is the edge+---- with vertices i and i+1 with respect to the current focus. All indices+-- -- modulo n.+-- --++-- | Given the vertices of the polygon. Produce a list of edges. The edges are+-- half-open.+toEdges    :: C.CSeq (Point 2 r :+ p) -> C.CSeq (LineSegment 2 p r)+toEdges vs = C.zipLWith (\p q -> LineSegment (Closed p) (Open q)) vs (C.rotateR vs)+  -- let vs' = F.toList vs in+  -- C.fromList $ zipWith (\p q -> LineSegment (Closed p) (Open q)) vs' (tail vs' ++ vs')+++-- | Test if q lies on the boundary of the polygon. Running time: O(n)+--+-- >>> Point2 1 1 `onBoundary` simplePoly+-- False+-- >>> Point2 0 0 `onBoundary` simplePoly+-- True+-- >>> Point2 10 0 `onBoundary` simplePoly+-- True+-- >>> Point2 5 13 `onBoundary` simplePoly+-- False+-- >>> Point2 5 10 `onBoundary` simplePoly+-- False+-- >>> Point2 10 5 `onBoundary` simplePoly+-- True+-- >>> Point2 20 5 `onBoundary` simplePoly+-- False+--+-- TODO: testcases multipolygon+onBoundary        :: (Fractional r, Ord r) => Point 2 r -> Polygon t p r -> Bool+q `onBoundary` pg = any (q `onSegment`) es+  where+    out = SimplePolygon $ pg^.outerBoundary+    es = concatMap (F.toList . outerBoundaryEdges) $ out : holeList pg++-- | Check if a point lies inside a polygon, on the boundary, or outside of the polygon.+-- Running time: O(n).+--+-- >>> Point2 1 1 `inPolygon` simplePoly+-- Inside+-- >>> Point2 0 0 `inPolygon` simplePoly+-- OnBoundary+-- >>> Point2 10 0 `inPolygon` simplePoly+-- OnBoundary+-- >>> Point2 5 13 `inPolygon` simplePoly+-- Inside+-- >>> Point2 5 10 `inPolygon` simplePoly+-- Inside+-- >>> Point2 10 5 `inPolygon` simplePoly+-- OnBoundary+-- >>> Point2 20 5 `inPolygon` simplePoly+-- Outside+--+-- TODO: Add some testcases with multiPolygons+-- TODO: Add some more onBoundary testcases+inPolygon                                :: forall t p r. (Fractional r, Ord r)+                                         => Point 2 r -> Polygon t p r+                                         -> PointLocationResult+q `inPolygon` pg+    | q `onBoundary` pg                             = OnBoundary+    | odd kl && odd kr && not (any (q `inHole`) hs) = Inside+    | otherwise                                     = Outside+  where+    l = horizontalLine $ q^.yCoord++    -- Given a line segment, compute the intersection point (if a point) with the+    -- line l+    intersectionPoint = asA @(Point 2 r) . (`intersect` l)++    -- Count the number of intersections that the horizontal line through q+    -- maxes with the polygon, that are strictly to the left and strictly to+    -- the right of q. If these numbers are both odd the point lies within the polygon.+    --+    --+    -- note that: - by the asA (Point 2 r) we ignore horizontal segments (as desired)+    --            - by the filtering, we effectively limit l to an open-half line, starting+    --               at the (open) point q.+    --            - by using half-open segments as edges we avoid double counting+    --               intersections that coincide with vertices.+    --            - If the point is outside, and on the same height as the+    --              minimum or maximum coordinate of the polygon. The number of+    --              intersections to the left or right may be one. Thus+    --              incorrectly classifying the point as inside. To avoid this,+    --              we count both the points to the left *and* to the right of+    --              p. Only if both are odd the point is inside.  so that if+    --              the point is outside, and on the same y-coordinate as one+    --              of the extermal vertices (one ofth)+    --+    -- See http://geomalgorithms.com/a03-_inclusion.html for more information.+    SP kl kr = count (\p -> (p^.xCoord) `compare` (q^.xCoord))+             . mapMaybe intersectionPoint . F.toList . outerBoundaryEdges $ pg++    -- For multi polygons we have to test if we do not lie in a hole .+    inHole = insidePolygon+    hs     = holeList pg++    count   :: (a -> Ordering) -> [a] -> SP Int Int+    count f = foldr (\x (SP lts gts) -> case f x of+                             LT -> SP (lts + 1) gts+                             EQ -> SP lts       gts+                             GT -> SP lts       (gts + 1)) (SP 0 0)+++-- | Test if a point lies strictly inside the polgyon.+insidePolygon        :: (Fractional r, Ord r) => Point 2 r -> Polygon t p r -> Bool+q `insidePolygon` pg = q `inPolygon` pg == Inside+++-- testQ = map (`inPolygon` testPoly) [ Point2 1 1    -- Inside+--                                    , Point2 0 0    -- OnBoundary+--                                    , Point2 5 14   -- Inside+--                                    , Point2 5 10   -- Inside+--                                    , Point2 10 5   -- OnBoundary+--                                    , Point2 20 5   -- Outside+--                                    ]++-- testPoly :: SimplePolygon () Rational+-- testPoly = SimplePolygon . C.fromList . map ext $ [ Point2 0 0+--                                                   , Point2 10 0+--                                                   , Point2 10 10+--                                                   , Point2 5 15+--                                                   , Point2 1 11+--                                                   ]++-- | Compute the area of a polygon+area                        :: Fractional r => Polygon t p r -> r+area poly@(SimplePolygon _) = abs $ signedArea poly+area (MultiPolygon vs hs)   = area (SimplePolygon vs) - sum [area h | h <- hs]+++-- | Compute the signed area of a simple polygon. The the vertices are in+-- clockwise order, the signed area will be negative, if the verices are given+-- in counter clockwise order, the area will be positive.+signedArea      :: Fractional r => SimplePolygon p r -> r+signedArea poly = x / 2+  where+    x = sum [ p^.core.xCoord * q^.core.yCoord - q^.core.xCoord * p^.core.yCoord+            | LineSegment' p q <- F.toList $ outerBoundaryEdges poly  ]+++-- | Compute the centroid of a simple polygon.+centroid      :: Fractional r => SimplePolygon p r -> Point 2 r+centroid poly = Point $ sum' xs ^/ (6 * signedArea poly)+  where+    xs = [ (toVec p ^+^ toVec q) ^* (p^.xCoord * q^.yCoord - q^.xCoord * p^.yCoord)+         | LineSegment' (p :+ _) (q :+ _) <- F.toList $ outerBoundaryEdges poly  ]++    sum' = F.foldl' (^+^) zero+++-- | Pick a  point that is inside the polygon.+--+-- (note: if the polygon is degenerate; i.e. has <3 vertices, we report a+-- vertex of the polygon instead.)+--+-- pre: the polygon is given in CCW order+--+-- running time: \(O(n)\)+pickPoint    :: (Ord r, Fractional r) => Polygon p t r -> Point 2 r+pickPoint pg | isTriangle pg = centroid . SimplePolygon $ pg^.outerBoundary+             | otherwise     = let LineSegment' (p :+ _) (q :+ _) = findDiagonal pg+                               in p .+^ (0.5 *^ (q .-. p))++-- | Test if the polygon is a triangle+--+-- running time: \(O(1)\)+isTriangle :: Polygon p t r -> Bool+isTriangle = \case+    SimplePolygon vs   -> go vs+    MultiPolygon vs [] -> go vs+    MultiPolygon _  _  -> False+  where+    go vs = case toNonEmpty vs of+              (_ :| [_,_]) -> True+              _            -> False++-- | Find a diagonal of the polygon.+--+-- pre: the polygon is given in CCW order+--+-- running time: \(O(n)\)+findDiagonal    :: (Ord r, Fractional r) => Polygon t p r -> LineSegment 2 p r+findDiagonal pg = List.head . catMaybes . F.toList $ diags+     -- note that a diagonal is guaranteed to exist, so the usage of head is safe.+  where+    vs      = pg^.outerBoundary+    diags   = C.zip3LWith f (C.rotateL vs) vs (C.rotateR vs)+    f u v w = case ccw (u^.core) (v^.core) (w^.core) of+                CCW      -> Just $ findDiag u v w+                            -- v is a convex vertex, so find a diagonal+                            -- (either uw) or from v to a point inside the+                            -- triangle+                CW       -> Nothing -- v is a reflex vertex+                CoLinear -> Nothing -- colinear vertex!?++    -- we test if uw is a diagonal by figuring out if there is a vertex+    -- strictly inside the triangle t. If there is no such vertex then uw must+    -- be a diagonal (i.e. uw intersects the polygon boundary iff there is a+    -- vtx inside t).  If there are vertices inside the triangle, we find the+    -- one z furthest from the line(segment) uw. It then follows that vz is a+    -- diagonal. Indeed this is pretty much the argument used to prove that any+    -- polygon can be triangulated. See BKOS Chapter 3 for details.+    findDiag u v w = let t  = Triangle u v w+                         uw = ClosedLineSegment u w+                     in maybe uw (ClosedLineSegment v)+                      . safeMaximumOn (distTo $ supportingLine uw)+                      . filter (\(z :+ _) -> z `inTriangle` t == Inside)+                      . F.toList . polygonVertices+                      $ pg++    distTo l (z :+ _) = sqDistanceTo z l+++safeMaximumOn   :: Ord b => (a -> b) -> [a] -> Maybe a+safeMaximumOn f = \case+  [] -> Nothing+  xs -> Just $ List.maximumBy (comparing f) xs+++-- | Test if the outer boundary of the polygon is in clockwise or counter+-- clockwise order.+--+-- running time: \(O(n)\)+--+isCounterClockwise :: (Eq r, Fractional r) => Polygon t p r -> Bool+isCounterClockwise = (\x -> x == abs x) . signedArea+                   . fromPoints . F.toList . (^.outerBoundary)+++-- | Make sure that every edge has the polygon's interior on its+-- right, by orienting the outer boundary into clockwise order, and+-- the inner borders (i.e. any holes, if they exist) into+-- counter-clockwise order.+--+-- running time: \(O(n)\)+-- | Orient the outer boundary of the polygon to clockwise order+toClockwiseOrder   :: (Eq r, Fractional r) => Polygon t p r -> Polygon t p r+toClockwiseOrder p = (toClockwiseOrder' p)&polygonHoles'.traverse %~ toCounterClockWiseOrder'++-- | Orient the outer boundary into clockwise order. Leaves any holes+-- as they are.+--+toClockwiseOrder'   :: (Eq r, Fractional r) => Polygon t p r -> Polygon t p r+toClockwiseOrder' pg+      | isCounterClockwise pg = reverseOuterBoundary pg+      | otherwise             = pg++-- | Make sure that every edge has the polygon's interior on its left,+-- by orienting the outer boundary into counter-clockwise order, and+-- the inner borders (i.e. any holes, if they exist) into clockwise order.+--+-- running time: \(O(n)\)+toCounterClockWiseOrder   :: (Eq r, Fractional r) => Polygon t p r -> Polygon t p r+toCounterClockWiseOrder p =+  (toCounterClockWiseOrder' p)&polygonHoles'.traverse %~ toClockwiseOrder'++-- | Orient the outer boundary into counter-clockwise order. Leaves+-- any holes as they are.+--+toCounterClockWiseOrder'   :: (Eq r, Fractional r) => Polygon t p r -> Polygon t p r+toCounterClockWiseOrder' p+      | not $ isCounterClockwise p = reverseOuterBoundary p+      | otherwise                  = p++reverseOuterBoundary   :: Polygon t p r -> Polygon t p r+reverseOuterBoundary p = p&outerBoundary %~ C.reverseDirection+++-- | Convert a Polygon to a simple polygon by forgetting about any holes.+asSimplePolygon                        :: Polygon t p r -> SimplePolygon p r+asSimplePolygon poly@(SimplePolygon _) = poly+asSimplePolygon (MultiPolygon vs _)    = SimplePolygon vs+++-- | assigns unique integer numbers to all vertices. Numbers start from 0, and+-- are increasing along the outer boundary. The vertices of holes+-- will be numbered last, in the same order.+--+-- >>> numberVertices simplePoly+-- SimplePolygon (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 ()])+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'

@@ -0,0 +1,39 @@+{-# LANGUAGE TemplateHaskell #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Data.Geometry.Polygon.Extremes+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- Finding the Extremal vertex of a polygon in a given direction.+--+--------------------------------------------------------------------------------+module Data.Geometry.Polygon.Extremes( cmpExtreme+                                     , extremesLinear+                                     ) where++import           Control.Lens hiding (Simple)+import           Data.Ext+import qualified Data.Foldable as F+import           Data.Geometry.Point+import           Data.Geometry.Polygon.Core+import           Data.Geometry.Vector++--------------------------------------------------------------------------------++-- | Comparison that compares which point is 'larger' in the direction given by+-- the vector u.+cmpExtreme       :: (Num r, Ord r)+                 => Vector 2 r -> Point 2 r :+ p -> Point 2 r :+ q -> Ordering+cmpExtreme u p q = u `dot` (p^.core .-. q^.core) `compare` 0+++-- | Finds the extreme points, minimum and maximum, in a given direction+--+-- running time: \(O(n)\)+extremesLinear     :: (Ord r, Num r) => Vector 2 r -> Polygon t p r+                   -> (Point 2 r :+ p, Point 2 r :+ p)+extremesLinear u p = let vs = p^.outerBoundary+                         f  = cmpExtreme u+                     in (F.minimumBy f vs, F.maximumBy f vs)

@@ -19,7 +19,7 @@                                        , queryRange                                        ) where -import           Algorithms.Geometry.ClosestPair.DivideAndConquer (mergeSortedListsBy)+import           Algorithms.DivideAndConquer (mergeSortedListsBy) import           Control.Lens import           Data.BinaryTree import           Data.Ext

@@ -69,8 +69,8 @@    (Slab h) `intersect` (Slab v) = coRec $ box low high     where-      low  = point2 (v^.start.core) (h^.start.core) :+ (v^.start.extra, h^.start.extra)-      high = point2 (v^.end.core)   (h^.end.core)   :+ (v^.end.extra,   h^.end.extra)+      low  = Point2 (v^.start.core) (h^.start.core) :+ (v^.start.extra, h^.start.extra)+      high = Point2 (v^.end.core)   (h^.end.core)   :+ (v^.end.extra,   h^.end.extra)   

@@ -1,10 +1,4 @@-{-# LANGUAGE FlexibleContexts #-}-{-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE UndecidableInstances #-}-{-# LANGUAGE StandaloneDeriving #-}-{-# LANGUAGE GeneralizedNewtypeDeriving  #-}-{-# LANGUAGE ScopedTypeVariables  #-}-{-# LANGUAGE DeriveFunctor  #-} module Data.Geometry.Transformation where  import           Control.Lens (lens,Lens',set)@@ -130,8 +124,6 @@  -------------------------------------------------------------------------------- -- * Functions that execute transformations---- type AlwaysTrueTransformation d = (Arity (1 + d), AlwaysTrueSnoc d, Arity d, Index' (1+d-1) (1+d))  translateBy :: ( IsTransformable g, Num (NumType g)                , Arity (Dimension g), Arity (Dimension g + 1)

@@ -26,6 +26,7 @@  -------------------------------------------------------------------------------- +-- | Triangles in \(d\)-dimensional space. data Triangle d p r = Triangle (Point d r :+ p)                                (Point d r :+ p)                                (Point d r :+ p)@@ -47,10 +48,11 @@ instance (Fractional r, Arity d, Arity (d + 1)) => IsTransformable (Triangle d p r) where   transformBy = transformPointFunctor - -- | convenience function to construct a triangle without associated data.-triangle'       :: Point d r -> Point d r -> Point d r -> Triangle d () r-triangle' p q r = Triangle (ext p) (ext q) (ext r)+pattern Triangle' :: Point d r -> Point d r -> Point d r -> Triangle d () r+pattern Triangle' p q r <- Triangle (p :+ ()) (q :+ ()) (r :+ ())+  where+    Triangle' p q r = Triangle (ext p) (ext q) (ext r)   sideSegments                  :: Triangle d p r -> [LineSegment d p r]

@@ -41,13 +41,13 @@ --------------------------------------------------------------------------------  type instance Dimension (Vector d r) = d-type instance NumType (Vector d r) = r+type instance NumType   (Vector d r) = r  instance (Arbitrary r, Arity d) => Arbitrary (Vector d r) where   arbitrary = vectorFromListUnsafe <$> infiniteList  --- | Test if v is a scalar multiple of u.+-- | 'isScalarmultipleof u v' test if v is a scalar multiple of u. -- -- >>> Vector2 1 1 `isScalarMultipleOf` Vector2 10 10 -- True@@ -69,7 +69,7 @@ {-# SPECIALIZE     isScalarMultipleOf :: (Eq r, Fractional r) => Vector 2 r -> Vector 2 r -> Bool  #-} --- | Get the scalar labmda s.t. v = lambda * u (if it exists)+-- | scalarMultiple u v computes the scalar labmda s.t. v = lambda * u (if it exists) scalarMultiple     :: (Eq r, Fractional r, Arity d)                    => Vector d r -> Vector d r -> Maybe r scalarMultiple u v

@@ -282,6 +282,6 @@ destruct v = (head $ F.toList v, vectorFromListUnsafe . tail $ F.toList v)   -- FIXME: this implementaion of tail is not particularly nice -snoc     :: (ImplicitArity d, ImplicitArity (S d), (1 + FromPeano d) ~ (FromPeano d + 1))-         => VectorFamily d r -> r -> VectorFamily (S d) r-snoc = flip V.snoc+-- snoc     :: (ImplicitArity d, ImplicitArity (S d))+--          => VectorFamily d r -> r -> VectorFamily (S d) r+-- snoc = flip V.snoc

@@ -223,7 +223,7 @@      sing x = x NonEmpty.:| [] -    vts    = map (\(p,sp) -> (p,map (^.extra) . sortArround (ext p) <$> sp))+    vts    = map (\(p,sp) -> (p,map (^.extra) . sortAround (ext p) <$> sp))            . M.assocs $ pts     -- vertex Data     vxData = V.fromList . map (\(p,sp) -> VertexData p (sp^._1)) $ vts@@ -270,6 +270,7 @@ vertices' = PG.vertices' . _graph  -- | Enumerate all vertices, together with their vertex data+-- -- >>> mapM_ print $ vertices smallG -- (VertexId 0,VertexData {_location = Point2 [0,0], _vData = 0}) -- (VertexId 1,VertexData {_location = Point2 [2,2], _vData = 1})

@@ -91,9 +91,8 @@  -------------------------------------------------------------------------------- --- | Transforms the planar graph into a format taht can be easily converted--- into JSON format. For every vertex, the adjacent vertices are given in--- counter clockwise order.+-- | Transforms the plane graph into adjacency lists. For every+-- vertex, the adjacent vertices are given in counter clockwise order. -- -- See 'toAdjacencyLists' for notes on how we handle self-loops. --@@ -102,6 +101,11 @@ toAdjRep = first (\(PGA.Vtx v aj (VertexData p x)) -> Vtx v p aj x) . PGIO.toAdjRep          .  view graph +-- | Given the AdjacencyList representation of a plane graph,+-- construct the plane graph representing it. All the adjacencylists+-- should be in counter clockwise order.+--+-- running time: \(O(n)\) fromAdjRep    :: proxy s -> Gr (Vtx v e r) (Face f) -> PlaneGraph s v e f r fromAdjRep px = PlaneGraph . PGIO.fromAdjRep px               . first (\(Vtx v p aj x) -> PGA.Vtx v aj $ VertexData p x)@@ -124,5 +128,8 @@     -- sort the adjacencies around every vertex v     sort' (Vtx v p ajs x) = Vtx v p (List.sortBy (around p) ajs) x     around p (a,_) (b,_) = ccwCmpAround (ext p) (location' V.! a) (location' V.! b)+                           -- note: since the graph is planar, there should not be+                           -- any pairs of points for which ccwCmpAround returns EQ+                           -- hence, no need to pick a secondary comparison  --------------------------------------------------------------------------------

@@ -1,21 +1,4 @@-faces:-- incidentEdge:-  - 0-  - 4-  fData: []-- incidentEdge:-  - 0-  - 2-  fData: []-- incidentEdge:-  - 0-  - 1-  fData: []-- incidentEdge:-  - 0-  - 3-  fData: []-ajacencies:+adjacencies: - adj:   - - 4     - []@@ -88,3 +71,20 @@   vData:   - []   - []+faces:+- fData: []+  incidentEdge:+  - 0+  - 4+- fData: []+  incidentEdge:+  - 0+  - 2+- fData: []+  incidentEdge:+  - 0+  - 1+- fData: []+  incidentEdge:+  - 0+  - 3

@@ -1,16 +1,3 @@-faces:-- incidentEdge:-  - 0-  - 2-  fData: OuterFace-- incidentEdge:-  - 0-  - 1-  fData: A-- incidentEdge:-  - 0-  - 3-  fData: B ajacencies: - adj:   - - 2@@ -56,3 +43,16 @@   - -1   - 4   vData: 3+faces:+- fData: OuterFace+  incidentEdge:+  - 0+  - 2+- fData: A+  incidentEdge:+  - 0+  - 1+- fData: B+  incidentEdge:+  - 0+  - 3