hgeometry 0.13 → 0.14

raw patch · 66 files changed

+2306/−954 lines, 66 filesdep ~hgeometry-combinatorialdep ~vector-circularPVP ok

version bump matches the API change (PVP)

Dependency ranges changed: hgeometry-combinatorial, vector-circular

API changes (from Hackage documentation)

Files

• benchmark/Algorithms/Geometry/ConvexHull/GrahamFamPeano.hs +0/−107
• benchmark/Algorithms/Geometry/LineSegmentIntersection/Bench.hs +10/−8
• benchmark/Algorithms/Geometry/LineSegmentIntersection/BentleyOttmannNoExt.hs +440/−0
• benchmark/Algorithms/Geometry/LineSegmentIntersection/BentleyOttmannOld.hs +15/−26
• benchmark/Algorithms/Geometry/LineSegmentIntersection/TypesNoExt.hs +200/−0
• benchmark/Algorithms/Geometry/PolygonTriangulation/MakeMonotoneOld.hs +6/−7
• changelog +21/−2
• changelog.org +21/−2
• doctests.hs +2/−0
• hgeometry.cabal +24/−10
• src/Algorithms/Geometry/DelaunayTriangulation/Types.hs +5/−6
• src/Algorithms/Geometry/EuclideanMST.hs +1/−2
• src/Algorithms/Geometry/LineSegmentIntersection.hs +16/−14
• src/Algorithms/Geometry/LineSegmentIntersection/BentleyOttmann.hs +205/−75
• src/Algorithms/Geometry/LineSegmentIntersection/BooleanSweep.hs +65/−46
• src/Algorithms/Geometry/LineSegmentIntersection/Naive.hs +28/−26
• src/Algorithms/Geometry/LineSegmentIntersection/Types.hs +162/−46
• src/Algorithms/Geometry/LowerEnvelope/DualCH.hs +1/−1
• src/Algorithms/Geometry/PolygonTriangulation/EarClip.hs +1/−1
• src/Algorithms/Geometry/PolygonTriangulation/MakeMonotone.hs +25/−26
• src/Algorithms/Geometry/PolygonTriangulation/Triangulate.hs +7/−9
• src/Algorithms/Geometry/PolygonTriangulation/TriangulateMonotone.hs +6/−8
• src/Algorithms/Geometry/PolygonTriangulation/Types.hs +13/−15
• src/Algorithms/Geometry/RayShooting/Naive.hs +88/−0
• src/Algorithms/Geometry/SSSP.hs +3/−3
• src/Algorithms/Geometry/VisibilityPolygon/Lee.hs +2/−8
• src/Algorithms/Geometry/WSPD.hs +4/−4
• src/Data/Geometry/Arrangement/Internal.hs +17/−20
• src/Data/Geometry/Ball.hs +1/−1
• src/Data/Geometry/BezierSpline.hs +6/−16
• src/Data/Geometry/Box.hs +16/−4
• src/Data/Geometry/Box/Internal.hs +31/−8
• src/Data/Geometry/Box/Sides.hs +1/−1
• src/Data/Geometry/Duality.hs +2/−2
• src/Data/Geometry/HyperPlane.hs +5/−3
• src/Data/Geometry/Interval.hs +42/−14
• src/Data/Geometry/KDTree.hs +1/−1
• src/Data/Geometry/Line.hs +4/−3
• src/Data/Geometry/Line/Internal.hs +27/−27
• src/Data/Geometry/LineSegment.hs +109/−2
• src/Data/Geometry/LineSegment/Internal.hs +155/−21
• src/Data/Geometry/PlanarSubdivision.hs +17/−22
• src/Data/Geometry/PlanarSubdivision/Basic.hs +9/−11
• src/Data/Geometry/PlanarSubdivision/Dynamic.hs +113/−117
• src/Data/Geometry/PlanarSubdivision/Merge.hs +5/−7
• src/Data/Geometry/Point.hs +26/−1
• src/Data/Geometry/Point/Class.hs +9/−8
• src/Data/Geometry/Point/Internal.hs +44/−7
• src/Data/Geometry/PointLocation/PersistentSweep.hs +1/−2
• src/Data/Geometry/PolyLine.hs +8/−1
• src/Data/Geometry/Polygon.hs +12/−1
• src/Data/Geometry/Polygon/Core.hs +13/−15
• src/Data/Geometry/RangeTree.hs +3/−4
• src/Data/Geometry/RangeTree/Measure.hs +4/−5
• src/Data/Geometry/Slab.hs +18/−17
• src/Data/Geometry/SubLine.hs +142/−67
• src/Data/Geometry/Transformation.hs +2/−2
• src/Data/Geometry/Transformation/Internal.hs +7/−2
• src/Data/Geometry/Triangle.hs +1/−1
• src/Data/Geometry/Vector.hs +3/−3
• src/Data/Geometry/Vector/VectorFamily.hs +24/−17
• src/Data/Geometry/Vector/VectorFamilyPeano.hs +4/−3
• src/Data/Geometry/Vector/VectorFixed.hs +5/−4
• src/Data/Geometry/VerticalRayShooting/PersistentSweep.hs +9/−26
• src/Data/PlaneGraph/Core.hs +33/−29
• src/Data/PlaneGraph/IO.hs +6/−7

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

@@ -1,7 +1,8 @@ module Algorithms.Geometry.LineSegmentIntersection.Bench (benchmark) where -import qualified Algorithms.Geometry.LineSegmentIntersection.BentleyOttmann    as BONew-import qualified Algorithms.Geometry.LineSegmentIntersection.BentleyOttmannOld as BOOld+import qualified Algorithms.Geometry.LineSegmentIntersection.BentleyOttmann     as BONew+import qualified Algorithms.Geometry.LineSegmentIntersection.BentleyOttmannNoExt as BONoExt+import qualified Algorithms.Geometry.LineSegmentIntersection.BentleyOttmannOld   as BOOld  import           Control.DeepSeq import           Control.Lens@@ -29,11 +30,11 @@ --------------------------------------------------------------------------------  genPts                 :: (Ord r, Random r, RandomGen g)-                       => Int -> Rand g [LineSegment 2 () r]-genPts n = replicateM n sampleLineSegment+                       => Int -> Rand g [LineSegment 2 () r :+ ()]+genPts n = map ext <$> replicateM n sampleLineSegment  -- | Benchmark computing the closest pair-benchBuild    :: (Ord r, Fractional r, NFData r) => [LineSegment 2 () r] -> Benchmark+benchBuild    :: (Ord r, Fractional r, NFData r) => [LineSegment 2 () r :+ ()] -> Benchmark benchBuild ss = bgroup "LineSegs" [ bgroup (show n) (build $ take n ss)                                   | n <- sizes' ss                                   ]@@ -42,9 +43,10 @@       -- let n = length pts in [ n*i `div` 100 | i <- [10,20,25,50,75,100]]      build segs = [ bench "sort"     $ nf sort' segs-                 , bench "Old"      $ nf BOOld.intersections segs+                 , bench "Old"      $ nf BOOld.intersections (map (^.core) segs)+                 , bench "NoExt"    $ nf BONoExt.intersections (map (^.core) segs)                  , bench "New"      $ nf BONew.intersections segs                  ] -sort' :: Ord r => [LineSegment 2 () r] -> [Point 2 r]-sort' = List.sort . concatMap (\s -> s^..endPoints.core)+sort' :: Ord r => [LineSegment 2 () r :+ ()] -> [Point 2 r]+sort' = List.sort . concatMap (\s -> s^..core.endPoints.core)

@@ -0,0 +1,440 @@+{-# LANGUAGE ScopedTypeVariables #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Algorithms.Geometry.LineSegmentIntersection.BentleyOttmann+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--+-- The \(O((n+k)\log n)\) time line segment intersection algorithm by Bentley+-- and Ottmann.+--+--------------------------------------------------------------------------------+module Algorithms.Geometry.LineSegmentIntersection.BentleyOttmannNoExt+  ( intersections+  , interiorIntersections+  ) where++import           Algorithms.Geometry.LineSegmentIntersection.TypesNoExt+import           Control.Lens hiding (contains)+import           Data.Ext+import qualified Data.Foldable as F+import           Data.Function (on)+import           Data.Geometry.Interval+import           Data.Geometry.LineSegment+import           Data.Geometry.Point+import           Data.Geometry.Properties+import qualified Data.List as L+import           Data.List.NonEmpty (NonEmpty(..))+import qualified Data.List.NonEmpty as NonEmpty+import qualified Data.Map as M+import           Data.Maybe+import           Data.Ord (Down(..), comparing)+import qualified Data.Set as EQ -- event queue+import qualified Data.Set as SS -- status struct+import qualified Data.Set.Util as SS -- status struct+import           Data.Vinyl+import           Data.Vinyl.CoRec++--------------------------------------------------------------------------------++-- | Compute all intersections+--+-- \(O((n+k)\log n)\), where \(k\) is the number of intersections.+intersections    :: (Ord r, Fractional r)+                 => [LineSegment 2 p r] -> Intersections p r+intersections ss = merge $ sweep pts SS.empty+  where+    pts = EQ.fromAscList . groupStarts . L.sort . concatMap asEventPts $ ss++-- | Computes all intersection points p s.t. p lies in the interior of at least+-- one of the segments.+--+--  \(O((n+k)\log n)\), where \(k\) is the number of intersections.+interiorIntersections :: (Ord r, Fractional r)+                       => [LineSegment 2 p r] -> Intersections p r+interiorIntersections = M.filter isInteriorIntersection . intersections++-- | Computes the event points for a given line segment+asEventPts   :: Ord r => LineSegment 2 p r -> [Event p r]+asEventPts s = let [p,q] = L.sortBy ordPoints [s^.start.core,s^.end.core]+               in [Event p (Start $ s :| []), Event q (End s)]++-- | Group the segments with the intersection points+merge :: (Ord r, Fractional r) =>  [IntersectionPoint p r] -> Intersections p r+merge = foldr (\(IntersectionPoint p a) -> M.insertWith (<>) p a) M.empty++-- | Group the startpoints such that segments with the same start point+-- correspond to one event.+groupStarts                          :: Eq r => [Event p r] -> [Event p r]+groupStarts []                       = []+groupStarts (Event p (Start s) : es) = Event p (Start ss) : groupStarts rest+  where+    (ss',rest) = L.span sameStart es+    -- FIXME: this seems to keep the segments on decreasing y, increasing x. shouldn't we+    -- sort them cyclically around p instead?+    ss         = let (x:|xs) = s+                 in x :| (xs ++ concatMap startSegs ss')++    sameStart (Event q (Start _)) = p == q+    sameStart _                   = False+groupStarts (e : es)                 = e : groupStarts es++--------------------------------------------------------------------------------+-- * Data type for Events++-- | Type of segment+data EventType s = Start !(NonEmpty s)| Intersection | End !s deriving (Show)++instance Eq (EventType s) where+  a == b = a `compare` b == EQ++instance Ord (EventType s) where+  (Start _)    `compare` (Start _)    = EQ+  (Start _)    `compare` _            = LT+  Intersection `compare` (Start _)    = GT+  Intersection `compare` Intersection = EQ+  Intersection `compare` (End _)      = LT+  (End _)      `compare` (End _)      = EQ+  (End _)      `compare` _            = GT++-- | The actual event consists of a point and its type+data Event p r = Event { eventPoint :: !(Point 2 r)+                       , eventType  :: !(EventType (LineSegment 2 p r))+                       } deriving (Show,Eq)++instance Ord r => Ord (Event p r) where+  -- decreasing on the y-coord, then increasing on x-coord, and increasing on event-type+  (Event p s) `compare` (Event q t) = case ordPoints p q of+                                        EQ -> s `compare` t+                                        x  -> x++-- | Get the segments that start at the given event point+startSegs   :: Event p r -> [LineSegment 2 p r]+startSegs e = case eventType e of+                Start ss -> NonEmpty.toList ss+                _        -> []++--------------------------------------------------------------------------------+++--------------------------------------------------------------------------------+-- * The Main Sweep++type EventQueue      p r = EQ.Set (Event p r)+type StatusStructure p r = SS.Set (LineSegment 2 p r)++-- | Run the sweep handling all events+sweep       :: (Ord r, Fractional r)+            => EventQueue p r -> StatusStructure p r -> [IntersectionPoint p r]+sweep eq ss = case EQ.minView eq of+    Nothing      -> []+    Just (e,eq') -> handle e eq' ss++isClosedStart                     :: Eq r => Point 2 r -> LineSegment 2 p r -> Bool+isClosedStart p (LineSegment s e)+  | p == s^.unEndPoint.core       = isClosed s+  | otherwise                     = isClosed e+++-- data AssocKind b a = Start b a | End b a | Neighter a++-- -- | test if the given segment has p as its endpoint, an construct the+-- -- appropriate associated representing that.+-- mkAssociated                :: Point 2 r -> LineSegment 2 p r -> AssocKind (LineSegment 2 p r)+-- mkAssociated p s@(LineSegment a b)+--   | p == a^.unEndPoint.core = Start a s+--   | p == b^.unEndPoint.core = End b s+--   | otherwise               = Neighter s++-- -- -- | We need to report a segment as an segment for starting point p if+-- -- -- it is a closed segment starting at p, or an open segment starting+-- -- -- at p that intersects with some other segment.  since the segments+-- -- -- are given in sorted order around s, we can just look at the next+-- -- -- segment to see if we should report such an open-ended segment.+-- -- shouldReportStart   :: Point 2 r -> [LineSegment 2 p r] -> Associated p r+-- -- shouldReportStart p = go . map (categorize p)+-- --   where+-- --     go []     = mempty+-- --     go (s:ss) = let (xs,ys) = List.span overlapsWith s ss+-- --                 in case s of+-- --                      Start (Closed _) s' -> Asso+++++++-- --     (s@(LineSegment a b):ss)+-- --         | p == a^.unEndPoint.core =+++-- --           if isClosed a || overlapsWithNext ss+-- --                                     then Associated [s] [] [] <> go ss+-- --         -- | p == b^.unEndPoint.core = Associated [] [s] []++++++++--     _  []                  = mempty+--     go certainlyReport (s:ss) = let x  = mkAssociated p s+--                                     x' = then x else mempty+--                                 in++++--       case shouldReport mp s of++++++--       mkAssociated mp s <> go (Just s) ss+++--     mkAsscoiated _ s@(LineSegment a b)+--       | p == a^.unEndPoint.core = if isClosed a ||++++--       = Associated [s] [] []+--       | p == b^.unEndPoint.core = Associated [] [s] []+--       | otherwise               = mempty++-- _ []     = []++++-- shouldReportStart _ []     = []+-- shouldReportStart p (s:ss) = case hasStartingPoint p s of+--                                Nothing            -> shouldReportStart ss -- don't report the seg+--                                Just (Closed _, s) -> s : shouldReportStart ss+--                                Just (Open _, )+++-- -- [s] | isClosedStart p s = [s]+-- --                         | otherwise         = []+-- -- shouldReportStart p (s:s':ss) | isStart p s =++++-- (s:ss) = isClosedStart p s ||+++-- shouldReport   :: Eq r => Point 2 r -> [LineSegment 2 p r] -> Associated p r+-- shouldReport p = foldMap (\(s,c) -> case c of+--                                       Start'   -> Associated [s] [] []+--                                       End'     -> Associated [] [s] []+--                                       Neighter -> Associated [] [] [s]+--                          )+--                . overlapsOr (\(LineSegment a b,c) -> case c of+--                                              Start'   -> isClosed a+--                                              End'     -> isClosed b+--                                              Neighter -> False+--                               ) (overlap p)+--                . map (\s -> (s, categorize p s))++overlap :: Point 2 r -> (LineSegment 2 q r, Cat) -> (LineSegment 2 q r, Cat) -> Bool+overlap p s1 s2 = go (toStart s1) (toStart s2)+  where+    toStart (s@(LineSegment a b),c) = case c of+                                        Start' -> (s,False)+                                        End'   -> (LineSegment b a,False) -- flip to start+                                        Neighter -> (s, True)+    go = undefined+++++data Cat = Start' | End' | Neighter++categorize p (LineSegment a b)+  | p == a^.unEndPoint.core = Start'+  | p == b^.unEndPoint.core = End'+  | otherwise               = Neighter++++overlapsOr     :: (a -> Bool)+               -> (a -> a -> Bool)+               -> [a]+               -> [a]+overlapsOr p q = map fst . filter snd . map (\((a,b),b') -> (a, b || b'))+               . overlapsWithNeighbour (q `on` fst)+               . map (\x -> (x, p x))++overlapsWithNeighbour   :: (a -> a -> Bool) -> [a] -> [(a,Bool)]+overlapsWithNeighbour p = go0+  where+    go0 = \case+      []     -> []+      (x:xs) -> go x False xs++    go x b = \case+      []     -> []+      (y:ys) -> let b' = p x y+                in (x,b || b') : go y b' ys++++++++++annotateReport   :: (a -> Bool) -> [a] -> [(a,Bool)]+annotateReport p = map (\x -> (x, p x))+++overlapsWithNext'   :: (a -> a -> Bool) -> [a] -> [(a,Bool)]+overlapsWithNext' p = go+  where+    go = \case+      []           -> []+      [x]          -> [(x,False)]+      (x:xs@(y:_)) -> (x,p x y) : go xs++overlapsWithPrev'   :: (a -> a -> Bool) -> [a] -> [(a,Bool)]+overlapsWithPrev' p = go0+  where+    go0 = \case+      []     -> []+      (x:xs) -> (x,False) : go x xs++    go x = \case+      []     -> []+      (y:ys) -> (y,p x y) : go y ys+++++++overlapsWithNeighbour2 p = map (\((a,b),b') -> (a, b || b'))+                         . overlapsWithNext' (p `on` fst)+                         . overlapsWithPrev' p++shouldBe :: Eq a => a -> a -> Bool+shouldBe = (==)++propSameAsSeparate p xs = overlapsWithNeighbour p xs `shouldBe` overlapsWithNeighbour2 p xs++test' = overlapsWithNeighbour (==) testOverlapNext+testOverlapNext = [1,2,3,3,3,5,6,6,8,10,11,34,2,2,3]++-- reportOverlappingBy :: Eq a => (a -> Bool) -> [a] -> [a]+-- reportOverlappingBy p = \case+--   []     -> []+--   (x:xs) -> L.span+++-- | Handle an event point+handle                           :: forall r p. (Ord r, Fractional r)+                                 => Event p r -> EventQueue p r -> StatusStructure p r+                                 -> [IntersectionPoint p r]+handle e@(eventPoint -> p) eq ss = toReport <> sweep eq' ss'+  where+    starts                   = startSegs e+    (before,contains',after) = extractContains p ss+    (ends,contains)          = L.partition (endsAt p) contains'+    -- starting segments, exluding those that have an open starting point+    starts' = filter (isClosedStart p) starts+++    -- starts'' = shouldReport p . SS.toAscList $ newSegs+    -- FIXME: we should look at the starts in-order (around p).+    -- closed endpoints we should report anyway. For an open endpoint+    -- we should check if it overlaps with a sucessor or predecessor+    -- to see if we have to report it.++    -- I think we could get those from the 'toStatusStruct' structure below++    -- any (closed) ending segments at this event point.+    closedEnds = filter (isClosedStart p) ends++    toReport = case starts' <> contains' of+                 (_:_:_) -> [mkIntersectionPoint p (starts' <> closedEnds) contains]+                 _       -> []++    -- new status structure+    ss' = before `SS.join` newSegs `SS.join` after+    newSegs = toStatusStruct p $ starts ++ contains+++    -- the new eeventqueue+    eq' = foldr EQ.insert eq es+    -- the new events:+    es | F.null newSegs  = maybeToList $ app (findNewEvent p) sl sr+       | otherwise       = let s'  = SS.lookupMin newSegs+                               s'' = SS.lookupMax newSegs+                           in catMaybes [ app (findNewEvent p) sl  s'+                                        , app (findNewEvent p) s'' sr+                                        ]+    sl = SS.lookupMax before+    sr = SS.lookupMin after++    app f x y = do { x' <- x ; y' <- y ; f x' y'}++-- | split the status structure, extracting the segments that contain p.+-- the result is (before,contains,after)+extractContains      :: (Fractional r, Ord r)+                     => Point 2 r -> StatusStructure p r+                     -> (StatusStructure p r, [LineSegment 2 p r], StatusStructure p r)+extractContains p ss = (before, F.toList mid1 <> F.toList mid2, after)+  where+    (before, mid1, after') = SS.splitOn (xCoordAt $ p^.yCoord) (p^.xCoord) ss+    -- Make sure to also select the horizontal segments containing p+    (mid2, after) = SS.spanAntitone (intersects p) after'+++-- | Given a point and the linesegements that contain it. Create a piece of+-- status structure for it.+toStatusStruct      :: (Fractional r, Ord r)+                    => Point 2 r -> [LineSegment 2 p r] -> StatusStructure p r+toStatusStruct p xs = ss `SS.join` hors+  -- ss { SS.nav = ordAtNav $ p^.yCoord } `SS.join` hors+  where+    (hors',rest) = L.partition isHorizontal xs+    ss           = SS.fromListBy (ordAtY $ maxY xs) rest+    hors         = SS.fromListBy (comparing rightEndpoint) hors'++    isHorizontal s  = s^.start.core.yCoord == s^.end.core.yCoord++    -- find the y coord of the first interesting thing below the sweep at y+    maxY = maximum . filter (< p^.yCoord)+         . concatMap (\s -> [s^.start.core.yCoord,s^.end.core.yCoord])++-- | Get the right endpoint of a segment+rightEndpoint   :: Ord r => LineSegment 2 p r -> r+rightEndpoint s = (s^.start.core.xCoord) `max` (s^.end.core.xCoord)++-- | Test if a segment ends at p+endsAt                      :: Ord r => Point 2 r -> LineSegment 2 p r -> Bool+endsAt p (LineSegment' a b) = all (\q -> ordPoints (q^.core) p /= GT) [a,b]++--------------------------------------------------------------------------------+-- * Finding New events++-- | Find all events+findNewEvent       :: (Ord r, Fractional r)+                   => Point 2 r -> LineSegment 2 p r -> LineSegment 2 p r+                   -> Maybe (Event p r)+findNewEvent p l r = match (l `intersect` r) $+     H (const Nothing) -- NoIntersection+  :& H (\q -> if ordPoints q p == GT then Just (Event q Intersection)+                                     else Nothing)+  :& H (const Nothing) -- full segment intersectsions are handled+                       -- at insertion time+  :& RNil++++type R = Rational++seg1, seg2 :: LineSegment 2 () R+seg1 = ClosedLineSegment (ext $ Point2 0 0) (ext $ Point2 0 10)+seg2 = ClosedLineSegment (ext $ Point2 0 1) (ext $ Point2 0 5)

@@ -12,7 +12,11 @@ -------------------------------------------------------------------------------- module Algorithms.Geometry.LineSegmentIntersection.BentleyOttmannOld where -import           Algorithms.Geometry.LineSegmentIntersection+import           Algorithms.Geometry.LineSegmentIntersection.TypesNoExt( Intersections+                                                            , IntersectionPoint(..)+                                                            , Associated(..)+                                                            , mkIntersectionPoint+                                                            ) import           Control.Lens hiding (contains) import           Data.Ext import qualified Data.Foldable as F@@ -33,6 +37,8 @@  -------------------------------------------------------------------------------- +-- todo; use an old copy of the imports as well.+ -- | Compute all intersections -- -- \(O((n+k)\log n)\), where \(k\) is the number of intersections.@@ -48,15 +54,19 @@ --  \(O((n+k)\log n)\), where \(k\) is the number of intersections. interiorIntersections :: (Ord r, Fractional r)                        => [LineSegment 2 p r] -> Intersections p r-interiorIntersections = M.filter (not . isEndPointIntersection) . intersections+interiorIntersections = M.filter isInteriorIntersection . intersections +isInteriorIntersection :: Associated p r -> Bool+isInteriorIntersection = not . null . _interiorTo++ -- | Computes the event points for a given line segment asEventPts   :: Ord r => LineSegment 2 p r -> [Event p r] asEventPts s = let [p,q] = L.sortBy ordPoints [s^.start.core,s^.end.core]                in [Event p (Start $ s :| []), Event q (End s)]  -- | Group the segments with the intersection points-merge :: Ord r =>  [IntersectionPoint p r] -> Intersections p r+merge :: (Ord r, Fractional r) =>  [IntersectionPoint p r] -> Intersections p r merge = foldr (\(IntersectionPoint p a) -> M.insertWith (<>) p a) M.empty  -- | Group the startpoints such that segments with the same start point@@ -113,27 +123,6 @@                 _        -> []  ------------------------------------------------------------------------------------ | Compare based on the x-coordinate of the intersection with the horizontal--- line through y-ordAt   :: (Fractional r, Ord r) => r -> Compare (LineSegment 2 p r)-ordAt y = comparing (xCoordAt y)---- | Given a y coord and a line segment that intersects the horizontal line--- through y, compute the x-coordinate of this intersection point.------ note that we will pretend that the line segment is closed, even if it is not-xCoordAt             :: (Fractional r, Ord r) => r -> LineSegment 2 p r -> r-xCoordAt y (LineSegment' (Point2 px py :+ _) (Point2 qx qy :+ _))-      | py == qy     = px `max` qx  -- s is horizontal, and since it by the-                                    -- precondition it intersects the sweep-                                    -- line, we return the x-coord of the-                                    -- rightmost endpoint.-      | otherwise    = px + alpha * (qx - px)-  where-    alpha = (y - py) / (qy - py)---------------------------------------------------------------------------------- -- * The Main Sweep  type EventQueue      p r = EQ.Set (Event p r)@@ -163,7 +152,7 @@     -- starting segments, exluding those that have an open starting point     starts'  = filter (isClosedStart p) starts     toReport = case starts' ++ contains' of-                 (_:_:_) -> [IntersectionPoint p $ associated (starts' <> ends) contains]+                 (_:_:_) -> [mkIntersectionPoint p (starts' <> ends) contains]                  _       -> []      -- new status structure@@ -203,7 +192,7 @@   -- ss { SS.nav = ordAtNav $ p^.yCoord } `SS.join` hors   where     (hors',rest) = L.partition isHorizontal xs-    ss           = SS.fromListBy (ordAt $ maxY xs) rest+    ss           = SS.fromListBy (ordAtY $ maxY xs) rest     hors         = SS.fromListBy (comparing rightEndpoint) hors'      isHorizontal s  = s^.start.core.yCoord == s^.end.core.yCoord

@@ -0,0 +1,200 @@+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE TemplateHaskell #-}+--------------------------------------------------------------------------------+-- |+-- Module      :  Algorithms.Geometry.LineSegmentIntersection.Types+-- Copyright   :  (C) Frank Staals+-- License     :  see the LICENSE file+-- Maintainer  :  Frank Staals+--------------------------------------------------------------------------------+module Algorithms.Geometry.LineSegmentIntersection.TypesNoExt where++-- import           Algorithms.DivideAndConquer+import           Control.DeepSeq+import           Control.Lens+import           Data.Ext+import           Data.Bifunctor+import           Data.Geometry.Interval+import           Data.Geometry.LineSegment+import           Data.Geometry.Point+import qualified Data.Map as Map+import qualified Data.Set as Set+import           Data.Ord (comparing, Down(..))+import           GHC.Generics+import           Data.Vinyl.CoRec+import           Data.Vinyl+import           Data.Intersection+++----------------------------------------------------------------------------------+++-- FIXME: What do we do when one segmnet lies *on* the other one. For+-- the short segment it should be an "around start", but then the+-- startpoints do not match.+--+-- for the long one it's an "on" segment, but they do not intersect+++-- | Assumes that two segments have the same start point+newtype AroundStart a = AroundStart a deriving (Show,Read,NFData)++instance Eq r => Eq (AroundStart (LineSegment 2 p r)) where+  -- | equality on endpoint+  (AroundStart s) == (AroundStart s') = s^.end.core == s'^.end.core++instance (Ord r, Num r) => Ord (AroundStart (LineSegment 2 p r)) where+  -- | ccw ordered around their suposed common startpoint+  (AroundStart s) `compare` (AroundStart s') =+    ccwCmpAround (s^.start.core) (s^.end.core)  (s'^.end.core)++----------------------------------------++-- | Assumes that two segments have the same end point+newtype AroundEnd a = AroundEnd a deriving (Show,Read,NFData)++instance Eq r => Eq (AroundEnd (LineSegment 2 p r)) where+  -- | equality on endpoint+  (AroundEnd s) == (AroundEnd s') = s^.start.core == s'^.start.core++instance (Ord r, Num r) => Ord (AroundEnd (LineSegment 2 p r)) where+  -- | ccw ordered around their suposed common end point+  (AroundEnd s) `compare` (AroundEnd s') =+    ccwCmpAround (s^.end.core) (s^.start.core)  (s'^.start.core)++--------------------------------------------------------------------------------++-- | Assumes that two segments intersect in a single point.+newtype AroundIntersection a = AroundIntersection a deriving (Show,Read,NFData)++instance Eq r => Eq (AroundIntersection (LineSegment 2 p r)) where+  -- | equality ignores the p type+  (AroundIntersection s) == (AroundIntersection s') = first (const ()) s == first (const ()) s'++instance (Ord r, Fractional r) => Ord (AroundIntersection (LineSegment 2 p r)) where+  -- | ccw ordered around their common intersection point.+  l@(AroundIntersection s) `compare` r@(AroundIntersection s') = match (s `intersect` s') $+        H (\NoIntersection     -> error "AroundIntersection: segments do not intersect!")+     :& H (\p                  -> cmpAroundP p s s')+     :& H (\_                  -> (squaredLength s) `compare` (squaredLength s'))+                                 -- if s and s' just happen to be the same length but+                                 -- intersect in different behaviour from using (==).+                                 -- but that situation doese not satisfy the precondition+                                 -- of aroundIntersection anyway.+     :& RNil+    where+      squaredLength (LineSegment' a b) = squaredEuclideanDist (a^.core) (b^.core)++-- | compare around p+cmpAroundP        :: (Ord r, Num r) => Point 2 r -> LineSegment 2 p r -> LineSegment 2 p r -> Ordering+cmpAroundP p s s' = ccwCmpAround p (s^.start.core)  (s'^.start.core)+++-- seg1 = ClosedLineSegment (ext $ Point2 0 0) (ext $ Point2 0 10)+-- seg2 = ClosedLineSegment (ext $ Point2 0 0) (ext $ Point2 0 10)+++--------------------------------------------------------------------------------++++-- | The line segments that contain a given point p may either have p+-- as the endpoint or have p in their interior.+--+-- if somehow the segment is degenerate, and p is both the start and+-- end it is reported only as the start point.+data Associated p r =+  Associated { _startPointOf :: Set.Set (AroundEnd (LineSegment 2 p r))+             -- ^ segments for which the intersection point is the+             -- start point (i.e. s^.start.core == p)+             , _endPointOf   :: Set.Set (AroundStart (LineSegment 2 p r))+             -- ^ segments for which the intersection point is the end+             -- point (i.e. s^.end.core == p)+             , _interiorTo   :: Set.Set (AroundIntersection (LineSegment 2 p r))+             } deriving stock (Show, Read, Generic, Eq)++makeLenses ''Associated++++-- | Reports whether this associated has any interior intersections+--+-- \(O(1)\)+isInteriorIntersection :: Associated p r -> Bool+isInteriorIntersection = not . null . _interiorTo+++-- | test if the given segment has p as its endpoint, an construct the+-- appropriate associated representing that.+--+-- pre: p intersects the segment+mkAssociated                :: (Ord r, Fractional r)+                            => Point 2 r -> LineSegment 2 p r -> Associated p r+mkAssociated p s@(LineSegment a b)+  | p == a^.unEndPoint.core = mempty&startPointOf .~  Set.singleton (AroundEnd s)+  | p == b^.unEndPoint.core = mempty&endPointOf   .~  Set.singleton (AroundStart s)+  | otherwise               = mempty&interiorTo   .~  Set.singleton (AroundIntersection s)+++-- | test if the given segment has p as its endpoint, an construct the+-- appropriate associated representing that.+--+-- If p is not one of the endpoints we concstruct an empty Associated!+--+mkAssociated'     :: (Ord r, Fractional r) => Point 2 r -> LineSegment 2 p r -> Associated p r+mkAssociated' p s = (mkAssociated p s)&interiorTo .~ mempty++instance (Ord r, Fractional r) => Semigroup (Associated p r) where+  (Associated ss es is) <> (Associated ss' es' is') =+    Associated (ss <> ss') (es <> es') (is <> is')++instance (Ord r, Fractional r) => Monoid (Associated p r) where+  mempty = Associated mempty mempty mempty++instance (NFData p, NFData r) => NFData (Associated p r)++-- | For each intersection point the segments intersecting there.+type Intersections p r = Map.Map (Point 2 r) (Associated p r)++-- | An intersection point together with all segments intersecting at+-- this point.+data IntersectionPoint p r =+  IntersectionPoint { _intersectionPoint :: !(Point 2 r)+                    , _associatedSegs    :: !(Associated p r)+                    } deriving (Show,Read,Eq,Generic)+makeLenses ''IntersectionPoint++instance (NFData p, NFData r) => NFData (IntersectionPoint p r)+++-- sameOrder           :: (Ord r, Num r, Eq p) => Point 2 r+--                     -> [LineSegment 2 p r] -> [LineSegment 2 p r] -> Bool+-- sameOrder c ss ss' = f ss == f ss'+--   where+--     f = map (^.extra) . sortAround' (ext c) . map (\s -> s^.end.core :+ s)+++++-- | Given a point p, and a bunch of segments that suposedly intersect+-- at p, correctly categorize them.+mkIntersectionPoint         :: (Ord r, Fractional r)+                            => Point 2 r+                            -> [LineSegment 2 p r] -- ^ uncategorized+                            -> [LineSegment 2 p r] -- ^ segments we know contain p,+                            -> IntersectionPoint p r+mkIntersectionPoint p as cs = IntersectionPoint p $ foldMap (mkAssociated p) $ as <> cs++  -- IntersectionPoint p+  --                           $ Associated mempty mempty (Set.fromAscList cs')+  --                           <> foldMap (mkAssociated p) as+  -- where+  --   cs' = map AroundIntersection . List.sortBy (cmpAroundP p) $ cs+  -- -- TODO: In the bentley ottman algo we already know the sorted order of the segments+  -- -- so we can likely save the additional sort++++-- | An ordering that is decreasing on y, increasing on x+ordPoints     :: Ord r => Point 2 r -> Point 2 r -> Ordering+ordPoints a b = let f p = (Down $ p^.yCoord, p^.xCoord) in comparing f a b

@@ -2,7 +2,6 @@ {-# LANGUAGE TemplateHaskell     #-} module Algorithms.Geometry.PolygonTriangulation.MakeMonotoneOld where -import Algorithms.Geometry.LineSegmentIntersection.BentleyOttmann (ordAt, xCoordAt) import Algorithms.Geometry.PolygonTriangulation.Types  import           Control.Lens@@ -16,7 +15,7 @@ import           Data.Ext import qualified Data.Foldable                         as F import           Data.Geometry.LineSegment-import           Data.Geometry.PlanarSubdivision.Basic+import           Data.Geometry.PlanarSubdivision import           Data.Geometry.Point import           Data.Geometry.Polygon import qualified Data.IntMap                           as IntMap@@ -144,11 +143,11 @@ -- pieces. -- -- running time: \(O(n\log n)\)-makeMonotone      :: (Fractional r, Ord r)+makeMonotone      :: forall proxy s t p r. (Fractional r, Ord r)                   => proxy s -> Polygon t p r                   -> PlanarSubdivision s p PolygonEdgeType PolygonFaceData r-makeMonotone px pg = let (e:es) = listEdges pg-                     in constructSubdivision px e es (computeDiagonals pg)+makeMonotone _ pg = let (e:es) = listEdges pg+                    in constructSubdivision @s e es (computeDiagonals pg)  type Sweep p r = WriterT (DList.DList (LineSegment 2 Int r))                    (StateT (StatusStruct r)@@ -181,11 +180,11 @@  insertAt   :: (Ord r, Fractional r) => Point 2 r -> LineSegment 2 q r            -> OrdSeq (LineSegment 2 q r) -> OrdSeq (LineSegment 2 q r)-insertAt v = SS.insertBy (ordAt $ v^.yCoord)+insertAt v = SS.insertBy (ordAtY $ v^.yCoord)  deleteAt   :: (Fractional r, Ord r) => Point 2 r -> LineSegment 2 p r            -> OrdSeq (LineSegment 2 p r) -> OrdSeq (LineSegment 2 p r)-deleteAt v = SS.deleteAllBy (ordAt $ v^.yCoord)+deleteAt v = SS.deleteAllBy (ordAtY $ v^.yCoord)   handleStart              :: (Fractional r, Ord r)

@@ -2,10 +2,29 @@  * Changelog +** 0.14++- Allow the associated/extra data of linesegments and intervals to+  differ when testing for intersections.+- Intersection testing between line segments and rectangles+- Testing if lines and/or line segments intersect no longer requires a+  Fractional constraint; Num is sufficient. However, in turn we now do+  need Ord rather than just Eq. That seemed a worthwile tradeoff though.+- Cleaning up the public API by hiding several internal modules.+- Introduced the 'HasSquaredEuclideanDistance' class describing+  geometry types for which we can compute the squared distance from a+  point to a geometry, and added instances for some of the basic+  geometries.+- Fixed a bug in computing lengths to open line segments.+- Removed some proxy arguments, in e.g. Data.Geometry.Point.coord,+  rather than take a Proxy to specify which coordinate we want, use+  type applications.+- Support for GHC 9.0 and 9.2+- Better support for open-ended line segments in the Bentley Ottmann+  line segment intersection algorithm.+ ** 0.13 -- Implementation of Logaritmic Method, wich allows us to transform a-  static data structure into an insertion only data structure - Moved 'intersects' from the HasIntersectionWith class into a new   class IsIntersectableWith. This allows separate (weaker) constraints   for checking *if* geometries intersect rather than computing exact

@@ -2,10 +2,29 @@  * Changelog +** 0.14++- Allow the associated/extra data of linesegments and intervals to+  differ when testing for intersections.+- Intersection testing between line segments and rectangles+- Testing if lines and/or line segments intersect no longer requires a+  Fractional constraint; Num is sufficient. However, in turn we now do+  need Ord rather than just Eq. That seemed a worthwile tradeoff though.+- Cleaning up the public API by hiding several internal modules.+- Introduced the 'HasSquaredEuclideanDistance' class describing+  geometry types for which we can compute the squared distance from a+  point to a geometry, and added instances for some of the basic+  geometries.+- Fixed a bug in computing lengths to open line segments.+- Removed some proxy arguments, in e.g. Data.Geometry.Point.coord,+  rather than take a Proxy to specify which coordinate we want, use+  type applications.+- Support for GHC 9.0 and 9.2+- Better support for open-ended line segments in the Bentley Ottmann+  line segment intersection algorithm.+ ** 0.13 -- Implementation of Logaritmic Method, wich allows us to transform a-  static data structure into an insertion only data structure - Moved 'intersects' from the HasIntersectionWith class into a new   class IsIntersectableWith. This allows separate (weaker) constraints   for checking *if* geometries intersect rather than computing exact

@@ -32,6 +32,8 @@           , "DeriveGeneric"           , "FlexibleInstances"           , "FlexibleContexts"+          , "DerivingStrategies"+          , "DerivingVia"           ]  files :: [String]

@@ -1,6 +1,6 @@ cabal-version:       2.4 name:                hgeometry-version:             0.13+version:             0.14 synopsis:            Geometric Algorithms, Data structures, and Data types. description:   HGeometry provides some basic geometry types, and geometric algorithms and@@ -50,17 +50,14 @@                     Data.Geometry.Vector                     Data.Geometry.Vector.VectorFixed                     Data.Geometry.Vector.VectorFamily-                    Data.Geometry.Vector.VectorFamilyPeano                      Data.Geometry.Matrix                      -- Data.Geometry.Vector.Vinyl                     Data.Geometry.Interval-                    Data.Geometry.Interval.Util                     Data.Geometry.Point                      Data.Geometry.Line-                    Data.Geometry.Line.Internal                     Data.Geometry.LineSegment                     Data.Geometry.LineSegment.Internal                     Data.Geometry.SubLine@@ -96,13 +93,10 @@                      Data.Geometry.PlanarSubdivision                     Data.Geometry.PlanarSubdivision.Raw-                    Data.Geometry.PlanarSubdivision.Basic-                    Data.Geometry.PlanarSubdivision.Merge                     Data.Geometry.PlanarSubdivision.Dynamic                     Data.Geometry.PlanarSubdivision.TreeRep                      Data.Geometry.Arrangement-                    Data.Geometry.Arrangement.Internal                      Data.Geometry.RangeTree                     Data.Geometry.RangeTree.Measure@@ -186,9 +180,10 @@                     Algorithms.Geometry.SSSP                     Algorithms.Geometry.SSSP.Naive +                    Algorithms.Geometry.RayShooting.Naive+                     -- * Embedded Planar Graphs                     Data.PlaneGraph-                    Data.PlaneGraph.Core                     Data.PlaneGraph.AdjRep                     Data.PlaneGraph.IO @@ -210,12 +205,19 @@                      Algorithms.Geometry.WSPD.Types +                    Data.Geometry.Vector.VectorFamilyPeano+                     Data.Geometry.Point.Internal                     Data.Geometry.Point.Orientation                     Data.Geometry.Point.Quadrants                     Data.Geometry.Point.Orientation.Degenerate                     Data.Geometry.Point.Class +                    Data.Geometry.Line.Internal+++                    Data.Geometry.Interval.Util+                     Algorithms.Geometry.SoS.Expr                     Algorithms.Geometry.SoS.AsPoint                     Algorithms.Geometry.SoS.Internal@@ -223,8 +225,13 @@                     Algorithms.Geometry.SoS.Determinant                     Algorithms.Geometry.SoS.Sign +                    Data.PlaneGraph.Core +                    Data.Geometry.Arrangement.Internal +                    Data.Geometry.PlanarSubdivision.Basic+                    Data.Geometry.PlanarSubdivision.Merge+   -- other-extensions:   build-depends:                 base                    >= 4.11      &&     < 5@@ -233,6 +240,7 @@               , bifunctors              >= 4.1               , bytestring              >= 0.10               , containers              >= 0.5.9+              -- , multi-containers        >= 0.2               , dlist                   >= 0.7               , lens                    >= 4.2               , semigroupoids           >= 5@@ -259,7 +267,7 @@                , vector                  >= 0.11               , data-clist              >= 0.1.2.3-              , vector-circular         >= 0.1.2+              , vector-circular         >= 0.1.4               , nonempty-vector         >= 0.2.0.0               , text                    >= 1.1.1.0               , vector-algorithms@@ -300,6 +308,8 @@                     , DeriveFoldable                     , DeriveTraversable                     , DeriveGeneric+                    , DerivingStrategies+                    , DerivingVia                       , FlexibleInstances@@ -329,7 +339,7 @@                  Algorithms.Geometry.ConvexHull.Bench                  Algorithms.Geometry.ConvexHull.GrahamV2                  Algorithms.Geometry.ConvexHull.GrahamFam-                 Algorithms.Geometry.ConvexHull.GrahamFamPeano+                 -- Algorithms.Geometry.ConvexHull.GrahamFamPeano                  Algorithms.Geometry.ConvexHull.GrahamFixed                  Data.Geometry.Vector.VectorFamily6                  Algorithms.Geometry.ConvexHull.GrahamFam6@@ -341,6 +351,8 @@                   Algorithms.Geometry.LineSegmentIntersection.Bench                  Algorithms.Geometry.LineSegmentIntersection.BentleyOttmannOld+                 Algorithms.Geometry.LineSegmentIntersection.BentleyOttmannNoExt+                 Algorithms.Geometry.LineSegmentIntersection.TypesNoExt                   Algorithms.Geometry.PolygonTriangulation.Bench                  Algorithms.Geometry.PolygonTriangulation.MakeMonotoneOld@@ -401,3 +413,5 @@                     , FlexibleInstances                     , FlexibleContexts                     , MultiParamTypeClasses+                    , DerivingStrategies+                    , DeriveGeneric

@@ -115,16 +115,15 @@ -- | convert the triangulation into a planarsubdivision -- -- running time: \(O(n)\).-toPlanarSubdivision    :: (Ord r, Fractional r)-                       => proxy s -> Triangulation p r -> PlanarSubdivision s p () () r-toPlanarSubdivision px = fromPlaneGraph . toPlaneGraph px+toPlanarSubdivision :: forall s p r. (Ord r, Fractional r)+                    => Triangulation p r -> PlanarSubdivision s p () () r+toPlanarSubdivision = fromPlaneGraph . toPlaneGraph  -- | convert the triangulation into a plane graph -- -- running time: \(O(n)\).-toPlaneGraph    :: forall proxy s p r.-                   proxy s -> Triangulation p r -> PG.PlaneGraph s p () () r-toPlaneGraph _ tr = PG.PlaneGraph $ g&PPG.vertexData .~ vtxData+toPlaneGraph    :: forall s p r. Triangulation p r -> PG.PlaneGraph s p () () r+toPlaneGraph tr = PG.PlaneGraph $ g&PPG.vertexData .~ vtxData   where     g       = PPG.fromAdjacencyLists . V.toList . V.imap f $ tr^.neighbours     f i adj = (VertexId i, C.leftElements $ VertexId <$> adj) -- report in CCW order

@@ -19,7 +19,6 @@ import           Data.Geometry import qualified Data.List.NonEmpty as NonEmpty import           Data.PlaneGraph-import           Data.Proxy import           Data.Tree  --------------------------------------------------------------------------------@@ -39,7 +38,7 @@     -- since we care only about the relative order of the edges we can use the     -- squared Euclidean distance rather than the Euclidean distance, thus     -- avoiding the Floating constraint-    g = withEdgeDistances squaredEuclideanDist . toPlaneGraph (Proxy :: Proxy MSTW)+    g = withEdgeDistances squaredEuclideanDist . toPlaneGraph @MSTW       . delaunayTriangulation $ pts     t = mst $ g^.graph 

@@ -6,28 +6,30 @@ -- Maintainer  :  Frank Staals -------------------------------------------------------------------------------- module Algorithms.Geometry.LineSegmentIntersection-  ( hasInteriorIntersections-  , hasSelfIntersections+  ( BooleanSweep.hasIntersections+  , BO.intersections+  , BO.interiorIntersections   , Intersections   , Associated(..)-  , IntersectionPoint(..)-  , isEndPointIntersection-  , associated-  , Compare+  , IntersectionPoint(..), mkIntersectionPoint+  -- , isInteriorIntersection+  , hasSelfIntersections   ) where  import qualified Algorithms.Geometry.LineSegmentIntersection.BentleyOttmann as BO+import qualified Algorithms.Geometry.LineSegmentIntersection.BooleanSweep as BooleanSweep import           Algorithms.Geometry.LineSegmentIntersection.Types+import           Data.Ext (ext) import           Data.Geometry.LineSegment import           Data.Geometry.Polygon --- Tests if there are any interior intersections.------ | \(O(n \log n)\)-hasInteriorIntersections :: (Ord r, Fractional r)-                         => [LineSegment 2 p r] -> Bool-hasInteriorIntersections = not . null . BO.interiorIntersections+import qualified Data.Map as Map --- | \(O(n \log n)\)+-- | Test if the polygon has self intersections.+--+-- \(O(n \log n)\) hasSelfIntersections :: (Ord r, Fractional r) => Polygon t p r -> Bool-hasSelfIntersections = hasInteriorIntersections . listEdges+hasSelfIntersections = not . Map.null . BO.interiorIntersections . map ext . listEdges+-- hasSelfIntersections :: (Ord r, Num r) => Polygon t p r -> Bool+-- hasSelfIntersections = BooleanSweep.hasIntersections . listEdges+-- FIXME: fix the open/closed bug, then switch to a boolean sweep based version

@@ -13,15 +13,14 @@ module Algorithms.Geometry.LineSegmentIntersection.BentleyOttmann   ( intersections   , interiorIntersections-    -- FIXME: Move ordAt and xCoordAt to Data.Geometry.LineSegment?-  , ordAt-  , xCoordAt   ) where  import           Algorithms.Geometry.LineSegmentIntersection.Types import           Control.Lens hiding (contains)+import           Data.Coerce import           Data.Ext import qualified Data.Foldable as F+import           Data.Function (on) import           Data.Geometry.Interval import           Data.Geometry.LineSegment import           Data.Geometry.Point@@ -32,49 +31,105 @@ import qualified Data.Map as M import           Data.Maybe import           Data.Ord (Down(..), comparing)+import qualified Data.Set as EQ -- event queue import qualified Data.Set as SS -- status struct+import qualified Data.Set as Set import qualified Data.Set.Util as SS -- status struct-import qualified Data.Set as EQ -- event queue import           Data.Vinyl import           Data.Vinyl.CoRec- --------------------------------------------------------------------------------  -- | Compute all intersections -- -- \(O((n+k)\log n)\), where \(k\) is the number of intersections.-intersections    :: (Ord r, Fractional r)-                 => [LineSegment 2 p r] -> Intersections p r-intersections ss = merge $ sweep pts SS.empty+intersections    :: forall p r e. (Ord r, Fractional r)+                 => [LineSegment 2 p r :+ e] -> Intersections p r e+intersections ss = fmap unflipSegs . merge $ sweep pts SS.empty   where-    pts = EQ.fromAscList . groupStarts . L.sort . concatMap asEventPts $ ss+    pts = EQ.fromAscList . groupStarts . L.sort . concatMap (asEventPts . tagFlipped) $ ss  -- | Computes all intersection points p s.t. p lies in the interior of at least -- one of the segments. -- --  \(O((n+k)\log n)\), where \(k\) is the number of intersections. interiorIntersections :: (Ord r, Fractional r)-                       => [LineSegment 2 p r] -> Intersections p r-interiorIntersections = M.filter (not . isEndPointIntersection) . intersections+                       => [LineSegment 2 p r :+ e] -> Intersections p r e+interiorIntersections = M.filter isInteriorIntersection . intersections +--------------------------------------------------------------------------------+-- * Flipping and unflipping++data Flipped = NotFlipped | Flipped deriving (Show,Eq)++-- | Make sure the 'start' endpoint occurs before the end-endpoints in+-- terms of the sweep order.+tagFlipped   :: Ord r => LineSegment 2 p r :+ e -> LineSegment 2 p r :+ (e :+ Flipped)+tagFlipped s = case (s^.core.start.core) `ordPoints` (s^.core.end.core) of+                 GT -> s&core  %~ flipSeg+                        &extra %~ (:+ Flipped)+                 _  -> s&extra %~ (:+ NotFlipped)++-- | Flips the segment+flipSeg     :: LineSegment d p r -> LineSegment d p r+flipSeg seg = seg&start .~ (seg^.end)+                 &end   .~ (seg^.start)++-- | Unflips the segments in an associated.+unflipSegs                       :: (Fractional r, Ord r)+                                 => Associated p r (e :+ Flipped) -> Associated p r e+unflipSegs (Associated ss es is) =+    Associated (dropFlipped ss1 <> unflipSegs' es')+               (dropFlipped es1 <> unflipSegs' ss')+               (dropFlipped is1 <> unflipSegs' is')+  where+    (ss',ss1) = Set.partition (\(AroundEnd          s) -> isFlipped s) ss+    (es',es1) = Set.partition (\(AroundStart        s) -> isFlipped s) es+    (is',is1) = Set.partition (\(AroundIntersection s) -> isFlipped s) is++    isFlipped s = Flipped == s^.extra.extra++    -- | For segments that are not acutally flipped, we can just drop the flipped bit+    dropFlipped :: Functor f+                => Set.Set (f (LineSegment 2 p r :+ (e :+ Flipped)))+                -> Set.Set (f (LineSegment 2 p r :+ e))+    dropFlipped = Set.mapMonotonic (fmap dropFlip)++    -- For flipped segs we unflip them (and appropriately coerce the+    -- so that they remain in the same order. I.e. if they were sorted+    -- around the start point they are now sorted around the endpoint.+    unflipSegs' :: ( Functor f+                   , Coercible (f (LineSegment 2 p r :+ e)) (g (LineSegment 2 p r :+ e))+                   )+                => Set.Set (f (LineSegment 2 p r :+ (e :+ Flipped)))+                -> Set.Set (g (LineSegment 2 p r :+ e))+    unflipSegs' = Set.mapMonotonic (coerce . fmap unflip)++    unflip   (s :+ (e :+ _)) = flipSeg s :+ e+    dropFlip (s :+ (e :+ _)) = s :+ e++--------------------------------------------------------------------------------+ -- | Computes the event points for a given line segment-asEventPts   :: Ord r => LineSegment 2 p r -> [Event p r]-asEventPts s = let [p,q] = L.sortBy ordPoints [s^.start.core,s^.end.core]-               in [Event p (Start $ s :| []), Event q (End s)]+asEventPts   :: LineSegment 2 p r :+ e -> [Event p r e]+asEventPts s = [ Event (s^.core.start.core) (Start $ s :| [])+               , Event (s^.core.end.core)   (End s)+               ]  -- | Group the segments with the intersection points-merge :: Ord r =>  [IntersectionPoint p r] -> Intersections p r+merge :: (Ord r, Fractional r) =>  [IntersectionPoint p r e] -> Intersections p r e merge = foldr (\(IntersectionPoint p a) -> M.insertWith (<>) p a) M.empty  -- | Group the startpoints such that segments with the same start point -- correspond to one event.-groupStarts                          :: Eq r => [Event p r] -> [Event p r]+groupStarts                          :: Eq r => [Event p r e] -> [Event p r e] groupStarts []                       = [] groupStarts (Event p (Start s) : es) = Event p (Start ss) : groupStarts rest   where     (ss',rest) = L.span sameStart es-    -- sort the segs on lower endpoint-    ss         = let (x:|xs) = s in x :| (xs ++ concatMap startSegs ss')+    -- FIXME: this seems to keep the segments on decreasing y, increasing x. shouldn't we+    -- sort them cyclically around p instead?+    ss         = let (x:|xs) = s+                 in x :| (xs ++ concatMap startSegs ss')      sameStart (Event q (Start _)) = p == q     sameStart _                   = False@@ -99,84 +154,68 @@   (End _)      `compare` _            = GT  -- | The actual event consists of a point and its type-data Event p r = Event { eventPoint :: !(Point 2 r)-                       , eventType  :: !(EventType (LineSegment 2 p r))-                       } deriving (Show,Eq)+data Event p r e = Event { eventPoint :: !(Point 2 r)+                         , eventType  :: !(EventType (LineSegment 2 p r :+ e))+                         } deriving (Show,Eq) -instance Ord r => Ord (Event p r) where+instance Ord r => Ord (Event p r e) where   -- decreasing on the y-coord, then increasing on x-coord, and increasing on event-type   (Event p s) `compare` (Event q t) = case ordPoints p q of                                         EQ -> s `compare` t                                         x  -> x --- | An ordering that is decreasing on y, increasing on x-ordPoints     :: Ord r => Point 2 r -> Point 2 r -> Ordering-ordPoints a b = let f p = (Down $ p^.yCoord, p^.xCoord) in comparing f a b- -- | Get the segments that start at the given event point-startSegs   :: Event p r -> [LineSegment 2 p r]+startSegs   :: Event p r e -> [LineSegment 2 p r :+ e] startSegs e = case eventType e of                 Start ss -> NonEmpty.toList ss                 _        -> []  -------------------------------------------------------------------------------- --- | Compare based on the x-coordinate of the intersection with the horizontal--- line through y-ordAt   :: (Fractional r, Ord r) => r -> Compare (LineSegment 2 p r)-ordAt y = comparing (xCoordAt y) --- | Given a y coord and a line segment that intersects the horizontal line--- through y, compute the x-coordinate of this intersection point.------ note that we will pretend that the line segment is closed, even if it is not-xCoordAt             :: (Fractional r, Ord r) => r -> LineSegment 2 p r -> r-xCoordAt y (LineSegment' (Point2 px py :+ _) (Point2 qx qy :+ _))-      | py == qy     = px `max` qx  -- s is horizontal, and since it by the-                                    -- precondition it intersects the sweep-                                    -- line, we return the x-coord of the-                                    -- rightmost endpoint.-      | otherwise    = px + alpha * (qx - px)-  where-    alpha = (y - py) / (qy - py)- -------------------------------------------------------------------------------- -- * The Main Sweep -type EventQueue      p r = EQ.Set (Event p r)-type StatusStructure p r = SS.Set (LineSegment 2 p r)+type EventQueue      p r e = EQ.Set (Event p r e)+type StatusStructure p r e = SS.Set (LineSegment 2 p r :+ e)  -- | Run the sweep handling all events sweep       :: (Ord r, Fractional r)-            => EventQueue p r -> StatusStructure p r -> [IntersectionPoint p r]+            => EventQueue p r e -> StatusStructure p r e -> [IntersectionPoint p r e] sweep eq ss = case EQ.minView eq of     Nothing      -> []     Just (e,eq') -> handle e eq' ss -isClosedStart                     :: Eq r => Point 2 r -> LineSegment 2 p r -> Bool-isClosedStart p (LineSegment s e)-  | p == s^.unEndPoint.core       = isClosed s-  | otherwise                     = isClosed e- -- | Handle an event point-handle                           :: forall r p. (Ord r, Fractional r)-                                 => Event p r -> EventQueue p r -> StatusStructure p r-                                 -> [IntersectionPoint p r]+handle                           :: forall r p e. (Ord r, Fractional r)+                                 => Event p r e -> EventQueue p r e -> StatusStructure p r e+                                 -> [IntersectionPoint p r e] handle e@(eventPoint -> p) eq ss = toReport <> sweep eq' ss'   where     starts                   = startSegs e     (before,contains',after) = extractContains p ss     (ends,contains)          = L.partition (endsAt p) contains'     -- starting segments, exluding those that have an open starting point-    starts'  = filter (isClosedStart p) starts-    toReport = case starts' ++ contains' of-                 (_:_:_) -> [IntersectionPoint p $ associated (starts' <> ends) contains]+    -- starts' = filter (isClosedStart p) starts+    starts' = shouldReport p $ SS.toAscList newSegs++    -- If we just inserted open-ended segments that start here, then+    -- don't consider them to be "contained" segments.+    pureContains = filter (\(LineSegment s _ :+ _) ->+                              not $ isOpen s && p == s^.unEndPoint.core) contains++    -- any (closed) ending segments at this event point.+    closedEnds = filter (\(LineSegment _ e' :+ _) -> isClosed e') ends++    toReport = case starts' <> closedEnds <> pureContains of+                 (_:_:_) -> [mkIntersectionPoint p (starts' <> closedEnds) pureContains]                  _       -> []      -- new status structure     ss' = before `SS.join` newSegs `SS.join` after     newSegs = toStatusStruct p $ starts ++ contains +     -- the new eeventqueue     eq' = foldr EQ.insert eq es     -- the new events:@@ -191,54 +230,145 @@      app f x y = do { x' <- x ; y' <- y ; f x' y'} +-- | given the starting point p, and the segments that either start in+-- p, or continue in p, in left to right order along a line just+-- epsilon below p, figure out which segments we should report as+-- intersecting at p.+--+-- in partcular; those that:+-- - have a closed endpoint at p+-- - those that have an open endpoint at p and have an intersection+--   with a segment eps below p. Those segments thus overlap wtih+--   their predecessor or successor in the cyclic order.+shouldReport   :: (Ord r, Num r)+               => Point 2 r -> [LineSegment 2 p r :+ e] -> [LineSegment 2 p r :+ e]+shouldReport _ = overlapsOr (\(LineSegment s _ :+ _) -> isClosed s)+                            (\(s :+ _) (s2 :+ _) -> s `intersects` s2)+ -- | split the status structure, extracting the segments that contain p. -- the result is (before,contains,after) extractContains      :: (Fractional r, Ord r)-                     => Point 2 r -> StatusStructure p r-                     -> (StatusStructure p r, [LineSegment 2 p r], StatusStructure p r)+                     => Point 2 r -> StatusStructure p r e+                     -> (StatusStructure p r e, [LineSegment 2 p r :+ e], StatusStructure p r e) extractContains p ss = (before, F.toList mid1 <> F.toList mid2, after)   where-    (before, mid1, after') = SS.splitOn (xCoordAt $ p^.yCoord) (p^.xCoord) ss+    (before, mid1, after') = SS.splitOn (xCoordAt' $ p^.yCoord) (p^.xCoord) ss     -- Make sure to also select the horizontal segments containing p-    (mid2, after) = SS.spanAntitone (intersects p) after'-+    (mid2, after) = SS.spanAntitone (intersects p . view core) after'+    xCoordAt' y sa = xCoordAt y (sa^.core)  -- | Given a point and the linesegements that contain it. Create a piece of -- status structure for it. toStatusStruct      :: (Fractional r, Ord r)-                    => Point 2 r -> [LineSegment 2 p r] -> StatusStructure p r+                    => Point 2 r -> [LineSegment 2 p r :+ e] -> StatusStructure p r e toStatusStruct p xs = ss `SS.join` hors   -- ss { SS.nav = ordAtNav $ p^.yCoord } `SS.join` hors   where     (hors',rest) = L.partition isHorizontal xs-    ss           = SS.fromListBy (ordAt $ maxY xs) rest+    ss           = SS.fromListBy (ordAtY' $ maxY xs) rest     hors         = SS.fromListBy (comparing rightEndpoint) hors' -    isHorizontal s  = s^.start.core.yCoord == s^.end.core.yCoord+    isHorizontal s  = s^.core.start.core.yCoord == s^.core.end.core.yCoord +    ordAtY' q sa sb = ordAtY q (sa^.core) (sb^.core)+     -- find the y coord of the first interesting thing below the sweep at y     maxY = maximum . filter (< p^.yCoord)-         . concatMap (\s -> [s^.start.core.yCoord,s^.end.core.yCoord])+         . concatMap (\s -> [s^.core.start.core.yCoord,s^.core.end.core.yCoord])  -- | Get the right endpoint of a segment-rightEndpoint   :: Ord r => LineSegment 2 p r -> r-rightEndpoint s = (s^.start.core.xCoord) `max` (s^.end.core.xCoord)+rightEndpoint   :: Ord r => LineSegment 2 p r :+ e -> r+rightEndpoint s = (s^.core.start.core.xCoord) `max` (s^.core.end.core.xCoord)  -- | Test if a segment ends at p-endsAt                      :: Ord r => Point 2 r -> LineSegment 2 p r -> Bool-endsAt p (LineSegment' a b) = all (\q -> ordPoints (q^.core) p /= GT) [a,b]+endsAt                                  :: Eq r => Point 2 r -> LineSegment 2 p r :+ e -> Bool+endsAt p (LineSegment' _ (b :+ _) :+ _) = p == b+  -- all (\q -> ordPoints (q^.core) p /= GT) [a,b]  -------------------------------------------------------------------------------- -- * Finding New events  -- | Find all events findNewEvent       :: (Ord r, Fractional r)-                   => Point 2 r -> LineSegment 2 p r -> LineSegment 2 p r-                   -> Maybe (Event p r)-findNewEvent p l r = match (l `intersect` r) $+                   => Point 2 r -> LineSegment 2 p r :+ e -> LineSegment 2 p r :+ e+                   -> Maybe (Event p r e)+findNewEvent p l r = match ((l^.core) `intersect` (r^.core)) $      H (const Nothing) -- NoIntersection   :& H (\q -> if ordPoints q p == GT then Just (Event q Intersection)                                      else Nothing)   :& H (const Nothing) -- full segment intersectsions are handled                        -- at insertion time   :& RNil++++type R = Rational++seg1, seg2 :: LineSegment 2 () R+seg1 = ClosedLineSegment (ext $ Point2 0 0) (ext $ Point2 0 10)+seg2 = ClosedLineSegment (ext $ Point2 0 1) (ext $ Point2 0 5)++++--------------------------------------------------------------------------------+-- *++-- | Given a predicate p on elements, and a predicate q on+-- (neighbouring) pairs of elements, filter the elements that satisfy+-- p, or together with one of their neighbours satisfy q.+overlapsOr     :: (a -> Bool)+               -> (a -> a -> Bool)+               -> [a]+               -> [a]+overlapsOr p q = map fst . filter snd . map (\((a,b),b') -> (a, b || b'))+               . overlapsWithNeighbour (q `on` fst)+               . map (\x -> (x, p x))++-- | Given a predicate, test and a list, annotate each element whether+-- it, together with one of its neighbors satisifies the predicate.+overlapsWithNeighbour   :: (a -> a -> Bool) -> [a] -> [(a,Bool)]+overlapsWithNeighbour p = go0+  where+    go0 = \case+      []     -> []+      (x:xs) -> go x False xs++    go x b = \case+      []     -> []+      (y:ys) -> let b' = p x y+                in (x,b || b') : go y b' ys++-- annotateReport   :: (a -> Bool) -> [a] -> [(a,Bool)]+-- annotateReport p = map (\x -> (x, p x))++overlapsWithNext'   :: (a -> a -> Bool) -> [a] -> [(a,Bool)]+overlapsWithNext' p = go+  where+    go = \case+      []           -> []+      [x]          -> [(x,False)]+      (x:xs@(y:_)) -> (x,p x y) : go xs++overlapsWithPrev'   :: (a -> a -> Bool) -> [a] -> [(a,Bool)]+overlapsWithPrev' p = go0+  where+    go0 = \case+      []     -> []+      (x:xs) -> (x,False) : go x xs++    go x = \case+      []     -> []+      (y:ys) -> (y,p x y) : go y ys+++overlapsWithNeighbour2 p = map (\((a,b),b') -> (a, b || b'))+                         . overlapsWithNext' (p `on` fst)+                         . overlapsWithPrev' p++shouldBe :: Eq a => a -> a -> Bool+shouldBe = (==)++propSameAsSeparate p xs = overlapsWithNeighbour p xs `shouldBe` overlapsWithNeighbour2 p xs++test' = overlapsWithNeighbour (==) testOverlapNext+testOverlapNext = [1,2,3,3,3,5,6,6,8,10,11,34,2,2,3]

@@ -6,31 +6,31 @@ -- License     :  see the LICENSE file -- Maintainer  :  David Himmelstrup ----- \( O(n \log n) \) algorithm for determining if any two line segments overlap.+-- \( O(n \log n) \) algorithm for determining if any two sets of line segments intersect. -- -- Shamos and Hoey. -- -------------------------------------------------------------------------------- module Algorithms.Geometry.LineSegmentIntersection.BooleanSweep   ( hasIntersections-  , segmentsOverlap   ) where -import           Control.Lens              hiding (contains)+import           Control.Lens hiding (contains) import           Data.Ext import           Data.Geometry.Interval-import           Data.Geometry.Line import           Data.Geometry.LineSegment import           Data.Geometry.Point-import           Data.Geometry.Triangle-import qualified Data.List                 as L++import           Data.Intersection+import qualified Data.List as L import           Data.Maybe-import           Data.Ord                  (Down (..), comparing)-import qualified Data.Set                  as SS-import qualified Data.Set.Util             as SS+import           Data.Ord (Down (..), comparing)+import qualified Data.Set as SS+import qualified Data.Set.Util as SS  -- import           Data.RealNumber.Rational--- import Debug.Trace+import Debug.Trace+import Data.Geometry.Polygon  -------------------------------------------------------------------------------- @@ -38,14 +38,14 @@ -- -- \(O(n\log n)\) hasIntersections    :: (Ord r, Num r)-                 => [LineSegment 2 p r] -> Bool+                 => [LineSegment 2 p r :+ e] -> Bool hasIntersections ss = sweep pts SS.empty   where     pts = L.sortBy ordEvents . concatMap asEventPts $ ss  -- | Computes the event points for a given line segment-asEventPts   :: Ord r => LineSegment 2 p r -> [Event p r]-asEventPts s =+asEventPts          :: Ord r => LineSegment 2 p r :+ e -> [Event p r]+asEventPts (s :+ _) =   case ordPoints (s^.start.core) (s^.end.core) of     LT -> [Insert s, Delete s]     _  -> let LineSegment a b = s@@ -57,6 +57,7 @@  -- | The actual event consists of a point and its type data Event p r = Insert (LineSegment 2 p r) | Delete (LineSegment 2 p r)+               deriving (Show)  eventPoint :: Event p r -> Point 2 r eventPoint (Insert l) = l^.start.core@@ -92,7 +93,7 @@   where     p = l^.end.core     (before,_contains,after) = splitBeforeAfter p ss-    overlaps = fromMaybe False (segmentsOverlap <$> sl <*> sr)+    overlaps = fromMaybe False (intersects <$> sl <*> sr)     sl = SS.lookupMax before     sr = SS.lookupMin after     ss' = before `SS.join` after@@ -101,14 +102,22 @@   where     p = l^.start.core     (before,contains,after) = splitBeforeAfter p ss-    endOverlap =-      (not (null contains) && isClosed startPoint)-    overlaps = or [ fromMaybe False (segmentsOverlap l <$> sl)-                  , fromMaybe False (segmentsOverlap l <$> sr) ]++    -- Check whether the endpoint is contained in one of the existing+    -- segments. The only segments that could qualify are the ones in+    -- 'contains'. Hence check only those. Note that it is not+    -- sufficient just to check whether 'contains' is empty or not,+    -- since there may be segments whose endpoint is open and coincides with p.+    endOverlap = isClosed startPoint && any (p `intersects`) contains++    overlaps =+      or [ fromMaybe False (intersects l <$> sl)+                  , fromMaybe False (intersects l <$> sr) ]     sl = SS.lookupMax before     sr = SS.lookupMin after     ss' = before `SS.join` SS.singleton l `SS.join` after + -- | split the status structure around p. -- the result is (before,contains,after) splitBeforeAfter      :: (Num r, Ord r)@@ -139,33 +148,43 @@ -------------------------------------------------------------------------------- -- * Finding New events -segmentsOverlap :: (Num r, Ord r) => LineSegment 2 p r -> LineSegment 2 p r -> Bool-segmentsOverlap a@(LineSegment aStart aEnd) b =-    (isClosed aStart && (aStart^.unEndPoint.core) `onSegment2` b) ||-    (isClosed aEnd && (aEnd^.unEndPoint.core) `onSegment2` b) ||-    (opposite (ccw' (a^.start) (b^.start) (a^.end)) (ccw' (a^.start) (b^.end) (a^.end)) &&-    not (onTriangleRelaxed (a^.end.core) t1) &&-    not (onTriangleRelaxed (a^.start.core) t2))-  where-    opposite CW CCW = True-    opposite CCW CW = True-    opposite _ _    = False-    t1 = Triangle (a^.start) (b^.start) (b^.end)-    t2 = Triangle (a^.end) (b^.start) (b^.end)+-- -- | Given two segments test if they intersect. Why don't we simply use 'intersect'+-- segmentsOverlap :: (Num r, Ord r) => LineSegment 2 p r -> LineSegment 2 p r -> Bool+-- segmentsOverlap a@(LineSegment aStart aEnd) b =+--     (isClosed aStart && (aStart^.unEndPoint.core) `intersects` b) ||+--     (isClosed aEnd && (aEnd^.unEndPoint.core) `intersects` b) ||+--     (opposite (ccw' (a^.start) (b^.start) (a^.end)) (ccw' (a^.start) (b^.end) (a^.end)) &&+--     not (onTriangleRelaxed (a^.end.core) t1) &&+--     not (onTriangleRelaxed (a^.start.core) t2))+--   where+--     opposite CW CCW = True+--     opposite CCW CW = True+--     opposite _ _    = False+--     t1 = Triangle (a^.start) (b^.start) (b^.end)+--     t2 = Triangle (a^.end) (b^.start) (b^.end) --- Copied from Data.Geometry.LineSegment.Internal. Delete when PR#62 is merged.-onSegment2                          :: (Ord r, Num r)-                                    => Point 2 r -> LineSegment 2 p r -> Bool-p `onSegment2` s@(LineSegment u v) = case ccw' (ext p) (u^.unEndPoint) (v^.unEndPoint) of-    CoLinear -> let su = p `onSide` lu-                    sv = p `onSide` lv-                in su /= sv-                && ((su == OnLine) `implies` isClosed u)-                && ((sv == OnLine) `implies` isClosed v)-    _        -> False-  where-    (Line _ w) = perpendicularTo $ supportingLine s-    lu = Line (u^.unEndPoint.core) w-    lv = Line (v^.unEndPoint.core) w -    a `implies` b = b || not a+bug' = hasIntersections $ map ext $ listEdges bug++bug :: SimplePolygon () Int+bug = fromPoints . map ext $ [+  Point2 144 592+  , Point2 336 624+  , Point2 320 544+  , Point2 240 624+  ]++s1, s2 :: LineSegment 2 () Int+s1 = read "LineSegment (Closed (Point2 240 620 :+ ())) (Open (Point2 320 544 :+ ()))"+s2 = read "LineSegment (Closed (Point2 144 592 :+ ())) (Open (Point2 336 624 :+ ()))"++tr s x = traceShow (s <> " : ", x) x++edges' :: [LineSegment 2 () Int]+edges' = [ LineSegment (Closed (Point2 240 624 :+ ())) (Open (Point2 320 544 :+ ()))+--         , LineSegment (Closed (Point2 320 544 :+ ())) (Open (Point2 336 624 :+ ()))+         , LineSegment (Closed (Point2 336 624 :+ ())) (Open (Point2 144 592 :+ ()))+         , LineSegment (Closed (Point2 144 592 :+ ())) (Open (Point2 240 624 :+ ()))+         ]++-- ah, I guess it selects the wrong predecessor/successor seg, since they overlap at the endpoint.

@@ -6,51 +6,53 @@   ) where  import           Algorithms.Geometry.LineSegmentIntersection.Types-import           Control.Lens+import           Control.Lens((^.)) import           Data.Ext-import           Data.Geometry.Interval+-- import           Data.Geometry.Interval import           Data.Geometry.LineSegment import           Data.Geometry.Point import           Data.Geometry.Properties import qualified Data.Map as M+import           Data.Util import           Data.Vinyl import           Data.Vinyl.CoRec+import qualified Data.List as List +--------------------------------------------------------------------------------  -- | Compute all intersections (naively) -- -- \(O(n^2)\)-intersections :: forall r p. (Ord r, Fractional r)-              => [LineSegment 2 p r] -> Intersections p r-intersections = foldr collect mempty . pairs+intersections :: forall r p e. (Ord r, Fractional r)+              => [LineSegment 2 p r :+ e] -> Intersections p r e+intersections = foldr collect mempty . uniquePairs  -- | Test if the two segments intersect, and if so add the segment to the map-collect          :: (Ord r, Fractional r)-                 => (LineSegment 2 p r, LineSegment 2 p r)-                 -> Intersections p r -> Intersections p r-collect (s,s') m = match (s `intersect` s') $+collect              :: (Ord r, Fractional r)+                     => Two (LineSegment 2 p r :+ e)+                     -> Intersections p r e -> Intersections p r e+collect (Two s s') m = match ((s^.core) `intersect` (s'^.core)) $      H (\NoIntersection -> m)   :& H (\p              -> handlePoint s s' p m)-  :& H (\s''            -> foldr (handlePoint s s') m [s''^.start.core, s''^.end.core])+  :& H (\s''            -> handlePoint s s' (topEndPoint s'') m)   :& RNil --- | Add s and s' to the map with key p-handlePoint        :: Ord r-                   => LineSegment 2 p r -> LineSegment 2 p r -> Point 2 r-                   -> Intersections p r -> Intersections p r-handlePoint s s' p = addTo p s . addTo p s' --- | figure out which map to add the point to-addTo                  :: Ord r => Point 2 r -> LineSegment 2 p r-                       -> Intersections p r -> Intersections p r-addTo p s-  | p `isEndPointOf` s = M.insertWith (<>) p (associated [s] [])-  | otherwise          = M.insertWith (<>) p (associated [] [s])+topEndPoint :: Ord r => LineSegment 2 p r -> Point 2 r+topEndPoint (LineSegment' (a :+ _) (b :+ _)) = List.minimumBy ordPoints [a,b] -isEndPointOf       :: Eq r => Point 2 r -> LineSegment 2 p r -> Bool-p `isEndPointOf` s = p == s^.start.core || p == s^.end.core +-- | Add s and s' to the map with key p+handlePoint        :: (Ord r, Fractional r)+                   => LineSegment 2 p r :+ e+                   -> LineSegment 2 p r :+ e+                   -> Point 2 r+                   -> Intersections p r e -> Intersections p r e+handlePoint s s' p = M.insertWith (<>) p (mkAssociated p s <> mkAssociated p s') -pairs        :: [a] -> [(a, a)]-pairs []     = []-pairs (x:xs) = map (x,) xs ++ pairs xs++type R = Rational++seg1, seg2 :: LineSegment 2 () R+seg1 = ClosedLineSegment (ext $ Point2 0 0) (ext $ Point2 0 10)+seg2 = ClosedLineSegment (ext $ Point2 0 1) (ext $ Point2 0 5)

@@ -1,3 +1,4 @@+{-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE TemplateHaskell #-} -------------------------------------------------------------------------------- -- |@@ -8,87 +9,202 @@ -------------------------------------------------------------------------------- module Algorithms.Geometry.LineSegmentIntersection.Types where +-- import           Algorithms.DivideAndConquer import           Control.DeepSeq import           Control.Lens import           Data.Ext+import           Data.Bifunctor import           Data.Geometry.Interval import           Data.Geometry.LineSegment import           Data.Geometry.Point-import           Data.Geometry.Properties-import qualified Data.List as L-import           Data.List.NonEmpty (NonEmpty(..))-import qualified Data.List.NonEmpty as NonEmpty import qualified Data.Map as Map+import qualified Data.Set as Set+import           Data.Ord (comparing, Down(..)) import           GHC.Generics+import           Data.Vinyl.CoRec+import           Data.Vinyl+import           Data.Intersection --------------------------------------------------------------------------------- -type Compare a = a -> a -> Ordering+---------------------------------------------------------------------------------- --- get the endpoints of a line segment-endPoints'   :: (HasEnd s, HasStart s) => s -> (StartCore s, EndCore s)-endPoints' s = (s^.start.core,s^.end.core) +-- FIXME: What do we do when one segmnet lies *on* the other one. For+-- the short segment it should be an "around start", but then the+-- startpoints do not match.+--+-- for the long one it's an "on" segment, but they do not intersect -type Set' l =-  Map.Map (Point (Dimension l) (NumType l), Point (Dimension l) (NumType l)) (NonEmpty l) -data Associated p r = Associated { _endPointOf        :: Set' (LineSegment 2 p r)-                                 , _interiorTo        :: Set' (LineSegment 2 p r)-                                 } deriving (Show, Generic)+-- | Assumes that two segments have the same start point+newtype AroundStart a = AroundStart a deriving (Show,Read,NFData,Functor) +instance Eq r => Eq (AroundStart (LineSegment 2 p r :+ e)) where+  -- | equality on endpoint+  (AroundStart s) == (AroundStart s') = s^.core.end.core == s'^.core.end.core -instance (Eq p, Eq r) => Eq (Associated p r) where-  (Associated es is) == (Associated es' is') = f es es' && f is is'+instance (Ord r, Num r) => Ord (AroundStart (LineSegment 2 p r :+ e)) where+  -- | ccw ordered around their suposed common startpoint+  (AroundStart s) `compare` (AroundStart s') =+    ccwCmpAround (s^.core.start.core) (s^.core.end.core)  (s'^.core.end.core)++----------------------------------------++-- | Assumes that two segments have the same end point+newtype AroundEnd a = AroundEnd a deriving (Show,Read,NFData,Functor)++instance Eq r => Eq (AroundEnd (LineSegment 2 p r :+ e)) where+  -- | equality on endpoint+  (AroundEnd s) == (AroundEnd s') = s^.core.start.core == s'^.core.start.core++instance (Ord r, Num r) => Ord (AroundEnd (LineSegment 2 p r :+ e)) where+  -- | ccw ordered around their suposed common end point+  (AroundEnd s) `compare` (AroundEnd s') =+    ccwCmpAround (s^.core.end.core) (s^.core.start.core)  (s'^.core.start.core)++--------------------------------------------------------------------------------++-- | Assumes that two segments intersect in a single point.+newtype AroundIntersection a = AroundIntersection a deriving (Show,Read,NFData,Functor)++instance Eq r => Eq (AroundIntersection (LineSegment 2 p r :+ e)) where+  -- | equality ignores the p and the e types+  (AroundIntersection (s :+ _)) == (AroundIntersection (s' :+ _))+    = first (const ()) s == first (const ()) s'++instance (Ord r, Fractional r) => Ord (AroundIntersection (LineSegment 2 p r :+ e)) where+  -- | ccw ordered around their common intersection point.+  (AroundIntersection (s :+ _)) `compare` (AroundIntersection (s' :+ _)) =+    match (s `intersect` s') $+        H (\NoIntersection     -> error "AroundIntersection: segments do not intersect!")+     :& H (\p                  -> cmpAroundP p s s')+     :& H (\_                  -> (squaredLength s) `compare` (squaredLength s'))+                                 -- if s and s' just happen to be the same length but+                                 -- intersect in different behaviour from using (==).+                                 -- but that situation doese not satisfy the precondition+                                 -- of aroundIntersection anyway.+     :& RNil     where-      f xs ys = and $ zipWith (\(p,pa) (q,qa) -> p == q && pa `sameElements` qa)-                        (Map.toAscList xs) (Map.toAscList ys)+      squaredLength (LineSegment' a b) = squaredEuclideanDist (a^.core) (b^.core) -      g = L.nub . NonEmpty.toList-      sameElements (g -> xs) (g -> ys) = L.null $ (xs L.\\ ys) ++ (ys L.\\ xs)+-- | compare around p+cmpAroundP        :: (Ord r, Num r) => Point 2 r -> LineSegment 2 p r -> LineSegment 2 p r -> Ordering+cmpAroundP p s s' = ccwCmpAround p (s^.start.core)  (s'^.start.core)  -instance (NFData p, NFData r) => NFData (Associated p r)+-- seg1 = ClosedLineSegment (ext $ Point2 0 0) (ext $ Point2 0 10)+-- seg2 = ClosedLineSegment (ext $ Point2 0 0) (ext $ Point2 0 10)  +--------------------------------------------------------------------------------  -associated       :: Ord r-                 => [LineSegment 2 p r] -> [LineSegment 2 p r] -> Associated p r-associated es is = Associated (f es) (f is)-  where-    f = foldr (\s -> Map.insertWith (<>) (endPoints' s) (s :| [])) mempty +-- | The line segments that contain a given point p may either have p+-- as the endpoint or have p in their interior.+--+-- if somehow the segment is degenerate, and p is both the start and+-- end it is reported only as the start point.+data Associated p r e =+  Associated { _startPointOf :: Set.Set (AroundEnd (LineSegment 2 p r :+ e))+             -- ^ segments for which the intersection point is the+             -- start point (i.e. s^.start.core == p)+             , _endPointOf   :: Set.Set (AroundStart (LineSegment 2 p r :+ e))+             -- ^ segments for which the intersection point is the end+             -- point (i.e. s^.end.core == p)+             , _interiorTo   :: Set.Set (AroundIntersection (LineSegment 2 p r :+ e))+             } deriving stock (Show, Read, Generic, Eq) -endPointOf :: Associated p r -> [LineSegment 2 p r]-endPointOf = concatMap NonEmpty.toList . Map.elems . _endPointOf+makeLenses ''Associated -interiorTo :: Associated p r -> [LineSegment 2 p r]-interiorTo = concatMap NonEmpty.toList . Map.elems . _interiorTo+instance Functor (Associated p r) where+  fmap f (Associated ss es is) = Associated (Set.mapMonotonic (g f) ss)+                                            (Set.mapMonotonic (g f) es)+                                            (Set.mapMonotonic (g f) is)+    where+      g   :: forall f c e b. Functor f => (e -> b) -> f (c :+ e) -> f (c :+ b)+      g f' = fmap (&extra %~ f')  -instance Ord r => Semigroup (Associated p r) where-  (Associated es is) <> (Associated es' is') = Associated (es <> es') (is <> is')+-- | Reports whether this associated has any interior intersections+--+-- \(O(1)\)+isInteriorIntersection :: Associated p r e -> Bool+isInteriorIntersection = not . null . _interiorTo -instance Ord r => Monoid (Associated p r) where-  mempty = Associated mempty mempty-  mappend = (<>) -type Intersections p r = Map.Map (Point 2 r) (Associated p r)+-- | test if the given segment has p as its endpoint, an construct the+-- appropriate associated representing that.+--+-- pre: p intersects the segment+mkAssociated                :: (Ord r, Fractional r)+                            => Point 2 r -> LineSegment 2 p r :+ e-> Associated p r e+mkAssociated p s@(LineSegment a b :+ _)+  | p == a^.unEndPoint.core = mempty&startPointOf .~  Set.singleton (AroundEnd s)+  | p == b^.unEndPoint.core = mempty&endPointOf   .~  Set.singleton (AroundStart s)+  | otherwise               = mempty&interiorTo   .~  Set.singleton (AroundIntersection s) -data IntersectionPoint p r =++-- | test if the given segment has p as its endpoint, an construct the+-- appropriate associated representing that.+--+-- If p is not one of the endpoints we concstruct an empty Associated!+--+mkAssociated'     :: (Ord r, Fractional r)+                  => Point 2 r -> LineSegment 2 p r :+ e -> Associated p r e+mkAssociated' p s = (mkAssociated p s)&interiorTo .~ mempty++instance (Ord r, Fractional r) => Semigroup (Associated p r e) where+  (Associated ss es is) <> (Associated ss' es' is') =+    Associated (ss <> ss') (es <> es') (is <> is')++instance (Ord r, Fractional r) => Monoid (Associated p r e) where+  mempty = Associated mempty mempty mempty++instance (NFData p, NFData r, NFData e) => NFData (Associated p r e)++-- | For each intersection point the segments intersecting there.+type Intersections p r e = Map.Map (Point 2 r) (Associated p r e)++-- | An intersection point together with all segments intersecting at+-- this point.+data IntersectionPoint p r e =   IntersectionPoint { _intersectionPoint :: !(Point 2 r)-                    , _associatedSegs    :: !(Associated p r)-                    } deriving (Show,Eq)+                    , _associatedSegs    :: !(Associated p r e)+                    } deriving (Show,Read,Eq,Generic,Functor) makeLenses ''IntersectionPoint +instance (NFData p, NFData r, NFData e) => NFData (IntersectionPoint p r e) --- | reports true if there is at least one segment for which this intersection--- point is interior.------ \(O(1)\)-isEndPointIntersection :: Associated p r -> Bool-isEndPointIntersection = Map.null . _interiorTo +-- sameOrder           :: (Ord r, Num r, Eq p) => Point 2 r+--                     -> [LineSegment 2 p r] -> [LineSegment 2 p r] -> Bool+-- sameOrder c ss ss' = f ss == f ss'+--   where+--     f = map (^.extra) . sortAround' (ext c) . map (\s -> s^.end.core :+ s) --- newtype E a b = E (a -> b)++++-- | Given a point p, and a bunch of segments that suposedly intersect+-- at p, correctly categorize them.+mkIntersectionPoint         :: (Ord r, Fractional r)+                            => Point 2 r+                            -> [LineSegment 2 p r :+ e] -- ^ uncategorized+                            -> [LineSegment 2 p r :+ e] -- ^ segments we know contain p,+                            -> IntersectionPoint p r e+mkIntersectionPoint p as cs = IntersectionPoint p $ foldMap (mkAssociated p) $ as <> cs++  -- IntersectionPoint p+  --                           $ Associated mempty mempty (Set.fromAscList cs')+  --                           <> foldMap (mkAssociated p) as+  -- where+  --   cs' = map AroundIntersection . List.sortBy (cmpAroundP p) $ cs+  -- -- TODO: In the bentley ottman algo we already know the sorted order of the segments+  -- -- so we can likely save the additional sort++++-- | An ordering that is decreasing on y, increasing on x+ordPoints     :: Ord r => Point 2 r -> Point 2 r -> Ordering+ordPoints a b = let f p = (Down $ p^.yCoord, p^.xCoord) in comparing f a b

@@ -36,7 +36,7 @@ -- the upper convex hull. It uses the given algorithm to do so -- -- running time: O(time required by the given upper hull algorithm)-lowerEnvelopeWith        :: (Fractional r, Eq r)+lowerEnvelopeWith        :: (Fractional r, Ord r)                          => UpperHullAlgorithm (Line 2 r :+ a) r                          -> NonEmpty (Line 2 r :+ a) -> Envelope a r lowerEnvelopeWith chAlgo = fromPts . chAlgo . toPts

@@ -9,7 +9,7 @@ -- -- Ear clipping triangulation algorithms. The baseline algorithm runs in \( O(n^2) \) -- but has a low constant factor overhead. The z-order hashed variant runs in--- \( O(n \log n) \).+-- \( O(n \log n) \) time. -- -- References: --

@@ -7,38 +7,37 @@ -- License     :  see the LICENSE file -- Maintainer  :  Frank Staals ---------------------------------------------------------------------------------module Algorithms.Geometry.PolygonTriangulation.MakeMonotone( makeMonotone-                                                            , computeDiagonals+module Algorithms.Geometry.PolygonTriangulation.MakeMonotone+  ( makeMonotone+  , computeDiagonals  -                                                            , VertexType(..)-                                                            , classifyVertices-                                                            ) where--import Algorithms.Geometry.LineSegmentIntersection.BentleyOttmann (ordAt, xCoordAt)-import Algorithms.Geometry.PolygonTriangulation.Types+  , VertexType(..)+  , classifyVertices+  ) where +import           Algorithms.Geometry.PolygonTriangulation.Types import           Control.Lens import           Control.Monad.Reader import           Control.Monad.State.Strict-import           Control.Monad.Writer                  (WriterT, execWriterT, tell)+import           Control.Monad.Writer (WriterT, execWriterT, tell) import           Data.Bifunctor-import qualified Data.DList                            as DList+import qualified Data.DList as DList import           Data.Ext-import qualified Data.Foldable                         as F+import qualified Data.Foldable as F import           Data.Geometry.LineSegment import           Data.Geometry.PlanarSubdivision.Basic import           Data.Geometry.Point import           Data.Geometry.Polygon-import qualified Data.IntMap                           as IntMap-import qualified Data.List.NonEmpty                    as NonEmpty-import           Data.Ord                              (Down (..), comparing)-import qualified Data.Set                              as SS-import qualified Data.Set.Util                         as SS+import qualified Data.IntMap as IntMap+import qualified Data.List.NonEmpty as NonEmpty+import           Data.Ord (Down (..), comparing)+import qualified Data.Set as SS+import qualified Data.Set.Util as SS import           Data.Util-import qualified Data.Vector                           as V-import qualified Data.Vector.Circular                  as CV-import qualified Data.Vector.Mutable                   as MV+import qualified Data.Vector as V+import qualified Data.Vector.Circular as CV+import qualified Data.Vector.Mutable as MV   -- import Debug.Trace@@ -162,11 +161,11 @@ -- 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-                  -> PlanarSubdivision s p PolygonEdgeType PolygonFaceData r-makeMonotone px pg = let (e:es) = listEdges pg-                     in constructSubdivision px e es (computeDiagonals pg)+makeMonotone    :: forall s t p r. (Fractional r, Ord r)+                => Polygon t p r+                -> PlanarSubdivision s p PolygonEdgeType PolygonFaceData r+makeMonotone pg = let (e:es) = listEdges pg+                  in constructSubdivision e es (computeDiagonals pg)  type Sweep p r = WriterT (DList.DList (LineSegment 2 Int r))                    (StateT (StatusStruct r)@@ -199,11 +198,11 @@  insertAt   :: (Ord r, Fractional r) => Point 2 r -> LineSegment 2 q r            -> SS.Set (LineSegment 2 q r) -> SS.Set (LineSegment 2 q r)-insertAt v = SS.insertBy (ordAt $ v^.yCoord)+insertAt v = SS.insertBy (ordAtY $ v^.yCoord)  deleteAt   :: (Fractional r, Ord r) => Point 2 r -> LineSegment 2 p r            -> SS.Set (LineSegment 2 p r) -> SS.Set (LineSegment 2 p r)-deleteAt v = SS.deleteAllBy (ordAt $ v^.yCoord)+deleteAt v = SS.deleteAllBy (ordAtY $ v^.yCoord)   handleStart              :: (Fractional r, Ord r)

@@ -24,10 +24,9 @@ -- | Triangulates a polygon of \(n\) vertices -- -- running time: \(O(n \log n)\)-triangulate        :: (Ord r, Fractional r)-                   => proxy s -> Polygon t p r-                   -> PlanarSubdivision s p PolygonEdgeType PolygonFaceData r-triangulate px pg' = constructSubdivision px e es diags+triangulate     :: forall s t p r. (Ord r, Fractional r)+                => Polygon t p r -> PlanarSubdivision s p PolygonEdgeType PolygonFaceData r+triangulate pg' = constructSubdivision e es diags   where     (pg, diags)   = computeDiagonals' pg'     (e:es)        = listEdges pg@@ -36,10 +35,9 @@ -- | Triangulates a polygon of \(n\) vertices -- -- running time: \(O(n \log n)\)-triangulate'        :: (Ord r, Fractional r)-                    => proxy s -> Polygon t p r-                    -> PlaneGraph s p PolygonEdgeType PolygonFaceData r-triangulate' px pg' = constructGraph px e es diags+triangulate'     :: forall s t p r. (Ord r, Fractional r)+                 => Polygon t p r -> PlaneGraph s p PolygonEdgeType PolygonFaceData r+triangulate' pg' = constructGraph e es diags   where     (pg, diags)   = computeDiagonals' pg'     (e:es)        = listEdges pg@@ -62,7 +60,7 @@ computeDiagonals' pg' = (pg, monotoneDiags <> extraDiags)   where     pg            = toCounterClockWiseOrder pg'-    monotoneP     = MM.makeMonotone (Identity pg) pg -- use some arbitrary proxy type+    monotoneP     = MM.makeMonotone @() pg -- use some arbitrary proxy type     -- outerFaceId'  = outerFaceId monotoneP      monotoneDiags = map (^._2.core) . filter (\e' -> e'^._2.extra == Diagonal)

@@ -47,10 +47,9 @@ -- | Triangulates a polygon of \(n\) vertices -- -- running time: \(O(n \log n)\)-triangulate        :: (Ord r, Fractional r)-                   => proxy s -> MonotonePolygon p r-                   -> PlanarSubdivision s p PolygonEdgeType PolygonFaceData r-triangulate px pg' = constructSubdivision px e es (computeDiagonals pg)+triangulate     :: forall s p r. (Ord r, Fractional r)+                => MonotonePolygon p r -> PlanarSubdivision s p PolygonEdgeType PolygonFaceData r+triangulate pg' = constructSubdivision e es (computeDiagonals pg)   where     pg     = toCounterClockWiseOrder pg'     (e:es) = listEdges pg@@ -59,10 +58,9 @@ -- | Triangulates a polygon of \(n\) vertices -- -- running time: \(O(n \log n)\)-triangulate'        :: (Ord r, Fractional r)-                    => proxy s -> MonotonePolygon p r-                    -> PlaneGraph s p PolygonEdgeType PolygonFaceData r-triangulate' px pg' = constructGraph px e es (computeDiagonals pg)+triangulate'     :: forall s p r. (Ord r, Fractional r)+                 => MonotonePolygon p r-> PlaneGraph s p PolygonEdgeType PolygonFaceData r+triangulate' pg' = constructGraph e es (computeDiagonals pg)   where     pg     = toCounterClockWiseOrder pg'     (e:es) = listEdges pg

@@ -31,16 +31,15 @@ -- -- -- running time: \(O(n\log n)\)-constructSubdivision                  :: forall proxy r s p. (Fractional r, Ord r)-                                      => proxy s-                                      -> LineSegment 2 p r -- ^ A counter-clockwise-                                                         -- edge along the outer-                                                         -- boundary-                                      -> [LineSegment 2 p r] -- ^ remaining original edges-                                      -> [LineSegment 2 p r] -- ^ diagonals-                                      -> PlanarSubdivision s-                                            p PolygonEdgeType PolygonFaceData r-constructSubdivision px e origs diags = fromPlaneGraph $ constructGraph px e origs diags+constructSubdivision               :: forall s r p. (Fractional r, Ord r)+                                   => LineSegment 2 p r -- ^ A counter-clockwise+                                                      -- edge along the outer+                                                      -- boundary+                                   -> [LineSegment 2 p r] -- ^ remaining original edges+                                   -> [LineSegment 2 p r] -- ^ diagonals+                                   -> PlanarSubdivision s+                                         p PolygonEdgeType PolygonFaceData r+constructSubdivision e origs diags = fromPlaneGraph $ constructGraph e origs diags  -- constructSubdivision px e origs diags = --     subdiv & rawVertexData.traverse.dataVal  %~ NonEmpty.head@@ -80,22 +79,21 @@ -- -- -- running time: \(O(n\log n)\)-constructGraph                  :: forall proxy r s p. (Fractional r, Ord r)-                                      => proxy s-                                      -> LineSegment 2 p r -- ^ A counter-clockwise+constructGraph                  :: forall s r p. (Fractional r, Ord r)+                                      => LineSegment 2 p r -- ^ A counter-clockwise                                                          -- edge along the outer                                                          -- boundary                                       -> [LineSegment 2 p r] -- ^ remaining original edges                                       -> [LineSegment 2 p r] -- ^ diagonals                                       -> PG.PlaneGraph s                                             p PolygonEdgeType PolygonFaceData r-constructGraph px e origs diags =+constructGraph e origs diags =     subdiv & PG.vertexData.traverse  %~ NonEmpty.head            & PG.faceData             .~ faceData'            & PG.rawDartData.traverse %~ snd   where     subdiv :: PG.PlaneGraph s (NonEmpty p) (Bool,PolygonEdgeType) () r-    subdiv = PG.fromConnectedSegments px $ e' : origs' <> diags'+    subdiv = PG.fromConnectedSegments $ e' : origs' <> diags'      diags' = (:+ (True, Diagonal)) <$> diags     origs' = (:+ (False,Original)) <$> origs

@@ -0,0 +1,88 @@+module Algorithms.Geometry.RayShooting.Naive+  ( firstHit+  , firstHit'+++  , firstHitSegments+  , intersectionDistance+  , labelWithDistances+  ) where++import           Control.Lens+import           Data.Bifunctor+import           Data.Ext+import           Data.Geometry.HalfLine+import           Data.Geometry.LineSegment+import           Data.Geometry.Point+import           Data.Geometry.Polygon+import           Data.Intersection+import qualified Data.List as List+import           Data.Maybe+import           Data.Ord (comparing)+import           Data.Vinyl.CoRec+import           Data.Vinyl++--------------------------------------------------------------------------------++-- |+--+-- pre: halfline should start in the interior+firstHit     :: (Fractional r, Ord r)+             => HalfLine 2 r+             -> Polygon t p r+             -> LineSegment 2 p r+firstHit ray = fromMaybe err . firstHit' ray+  where+    err = error "Algorithms.Geometry.RayShooting.Naive: no intersections; ray must have started outside the polygon"++-- | Compute the first edge hit by the ray, if it exists+firstHit'        :: (Fractional r, Ord r)+                 => HalfLine 2 r+                 -> Polygon t p r+                 -> Maybe (LineSegment 2 p r)+firstHit' ray pg = fmap (^.core) . firstHitSegments ray . map ext $ listEdges pg+++-- | Compute the first segment hit by the ray, if it exists+firstHitSegments     :: (Ord r, Fractional r)+                     => HalfLine 2 r+                     -> [LineSegment 2 p r :+ e]+                     -> Maybe (LineSegment 2 p r :+ e)+firstHitSegments ray = fmap (^.extra) . minimumOn (^.core)+                     . mapMaybe (\(s :+ (md, e)) -> (:+ (s :+ e)) <$> md)+                     . labelWithDistances (ray^.startPoint) ray++minimumOn   :: Ord b => (a -> b) -> [a] -> Maybe a+minimumOn f = \case+  [] -> Nothing+  xs -> Just . List.minimumBy (comparing f) $ xs+++-- | Given q, a ray, and a segment s, computes if the+-- segment intersects the initial, rightward ray starting in q, and if+-- so returns the (squared) distance from q to that point together+-- with the segment.+intersectionDistance         :: forall r p. (Ord r, Fractional r)+                             => Point 2 r -> HalfLine 2 r -> LineSegment 2 p r+                             -> Maybe r+intersectionDistance q ray s = match (seg `intersect` ray) $+       H (\NoIntersection                   -> Nothing)+    :& H (\p                                -> Just $ d p)+    :& H (\(LineSegment' (a :+ _) (b :+ _)) -> Just $ d a `min` d b)+    :& RNil+    -- TODO: there is some slight subtility if the segment is open.+  where+    d = squaredEuclideanDist q+    seg = first (const ()) s+++-- | Labels the segments with the distance from q to their+-- intersection point with the ray.+labelWithDistances       :: (Ord r, Fractional r)+                         => Point 2 r -> HalfLine 2 r -> [LineSegment 2 p r :+ b]+                         -> [LineSegment 2 p r :+ (Maybe r, b)]+labelWithDistances q ray = map (\(s :+ e) -> s :+ (intersectionDistance q ray s, e))++++--------------------------------------------------------------------------------

@@ -38,7 +38,6 @@                                                   VertexId, VertexId', dual, graph, incidentEdges,                                                   leftFace, vertices) import qualified Data.PlaneGraph                 as PlaneGraph-import           Data.Proxy import           Data.Tree                       (Tree (Node)) import qualified Data.Vector                     as V import qualified Data.Vector.Circular            as CV@@ -69,10 +68,11 @@ --        we're running 'unsafeFromPoints . toPoints' to reset the polygon. --        Super silly. Please fix. -- | \( O(n \log n) \)-triangulate :: (Ord r, Fractional r) => SimplePolygon p r -> PlaneGraph s Int PolygonEdgeType PolygonFaceData r+triangulate   :: forall s p r. (Ord r, Fractional r)+              => SimplePolygon p r -> PlaneGraph s Int PolygonEdgeType PolygonFaceData r triangulate p =   let poly' = snd $ bimapAccumL (\a _ -> (a+1,a)) (,) 0 $ unsafeFromPoints $ toPoints p-  in triangulate' Proxy poly'+  in triangulate' @s poly'  -- | \( O(n) \) Single-Source shortest path. sssp :: (Ord r, Fractional r)

@@ -24,6 +24,7 @@   , compareAroundEndPoint   ) where +import           Algorithms.Geometry.RayShooting.Naive import           Control.Lens import           Control.Monad ((<=<)) import           Data.Bifunctor (first)@@ -38,9 +39,9 @@ import           Data.Geometry.Vector import           Data.Intersection import qualified Data.List as List-import qualified Data.List.Util as List import           Data.List.NonEmpty (NonEmpty(..)) import qualified Data.List.NonEmpty as NonEmpty+import qualified Data.List.Util as List import           Data.Maybe (mapMaybe, isJust) import           Data.Ord (comparing) import           Data.RealNumber.Rational@@ -416,13 +417,6 @@       Just z  -> Just $ squaredEuclideanDist q z   where     seg = first (const ()) s---- | Labels the segments with the distance from q to their--- intersection point with the ray.-labelWithDistances       :: (Ord r, Fractional r)-                         => Point 2 r -> HalfLine 2 r -> [LineSegment 2 p r :+ b]-                         -> [LineSegment 2 p r :+ (Maybe r, b)]-labelWithDistances q ray = map (\(s :+ e) -> s :+ (initialIntersection q ray s, e))  -------------------------------------------------------------------------------- -- * Comparators for the rotating ray

@@ -62,7 +62,7 @@ fairSplitTree pts = foldUp node' Leaf $ fairSplitTree' n pts'   where     pts' = imap sortOn . pure . g $ pts-    n    = length $ pts'^.GV.element (C :: C 0)+    n    = length $ pts'^.GV.element @0      sortOn' i = NonEmpty.sortWith (^.core.unsafeCoord i)     sortOn  i = LSeq.fromNonEmpty . sortOn' (i + 1)@@ -130,7 +130,7 @@                      => Int -> GV.Vector d (PointSeq d (Idx :+ p) r)                      -> BinLeafTree Int (Point d r :+ p) fairSplitTree' n pts-    | n <= 1    = let p = LSeq.head $ pts^.GV.element (C :: C 0) in Leaf (dropIdx p)+    | n <= 1    = let p = LSeq.head $ pts^.GV.element @0 in Leaf (dropIdx p)     | otherwise = foldr node' (V.last path) $ V.zip nodeLevels (V.init path)   where     -- note that points may also be assigned level 'Nothing'.@@ -156,7 +156,7 @@     node' (lvl,lc) rc = case lvl^?widestDim._Just of                           Nothing -> error "Unknown widest dimension"                           Just j  -> Node lc j rc-    recurse pts' = fairSplitTree' (length $ pts'^.GV.element (C :: C 0))+    recurse pts' = fairSplitTree' (length $ pts'^.GV.element @0)                                   (reIndexPoints pts')  -- | Assign the points to their the correct class. The 'Nothing' class is@@ -209,7 +209,7 @@                    -> GV.Vector d (PointSeq d (Idx :+ p) r) reIndexPoints ptsV = fmap reIndex ptsV   where-    pts = ptsV^.GV.element (C :: C 0)+    pts = ptsV^.GV.element @0      reIndex = fmap (\p -> p&extra.core %~ fromJust . flip IntMap.lookup mapping')     mapping' = IntMap.fromAscList $ zip (map (^.extra.core) . F.toList $ pts) [0..]

@@ -53,39 +53,36 @@ -- | Builds an arrangement of \(n\) lines -- -- running time: \(O(n^2\log n\)-constructArrangement       :: (Ord r, Fractional r)-                           => proxy s-                           -> [Line 2 r :+ l]-                           -> Arrangement s l () (Maybe l) () r-constructArrangement px ls = let b  = makeBoundingBox ls-                             in constructArrangementInBox' px b ls+constructArrangement    :: forall s l r. (Ord r, Fractional r)+                        => [Line 2 r :+ l]+                        -> Arrangement s l () (Maybe l) () r+constructArrangement ls = let b  = makeBoundingBox ls+                          in constructArrangementInBox' b ls  -- | Constructs the arrangemnet inside the box.  note that the resulting box -- may be larger than the given box to make sure that all vertices of the -- arrangement actually fit. -- -- running time: \(O(n^2\log n\)-constructArrangementInBox            :: (Ord r, Fractional r)-                                     => proxy s-                                     -> Rectangle () r-                                     -> [Line 2 r :+ l]-                                     -> Arrangement s l () (Maybe l) () r-constructArrangementInBox px rect ls = let b  = makeBoundingBox ls-                                       in constructArrangementInBox' px (b <> rect) ls+constructArrangementInBox         :: forall s l r. (Ord r, Fractional r)+                                  => Rectangle () r+                                  -> [Line 2 r :+ l]+                                  -> Arrangement s l () (Maybe l) () r+constructArrangementInBox rect ls = let b  = makeBoundingBox ls+                                    in constructArrangementInBox' (b <> rect) ls   -- | Constructs the arrangemnet inside the box. (for parts to be useful, it is -- assumed this boxfits at least the boundingbox of the intersections in the -- Arrangement)-constructArrangementInBox'            :: (Ord r, Fractional r)-                                      => proxy s-                                      -> Rectangle () r-                                      -> [Line 2 r :+ l]-                                      -> Arrangement s l () (Maybe l) () r-constructArrangementInBox' px rect ls =+constructArrangementInBox'         :: forall s l r. (Ord r, Fractional r)+                                   => Rectangle () r+                                   -> [Line 2 r :+ l]+                                   -> Arrangement s l () (Maybe l) () r+constructArrangementInBox' rect ls =     Arrangement (V.fromList ls) subdiv rect (link parts' subdiv)   where-    subdiv = fromConnectedSegments px segs+    subdiv = fromConnectedSegments segs                 & rawVertexData.traverse.dataVal .~ ()     (segs,parts') = computeSegsAndParts rect ls 

@@ -156,7 +156,7 @@ -- -- >>> 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)+disk       :: (Ord r, Fractional r)            => Point 2 r -> Point 2 r -> Point 2 r -> Maybe (Disk () r) disk p q r = match (f p `intersect` f q) $        H (\NoIntersection -> Nothing)

@@ -362,7 +362,7 @@ flat r b = let p = fst $ endPoints b                q = snd $ endPoints b                s = ClosedLineSegment (p :+ ()) (q :+ ())-               e t = sqDistanceToSeg (evaluate b t) s < r ^ 2+               e t = squaredEuclideanDistTo (evaluate b t) s < r ^ 2            in qdA p q < r ^ 2 || all e [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]  -- seems this is now covered by approximate@@ -394,7 +394,7 @@ --   twice the given tolerance. May return false negatives but not false positives. colinear :: (Ord r, Fractional r) => r -> BezierSpline 3 2 r -> Bool colinear eps (Bezier3 !a !b !c !d) = sqBound < eps*eps-  where ld = flip sqDistanceTo (lineThrough a d)+  where ld = flip squaredEuclideanDistTo (lineThrough a d)         sameSide = ccw a d b == ccw a d c         maxDist = max (ld b) (ld c)         sqBound@@ -426,13 +426,13 @@       recurse2 = 0.5 + 0.5 * parameterInterior treshold b2 p       chb1     = _simplePolygon $ convexHullB b1       chb2     = _simplePolygon $ convexHullB b2-      in1      = sqDistanceToPolygon p chb1 < treshold^2-      in2      = sqDistanceToPolygon p chb2 < treshold^2+      in1      = squaredEuclideanDistTo p chb1 < treshold^2+      in2      = squaredEuclideanDistTo p chb2 < treshold^2       result |     in1 &&     in2 = betterFit b p recurse1 recurse2              |     in2 && not in2 = recurse1              | not in2 &&     in2 = recurse2-             | sqDistanceToPolygon p chb1 < sqDistanceToPolygon p chb2 = recurse1-             | otherwise                                               = recurse2+             | squaredEuclideanDistTo p chb1 < squaredEuclideanDistTo p chb2 = recurse1+             | otherwise                                                     = recurse2   in result  -- | Given a point on (or close to) the extension of a Bezier curve, return the corresponding@@ -458,16 +458,6 @@   let q = evaluate b t       r = evaluate b u   in if qdA q p < qdA r p then t else u--sqDistanceToPolygon :: (Ord r, Fractional r) => Point 2 r -> SimplePolygon p r -> r-sqDistanceToPolygon point poly | insidePolygon point poly = 0-                               | otherwise = minimum $ map (sqDistanceToSeg point) $ listEdges poly--------------------------------------------------------------------------------------  -------------------------------------------------------------------------------- 

@@ -11,17 +11,29 @@ -- Orthogonal \(d\)-dimensiontal boxes (e.g. rectangles) -- ---------------------------------------------------------------------------------module Data.Geometry.Box( module Data.Geometry.Box.Internal-                        , module Data.Geometry.Box.Corners-                        , module Data.Geometry.Box.Sides-                        ) where+module Data.Geometry.Box+  ( module Data.Geometry.Box.Internal+  , module Data.Geometry.Box.Corners+  , module Data.Geometry.Box.Sides+  , inBox'+  ) where  import Control.DeepSeq import Data.Geometry.Box.Corners import Data.Geometry.Box.Internal import Data.Geometry.Box.Sides import Data.Geometry.Vector+import Data.Geometry.Point+import Data.Geometry.Boundary  --------------------------------------------------------------------------------  deriving instance (NFData p, NFData r, Arity d) => NFData (Box d p r)+++-- | Compute whether the point lies inside, on the boundary of, or+-- outside the box.+inBox' :: (Arity d, Ord r) => Point d r -> Box d p r -> PointLocationResult+q `inBox'` b | q `insideBox` b = Inside+             | q `inBox`     b = OnBoundary+             | otherwise       = Outside

@@ -1,6 +1,7 @@ {-# LANGUAGE TemplateHaskell  #-} {-# LANGUAGE UndecidableInstances  #-} {-# LANGUAGE InstanceSigs  #-}+{-# LANGUAGE AllowAmbiguousTypes #-} -------------------------------------------------------------------------------- -- | -- Module      :  Data.Geometry.Box.Internal@@ -205,7 +206,6 @@   where     toOpenRange (R.Range' l r) = R.OpenRange l r - -- | Get a vector with the extent of the box in each dimension. Note that the -- resulting vector is 0 indexed whereas one would normally count dimensions -- starting at zero.@@ -226,13 +226,13 @@  -- | 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 @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 @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-widthIn _ = view (V.element (C :: C (i - 1))) . size+widthIn :: forall i p d r. (Arity d, Arity (i - 1), Num r, ((i-1)+1) <= d)+        => Box d p r -> r+widthIn = view (V.element @(i-1)) . size   -- | Same as 'widthIn' but with a runtime int instead of a static dimension.@@ -258,7 +258,7 @@ -- >>> width (boundingBoxList' [origin] :: Rectangle () Int) -- 0 width :: Num r => Rectangle p r -> r-width = widthIn (C :: C 1)+width = widthIn @1  -- | -- >>> height (boundingBoxList' [origin, Point2 1 2] :: Rectangle () Int)@@ -266,7 +266,7 @@ -- >>> height (boundingBoxList' [origin] :: Rectangle () Int) -- 0 height :: Num r => Rectangle p r -> r-height = widthIn (C :: C 2)+height = widthIn @2   --------------------------------------------------------------------------------@@ -298,3 +298,26 @@  instance IsBoxable c => IsBoxable (c :+ e) where   boundingBox = boundingBox . view core++--------------------------------------------------------------------------------+-- * Distances++instance (Num r, Ord r) => HasSquaredEuclideanDistance (Box 2 p r) where+  pointClosestToWithDistance q bx =+      case ((q^.xCoord) `R.inRange` hor, (q^.yCoord) `R.inRange` ver) of+                      (False,False) -> if q^.yCoord < b+                                       then closest (Point2 l b) (Point2 r b)+                                       else closest (Point2 l t) (Point2 r t)+                      (True, False) -> if q^.yCoord < b+                                       then (q&yCoord .~ b, sq $ q^.yCoord - b)+                                       else (q&yCoord .~ t, sq $ q^.yCoord - t)+                      (False, True) -> if q^.xCoord < l+                                       then (q&yCoord .~ l, sq $ q^.xCoord - l)+                                       else (q&yCoord .~ r, sq $ q^.xCoord - r)+                      (True, True)  -> (q, 0) -- point lies inside the box+    where+      Vector2 hor@(R.Range' l r) ver@(R.Range' b t) = extent bx+      sq x = x*x+      closest p1 p2 = let d1 = squaredEuclideanDist q p1+                          d2 = squaredEuclideanDist q p2+                      in if d1 < d2 then (p1, d1) else (p2, d2)

@@ -16,7 +16,7 @@ import Data.Geometry.Directions import Data.Geometry.Box.Internal import Data.Geometry.Box.Corners-import Data.Geometry.LineSegment+import Data.Geometry.LineSegment.Internal import Data.Functor.Apply import Data.Semigroup.Foldable.Class import Data.Semigroup.Traversable.Class

@@ -19,9 +19,9 @@  -- | 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   :: (Fractional r, Ord r) => Line 2 r -> Maybe (Point 2 r) 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+dualPoint' :: (Fractional r, Ord r) => Line 2 r -> Point 2 r dualPoint' = fromJust . dualPoint

@@ -18,14 +18,16 @@ import Data.Geometry.Transformation import Data.Geometry.Vector import GHC.Generics (Generic)+import Data.Kind import GHC.TypeLits  --------------------------------------------------------------------------------  -- | Hyperplanes embedded in a \(d\) dimensional space.-data HyperPlane (d :: Nat) (r :: *) = HyperPlane { _inPlane   :: !(Point d r)-                                                 , _normalVec :: !(Vector d r)-                                                 } deriving Generic+data HyperPlane (d :: Nat) (r :: Type) =+  HyperPlane { _inPlane   :: !(Point d r)+             , _normalVec :: !(Vector d r)+             } deriving Generic makeLenses ''HyperPlane  type instance Dimension (HyperPlane d r) = d

@@ -14,7 +14,7 @@                                -- * querying the start and end of intervals                              , HasStart(..), HasEnd(..)                              -- * Working with intervals-                             , inInterval+                             , intersectsInterval, inInterval                              , shiftLeft'                               , asProperInterval, flipInterval@@ -23,18 +23,19 @@                              ) where  import           Control.DeepSeq-import           Control.Lens             (Iso', Lens', iso, (%~), (&), (^.))+import           Control.Lens (Iso', Lens', iso, (%~), (&), (^.)) import           Data.Bifunctor import           Data.Bitraversable import           Data.Ext-import qualified Data.Foldable            as F+import qualified Data.Foldable as F+import           Data.Geometry.Boundary import           Data.Geometry.Properties import           Data.Range-import           Data.Semigroup           (Arg (..))-import qualified Data.Traversable         as T+import           Data.Semigroup (Arg (..))+import qualified Data.Traversable as T import           Data.Vinyl import           Data.Vinyl.CoRec-import           GHC.Generics             (Generic)+import           GHC.Generics (Generic) import           Test.QuickCheck  --------------------------------------------------------------------------------@@ -80,14 +81,39 @@   bimap f g (GInterval r) = GInterval $ fmap (bimap g f) r  +-- type instance IntersectionOf r (Interval b r) = [NoIntersection, r]+-- -- somehow: GHC does not understand the r here cannot be 'Interval a r' itself :( +-- instance Ord r => r `HasIntersectionWith` Interval b r where+--   x `intersects` r = x `inRange` fmap (^.core) (r^._Range )+++-- instance Ord r => r `IsIntersectableWith` Interval b r where+--   x `intersect` r | x `intersects` r = coRec x+--                   | otherwise        = coRec NoIntersection+ -- | Test if a value lies in an interval. Note that the difference between --  inInterval and inRange is that the extra value is *not* used in the --  comparison with inInterval, whereas it is in inRange.-inInterval       :: Ord r => r -> Interval a r -> Bool-x `inInterval` r = x `inRange` fmap (^.core) (r^._Range )+intersectsInterval       :: Ord r => r -> Interval a r -> Bool+x `intersectsInterval` r = x `inRange` fmap (^.core) (r^._Range )  +-- | Compute where the given query value is with respect to the interval.+--+-- Note that even if the boundary of the interval is open we may+-- return "OnBoundary".+inInterval :: Ord r => r -> Interval a r -> PointLocationResult+x `inInterval` (Interval l r) =+  case x `compare` (l^.unEndPoint.core) of+    LT -> Outside+    EQ -> OnBoundary+    GT -> case x `compare` (r^.unEndPoint.core) of+            LT -> Inside+            EQ -> OnBoundary+            GT -> Outside++ pattern OpenInterval       :: (r :+ a) -> (r :+ a) -> Interval a r pattern OpenInterval   l u = GInterval (OpenRange   l u) @@ -127,10 +153,11 @@ type instance NumType   (Interval a r) = r  -type instance IntersectionOf (Interval a r) (Interval a r) = [NoIntersection, Interval a r]+type instance IntersectionOf (Interval a r) (Interval b r)+  = [NoIntersection, Interval (Either a b) r] -instance Ord r => Interval a r `HasIntersectionWith` Interval a r-instance Ord r => Interval a r `IsIntersectableWith` Interval a r where+instance Ord r => Interval a r `HasIntersectionWith` Interval b r+instance Ord r => Interval a r `IsIntersectableWith` Interval b r where    nonEmptyIntersection = defaultNonEmptyIntersection @@ -140,9 +167,10 @@                                                          (u&unEndPoint %~ g) )       :& RNil     where-      f x = Arg (x^.core) x-      r' = fmap f r-      s' = fmap f s+      r' :: Range (Arg r (r :+ Either a b))+      r' = fmap (\(x :+ a) -> Arg x (x :+ Left a))  r+      s' :: Range (Arg r (r :+ Either a b))+      s' = fmap (\(x :+ b) -> Arg x (x :+ Right b)) s        g (Arg _ x) = x 

@@ -189,6 +189,6 @@ asSingleton   :: (1 <= d, Arity d)               => PointSet (LSeq 1) d p r               -> Either (Point d r :+ p) (PointSet (LSeq 2) d p r)-asSingleton v = case v^.element (C :: C 0) of+asSingleton v = case v^.element @0 of                   (p :<| s) | null s -> Left p -- only one lement                   _                  -> Right $ unsafeCoerce v

@@ -16,12 +16,13 @@                          ) where  import           Control.Lens+import           Data.Bifunctor import           Data.Ext import           Data.Geometry.Boundary import           Data.Geometry.Box import           Data.Geometry.Line.Internal import           Data.Geometry.LineSegment-import           Data.Geometry.Point+import           Data.Geometry.Point.Internal import           Data.Geometry.Properties import           Data.Geometry.SubLine import           Data.Geometry.Transformation@@ -75,7 +76,7 @@   line' `intersect` (Boundary rect)  = case asAP segP of       [sl'] -> case fromUnbounded sl' of         Nothing   -> error "intersect: line x boundary rect; unbounded line? absurd"-        Just sl'' -> coRec $ sl''^.re _SubLine+        Just sl'' -> coRec $ first (either id id) $ sl''^.re _SubLine       []    -> case nub' $ asAP pointP of         [p]   -> coRec p         [p,q] -> coRec (p,q)@@ -97,7 +98,7 @@              => proxy t -> [t]       asAP _ = mapMaybe (asA @t) ints -      segP   = Proxy :: Proxy (SubLine 2 () (UnBounded r) r)+      segP   = Proxy :: Proxy (SubLine 2 (Either () ()) (UnBounded r) r)       pointP = Proxy :: Proxy (Point 2 r)  

@@ -15,7 +15,9 @@ import           Control.DeepSeq import           Control.Lens import qualified Data.Foldable as F-import           Data.Geometry.Point+import           Data.Geometry.Point.Internal+import           Data.Geometry.Point.Orientation.Degenerate+import           Data.Geometry.Point.Class import           Data.Geometry.Properties import           Data.Geometry.Vector import           Data.Ord (comparing)@@ -111,8 +113,17 @@ isParallelTo                         :: (Eq r, Fractional r, Arity d)                                      => Line d r -> Line d r -> Bool (Line _ u) `isParallelTo` (Line _ v) = u `isScalarMultipleOf` v-  -- TODO: Maybe use a specialize pragma for 2D (see intersect instance for two lines.)+{-# RULES+"isParallelTo/isParallelTo2" [3]+     forall (l1 :: forall r. Line 2 r) l2. isParallelTo l1 l2 = isParallelTo2 l1 l2+#-}+{-# INLINE[2] isParallelTo #-} +-- | Check whether two lines are parallel+isParallelTo2 :: (Eq r, Num r) => Line 2 r -> Line 2 r -> Bool+isParallelTo2 (Line _ (Vector2 ux uy)) (Line _ (Vector2 vx vy)) = denom == 0+    where+      denom       = vy * ux - vx * uy  -- | Test if point p lies on line l --@@ -130,9 +141,6 @@ p `onLine2` (Line q v) = ccw p q (q .+^ v) == CoLinear  --- -- | Get the point at the given position along line, where 0 corresponds to the -- anchorPoint of the line, and 1 to the point anchorPoint .+^ directionVector pointAt              :: (Num r, Arity d) => r -> Line d r -> Point d r@@ -154,7 +162,7 @@ -- >>> toOffset' (Point2 5 5) (lineThrough origin $ Point2 10 10) -- 0.5 ----- \<6,4\> is not on the line but we can still point closest to it.+-- The point (6,4) is not on the line but we can still point closest to it. -- >>> toOffset' (Point2 6 4) (lineThrough origin $ Point2 10 10) -- 0.5 toOffset'             :: (Eq r, Fractional r, Arity d) => Point d r -> Line d r -> r@@ -171,16 +179,14 @@                                                      , Line 2 r                                                      ] -instance (Eq r, Fractional r) => Line 2 r `HasIntersectionWith` Line 2 r--instance (Eq r, Fractional r) => Line 2 r `IsIntersectableWith` Line 2 r where-+instance (Ord r, Num r) => Line 2 r `HasIntersectionWith` Line 2 r where+  l1 `intersects` l2@(Line q _) = not (l1 `isParallelTo2` l2) || q `onLine2` l1 +instance (Ord r, Fractional r) => Line 2 r `IsIntersectableWith` Line 2 r where   nonEmptyIntersection = defaultNonEmptyIntersection-   l@(Line p ~(Vector2 ux uy)) `intersect` (Line q ~v@(Vector2 vx vy))-      | areParallel = if q `onLine` l then coRec l-                                      else coRec NoIntersection+      | areParallel = if q `onLine2` l then coRec l+                                       else coRec NoIntersection       | otherwise   = coRec r     where       r = q .+^ alpha *^ v@@ -229,7 +235,7 @@ {- HLINT ignore toLinearFunction -} -- | get values a,b s.t. the input line is described by y = ax + b. -- returns Nothing if the line is vertical-toLinearFunction                             :: forall r. (Fractional r, Eq r)+toLinearFunction                             :: forall r. (Fractional r, Ord r)                                              => Line 2 r -> Maybe (r,r) toLinearFunction l@(Line _ ~(Vector2 vx vy)) = match (l `intersect` verticalLine (0 :: r)) $        (H $ \NoIntersection -> Nothing)    -- l is a vertical line@@ -237,20 +243,14 @@     :& (H $ \_              -> Nothing)    -- l is a vertical line (through x=0)     :& RNil --- -- | get values a,b,c s.t. the input line is described by ax + by + c = 0--- toLinearFunction' :: Line 2 r -> (r,r,r)--- toLinearFunction' () --- | Given a point p and a line l, computes the point q on l closest to p.-pointClosestTo              :: (Fractional r, Arity d) => Point d r -> Line d r -> Point d r-pointClosestTo p (Line a m) = a .+^ (t0 *^ m)-  where-    -- see https://monkeyproofsolutions.nl/wordpress/how-to-calculate-the-shortest-distance-between-a-point-and-a-line/-    t0 = numerator / divisor-    numerator = (p .-. a) `dot` m-    divisor  = m `dot` m--+instance (Fractional r, Arity d) => HasSquaredEuclideanDistance (Line d r) where+  pointClosestTo p (Line a m) = a .+^ (t0 *^ m)+    where+      -- see https://monkeyproofsolutions.nl/wordpress/how-to-calculate-the-shortest-distance-between-a-point-and-a-line/+      t0 = numerator / divisor+      numerator = (p .-. a) `dot` m+      divisor  = m `dot` m   -- | Result of a side test

@@ -1,5 +1,6 @@ {-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE UndecidableInstances #-}+{-# OPTIONS_GHC -fno-warn-orphans #-} -------------------------------------------------------------------------------- -- | -- Module      :  Data.Geometry.LineSegment@@ -25,7 +26,113 @@   , flipSegment    , interpolate, sampleLineSegment+  , ordAtX, ordAtY, xCoordAt, yCoordAt   ) where -import Data.Geometry.Interval hiding (width, midPoint)-import Data.Geometry.LineSegment.Internal+-- import           Control.Lens+import           Data.Ext+-- import qualified Data.Foldable as F+import           Data.Geometry.Boundary+import           Data.Geometry.Box.Internal+import           Data.Geometry.Box.Sides+import           Data.Geometry.Interval hiding (width, midPoint)+import           Data.Geometry.LineSegment.Internal+import           Data.Geometry.Point+import           Data.Geometry.Properties+-- import           Data.Geometry.SubLine+import           Data.Util+-- import           Data.Vinyl.CoRec+-- import           Data.Bifunctor+-- import           Data.Either+-- import           Data.Maybe (mapMaybe)++++--------------------------------------------------------------------------------+++type instance IntersectionOf (LineSegment 2 p r) (Boundary (Rectangle q r)) =+  [ NoIntersection, Point 2 r, Two (Point 2 r) , LineSegment 2 () r ]+++type instance IntersectionOf (LineSegment 2 p r) (Rectangle q r) =+  [ NoIntersection, Point 2 r, LineSegment 2 (Maybe p) r ]++instance (Fractional r, Ord r)+         => LineSegment 2 p r `HasIntersectionWith` Boundary (Rectangle q r) where+  seg `intersects` (Boundary rect) = any (seg `intersects`) $ sides rect++instance (Fractional r, Ord r) => LineSegment 2 p r `HasIntersectionWith` Rectangle q r where+  seg@(LineSegment p q) `intersects` rect =+      inRect p || inRect q || any (seg `intersects`) (sides rect) || bothOpenAndOnBoundary seg+    where+      inRect = \case+        Open   (a :+ _) -> a `insideBox`  rect -- if strictly inside the seg intersects.+        Closed (a :+ _) -> a `inBox`      rect -- in or on the boundary is fine++      -- if somehow the segment is open, and both endpoints lie on+      -- different sides of the boundary, (so the segment crosses the+      -- interior) it also intersects. Handle that case.+      bothOpenAndOnBoundary (LineSegment (Open _) (Open _)) =+        interpolate (1/2) seg `intersects` rect+      bothOpenAndOnBoundary _                               = False++-- instance (Num r, Ord r)+--          => (LineSegment 2 p r) `IsIntersectableWith` (Boundary (Rectangle q r)) where+--   seg `intersect` (Boundary rect) = case partitionEithers res of+--     (s : _, _)    -> coRec s -- if we find a segment that should be the+--                              -- answer; we shouldn't fine more than one+--                              -- by the way.+--     ([], [])      -> coRec  NoIntersection+--     ([], [p])     -> coRec p+--     ([], (p:q:_)) -> coRec $ Two p q+--                      -- more than two points is impossible anwyay+--     where+--       res = mapMaybe (\side -> match (seg `intersect` side) $+--                        (H $ \NoIntersection            -> Nothing)+--                     :& (H $ \(p :: Point 2 r)          -> Just $ Right p)+--                     :& (H $ \(s :: LineSegment 2 () r) -> Just $ Left s)+--                     :& RNil+--              ) . F.toList $ sides rect++++-- -- instance (Num r, Ord r) => (LineSegment 2 p r) `IsIntersectableWith` (Rectangle q r) where+-- --   seg@(LineSegment' (p :+ _) (q :+ _)) `intersect` rect =+-- --       case (p `intersects` rect, q `intersects` rect) of+-- --         (True,True)   -> coRec seg'+-- --         (False,False) -> match boundaryIntersection $ -- both endpoints outside+-- --              (H $ \NoIntersection   -> coRec NoIntersection)+-- --           :& (H $ \(a :: Point 2 r) -> coRec a)+-- --           :& (H $ \(Two a b)        -> coRec $ ClosedLineSegment (ext a) (ext b))+-- --           :& (H $ \s                -> coRec s)+-- --           :& RNil+-- --         (True,False)  -> withInside p (\other -> LineSegment p' (closed other))+-- --         (False,True)  -> withInside q (\other -> LineSegment (closed other) q')+-- --     where+-- --       seg'@(LineSegment p' q') = first (const ()) seg++-- --       boundaryIntersection = seg `intersect` (Boundary rect)+-- --       closed :: Point 2 r -> EndPoint (Point 2 r :+ ())+-- --       closed = Closed . ext++-- --       -- the given endpoint endPt is inside the box [*], while the+-- --       -- other endpoint is not. The second arg is a function that+-- --       -- rebuilds the segment given the replacement endpoint, compute+-- --       -- the right segment that is inside the rectangle.+-- --       --+-- --       -- [*] We require that the *point* lies in or on the box. If the+-- --       -- endpoint was open, it may still be the case that we do not+-- --       -- actually intersect the rectangle (i.e. if the open endPoint+-- --       -- was on a corner of the rect).+-- --       -- withInside                      :: Point 2 r+-- --       --                                 -> (Point 2 r -> LineSegment 2 () r)+-- --       --                                 -> IntersectionOf ....+-- --       withInside endPt mkSeg = match boundaryIntersection $+-- --            (H $ \NoIntersection   -> coRec NoIntersection)+-- --            -- seems this should happen only if the endpoint that was+-- --            -- suposedly in/on the rect was open.+-- --         :& (H $ \(a :: Point 2 r) -> coRec . mkSeg $ a)+-- --         :& (H $ \(Two a b)        -> coRec . mkSeg $ if a == endPt then b else a)+-- --         :& (H $ \s                -> coRec s)+-- --         :& RNil

@@ -23,12 +23,14 @@   , orderedEndPoints   , segmentLength   , sqSegmentLength-  , sqDistanceToSeg, sqDistanceToSegArg+  , sqDistanceToSeg, sqDistanceToSegArg -- todo, at some point remove these. They are superfluous   , flipSegment    , interpolate   , validSegment   , sampleLineSegment++  , ordAtX, ordAtY, xCoordAt, yCoordAt   ) where  import           Control.Arrow ((&&&))@@ -46,12 +48,14 @@ import           Data.Geometry.Transformation.Internal import           Data.Geometry.Vector import           Data.Ord (comparing)+import           Data.Tuple (swap) import           Data.Vinyl import           Data.Vinyl.CoRec import           GHC.TypeLits import           Test.QuickCheck (Arbitrary(..), suchThatMap) import           Text.Read + -------------------------------------------------------------------------------- -- * d-dimensional LineSegments @@ -111,6 +115,7 @@  deriving instance (Arity d, NFData r, NFData p) => NFData (LineSegment d p r) +-- | Compute a random line segmeent sampleLineSegment :: (Arity d, RandomGen g, Random r) => Rand g (LineSegment d () r) sampleLineSegment = do   a <- ext <$> getRandom@@ -207,11 +212,14 @@ instance Arity d => Bifunctor (LineSegment d) where   bimap f g (GLineSegment i) = GLineSegment $ bimap f (fmap g) i +-- | Transform a segment into a closed line segment+toClosedSegment                    :: LineSegment d p r -> LineSegment d p r+toClosedSegment (LineSegment' s t) = ClosedLineSegment s t   -- ** Converting between Lines and LineSegments --- | Directly convert a line into a line segment.+-- | Directly convert a line into a Closed line segment. toLineSegment            :: (Monoid p, Num r, Arity d) => Line d r -> LineSegment d p r toLineSegment (Line p v) = ClosedLineSegment (p       :+ mempty)                                              (p .+^ v :+ mempty)@@ -222,11 +230,14 @@                                                                , Point d r                                                                ] -type instance IntersectionOf (LineSegment 2 p r) (LineSegment 2 p r) = [ NoIntersection-                                                                       , Point 2 r-                                                                       , LineSegment 2 p r-                                                                       ]+-- type instance IntersectionOf (LineSegment 2 p r) (LineSegment 2 p r) = [ NoIntersection+--                                                                        , Point 2 r+--                                                                        , LineSegment 2 p r+--                                                                        ] +type instance IntersectionOf (LineSegment 2 p r) (LineSegment 2 q r) =+  [ NoIntersection, Point 2 r, LineSegment 2 (Either p q) r]+ type instance IntersectionOf (LineSegment 2 p r) (Line 2 r) = [ NoIntersection                                                               , Point 2 r                                                               , LineSegment 2 p r@@ -272,11 +283,69 @@   -- work in higher dimensions that might allow us to drop the   -- Fractional constraint -instance (Ord r, Fractional r) =>-         LineSegment 2 p r `HasIntersectionWith` LineSegment 2 p r +-- | Orders the endpoints of the segments in the given direction.+withRank                                       :: forall p q r. (Ord r, Num r)+                                               => Vector 2 r+                                               -> LineSegment 2 p r  -> LineSegment 2 q r+                                               -> (Interval p Int, Interval q Int)+withRank v (LineSegment p q) (LineSegment a b) = (i1,i2)+  where+    -- let rank p = 3, rank q = 6+    i1 = Interval (p&unEndPoint.core .~ 3) (q&unEndPoint.core .~ 6)++    i2 = Interval (a&unEndPoint.core .~ assign' 1 a') (a&unEndPoint.core .~ assign' 2 b')++    -- make sure the intervals are in the same order, otherwise flip them.+    (a',b') = case cmp a b of+                LT -> (a,b)+                EQ -> (a,b)+                GT -> (b,a)++    assign' x c = case cmp c p of+                    LT -> x+                    EQ -> 3+                    GT -> case cmp c q of+                            LT -> 4 + x+                            EQ -> 6+                            GT -> 7 + x++    cmp     :: EndPoint (Point 2 r :+ a) -> EndPoint (Point 2 r :+ b) -> Ordering+    cmp c d = cmpInDirection v (c^.unEndPoint.core) (d^.unEndPoint.core)++instance (Ord r, Num r) =>+         LineSegment 2 p r `HasIntersectionWith` LineSegment 2 q r where+  s1@(LineSegment p _) `intersects` s2+    | l1 `isParallelTo2` l2 = parallelCase+    | otherwise             = s1 `intersects` l2  && s2 `intersects` l1+    where+      l1@(Line _ v) = supportingLine s1+      l2 = supportingLine s2++      parallelCase = (p^.unEndPoint.core) `onLine2` l2 && i1 `intersects` i2+      (i1,i2) = withRank v s1 s2++    -- correctness argument:+    -- if the segments share a supportingLine (l1 and l2 parallel, and point of l1 on l2)+    -- the segments intersect iff their intervals along the line intersect.++    -- if the supporting lines intersect in a point, say x the+    -- segments intersect iff s1 intersects the supporting line and+    -- vice versa:+    ---+    -- => direction: is trivial+    -- <= direction: s1 intersects l2 means x+    -- lies on s1. Symmetrically s2 intersects l1 means x lies on+    -- s2. Hence, x lies on both s1 and s2, and thus the segments+    -- intersect.++++++ instance (Ord r, Fractional r) =>-         LineSegment 2 p r `IsIntersectableWith` LineSegment 2 p r where+         LineSegment 2 p r `IsIntersectableWith` LineSegment 2 q r where   nonEmptyIntersection = defaultNonEmptyIntersection    a `intersect` b = match ((a^._SubLine) `intersect` (b^._SubLine)) $@@ -285,9 +354,17 @@       :& H (coRec . subLineToSegment)       :& RNil -instance (Ord r, Fractional r) =>+instance (Ord r, Num r) =>          LineSegment 2 p r `HasIntersectionWith` Line 2 r where+  (LineSegment p q) `intersects` l = case onSide (p^.unEndPoint.core) l of+    OnLine -> isClosed p || case onSide (q^.unEndPoint.core) l of+                              OnLine -> isClosed q || (p^.unEndPoint.core) /= (q^.unEndPoint.core)+                              _      -> False+    sp     -> case onSide (q^.unEndPoint.core) l of+                OnLine -> isClosed q+                sq     -> sp /= sq + instance (Ord r, Fractional r) =>          LineSegment 2 p r `IsIntersectableWith` Line 2 r where   nonEmptyIntersection = defaultNonEmptyIntersection@@ -348,27 +425,42 @@ segmentLength                     :: (Arity d, Floating r) => LineSegment d p r -> r segmentLength ~(LineSegment' p q) = distanceA (p^.core) (q^.core) +-- | Squared length of a line segment. sqSegmentLength                     :: (Arity d, Num r) => LineSegment d p r -> r sqSegmentLength ~(LineSegment' p q) = qdA (p^.core) (q^.core)  -- | Squared distance from the point to the Segment s. The same remark as for -- the 'sqDistanceToSegArg' applies here.+{-# DEPRECATED sqDistanceToSeg "use squaredEuclideanDistTo instead" #-} sqDistanceToSeg   :: (Arity d, Fractional r, Ord r) => Point d r -> LineSegment d p r -> r sqDistanceToSeg p = fst . sqDistanceToSegArg p - -- | Squared distance from the point to the Segment s, and the point on s--- realizing it.  Note that if the segment is *open*, the closest point--- returned may be one of the (open) end points, even though technically the--- end point does not lie on the segment. (The true closest point then lies--- arbitrarily close to the end point).-sqDistanceToSegArg     :: (Arity d, Fractional r, Ord r)-                       => Point d r -> LineSegment d p r -> (r, Point d r)-sqDistanceToSegArg p s = let m  = sqDistanceToArg p (supportingLine s)-                             xs = m : map (\(q :+ _) -> (qdA p q, q)) [s^.start, s^.end]-                         in   F.minimumBy (comparing fst)-                            . filter (flip onSegment s . snd) $ xs+-- realizing it.+--+-- Note that if the segment is *open*, the closest point returned may+-- be one of the (open) end points, even though technically the end+-- point does not lie on the segment. (The true closest point then+-- lies arbitrarily close to the end point).+--+-- >>> :{+-- let ls = OpenLineSegment (Point2 0 0 :+ ()) (Point2 1 0 :+ ())+--     p  = Point2 2 0+-- in  snd (sqDistanceToSegArg p ls) == Point2 1 0+-- :}+-- True+sqDistanceToSegArg                          :: (Arity d, Fractional r, Ord r)+                                            => Point d r -> LineSegment d p r -> (r, Point d r)+sqDistanceToSegArg p (toClosedSegment -> s) =+  let m  = sqDistanceToArg p (supportingLine s)+      xs = m : map (\(q :+ _) -> (qdA p q, q)) [s^.start, s^.end]+  in   F.minimumBy (comparing fst)+     . filter (flip onSegment s . snd) $ xs +instance (Fractional r, Arity d, Ord r) => HasSquaredEuclideanDistance (LineSegment d p r) where+  pointClosestToWithDistance q = swap . sqDistanceToSegArg q++ -- | flips the start and end point of the segment flipSegment   :: LineSegment d p r -> LineSegment d p r flipSegment s = let p = s^.start@@ -415,3 +507,45 @@                  -> Maybe (LineSegment d p r) validSegment u v = let s = LineSegment u v                    in if s^.start.core /= s^.end.core then Just s else Nothing++++-- | Given a y-coordinate, compare the segments based on the+-- x-coordinate of the intersection with the horizontal line through y+ordAtY   :: (Fractional r, Ord r) => r+         -> LineSegment 2 p r -> LineSegment 2 p r -> Ordering+ordAtY y = comparing (xCoordAt y)++-- | Given an x-coordinate, compare the segments based on the+-- y-coordinate of the intersection with the horizontal line through y+ordAtX   :: (Fractional r, Ord r) => r+         -> LineSegment 2 p r -> LineSegment 2 p r -> Ordering+ordAtX x = comparing (yCoordAt x)++-- | Given a y coord and a line segment that intersects the horizontal line+-- through y, compute the x-coordinate of this intersection point.+--+-- note that we will pretend that the line segment is closed, even if it is not+xCoordAt             :: (Fractional r, Ord r) => r -> LineSegment 2 p r -> r+xCoordAt y (LineSegment' (Point2 px py :+ _) (Point2 qx qy :+ _))+      | py == qy     = px `max` qx  -- s is horizontal, and since it by the+                                    -- precondition it intersects the sweep+                                    -- line, we return the x-coord of the+                                    -- rightmost endpoint.+      | otherwise    = px + alpha * (qx - px)+  where+    alpha = (y - py) / (qy - py)+++-- | Given an x-coordinate and a line segment that intersects the vertical line+-- through x, compute the y-coordinate of this intersection point.+--+-- note that we will pretend that the line segment is closed, even if it is not+yCoordAt :: (Fractional r, Ord r) => r -> LineSegment 2 p r -> r+yCoordAt x (LineSegment' (Point2 px py :+ _) (Point2 qx qy :+ _))+    | px == qx  = py `max` qy -- s is vertical, since by the precondition it+                              -- intersects we return the y-coord of the topmost+                              -- endpoint.+    | otherwise = py + alpha * (qy - py)+  where+    alpha = (x - px) / (qx - px)

@@ -24,7 +24,6 @@ import           Data.Geometry.PlanarSubdivision.Merge import           Data.Geometry.PlanarSubdivision.TreeRep import           Data.Geometry.Polygon-import           Data.Proxy   -- import Data.Geometry.Point@@ -42,25 +41,23 @@ -- -- runningtime: \(O(n\log n\log k)\) in case of polygons with holes, -- and \(O(n\log k)\) in case of simple polygons.-fromPolygons       :: (Foldable1 c, Ord r, Fractional r)-                   => proxy s-                   -> f -- ^ outer face data-                   -> c (Polygon t p r :+ f) -- ^ the disjoint polygons-                   -> PlanarSubdivision s p () f r-fromPolygons px oD = mergeAllWith const-                   . fmap (\(pg :+ iD) -> fromPolygon px pg iD oD) . toNonEmpty+fromPolygons    :: forall s c t p r f. (Foldable1 c, Ord r, Num r)+                => f -- ^ outer face data+                -> c (Polygon t p r :+ f) -- ^ the disjoint polygons+                -> PlanarSubdivision s p () f r+fromPolygons oD = mergeAllWith const+                . fmap (\(pg :+ iD) -> fromPolygon pg iD oD) . toNonEmpty  -- | Version of 'fromPolygons' that accepts 'SomePolygon's as input.-fromPolygons'      :: forall proxy c s p r f. (Foldable1 c, Ord r, Fractional r)-                   => proxy s-                   -> f -- ^ outer face data+fromPolygons'      :: forall s c p r f. (Foldable1 c, Ord r, Num r)+                   => f -- ^ outer face data                    -> c (SomePolygon p r :+ f) -- ^ the disjoint polygons                    -> PlanarSubdivision s p () f r-fromPolygons' px oD =+fromPolygons' oD =     mergeAllWith const . fmap (\(pg :+ iD) -> either (build iD) (build iD) pg) . toNonEmpty   where     build       :: f -> Polygon t p r -> PlanarSubdivision s p () f r-    build iD pg = fromPolygon px pg iD oD+    build iD pg = fromPolygon pg iD oD  -- | Construct a planar subdivision from a polygon. Since our PlanarSubdivision -- models only connected planar subdivisions, this may add dummy/invisible@@ -70,21 +67,19 @@ -- -- running time: \(O(n)\) for a simple polygon, \(O(n\log n)\) for a -- polygon with holes.-fromPolygon                              :: forall proxy t p f r s. (Ord r, Fractional r)-                                         => proxy s-                                         -> Polygon t p r+fromPolygon                              :: forall s t p f r. (Ord r, Num r)+                                         => Polygon t p r                                          -> f -- ^ data inside                                          -> f -- ^ data outside the polygon                                          -> PlanarSubdivision s p () f r-fromPolygon p pg@SimplePolygon{} iD oD   = fromSimplePolygon p pg iD oD-fromPolygon p (MultiPolygon vs hs) iD oD = case NonEmpty.nonEmpty hs of+fromPolygon pg@SimplePolygon{} iD oD   = fromSimplePolygon @s pg iD oD+fromPolygon (MultiPolygon vs hs) iD oD = case NonEmpty.nonEmpty hs of     Nothing  -> outerPG-    Just hs' -> let hs'' = (\pg -> fromSimplePolygon wp (toCounterClockWiseOrder pg) oD iD) <$> hs'+    Just hs' -> let hs'' = (\pg -> fromSimplePolygon @(Wrap s)+                                   (toCounterClockWiseOrder pg) oD iD) <$> hs'                 in embedAsHolesIn hs'' (\_ x -> x) i outerPG   where-    wp = Proxy :: Proxy (Wrap s)--    outerPG = fromSimplePolygon p vs iD oD+    outerPG = fromSimplePolygon @s vs iD oD     i = V.last $ faces' outerPG  

@@ -160,7 +160,7 @@ -- | Constructs a planarsubdivision from a PlaneGraph -- -- runningTime: \(O(n)\)-fromPlaneGraph   :: forall s v e f r. (Ord r, Fractional r)+fromPlaneGraph   :: forall s v e f r. (Ord r, Num r)                       => PlaneGraph s v e f r -> PlanarSubdivision s v e f r fromPlaneGraph g = fromPlaneGraph' g (PG.outerFaceDart g) @@ -212,24 +212,22 @@ -- -- pre: the input polygon is given in counterclockwise order -- running time: \(O(n)\).-fromSimplePolygon            :: (Ord r, Fractional r)-                             => proxy s-                             -> SimplePolygon p r+fromSimplePolygon            :: forall s p f r. (Ord r, Num r)+                             => SimplePolygon p r                              -> f -- ^ data inside                              -> f -- ^ data outside the polygon                              -> PlanarSubdivision s p () f r-fromSimplePolygon p pg iD oD =-  fromPlaneGraph (PG.fromSimplePolygon p pg iD oD)+fromSimplePolygon pg iD oD =+  fromPlaneGraph (PG.fromSimplePolygon pg iD oD)  -- | Constructs a connected planar subdivision. -- -- pre: the segments form a single connected component -- running time: \(O(n\log n)\)-fromConnectedSegments    :: (Foldable f, Ord r, Fractional r)-                         => proxy s-                         -> f (LineSegment 2 p r :+ e)-                         -> PlanarSubdivision s (NonEmpty p) e () r-fromConnectedSegments px = fromPlaneGraph . PG.fromConnectedSegments px+fromConnectedSegments :: forall s p e r f. (Foldable f, Ord r, Num r)+                      => f (LineSegment 2 p r :+ e)+                      -> PlanarSubdivision s (NonEmpty p) e () r+fromConnectedSegments = fromPlaneGraph . PG.fromConnectedSegments  -- g1 = PG.fromConnectedSegments (Identity Test1) testSegs -- ps1 = fromConnectedSegments (Identity Test1) testSegs

@@ -1,34 +1,30 @@-module Data.Geometry.PlanarSubdivision.Dynamic +module Data.Geometry.PlanarSubdivision.Dynamic   ( splitEdge, unSplitEdge   , sproutIntoFace   , splitFace   ) where -import Control.Lens--import Data.Vector (Vector, toList, (//), empty)-import qualified Data.Vector as V-import Data.List (sort, sortOn, findIndex)--import Data.Functor.Identity-import Data.Ext-import Data.Geometry hiding (Vector, head, imap)-import Data.Geometry.PlanarSubdivision-import Data.Geometry.PlanarSubdivision.Raw--import Data.PlanarGraph (Dart (Dart), Arc (Arc), VertexId (VertexId), FaceId (FaceId), Direction (Positive, Negative))-import Data.PlaneGraph (PlaneGraph)+import           Control.Lens+import           Data.Ext+import           Data.Functor.Identity+import           Data.Geometry hiding (Vector, head, imap)+import           Data.Geometry.PlanarSubdivision+import           Data.Geometry.PlanarSubdivision.Basic+import           Data.Geometry.PlanarSubdivision.Raw+import           Data.List (sort, sortOn, findIndex)+import           Data.List.NonEmpty (NonEmpty(..))+import qualified Data.List.NonEmpty as NonEmpty+import           Data.PlanarGraph (Dart (Dart), Arc (Arc), VertexId (VertexId), FaceId (FaceId), Direction (Positive, Negative))+import           Data.PlaneGraph (PlaneGraph) import qualified Data.PlaneGraph as PG-import Data.PlaneGraph.AdjRep hiding (id, vData, faces) import qualified Data.PlaneGraph.AdjRep as AR (id, vData, fData, faces, Face (..))+import           Data.PlaneGraph.AdjRep hiding (id, vData, faces)+import           Data.Vector (Vector, toList, (//), empty)+import qualified Data.Vector as V -import           Data.List.NonEmpty (NonEmpty(..))-import qualified Data.List.NonEmpty as NonEmpty+import           Debug.Trace -import Debug.Trace -import Data.Geometry.PlanarSubdivision.Basic- tracingOn = False  tr :: Show a => String -> a -> a@@ -44,17 +40,17 @@  -- | Splits a given edge of a planar subdivision by inserting a new vertex on the edges. --   Increases #vertices and #edges by 1.-splitEdge +splitEdge   :: (Show v, Show e, Show f, Show r)-  => VertexId' s -  -> VertexId' s -  -> Point 2 r -  -> v -  -> (e -> (e, e)) -  -> PlanarSubdivision s v e f r +  => VertexId' s+  -> VertexId' s+  -> Point 2 r+  -> v+  -> (e -> (e, e))   -> PlanarSubdivision s v e f r+  -> PlanarSubdivision s v e f r -splitEdge a b p v f d = +splitEdge a b p v f d =   let (_, la, _) = asLocalV a d       (_, lb, _) = asLocalV b d       v' = (freeVertexId d, v)@@ -67,13 +63,13 @@ --   Increases #vertices and #edges by 1. sproutIntoFace   :: (Show v, Show e, Show f, Show r)-  => VertexId' s -  -> FaceId' s -  -> Point 2 r -  -> v                       +  => VertexId' s+  -> FaceId' s+  -> Point 2 r+  -> v   -> (e, e)-  -> PlanarSubdivision s v e f r    -> PlanarSubdivision s v e f r+  -> PlanarSubdivision s v e f r  sproutIntoFace a f p v (e1, e2) d =   let [ea] = tr "[ea]" $ filter (\e -> headOf e d == a && leftFace e d == f) $ commonDarts d a f@@ -89,12 +85,12 @@ --   Increases #edges and #faces by 1. splitFace   :: (Show v, Show e, Show f, Show r)-  => VertexId' s -  -> VertexId' s -  -> (e, e)                       -  -> (f -> (f, f)) -  -> PlanarSubdivision s v e f r +  => VertexId' s+  -> VertexId' s+  -> (e, e)+  -> (f -> (f, f))   -> PlanarSubdivision s v e f r+  -> PlanarSubdivision s v e f r  splitFace a b e g d =   let (ca, _, _) = asLocalV a d@@ -124,14 +120,14 @@  -- | Splits a given edge of a planar subdivision by inserting a new vertex on the edges. --   Increases #vertices and #edges by 1.-unSplitEdge +unSplitEdge   :: (Show v, Show e, Show f, Show r)-  => VertexId' s +  => VertexId' s   -> ((e, e) -> e)-  -> PlanarSubdivision s v e f r    -> PlanarSubdivision s v e f r+  -> PlanarSubdivision s v e f r -unSplitEdge b f d = +unSplitEdge b f d =   let [a, c] = tr "[a, c]" $ toList $ neighboursOf b d       (_, la, _) = asLocalV a d       (_, lb, _) = asLocalV b d@@ -177,7 +173,7 @@                         & rawDartData   .~ (tr "rawDartData"   . vectorise $ getRawEdgeData cs)                         & rawFaceData   .~ (tr "rawFaceData"   . vectorise $ getRawFaceData cs) -getRawVertexData :: Vector (Component' s v e f r) +getRawVertexData :: Vector (Component' s v e f r)                  -> [(VertexId' s, Raw s (VertexId' (Wrap s)) v)] getRawVertexData = concat . imap (\ci g -> map (\(li, VertexData _ (gi, v)) -> (gi, Raw (toEnum ci) li v)) $ toList $ PG.vertices g) . toList @@ -194,7 +190,7 @@   -- data RawFace	s f--- _faceIdx :: !(Maybe (ComponentId s, FaceId' (Wrap s)))	 +-- _faceIdx :: !(Maybe (ComponentId s, FaceId' (Wrap s))) -- _faceDataVal :: !(FaceData (Dart s) f)  -- | Something in this implementation is not right. It makes asLocalF produce an error.@@ -229,36 +225,36 @@ -- INSERTIONS --  -splitEdgeInPlaneGraph -  :: (Show v, Show e, Show f, Show r) -  => VertexId' s -  -> VertexId' s -  -> Point 2 r -  -> v -  -> (e -> (e, e)) -  -> PlaneGraph s v e f r +splitEdgeInPlaneGraph+  :: (Show v, Show e, Show f, Show r)+  => VertexId' s+  -> VertexId' s+  -> Point 2 r+  -> v+  -> (e -> (e, e))   -> PlaneGraph s v e f r+  -> PlaneGraph s v e f r -- LET OP! TEST OF a EN b WEL VOORKOMEN!-splitEdgeInPlaneGraph a b p v f -  = tr "splitEdgeInPlaneGraph" -  . PG.fromAdjRep undefined -  . splitEdgeInAdjRep (fromEnum a) (fromEnum b) p v f +splitEdgeInPlaneGraph a b p v f+  = tr "splitEdgeInPlaneGraph"+  . PG.fromAdjRep+  . splitEdgeInAdjRep (fromEnum a) (fromEnum b) p v f   . PG.toAdjRep  sproutIntoFaceInPlaneGraph-  :: (Show v, Show e, Show f, Show r) -  => VertexId' s -  -> VertexId' s -  -> Point 2 r -  -> v +  :: (Show v, Show e, Show f, Show r)+  => VertexId' s+  -> VertexId' s+  -> Point 2 r+  -> v   -> (e, e)-  -> PlaneGraph s v e f r    -> PlaneGraph s v e f r+  -> PlaneGraph s v e f r sproutIntoFaceInPlaneGraph a c p v e g =   let ai = fromEnum a       ci = fromEnum c-  in tr "splitEdgeInPlaneGraph" -   $ PG.fromAdjRep undefined +  in tr "splitEdgeInPlaneGraph"+   $ PG.fromAdjRep    $ sproutInAdjRep ai ci p v e    $ PG.toAdjRep g @@ -279,7 +275,7 @@   -> PlaneGraph s v e f r -- input graaf   -> PlaneGraph s v e f r -- output graaf -splitFaceInPlaneGraph a b c d f e h g = +splitFaceInPlaneGraph a b c d f e h g =   let ai = fromEnum a       bi = fromEnum b       ci = fromEnum c@@ -287,28 +283,28 @@       fi = fromEnum $ tr "fi" $ traceShow (g ^. dataOf f) $ PG.tailOf (PG.boundaryDart f g) g       fj = fromEnum $ tr "fj" $ PG.headOf (PG.boundaryDart f g) g       -- ^ boundaryDart seems not working either-  in tr "splitFaceInPlaneGraph" -   $ PG.fromAdjRep undefined -   $ splitFaceInAdjRep ai bi ci di fi fj e h +  in tr "splitFaceInPlaneGraph"+   $ PG.fromAdjRep+   $ splitFaceInAdjRep ai bi ci di fi fj e h    $ PG.toAdjRep g   -- DELETIONS --  -unSplitEdgeInPlaneGraph -  :: (Show v, Show e, Show f, Show r) -  => VertexId' s -  -> VertexId' s -  -> VertexId' s -  -> ((e, e) -> e) -  -> PlaneGraph s v e f r +unSplitEdgeInPlaneGraph+  :: (Show v, Show e, Show f, Show r)+  => VertexId' s+  -> VertexId' s+  -> VertexId' s+  -> ((e, e) -> e)   -> PlaneGraph s v e f r+  -> PlaneGraph s v e f r -unSplitEdgeInPlaneGraph a b c f -  = tr "unSplitEdgeInPlaneGraph" -  . PG.fromAdjRep undefined -  . unSplitEdgeInAdjRep (fromEnum a) (fromEnum b) (fromEnum c) f +unSplitEdgeInPlaneGraph a b c f+  = tr "unSplitEdgeInPlaneGraph"+  . PG.fromAdjRep+  . unSplitEdgeInAdjRep (fromEnum a) (fromEnum b) (fromEnum c) f   . PG.toAdjRep  @@ -316,18 +312,18 @@ -- ADJREPS -- ------------- --- Gr   --- adjacencies :: [v]  --- faces :: [f]   +-- Gr+-- adjacencies :: [v]+-- faces :: [f] --- Vtx  --- id :: Int  --- loc :: Point 2 r   --- adj :: [(Int, e)] --- vData :: v   +-- Vtx+-- id :: Int+-- loc :: Point 2 r+-- adj :: [(Int, e)]+-- vData :: v --- Face   --- incidentEdge :: (Int, Int)   +-- Face+-- incidentEdge :: (Int, Int) -- fData :: f  --deriving instance (Show v, Show f) => Show (Gr v f)@@ -337,7 +333,7 @@  -- instance {-# OVERLAPS #-} Show (VertexId s Primal) where show i = 'v' : show (fromEnum i) -- instance {-# OVERLAPS #-} Show (FaceId   s Primal) where show i = 'f' : show (fromEnum i)--- instance {-# OVERLAPS #-} Show (Dart s, v) where +-- instance {-# OVERLAPS #-} Show (Dart s, v) where --   show (Dart (Arc s) Positive, _) = 'd' : show (fromEnum s) ++ "+" --   show (Dart (Arc s) Negative, _) = 'd' : show (fromEnum s) ++ "-" @@ -346,7 +342,7 @@ -- instance (Show v, Show f) => Show (Gr v f) where show g = "Gr " ++ (show $ adjacencies g) ++ " " ++ (show $ AR.faces g)  -- ik heb:-splitEdgeInAdjRep +splitEdgeInAdjRep   :: (Show v, Show e, Show f, Show r)   => Int                     -- index van vertex a   -> Int                     -- index van vertex b@@ -356,7 +352,7 @@   -> Gr (Vtx v e r) (Face f) -- input graaf   -> Gr (Vtx v e r) (Face f) -- output graaf -splitEdgeInAdjRep a b p v f g = +splitEdgeInAdjRep a b p v f g =   let n  = length $ adjacencies g       -- first find vertices a and b       oa = headTrace "splitEdgeInAdjRep oa" $ filter ((== a) . AR.id) $ adjacencies g@@ -371,12 +367,12 @@       -- create new vertex c       nc = Vtx {AR.id = n, loc = p, adj = [(a, snd $ f e2), (b, snd $ f e1)], AR.vData = v}       -- update faces (only if incidentEdge happens to point to ab)-      nf = replace ((== (a, b)) . incidentEdge) (\f -> f {incidentEdge = (a, n)}) -         $ replace ((== (b, a)) . incidentEdge) (\f -> f {incidentEdge = (b, n)}) +      nf = replace ((== (a, b)) . incidentEdge) (\f -> f {incidentEdge = (a, n)})+         $ replace ((== (b, a)) . incidentEdge) (\f -> f {incidentEdge = (b, n)})          $ AR.faces g   in tr "splitEdgeInAdjRep" $ (tr "original" g) {adjacencies = sortOn AR.id $ na : nb : nc : os, AR.faces = nf}-   + sproutInAdjRep   :: (Show v, Show e, Show f, Show r)   => Int                     -- index van vertex a@@ -394,7 +390,7 @@       os = tr "os" $ filter ((/= a) . AR.id) $ adjacencies g       -- need to find index of c       fj (Just x) = x-      fj Nothing  = error "splitFaceInAdjRep got Nothing"      +      fj Nothing  = error "splitFaceInAdjRep got Nothing"       ci = tr "ci" $ fj $ findIndex ((== c) . fst) $ adj oa       -- create new adjacency to new vertex z in a       na = tr "na" $ oa {adj = take ci (adj oa) ++ (n, fst e) : drop ci (adj oa)}@@ -402,7 +398,7 @@       nz = Vtx {AR.id = n, loc = p, adj = [(a, snd e)], AR.vData = v}   in tr "splitFaceInAdjRep" $ (tr "original" g) {adjacencies = sortOn AR.id $ na : nz : os} -splitFaceInAdjRep +splitFaceInAdjRep   :: (Show v, Show e, Show f, Show r)   => Int                     -- index van vertex a   -> Int                     -- index van vertex b@@ -418,14 +414,14 @@ -- is it easier to split a vertex than a face?  splitFaceInAdjRep a b c d u v e f g =-  let +  let       -- first find vertices a and b       oa = tr "oa" $ headTrace "splitFaceInAdjRep oa" $ filter ((== a) . AR.id) $ adjacencies g       ob = tr "ob" $ headTrace "splitFaceInAdjRep ob" $ filter ((== b) . AR.id) $ adjacencies g       os = tr "os" $ filter ((lift (&&) (/= a) (/= b)) . AR.id) $ adjacencies g       -- insert new adjacency between a and b       fj (Just x) = x-      fj Nothing  = error "splitFaceInAdjRep got Nothing"      +      fj Nothing  = error "splitFaceInAdjRep got Nothing"       -- need to find indices c and d!       ci = tr "ci" $ fj $ findIndex ((== c) . fst) $ adj oa       di = tr "di" $ fj $ findIndex ((== d) . fst) $ adj ob@@ -439,12 +435,12 @@       f1 = tr "f1" $ AR.Face {incidentEdge = (a, b), AR.fData = fst $ f fd}       f2 = tr "f2" $ AR.Face {incidentEdge = (b, a), AR.fData = snd $ f fd}   in tr "splitFaceInAdjRep" $ (tr "original" g) {adjacencies = sortOn AR.id $ na : nb : os, AR.faces = ef ++ [f1, f2]}-      -unSplitEdgeInAdjRep ++unSplitEdgeInAdjRep   :: (Show v, Show e, Show f, Show r)   => Int                     -- index van vertex a   -> Int                     -- index van vertex b (te verwijderen)@@ -453,7 +449,7 @@   -> Gr (Vtx v e r) (Face f) -- input graaf   -> Gr (Vtx v e r) (Face f) -- output graaf -unSplitEdgeInAdjRep a b c f g = +unSplitEdgeInAdjRep a b c f g =   let n  = length $ adjacencies g       -- first find vertices a, b and c       oa = head $ filter ((== a) . AR.id) $ adjacencies g@@ -470,27 +466,27 @@       nc = oc {adj = replace ((== b) . fst) (const (a, f (ecb, eba))) $ adj oc}       nv = sortOn AR.id $ na : nc : os       -- update faces (only if incidentEdge happens to point to ab or bc)-      nf = replace ((== (a, b)) . incidentEdge) (\f -> f {incidentEdge = (a, c)}) -         $ replace ((== (b, a)) . incidentEdge) (\f -> f {incidentEdge = (c, a)}) -         $ replace ((== (b, c)) . incidentEdge) (\f -> f {incidentEdge = (a, c)}) -         $ replace ((== (c, b)) . incidentEdge) (\f -> f {incidentEdge = (c, a)}) +      nf = replace ((== (a, b)) . incidentEdge) (\f -> f {incidentEdge = (a, c)})+         $ replace ((== (b, a)) . incidentEdge) (\f -> f {incidentEdge = (c, a)})+         $ replace ((== (b, c)) . incidentEdge) (\f -> f {incidentEdge = (a, c)})+         $ replace ((== (c, b)) . incidentEdge) (\f -> f {incidentEdge = (c, a)})          $ AR.faces g       -- restore consecutive numbering: replace vertex n-1 by b       ng = replaceIndex (n - 1) b $ (tr "original" g) {adjacencies = nv, AR.faces = nf}   in tr "unSplitEdgeInAdjRep" $ ng --- Gr   --- adjacencies :: [v]  --- faces :: [f]   +-- Gr+-- adjacencies :: [v]+-- faces :: [f] --- Vtx  --- id :: Int  --- loc :: Point 2 r   --- adj :: [(Int, e)] --- vData :: v   +-- Vtx+-- id :: Int+-- loc :: Point 2 r+-- adj :: [(Int, e)]+-- vData :: v --- Face   --- incidentEdge :: (Int, Int)   +-- Face+-- incidentEdge :: (Int, Int) -- fData :: f  replaceIndex :: Int -> Int -> Gr (Vtx v e r) (Face f) -> Gr (Vtx v e r) (Face f)

@@ -211,19 +211,17 @@  -------------------------------------------------------------------------------- -data Test = Test-newtype Id a = Id a-+data Test  triangle1 :: PlanarSubdivision Test () () Int Rational-triangle1 = (\pg -> fromSimplePolygon (Id Test) pg 1 0)+triangle1 = (\pg -> fromSimplePolygon @Test pg 1 0)           trianglePG1 trianglePG1 :: SimplePolygon () Rational trianglePG1 = fromPoints . map ext $ [origin, Point2 200 0, Point2 200 200]   triangle2 :: PlanarSubdivision Test () () Int Rational-triangle2 = (\pg -> fromSimplePolygon (Id Test) pg 2 0)+triangle2 = (\pg -> fromSimplePolygon @Test pg 2 0)           trianglePG2 trianglePG2 :: SimplePolygon () Rational trianglePG2 = fromPoints . map ext $ [Point2 0 30, Point2 10 30, Point2 10 40]@@ -231,13 +229,13 @@   triangle4 :: PlanarSubdivision Test () () Int Rational-triangle4 = (\pg -> fromSimplePolygon (Id Test) pg 1 0)+triangle4 = (\pg -> fromSimplePolygon @Test pg 1 0)           trianglePG4 trianglePG4 :: SimplePolygon () Rational trianglePG4 = fromPoints . map ext $ [Point2 400 400, Point2 600 400, Point2 600 600]  triangle3 :: PlanarSubdivision Test () () Int Rational-triangle3 = (\pg -> fromSimplePolygon (Id Test) pg 3 0)+triangle3 = (\pg -> fromSimplePolygon @Test pg 3 0)           trianglePG3 trianglePG3 :: SimplePolygon () Rational trianglePG3 = fromPoints . map ext $ [Point2 401 530, Point2 410 530, Point2 410 540]

@@ -31,9 +31,10 @@                            , Quadrant(..), quadrantWith, quadrant, partitionIntoQuadrants -                          , cmpByDistanceTo, cmpByDistanceTo'+                          , cmpByDistanceTo, cmpByDistanceTo', cmpInDirection                            , squaredEuclideanDist, euclideanDist+                          , HasSquaredEuclideanDistance(..)                            , coord, unsafeCoord                           ) where@@ -42,3 +43,27 @@ import Data.Geometry.Point.Internal hiding (coord, unsafeCoord) import Data.Geometry.Point.Orientation.Degenerate import Data.Geometry.Point.Quadrants+import Data.Geometry.Line.Internal+import Data.Geometry.Vector++--------------------------------------------------------------------------------++-- | Compare the points with respect to the direction given by the+-- vector, i.e. by taking planes whose normal is the given vector.+--+-- >>> cmpInDirection (Vector2 1 0) (Point2 5 0) (Point2 10 0)+-- LT+-- >>> cmpInDirection (Vector2 1 1) (Point2 5 0) (Point2 10 0)+-- LT+-- >>> cmpInDirection (Vector2 1 1) (Point2 5 0) (Point2 10 10)+-- LT+-- >>> cmpInDirection (Vector2 1 1) (Point2 15 15) (Point2 10 10)+-- GT+-- >>> cmpInDirection (Vector2 1 0) (Point2 15 15) (Point2 15 10)+-- EQ+cmpInDirection       :: (Ord r, Num r) => Vector 2 r -> Point 2 r -> Point 2 r -> Ordering+cmpInDirection n p q = case p `onSide` perpendicularTo (Line q n) of+                         LeftSide  -> LT+                         OnLine    -> EQ+                         RightSide -> GT+  -- TODO: Generalize to arbitrary dimension

@@ -1,3 +1,4 @@+{-# LANGUAGE  AllowAmbiguousTypes  #-} module Data.Geometry.Point.Class where  import           Control.Lens@@ -28,14 +29,14 @@  -- | Get the coordinate in a given dimension ----- >>> Point3 1 2 3 ^. coord (C :: C 2)+-- >>> Point3 1 2 3 ^. coord @2 -- 2--- >>> Point3 1 2 3 & coord (C :: C 1) .~ 10+-- >>> Point3 1 2 3 & coord @1 .~ 10 -- Point3 10 2 3--- >>> Point3 1 2 3 & coord (C :: C 3) %~ (+1)+-- >>> Point3 1 2 3 & coord @3 %~ (+1) -- Point3 1 2 4-coord   :: (1 <= i, i <= d, KnownNat i, Arity d, AsAPoint p) => proxy i -> Lens' (p d r) r-coord i = asAPoint.Internal.coord i+coord :: forall i p d r. (1 <= i, i <= d, KnownNat i, Arity d, AsAPoint p) => Lens' (p d r) r+coord = asAPoint.Internal.coord @i  -- | Get the coordinate in a given dimension. This operation is unsafe in the -- sense that no bounds are checked. Consider using `coord` instead.@@ -60,7 +61,7 @@ -- >>> Point2 1 2 & xCoord .~ 10 -- Point2 10 2 xCoord :: (1 <= d, Arity d, AsAPoint point) => Lens' (point d r) r-xCoord = coord (C :: C 1)+xCoord = coord @1 {-# INLINABLE xCoord #-}  -- | Shorthand to access the second coordinate C 2@@ -70,7 +71,7 @@ -- >>> Point3 1 2 3 & yCoord %~ (+1) -- Point3 1 3 3 yCoord :: (2 <= d, Arity d, AsAPoint point) => Lens' (point d r) r-yCoord = coord (C :: C 2)+yCoord = coord @2 {-# INLINABLE yCoord #-}  -- | Shorthand to access the third coordinate C 3@@ -80,5 +81,5 @@ -- >>> Point3 1 2 3 & zCoord %~ (+1) -- Point3 1 2 4 zCoord :: (3 <= d, Arity d, AsAPoint point) => Lens' (point d r) r-zCoord = coord (C :: C 3)+zCoord = coord @3 {-# INLINABLE zCoord #-}

@@ -1,5 +1,6 @@ {-# LANGUAGE ScopedTypeVariables  #-} {-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE AllowAmbiguousTypes #-} -------------------------------------------------------------------------------- -- | -- Module      :  Data.Geometry.Point@@ -27,6 +28,7 @@   , cmpByDistanceTo   , cmpByDistanceTo'   , squaredEuclideanDist, euclideanDist+  , HasSquaredEuclideanDistance(..)   ) where  import           Control.DeepSeq@@ -175,15 +177,15 @@  -- | Get the coordinate in a given dimension ----- >>> Point3 1 2 3 ^. coord (C :: C 2)+-- >>> Point3 1 2 3 ^. coord @2 -- 2--- >>> Point3 1 2 3 & coord (C :: C 1) .~ 10+-- >>> Point3 1 2 3 & coord @1 .~ 10 -- Point3 10 2 3--- >>> Point3 1 2 3 & coord (C :: C 3) %~ (+1)+-- >>> Point3 1 2 3 & coord @3 %~ (+1) -- Point3 1 2 4-coord   :: forall proxy i d r. (1 <= i, i <= d, Arity d, KnownNat i)-        => proxy i -> Lens' (Point d r) r-coord _ = unsafeCoord $ fromIntegral (natVal $ C @i)+coord :: forall i d r. (1 <= i, i <= d, Arity d, KnownNat i)+      => Lens' (Point d r) r+coord = unsafeCoord $ fromIntegral (natVal $ C @i) {-# INLINABLE coord #-}   -- somehow these rules don't fire@@ -242,7 +244,13 @@   pmap f = f  + --------------------------------------------------------------------------------+++++-------------------------------------------------------------------------------- -- * Functions specific to Two Dimensional points  -- | Compare by distance to the first argument@@ -256,7 +264,6 @@ cmpByDistanceTo' c p q = cmpByDistanceTo (c^.core) (p^.core) (q^.core)  - -- | Squared Euclidean distance between two points squaredEuclideanDist :: (Num r, Arity d) => Point d r -> Point d r -> r squaredEuclideanDist = qdA@@ -264,3 +271,33 @@ -- | Euclidean distance between two points euclideanDist :: (Floating r, Arity d) => Point d r -> Point d r -> r euclideanDist = distanceA+++--------------------------------------------------------------------------------+-- * Distances++class HasSquaredEuclideanDistance g where+  -- | Given a point q and a geometry g, the squared Euclidean distance between q and g.+  squaredEuclideanDistTo   :: (Num (NumType g), Arity (Dimension g))+                           => Point (Dimension g) (NumType g) -> g -> NumType g+  squaredEuclideanDistTo q = snd . pointClosestToWithDistance q++  -- | Given q and g, computes the point p in g closest to q according+  -- to the Squared Euclidean distance.+  pointClosestTo   :: (Num (NumType g), Arity (Dimension g))+                   => Point (Dimension g) (NumType g) -> g+                   -> Point (Dimension g) (NumType g)+  pointClosestTo q = fst . pointClosestToWithDistance q++  -- | Given q and g, computes the point p in g closest to q according+  -- to the Squared Euclidean distance. Returns both the point and the+  -- distance realized by this point.+  pointClosestToWithDistance     :: (Num (NumType g), Arity (Dimension g))+                                 => Point (Dimension g) (NumType g) -> g+                                 -> (Point (Dimension g) (NumType g), NumType g)+  pointClosestToWithDistance q g = let p = pointClosestTo q g+                                   in (p, squaredEuclideanDist p q)+  {-# MINIMAL pointClosestToWithDistance | pointClosestTo #-}++instance (Num r, Arity d) => HasSquaredEuclideanDistance (Point d r) where+  pointClosestTo _ p = p

@@ -32,7 +32,6 @@ import           Data.Geometry.Point import           Data.Geometry.Polygon import qualified Data.List.NonEmpty as NonEmpty-import           Data.Proxy import           Data.Util (SP(..)) import qualified Data.Vector as V @@ -153,7 +152,7 @@ -- type Vertex v r = Int :+ (Point 2 r :+ v)  inPolygonDS    :: (Fractional r, Ord r) => SimplePolygon v r -> InPolygonDS v r-inPolygonDS pg = pointLocationDS $ fromSimplePolygon (Proxy @Dummy) (numberVertices pg) In Out+inPolygonDS pg = pointLocationDS $ fromSimplePolygon @Dummy (numberVertices pg) In Out  -- | Finds the edge on or above the query point, if it exists --

@@ -24,6 +24,7 @@ import           Data.LSeq (LSeq, pattern (:<|)) import qualified Data.LSeq as LSeq import qualified Data.List.NonEmpty as NE+import           Data.Ord (comparing) import           GHC.Generics (Generic) import           GHC.TypeLits @@ -88,13 +89,19 @@   type EndExtra (PolyLine d p r) = p   end = points.last1 +instance (Fractional r, Arity d, Ord r) => HasSquaredEuclideanDistance (PolyLine d p r) where+  pointClosestToWithDistance q = F.minimumBy (comparing snd)+                               . fmap (pointClosestToWithDistance q)+                               . edgeSegments++ -- | Builds a Polyline from a list of points, if there are sufficiently many points fromPoints :: [Point d r :+ p] -> Maybe (PolyLine d p r) fromPoints = fmap PolyLine . LSeq.eval @2 . LSeq.fromList  -- | pre: The input list contains at least two points fromPointsUnsafe :: [Point d r :+ p] -> PolyLine d p r-fromPointsUnsafe = PolyLine . LSeq.forceLSeq (C @ 2) . LSeq.fromList+fromPointsUnsafe = PolyLine . LSeq.forceLSeq (C @2) . LSeq.fromList  -- | pre: The input list contains at least two points. All extra vields are -- initialized with mempty.

@@ -85,14 +85,15 @@ import           Control.Monad.Random.Class import           Data.Ext import qualified Data.Foldable as F+import           Data.Geometry.Boundary import           Data.Geometry.HalfSpace (rightOf) import           Data.Geometry.Line import           Data.Geometry.LineSegment import           Data.Geometry.Point-import           Data.Geometry.Boundary import           Data.Geometry.Polygon.Core import           Data.Geometry.Polygon.Extremes import           Data.Geometry.Properties+import           Data.Ord (comparing) import qualified Data.Sequence as Seq  --------------------------------------------------------------------------------@@ -141,3 +142,13 @@   --   where   --     unpack (CoRec x) = x   --     f = undefined++instance (Fractional r, Ord r) => HasSquaredEuclideanDistance (Boundary (Polygon t p r)) where+  pointClosestToWithDistance q = F.minimumBy (comparing snd)+                               . fmap (pointClosestToWithDistance q)+                               . listEdges . review _Boundary++instance (Fractional r, Ord r) => HasSquaredEuclideanDistance (Polygon t p r) where+  pointClosestToWithDistance q pg+    | q `intersects` pg = (q, 0)+    | otherwise         = pointClosestToWithDistance q (Boundary pg)

@@ -84,7 +84,7 @@ import           Data.Bifunctor import           Data.Bitraversable import           Data.Ext-import qualified Data.Foldable                                              as F+import qualified Data.Foldable as F import           Data.Geometry.Boundary import           Data.Geometry.Box                                          (IsBoxable (..),                                                                              boundingBoxList')@@ -98,18 +98,18 @@ import           Data.Geometry.Vector                                       (Additive (zero, (^+^)),                                                                              Affine ((.+^), (.-.)),                                                                              (*^), (^*), (^/))-import qualified Data.List                                                  as List-import qualified Data.List.NonEmpty                                         as NonEmpty-import           Data.Maybe                                                 (catMaybes)-import           Data.Ord                                                   (comparing)-import           Data.Semigroup                                             (sconcat)+import qualified Data.List as List+import qualified Data.List.NonEmpty as NonEmpty+import           Data.Maybe (catMaybes)+import           Data.Ord (comparing)+import           Data.Semigroup (sconcat) import           Data.Semigroup.Foldable import           Data.Util-import           Data.Vector                                                (Vector)-import qualified Data.Vector                                                as V-import           Data.Vector.Circular                                       (CircularVector)-import qualified Data.Vector.Circular                                       as CV-import qualified Data.Vector.Circular.Util                                  as CV+import           Data.Vector (Vector)+import qualified Data.Vector as V+import           Data.Vector.Circular (CircularVector)+import qualified Data.Vector.Circular as CV+import qualified Data.Vector.Circular.Util as CV   -- import Data.RealNumber.Rational@@ -272,8 +272,6 @@       pMulti  o = (\vs hs -> MultiPolygon (fromPoints vs) (map fromPoints hs))                <$> o .: "outerBoundary" <*> o .: "holes" -- -- * Functions on Polygons  -- | Getter access to the outer boundary vector of a polygon.@@ -386,7 +384,7 @@ -- | \( O(n \log n) \) Check if a polygon has any holes, duplicate points, or --   self-intersections. isSimple :: (Ord r, Fractional r) => Polygon p t r -> Bool-isSimple p@SimplePolygon{}   = null . BO.interiorIntersections $ listEdges p+isSimple p@SimplePolygon{}   = null . BO.interiorIntersections . map ext $ listEdges p isSimple (MultiPolygon b []) = isSimple b isSimple MultiPolygon{}      = False @@ -429,7 +427,7 @@   => CircularVector (Point 2 r :+ p) -> SimplePolygon p r simpleFromCircularVector v =   let p = fromCircularVector v-      hasInteriorIntersections = not . null . BO.interiorIntersections+      hasInteriorIntersections = not . null . BO.interiorIntersections . map ext   in if hasInteriorIntersections (listEdges p)       then error "Data.Geometry.Polygon.simpleFromCircularVector: \                  \Found self-intersections or repeated vertices."

@@ -18,7 +18,6 @@ import           Data.List.NonEmpty (NonEmpty(..)) import qualified Data.List.NonEmpty as NonEmpty import           Data.Measured.Class-import           Data.Proxy import           Data.Range import           GHC.TypeLits import           Prelude hiding (last,init,head)@@ -115,7 +114,7 @@                                     , 1 <= d -- this one is kind of silly                  ) => NonEmpty (Point d r :+ p) -> RT 2 d v p r createRangeTree2 = RangeTree . GRT.createTree-                 . fmap (\p -> p^.core.coord (Proxy :: Proxy 2) :+ Leaf [p])+                 . fmap (\p -> p^.core.coord @2 :+ Leaf [p])  -------------------------------------------------------------------------------- -- * Querying@@ -132,11 +131,11 @@ instance (1 <= d, Arity d) => Query 1 d where   search' qr = map unAssoc . GRT.search' r . _unRangeTree     where-      r = qr^.element (Proxy :: Proxy 0)+      r = qr^.element @0  instance ( 1 <= d, i <= d, Query (i-1) d, Arity d          , i ~ 2          ) => Query 2 d where   search' qr = concatMap (maybe [] (search' qr) . unAssoc) . GRT.search' r . _unRangeTree     where-      r = qr^.element (Proxy :: Proxy (i-1))+      r = qr^.element @(i-1)

@@ -1,4 +1,3 @@-{-# OPTIONS_GHC -Wno-orphans #-} -------------------------------------------------------------------------------- -- | -- Module      :  Data.Geometry.RangeTree.Measure@@ -53,11 +52,11 @@ instance (LabeledMeasure l, LabeledMeasure r) => LabeledMeasure (l :*: r) where   labeledMeasure xs = Pair (labeledMeasure xs) (labeledMeasure xs) -instance (Semigroup (l a), Semigroup (r a)) => Semigroup ((l :*: r) a) where-  (Pair l r) <> (Pair l' r') = Pair (l <> l') (r <> r')+-- instance (Semigroup (l a), Semigroup (r a)) => Semigroup ((l :*: r) a) where+--   (Pair l r) <> (Pair l' r') = Pair (l <> l') (r <> r') -instance (Monoid (l a), Monoid (r a)) => Monoid ((l :*: r) a) where-  mempty = Pair mempty mempty+-- instance (Monoid (l a), Monoid (r a)) => Monoid ((l :*: r) a) where+--   mempty = Pair mempty mempty   

@@ -26,11 +26,13 @@  -------------------------------------------------------------------------------- +-- | Orthogonal directions data Orthogonal = Horizontal | Vertical                 deriving (Show,Eq,Read)  -+-- | An strip between two parallel lines. The lines can be either+-- horizontal or vertical. newtype Slab (o :: Orthogonal) a r = Slab { _unSlab :: Interval a r }                                      deriving (Show,Eq) makeLenses ''Slab@@ -57,53 +59,52 @@   bimap f g (Slab i) = Slab $ bimap f g i  -type instance IntersectionOf (Slab o a r)          (Slab o a r) =-  [NoIntersection, Slab o a r]-type instance IntersectionOf (Slab Horizontal a r) (Slab Vertical a r) =-  '[Rectangle (a,a) r]+type instance IntersectionOf (Slab o a r) (Slab o b r) =+  [NoIntersection, Slab o (Either a b) r]+type instance IntersectionOf (Slab Horizontal a r) (Slab Vertical b r) =+  '[Rectangle (a,b) r]  -instance Ord r => Slab o a r `HasIntersectionWith` Slab o a r+instance Ord r => Slab o a r `HasIntersectionWith` Slab o b r -instance Ord r => Slab o a r `IsIntersectableWith` Slab o a r where+instance Ord r => Slab o a r `IsIntersectableWith` Slab o b r where   nonEmptyIntersection = defaultNonEmptyIntersection    (Slab i) `intersect` (Slab i') = match (i `intersect` i') $-        H (\NoIntersection -> coRec NoIntersection)-     :& H (\i''            -> coRec (Slab i'' :: Slab o a r))+        H (\NoIntersection                   -> coRec NoIntersection)+     :& H (\i''                              -> coRec $ (Slab i'' :: Slab o (Either a b) r))      :& RNil -instance Slab Horizontal a r `HasIntersectionWith` Slab Vertical a r where+instance Slab Horizontal a r `HasIntersectionWith` Slab Vertical b r where   _ `intersects` _ = True -instance Slab Horizontal a r `IsIntersectableWith` Slab Vertical a r where+instance Slab Horizontal a r `IsIntersectableWith` Slab Vertical b r where   nonEmptyIntersection _ _ _ = True    (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) :+ (h^.start.extra, v^.start.extra)+      high = Point2 (v^.end.core)   (h^.end.core)   :+ (h^.end.extra, v^.end.extra)    class HasBoundingLines (o :: Orthogonal) where   -- | The two bounding lines of the slab, first the lower one, then the higher one:-  --   boundingLines :: Num r => Slab o a r -> (Line 2 r :+ a, Line 2 r :+ a)-+  -- | Test if a point lies inside a slab.   inSlab :: Ord r => Point 2 r -> Slab o a r -> Bool   instance HasBoundingLines Horizontal where   boundingLines (Slab i) = (i^.start, i^.end)&both.core %~ horizontalLine -  p `inSlab` (Slab i) = (p^.yCoord) `inInterval` i+  p `inSlab` (Slab i) = (p^.yCoord) `intersectsInterval` i   instance HasBoundingLines Vertical where   boundingLines (Slab i) = (i^.start, i^.end)&both.core %~ verticalLine -  p `inSlab` (Slab i) = (p^.xCoord) `inInterval` i+  p `inSlab` (Slab i) = (p^.xCoord) `intersectsInterval` i   type instance IntersectionOf (Line 2 r) (Slab o a r) =

@@ -15,15 +15,16 @@   , subRange   , fixEndPoints   , dropExtra-  , _unBounded-  , toUnbounded-  , fromUnbounded   , onSubLine   , onSubLineUB   , onSubLine2   , onSubLine2UB+  , reorient   , getEndPointsUnBounded   , fromLine+  , _unBounded+  , toUnbounded+  , fromUnbounded   ) where  import           Control.Lens@@ -57,6 +58,7 @@ subRange :: Lens (SubLine d p1 s1 r) (SubLine d p2 s2 r) (Interval p1 s1) (Interval p2 s2) subRange = lens _subRange (SubLine . _line) + type instance Dimension (SubLine d p s r) = d  @@ -71,6 +73,7 @@          => Arbitrary (SubLine d p s r) where   arbitrary = SubLine <$> arbitrary <*> arbitrary + -- | Annotate the subRange with the actual ending points fixEndPoints    :: (Num r, Arity d) => SubLine d p r r -> SubLine d (Point d r :+ p) r r fixEndPoints sl = sl&subRange %~ f@@ -84,17 +87,8 @@ dropExtra :: SubLine d p s r -> SubLine d () s r dropExtra = over subRange (first (const ())) --- | Prism for downcasting an unbounded subline to a subline.-_unBounded :: Prism' (SubLine d p (UnBounded r) r) (SubLine d p r r)-_unBounded = prism' toUnbounded fromUnbounded --- | Transform into an subline with a potentially unbounded interval-toUnbounded :: SubLine d p r r -> SubLine d p (UnBounded r) r-toUnbounded = over subRange (fmap Val) --- | Try to make a potentially unbounded subline into a bounded one.-fromUnbounded               :: SubLine d p (UnBounded r) r -> Maybe (SubLine d p r r)-fromUnbounded (SubLine l i) = SubLine l <$> mapM unBoundedToMaybe i  -- | given point p, and a Subline l r such that p lies on line l, test if it -- lies on the subline, i.e. in the interval r@@ -102,24 +96,13 @@                           => Point d r -> SubLine d p r r -> Bool onSubLine p (SubLine l r) = case toOffset p l of                               Nothing -> False-                              Just x  -> x `inInterval` r+                              Just x  -> x `intersectsInterval` r --- | given point p, and a Subline l r such that p lies on line l, test if it--- lies on the subline, i.e. in the interval r-onSubLineUB                   :: (Ord r, Fractional r)-                              => Point 2 r -> SubLine 2 p (UnBounded r) r -> Bool-p `onSubLineUB` (SubLine l r) =-  p `onLine2` l &&-  Val (toOffset' p l) `inInterval` r -inSubLineIntervalUB                   :: (Ord r, Fractional r)-                              => Point 2 r -> SubLine 2 p (UnBounded r) r -> Bool-p `inSubLineIntervalUB` (SubLine l r) = Val (toOffset' p l) `inInterval` r- -- | given point p, and a Subline l r such that p lies on line l, test if it -- lies on the subline, i.e. in the interval r onSubLine2        :: (Ord r, Num r) => Point 2 r -> SubLine 2 p r r -> Bool-p `onSubLine2` sl = d `inInterval` r+p `onSubLine2` sl = d `intersectsInterval` r   where     -- get the endpoints (a,b) of the subline     SubLine _ (Interval s e) = fixEndPoints sl@@ -130,23 +113,60 @@     r = Interval (s&unEndPoint.core .~ 0) (e&unEndPoint.core .~ squaredEuclideanDist b a)  --- | given point p, and a Subline l r such that p lies on line l, test if it--- lies on the subline, i.e. in the interval r-onSubLine2UB        :: (Ord r, Fractional r)-                    => Point 2 r -> SubLine 2 p (UnBounded r) r -> Bool-p `onSubLine2UB` sl = p `onSubLineUB` sl+type instance IntersectionOf (SubLine 2 p s r) (SubLine 2 q s r) =+  [ NoIntersection, Point 2 r, SubLine 2 (Either p q) s r]  -type instance IntersectionOf (SubLine 2 p s r) (SubLine 2 q s r) = [ NoIntersection-                                                                   , Point 2 r-                                                                   , SubLine 2 p s r-                                                                   ]++ instance (Ord r, Fractional r) =>-         SubLine 2 p r r `HasIntersectionWith` SubLine 2 p r r+         SubLine 2 p r r `HasIntersectionWith` SubLine 2 q r r ++-- -- | Given two sublines that supposedly have the same line (but+-- -- possibly represented differently), test if they intersect.+-- intersectsSLRange :: SubLine 2 p r r -> SubLine 2 q r r -> Bool+-- intersectsSLRange = undefined+++-- -- | Given two sublines of the s ame line (but possibly represented differently)+-- -- align the first one to the second one.+-- --+-- -- pre: the+-- alignTo :: (Eq r, Num r, Arity d) => SubLine d p r r -> SubLine d q r r -> SubLine d p r r+-- sl `alignTo` (SubLine l@(Line p v) i2) = SubLine l i'+--   where+--     SubLine (Line q u) i = reorient sl v+++--     i' = undefined++++++++++-- | Given a subline with vector u, and a vector v that is parallel to+-- u (but possibly pointing in the exact opposite direction). Make the+-- subline point in direction v as well (but keep the magnitude of the+-- original vector.)+--+-- pre: the lines are parallel.+reorient :: (Eq r,Num r, Arity d) => SubLine d p r r -> Vector d r -> SubLine d p r r+reorient sl@(SubLine (Line p u) i) v+  | sameDirection u v = sl+  | otherwise         = SubLine (Line p ((-1) *^ u)) (flipInterval i)+++++ {- HLINT ignore "Redundant bracket" -} instance (Ord r, Fractional r) =>-         SubLine 2 p r r `IsIntersectableWith` SubLine 2 p r r where+         SubLine 2 p r r `IsIntersectableWith` SubLine 2 q r r where    nonEmptyIntersection = defaultNonEmptyIntersection @@ -167,11 +187,95 @@           $ s'&start.core .~ toOffset' (s'^.start.extra.core) l               &end.core   .~ toOffset' (s'^.end.extra.core)   l +++++-- testL :: SubLine 2 () (UnBounded Rational)+-- testL = SubLine (horizontalLine 0) (Interval (Closed (only 0)) (Open $ only 10))++-- horL :: SubLine 2 () (UnBounded Rational)+-- horL = fromLine $ horizontalLine 0+++-- test = (testL^.subRange) `intersect` (horL^.subRange)++-- toOffset (Point2 minInfinity minInfinity) (horizontalLine 0)+-- testzz = let f  = bimap (fmap Val) (const ())+--          in++-- testz :: SubLine 2 () Rational Rational+-- testz = SubLine (Line (Point2 0 0) (Vector2 10 0))+--                 (Interval (Closed (0 % 1 :+ ())) (Closed (1 % 1 :+ ())))+++++--------------------------------------------------------------------------------+-- * Anything that deals with Unbounded intervals++-- | Create a SubLine that covers the original line from -infinity to +infinity.+fromLine   :: Arity d => Line d r -> SubLine d () (UnBounded r) r+fromLine l = SubLine l (ClosedInterval (ext MinInfinity) (ext MaxInfinity))+++-- | Prism for downcasting an unbounded subline to a subline.+_unBounded :: Prism' (SubLine d p (UnBounded r) r) (SubLine d p r r)+_unBounded = prism' toUnbounded fromUnbounded++-- | Transform into an subline with a potentially unbounded interval+toUnbounded :: SubLine d p r r -> SubLine d p (UnBounded r) r+toUnbounded = over subRange (fmap Val)++-- | Try to make a potentially unbounded subline into a bounded one.+fromUnbounded               :: SubLine d p (UnBounded r) r -> Maybe (SubLine d p r r)+fromUnbounded (SubLine l i) = SubLine l <$> mapM unBoundedToMaybe i+++-- | Get the endpoints of an unbounded interval+getEndPointsUnBounded    :: (Num r, Arity d) => SubLine d p (UnBounded r) r+                         -> Interval p (UnBounded (Point d r))+getEndPointsUnBounded sl = second (fmap f) $ sl^.subRange+  where+    f = flip pointAt (sl^.line)++++++-- | given point p, and a Subline l r such that p lies on line l, test if it+-- lies on the subline, i.e. in the interval r+onSubLineUB                   :: (Ord r, Fractional r)+                              => Point 2 r -> SubLine 2 p (UnBounded r) r -> Bool+p `onSubLineUB` (SubLine l r) =+  p `onLine2` l &&+  Val (toOffset' p l) `intersectsInterval` r++inSubLineIntervalUB                   :: (Ord r, Fractional r)+                              => Point 2 r -> SubLine 2 p (UnBounded r) r -> Bool+p `inSubLineIntervalUB` (SubLine l r) = Val (toOffset' p l) `intersectsInterval` r++++-- | given point p, and a Subline l r such that p lies on line l, test if it+-- lies on the subline, i.e. in the interval r+onSubLine2UB        :: (Ord r, Fractional r)+                    => Point 2 r -> SubLine 2 p (UnBounded r) r -> Bool+p `onSubLine2UB` sl = p `onSubLineUB` sl++++++++--------+ instance (Ord r, Fractional r) =>-         SubLine 2 p (UnBounded r) r `HasIntersectionWith` SubLine 2 p (UnBounded r) r+         SubLine 2 p (UnBounded r) r `HasIntersectionWith` SubLine 2 q (UnBounded r) r  instance (Ord r, Fractional r) =>-         SubLine 2 p (UnBounded r) r `IsIntersectableWith` SubLine 2 p (UnBounded r) r where+         SubLine 2 p (UnBounded r) r `IsIntersectableWith` SubLine 2 q (UnBounded r) r where   nonEmptyIntersection = defaultNonEmptyIntersection    sl@(SubLine l r) `intersect` sm@(SubLine m _) = match (l `intersect` m) $@@ -190,32 +294,3 @@       s'  = getEndPointsUnBounded sm       s'' = second (fmap f) s'       f = flip toOffset' l---- | Get the endpoints of an unbounded interval-getEndPointsUnBounded    :: (Num r, Arity d) => SubLine d p (UnBounded r) r-                         -> Interval p (UnBounded (Point d r))-getEndPointsUnBounded sl = second (fmap f) $ sl^.subRange-  where-    f = flip pointAt (sl^.line)---- | Create a SubLine that covers the original line from -infinity to +infinity.-fromLine   :: Arity d => Line d r -> SubLine d () (UnBounded r) r-fromLine l = SubLine l (ClosedInterval (ext MinInfinity) (ext MaxInfinity))----- testL :: SubLine 2 () (UnBounded Rational)--- testL = SubLine (horizontalLine 0) (Interval (Closed (only 0)) (Open $ only 10))---- horL :: SubLine 2 () (UnBounded Rational)--- horL = fromLine $ horizontalLine 0----- test = (testL^.subRange) `intersect` (horL^.subRange)---- toOffset (Point2 minInfinity minInfinity) (horizontalLine 0)--- testzz = let f  = bimap (fmap Val) (const ())---          in---- testz :: SubLine 2 () Rational Rational--- testz = SubLine (Line (Point2 0 0) (Vector2 10 0))---                 (Interval (Closed (0 % 1 :+ ())) (Closed (1 % 1 :+ ())))

@@ -28,8 +28,8 @@  import           Control.Lens import           Data.Ext-import           Data.Geometry.Box (Rectangle, IsBoxable)-import qualified Data.Geometry.Box as Box+import           Data.Geometry.Box.Internal (Rectangle, IsBoxable)+import qualified Data.Geometry.Box.Internal as Box import           Data.Geometry.Properties import           Data.Geometry.Point import           Data.Geometry.Transformation.Internal

@@ -15,7 +15,6 @@ import           Data.Geometry.Properties import           Data.Geometry.Vector import qualified Data.Geometry.Vector as V-import           Data.Proxy import           GHC.TypeLits  {- $setup@@ -51,6 +50,12 @@ identity :: (Num r, Arity (d + 1)) => Transformation d r identity = Transformation identityMatrix +instance (Num r, Arity (d+1)) => Semigroup (Transformation d r) where+  (<>) = (|.|)+instance (Num r, Arity (d+1)) => Monoid (Transformation d r) where+  mempty = identity++ -- if it exists?  -- | Compute the inverse transformation@@ -168,7 +173,7 @@  transRow     :: forall n r. (Arity n, Arity (n + 1), Num r)              => Int -> r -> Vector (n + 1) r-transRow i x = set (V.element (Proxy :: Proxy n)) x $ mkRow i 1+transRow i x = set (V.element @n) x $ mkRow i 1  -------------------------------------------------------------------------------- -- * 3D Rotations

@@ -105,7 +105,7 @@  -- | Get the inscribed disk. Returns Nothing if the triangle is degenerate, -- i.e. if the points are colinear.-inscribedDisk                  :: (Eq r, Fractional r)+inscribedDisk                  :: (Ord r, Fractional r)                                => Triangle 2 p r -> Maybe (Disk () r) inscribedDisk (Triangle p q r) = disk (p^.core) (q^.core) (r^.core) 

@@ -174,7 +174,7 @@ -- >>> Vector2 1 2 & xComponent .~ 10 -- Vector2 10 2 xComponent :: (1 <= d, Arity d) => Lens' (Vector d r) r-xComponent = element (C :: C 0)+xComponent = element @0 {-# INLINABLE xComponent #-}  -- | Shorthand to access the second component@@ -184,7 +184,7 @@ -- >>> Vector2 1 2 & yComponent .~ 10 -- Vector2 1 10 yComponent :: (2 <= d, Arity d) => Lens' (Vector d r) r-yComponent = element (C :: C 1)+yComponent = element @1 {-# INLINABLE yComponent #-}  -- | Shorthand to access the third component@@ -194,5 +194,5 @@ -- >>> Vector3 1 2 3 & zComponent .~ 10 -- Vector3 1 2 10 zComponent :: (3 <= d, Arity d) => Lens' (Vector d r) r-zComponent = element (C :: C 2)+zComponent = element @2 {-# INLINABLE zComponent #-}

@@ -1,5 +1,6 @@ {-# LANGUAGE ScopedTypeVariables  #-} {-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE AllowAmbiguousTypes #-} -------------------------------------------------------------------------------- -- | -- Module      :  Data.Geometry.Vector.VectorFamily@@ -15,29 +16,30 @@ module Data.Geometry.Vector.VectorFamily where  import           Control.DeepSeq-import           Control.Lens                           hiding (element)+import           Control.Lens hiding (element) import           Control.Monad import           Data.Aeson-import qualified Data.Foldable                          as F+import qualified Data.Foldable as F import           Data.Functor.Classes import           Data.Geometry.Vector.VectorFamilyPeano (ImplicitArity, VectorFamily (..),                                                          VectorFamilyF) import qualified Data.Geometry.Vector.VectorFamilyPeano as Fam-import           Data.Geometry.Vector.VectorFixed       (C (..))+import           Data.Geometry.Vector.VectorFixed (C (..)) import           Data.Hashable+import           Data.Kind import           Data.List-import qualified Data.List                              as L+import qualified Data.List as L import           Data.Proxy-import qualified Data.Vector.Fixed                      as V-import           Data.Vector.Fixed.Cont                 (Peano)+import qualified Data.Vector.Fixed as V+import           Data.Vector.Fixed.Cont (Peano) import           GHC.TypeLits-import           Linear.Affine                          (Affine (..))+import           Linear.Affine (Affine (..)) import           Linear.Metric-import qualified Linear.V2                              as L2-import qualified Linear.V3                              as L3-import qualified Linear.V4                              as L4+import qualified Linear.V2 as L2+import qualified Linear.V3 as L3+import qualified Linear.V4 as L4 import           Linear.Vector-import           Text.Read                              (Read (..), readListPrecDefault)+import           Text.Read (Read (..), readListPrecDefault)  -------------------------------------------------------------------------------- -- * d dimensional Vectors@@ -46,7 +48,7 @@ -- | Datatype representing d dimensional vectors. The default implementation is -- based n VectorFixed. However, for small vectors we automatically select a -- more efficient representation.-newtype Vector (d :: Nat) (r :: *) = MKVector { _unV :: VectorFamily (Peano d) r }+newtype Vector (d :: Nat) (r :: Type) = MKVector { _unV :: VectorFamily (Peano d) r }  type instance V.Dim   (Vector d)   = Fam.FromPeano (Peano d) -- the above definition is a bit convoluted, but it allows us to make Vector an instance of@@ -196,17 +198,22 @@  -- | \( O(1) \) First element. Since arity is at least 1, this function is total. head   :: (Arity d, 1 <= d) => Vector d r -> r-head = view $ element (C :: C 0)+head = view $ element @0  -------------------------------------------------------------------------------- -- * Indexing vectors  -- | Lens into the i th element-element   :: forall proxy i d r. (Arity d, KnownNat i, (i + 1) <= d)-          => proxy i -> Lens' (Vector d r) r-element _ = singular . element' . fromInteger $ natVal (C :: C i)+element :: forall i d r. (Arity d, KnownNat i, (i + 1) <= d)+        => Lens' (Vector d r) r+element = elementProxy (C @i) {-# INLINE element #-} +-- | Lens into the i th element+elementProxy   :: forall proxy i d r. (Arity d, KnownNat i, (i + 1) <= d)+               => proxy i -> Lens' (Vector d r) r+elementProxy _ = singular $ element' $ fromInteger . natVal $ C @i+{-# INLINE elementProxy #-}  -- | Similar to 'element' above. Except that we don't have a static guarantee -- that the index is in bounds. Hence, we can only return a Traversal@@ -235,7 +242,7 @@  -- | \( O(1) \) Last element. Since the vector is non-empty, runtime bounds checks are bypassed. last :: forall d r. (KnownNat d, Arity (d + 1)) => Vector (d + 1) r -> r-last = view $ element (C :: C d)+last = view $ element @d  -- | Get a prefix of i elements of a vector prefix :: forall i d r. (Arity d, Arity i, i <= d)

@@ -19,6 +19,7 @@ import           Control.DeepSeq import           Control.Lens hiding (element) import           Data.Aeson (FromJSON(..),ToJSON(..))+import           Data.Kind -- import           Data.Aeson (ToJSON(..),FromJSON(..)) import qualified Data.Foldable as F import qualified Data.Geometry.Vector.VectorFixed as FV@@ -67,11 +68,11 @@ -- | Datatype representing d dimensional vectors. The default implementation is -- based n VectorFixed. However, for small vectors we automatically select a -- more efficient representation.-newtype VectorFamily (d :: PeanoNum) (r :: *) =+newtype VectorFamily (d :: PeanoNum) (r :: Type) =   VectorFamily { _unVF :: VectorFamilyF d r }  -- | Mapping between the implementation type, and the actual implementation.-type family VectorFamilyF (d :: PeanoNum) :: * -> * where+type family VectorFamilyF (d :: PeanoNum) :: Type -> Type where   VectorFamilyF Z        = Const ()   VectorFamilyF One      = Identity   VectorFamilyF Two      = L2.V2@@ -215,7 +216,7 @@   {-# INLINE rnf #-}  -instance (ImplicitPeano d, Hashable r) => Hashable (VectorFamily d r) where+instance (ImplicitArity d, Hashable r) => Hashable (VectorFamily d r) where   hashWithSalt = case (implicitPeano :: SingPeano d) of                    SZ                         -> hashWithSalt                    (SS SZ)                    -> hashWithSalt

@@ -13,10 +13,11 @@ import           Control.Lens hiding (element) import           Data.Aeson import qualified Data.Foldable as F+import           Data.Functor.Classes+import           Data.Kind import           Data.Proxy-import qualified Data.Vector.Fixed as V import           Data.Vector.Fixed (Arity)-import           Data.Functor.Classes+import qualified Data.Vector.Fixed as V import           Data.Vector.Fixed.Boxed import           GHC.Generics (Generic) import           GHC.TypeLits@@ -36,8 +37,8 @@  -- | Datatype representing d dimensional vectors. Our implementation wraps the -- implementation provided by fixed-vector.-newtype Vector (d :: Nat)  (r :: *) = Vector { _unV :: Vec d r }-                                    deriving (Generic)+newtype Vector (d :: Nat)  (r :: Type) = Vector { _unV :: Vec d r }+                                       deriving (Generic)  unV :: Lens' (Vector d r) (Vec d r) unV = lens _unV (const Vector)

@@ -16,13 +16,13 @@   , segmentAbove, segmentAboveOrOn   , findSlab   , lookupAbove, lookupAboveOrOn, searchInSlab-  , ordAt, yCoordAt   ) where  import           Algorithms.BinarySearch (binarySearchIn) import           Control.Lens hiding (contains, below) import           Data.Ext import           Data.Foldable (toList)+import           Data.Function (on) import           Data.Geometry.Line import           Data.Geometry.LineSegment import           Data.Geometry.Point@@ -30,16 +30,15 @@ import           Data.List.NonEmpty (NonEmpty(..)) import qualified Data.List.NonEmpty as NonEmpty import           Data.Maybe (mapMaybe)-import           Data.Ord (comparing) import           Data.Semigroup.Foldable import qualified Data.Set as SS -- status struct import qualified Data.Set.Util as SS import qualified Data.Vector as V  -import           Data.RealNumber.Rational+-- import           Data.RealNumber.Rational -type R = RealNumber 5+-- type R = RealNumber 5 --------------------------------------------------------------------------------  -- | The vertical ray shooting data structure@@ -138,10 +137,15 @@                       -> r :+ StatusStructure p e r handle ss ( l :+ acts           , r :+ _)   = let mid               = (l+r)/2-                            runActionAt x act = interpret act (ordAt x)+                            runActionAt x act = interpret act (ordAtX' x)                         in r :+ foldr (runActionAt mid) ss (orderActs acts)                            -- run deletions first +-- | order by x coord+ordAtX'   :: (Ord r, Fractional r)+          => r -> LineSegment 2 p r :+ a -> LineSegment 2 p r :+ a -> Ordering+ordAtX' x = ordAtX x `on` view core+ -- | orders the actions to put insertions first and then all deletions orderActs      :: NonEmpty (Action a) -> NonEmpty (Action a) orderActs acts = let (dels,ins) = NonEmpty.partition (\case@@ -207,24 +211,3 @@   ------------------------------------------------------------------------------------type Compare a = a -> a -> Ordering---- | Compare based on the y-coordinate of the intersection with the horizontal--- line through y-ordAt   :: (Fractional r, Ord r) => r -> Compare (LineSegment 2 p r :+ e)-ordAt x = comparing (yCoordAt x)----- | Given an x-coordinate and a line segment that intersects the vertical line--- through x, compute the y-coordinate of this intersection point.------ note that we will pretend that the line segment is closed, even if it is not-yCoordAt :: (Fractional r, Ord r) => r -> LineSegment 2 p r :+ e -> r-yCoordAt x (LineSegment' (Point2 px py :+ _) (Point2 qx qy :+ _) :+ _)-    | px == qx  = py `max` qy -- s is vertical, since by the precondition it-                              -- intersects we return the y-coord of the topmost-                              -- endpoint.-    | otherwise = py + alpha * (qy - py)-  where-    alpha = (x - px) / (qx - px)

@@ -70,6 +70,7 @@ import           Data.Geometry.Line (cmpSlope, supportingLine) import           Data.Geometry.LineSegment hiding (endPoints) import           Data.Geometry.Point+import           Data.Geometry.Vector import           Data.Geometry.Polygon import           Data.Geometry.Properties import qualified Data.List.NonEmpty as NonEmpty@@ -88,7 +89,6 @@ --------------------------------------------------------------------------------  -- $setup--- >>> import Data.Proxy -- >>> import Data.PlaneGraph.AdjRep(Gr(Gr),Face(Face),Vtx(Vtx)) -- >>> import Data.PlaneGraph.IO(fromAdjRep) -- >>> import Data.PlanarGraph.Dart(Dart(..),Arc(..))@@ -114,7 +114,7 @@ --                , Face (0,1) "A" --                , Face (1,0) "B" --                ]---     smallG = fromAdjRep (Proxy :: Proxy ()) small+--     smallG = fromAdjRep @() small -- :} -- --@@ -131,7 +131,7 @@ -- >>> data MyWorld -- >>> :{ -- let myPlaneGraph :: PlaneGraph MyWorld Int () String (RealNumber 5)---     myPlaneGraph = fromAdjRep (Proxy @MyWorld) myPlaneGraphAdjrep+--     myPlaneGraph = fromAdjRep @MyWorld myPlaneGraphAdjrep --     myPlaneGraphAdjrep :: Gr (Vtx Int () (RealNumber 5)) (Face String) --     myPlaneGraphAdjrep = Gr [ vtx 0 (Point2 0   0   ) [e 9, e 5, e 1, e 2] --                             , vtx 1 (Point2 4   4   ) [e 0, e 5, e 12]@@ -215,23 +215,23 @@ -- -- pre: the input polygon is given in counterclockwise order -- running time: \(O(n)\).-fromSimplePolygon                            :: proxy s-                                             -> SimplePolygon p r+fromSimplePolygon                            :: forall s p r f.+                                                SimplePolygon p r                                              -> f -- ^ data inside                                              -> f -- ^ data outside the polygon                                              -> PlaneGraph s p () f r-fromSimplePolygon p poly iD oD = PlaneGraph g'+fromSimplePolygon poly iD oD = PlaneGraph g'   where     vs     = poly ^. outerBoundaryVector-    g      = fromVertices p vs+    g      = fromVertices vs     fData' = V.fromList [iD, oD]     g'     = g & PG.faceData .~ fData'  -- | Constructs a planar from the given vertices-fromVertices      :: proxy s-                  -> CircularVector (Point 2 r :+ p)-                  -> PlanarGraph s Primal (VertexData r p) () ()-fromVertices _ vs = g&PG.vertexData .~ vData'+fromVertices    :: forall s r p.+                   CircularVector (Point 2 r :+ p)+                -> PlanarGraph s Primal (VertexData r p) () ()+fromVertices vs = g&PG.vertexData .~ vData'   where     n = length vs     g = planarGraph [ [ (Dart (Arc i)               Positive, ())@@ -245,11 +245,10 @@ -- pre: The segments form a single connected component -- -- running time: \(O(n\log n)\)-fromConnectedSegments      :: (Foldable f, Ord r, Num r)-                           => proxy s-                           -> f (LineSegment 2 p r :+ e)-                           -> PlaneGraph s (NonEmpty.NonEmpty p) e () r-fromConnectedSegments _ ss = PlaneGraph $ planarGraph dts & PG.vertexData .~ vxData+fromConnectedSegments    :: forall s p r e f. (Foldable f, Ord r, Num r)+                         => f (LineSegment 2 p r :+ e)+                         -> PlaneGraph s (NonEmpty.NonEmpty p) e () r+fromConnectedSegments ss = PlaneGraph $ planarGraph dts & PG.vertexData .~ vxData   where     pts         = M.fromListWith (<>) . concatMap f . zipWith g [0..] . F.toList $ ss     f (s :+ e)  = [ ( s^.start.core@@ -381,7 +380,7 @@ -- | face Ids of all internal faces in the plane graph -- -- running time: \(O(n)\)-internalFaces'   :: (Ord r, Fractional r) => PlaneGraph s v e f r  -> V.Vector (FaceId' s)+internalFaces'   :: (Ord r, Num r) => PlaneGraph s v e f r  -> V.Vector (FaceId' s) internalFaces' g = let i = outerFaceId g in V.filter (/= i) $ faces' g  -- | All faces with their face data.@@ -403,14 +402,14 @@  -- | Reports the outerface and all internal faces separately. -- running time: \(O(n)\)-faces''   :: (Ord r, Fractional r)+faces''   :: (Ord r, Num r)           => PlaneGraph s v e f r -> ((FaceId' s, f), V.Vector (FaceId' s, f)) faces'' g = let i = outerFaceId g             in ((i,g^.dataOf i), V.filter (\(j,_) -> i /= j) $ faces g)  -- | Reports all internal faces. -- running time: \(O(n)\)-internalFaces :: (Ord r, Fractional r)+internalFaces :: (Ord r, Num r)               => PlaneGraph s v e f r -> V.Vector (FaceId' s, f) internalFaces = snd . faces'' @@ -469,8 +468,8 @@   --- | All edges incident to vertex v in incoming direction--- (i.e. pointing into v) in counterclockwise order around v.+-- | All edges incident to vertex v in outgoing direction+-- (i.e. pointing out of v) in counterclockwise order around v. -- -- running time: \(O(k)\), where \(k) is the total number of incident edges of v --@@ -763,7 +762,7 @@ -- -- running time: \(O(n)\) ---outerFaceId    :: (Ord r, Fractional r) => PlaneGraph s v e f r -> FaceId' s+outerFaceId    :: (Ord r, Num r) => PlaneGraph s v e f r -> FaceId' s outerFaceId ps = leftFace (outerFaceDart ps) ps  @@ -772,19 +771,24 @@ -- -- running time: \(O(n)\) ---outerFaceDart    :: (Ord r, Fractional r) => PlaneGraph s v e f r -> Dart s-outerFaceDart ps = d+outerFaceDart    :: (Ord r, Num r) => PlaneGraph s v e f r -> Dart s+outerFaceDart pg = d   where-    (v,_)  = V.minimumBy (comparing (^._2.location)) . vertices $ ps+    (v,_)  = V.minimumBy (comparing (^._2.location)) . vertices $ pg            -- compare lexicographically; i.e. if same x-coord prefer the one with the            -- smallest y-coord-    d :+ _ = V.maximumBy (cmpSlope `on` (^.extra))-           .  fmap (\d' -> d' :+ edgeSegment d' ps ^. core.to supportingLine)-           $ incidentEdges v ps++    (_ :+ d) = V.minimumBy (cwCmpAroundWith' (Vector2 (-1) 0) (pg^.locationOf v :+ ()))+             . fmap (\d' -> let u = headOf d' pg in (pg^.locationOf u) :+ d')+             $ outgoingEdges v pg     -- based on the approach sketched at https://cstheory.stackexchange.com/questions/27586/finding-outer-face-in-plane-graph-embedded-planar-graph     -- basically: find the leftmost vertex, find the incident edge with the largest slope     -- and take the face left of that edge. This is the outerface.     -- note that this requires that the edges are straight line segments+    --+    -- note that rather computing slopes we just ask for the first+    -- vertec cw vertex around v. First with respect to some direction+    -- pointing towards the left.   --------------------------------------------------------------------------------@@ -913,7 +917,7 @@  -- | lists all internal faces of the plane graph with their -- boundaries.-internalFacePolygons    :: (Ord r, Fractional r)+internalFacePolygons    :: (Ord r, Num r)                         => PlaneGraph s v e f r ->  V.Vector (FaceId' s, SimplePolygon v r :+ f) internalFacePolygons pg = facePolygons' (outerFaceId pg) pg 

@@ -22,7 +22,6 @@ import qualified Data.PlanarGraph.IO as PGIO import           Data.PlaneGraph.Core import           Data.PlaneGraph.AdjRep-import           Data.Proxy import qualified Data.Vector as V import qualified Data.Vector.Mutable as MV import           Data.Yaml (ParseException)@@ -60,7 +59,7 @@ --                , Face (0,1) "A" --                , Face (1,0) "B" --                ]---     smallG = fromAdjRep (Proxy :: Proxy ()) small+--     smallG = fromAdjRep @() small -- :} -- --@@ -91,7 +90,7 @@  instance (FromJSON v, FromJSON e, FromJSON f, FromJSON r)          => FromJSON (PlaneGraph s v e f r) where-  parseJSON v = fromAdjRep (Proxy :: Proxy s) <$> parseJSON v+  parseJSON v = fromAdjRep @s <$> parseJSON v  -------------------------------------------------------------------------------- @@ -110,9 +109,9 @@ -- 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)+fromAdjRep :: forall s v e f r. Gr (Vtx v e r) (Face f) -> PlaneGraph s v e f r+fromAdjRep = PlaneGraph . PGIO.fromAdjRep+           . first (\(Vtx v p aj x) -> PGA.Vtx v aj $ VertexData p x)  -------------------------------------------------------------------------------- @@ -167,7 +166,7 @@  -- ![myGraph](docs/Data/PlaneGraph/planegraph.png) myPlaneGraph :: PlaneGraph MyWorld Int () String (RealNumber 5)-myPlaneGraph = fromAdjRep (Proxy @MyWorld) myPlaneGraphAdjrep+myPlaneGraph = fromAdjRep @MyWorld myPlaneGraphAdjrep  myPlaneGraphAdjrep :: Gr (Vtx Int () (RealNumber 5)) (Face String) myPlaneGraphAdjrep = Gr [ vtx 0 (Point2 0   0   ) [e 9, e 5, e 1, e 2]