diff --git a/benchmark/Algorithms/Geometry/ConvexHull/GrahamFamPeano.hs b/benchmark/Algorithms/Geometry/ConvexHull/GrahamFamPeano.hs deleted file mode 100644 --- a/benchmark/Algorithms/Geometry/ConvexHull/GrahamFamPeano.hs +++ /dev/null @@ -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 diff --git a/benchmark/Algorithms/Geometry/LineSegmentIntersection/Bench.hs b/benchmark/Algorithms/Geometry/LineSegmentIntersection/Bench.hs --- a/benchmark/Algorithms/Geometry/LineSegmentIntersection/Bench.hs +++ b/benchmark/Algorithms/Geometry/LineSegmentIntersection/Bench.hs @@ -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) diff --git a/benchmark/Algorithms/Geometry/LineSegmentIntersection/BentleyOttmannNoExt.hs b/benchmark/Algorithms/Geometry/LineSegmentIntersection/BentleyOttmannNoExt.hs new file mode 100644 --- /dev/null +++ b/benchmark/Algorithms/Geometry/LineSegmentIntersection/BentleyOttmannNoExt.hs @@ -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) diff --git a/benchmark/Algorithms/Geometry/LineSegmentIntersection/BentleyOttmannOld.hs b/benchmark/Algorithms/Geometry/LineSegmentIntersection/BentleyOttmannOld.hs --- a/benchmark/Algorithms/Geometry/LineSegmentIntersection/BentleyOttmannOld.hs +++ b/benchmark/Algorithms/Geometry/LineSegmentIntersection/BentleyOttmannOld.hs @@ -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 diff --git a/benchmark/Algorithms/Geometry/LineSegmentIntersection/TypesNoExt.hs b/benchmark/Algorithms/Geometry/LineSegmentIntersection/TypesNoExt.hs new file mode 100644 --- /dev/null +++ b/benchmark/Algorithms/Geometry/LineSegmentIntersection/TypesNoExt.hs @@ -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 diff --git a/benchmark/Algorithms/Geometry/PolygonTriangulation/MakeMonotoneOld.hs b/benchmark/Algorithms/Geometry/PolygonTriangulation/MakeMonotoneOld.hs --- a/benchmark/Algorithms/Geometry/PolygonTriangulation/MakeMonotoneOld.hs +++ b/benchmark/Algorithms/Geometry/PolygonTriangulation/MakeMonotoneOld.hs @@ -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) diff --git a/changelog b/changelog --- a/changelog +++ b/changelog @@ -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 diff --git a/changelog.org b/changelog.org --- a/changelog.org +++ b/changelog.org @@ -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 diff --git a/doctests.hs b/doctests.hs --- a/doctests.hs +++ b/doctests.hs @@ -32,6 +32,8 @@ , "DeriveGeneric" , "FlexibleInstances" , "FlexibleContexts" + , "DerivingStrategies" + , "DerivingVia" ] files :: [String] diff --git a/hgeometry.cabal b/hgeometry.cabal --- a/hgeometry.cabal +++ b/hgeometry.cabal @@ -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 diff --git a/src/Algorithms/Geometry/DelaunayTriangulation/Types.hs b/src/Algorithms/Geometry/DelaunayTriangulation/Types.hs --- a/src/Algorithms/Geometry/DelaunayTriangulation/Types.hs +++ b/src/Algorithms/Geometry/DelaunayTriangulation/Types.hs @@ -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 diff --git a/src/Algorithms/Geometry/EuclideanMST.hs b/src/Algorithms/Geometry/EuclideanMST.hs --- a/src/Algorithms/Geometry/EuclideanMST.hs +++ b/src/Algorithms/Geometry/EuclideanMST.hs @@ -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 diff --git a/src/Algorithms/Geometry/LineSegmentIntersection.hs b/src/Algorithms/Geometry/LineSegmentIntersection.hs --- a/src/Algorithms/Geometry/LineSegmentIntersection.hs +++ b/src/Algorithms/Geometry/LineSegmentIntersection.hs @@ -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 diff --git a/src/Algorithms/Geometry/LineSegmentIntersection/BentleyOttmann.hs b/src/Algorithms/Geometry/LineSegmentIntersection/BentleyOttmann.hs --- a/src/Algorithms/Geometry/LineSegmentIntersection/BentleyOttmann.hs +++ b/src/Algorithms/Geometry/LineSegmentIntersection/BentleyOttmann.hs @@ -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] diff --git a/src/Algorithms/Geometry/LineSegmentIntersection/BooleanSweep.hs b/src/Algorithms/Geometry/LineSegmentIntersection/BooleanSweep.hs --- a/src/Algorithms/Geometry/LineSegmentIntersection/BooleanSweep.hs +++ b/src/Algorithms/Geometry/LineSegmentIntersection/BooleanSweep.hs @@ -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. diff --git a/src/Algorithms/Geometry/LineSegmentIntersection/Naive.hs b/src/Algorithms/Geometry/LineSegmentIntersection/Naive.hs --- a/src/Algorithms/Geometry/LineSegmentIntersection/Naive.hs +++ b/src/Algorithms/Geometry/LineSegmentIntersection/Naive.hs @@ -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) diff --git a/src/Algorithms/Geometry/LineSegmentIntersection/Types.hs b/src/Algorithms/Geometry/LineSegmentIntersection/Types.hs --- a/src/Algorithms/Geometry/LineSegmentIntersection/Types.hs +++ b/src/Algorithms/Geometry/LineSegmentIntersection/Types.hs @@ -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 diff --git a/src/Algorithms/Geometry/LowerEnvelope/DualCH.hs b/src/Algorithms/Geometry/LowerEnvelope/DualCH.hs --- a/src/Algorithms/Geometry/LowerEnvelope/DualCH.hs +++ b/src/Algorithms/Geometry/LowerEnvelope/DualCH.hs @@ -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 diff --git a/src/Algorithms/Geometry/PolygonTriangulation/EarClip.hs b/src/Algorithms/Geometry/PolygonTriangulation/EarClip.hs --- a/src/Algorithms/Geometry/PolygonTriangulation/EarClip.hs +++ b/src/Algorithms/Geometry/PolygonTriangulation/EarClip.hs @@ -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: -- diff --git a/src/Algorithms/Geometry/PolygonTriangulation/MakeMonotone.hs b/src/Algorithms/Geometry/PolygonTriangulation/MakeMonotone.hs --- a/src/Algorithms/Geometry/PolygonTriangulation/MakeMonotone.hs +++ b/src/Algorithms/Geometry/PolygonTriangulation/MakeMonotone.hs @@ -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) diff --git a/src/Algorithms/Geometry/PolygonTriangulation/Triangulate.hs b/src/Algorithms/Geometry/PolygonTriangulation/Triangulate.hs --- a/src/Algorithms/Geometry/PolygonTriangulation/Triangulate.hs +++ b/src/Algorithms/Geometry/PolygonTriangulation/Triangulate.hs @@ -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) diff --git a/src/Algorithms/Geometry/PolygonTriangulation/TriangulateMonotone.hs b/src/Algorithms/Geometry/PolygonTriangulation/TriangulateMonotone.hs --- a/src/Algorithms/Geometry/PolygonTriangulation/TriangulateMonotone.hs +++ b/src/Algorithms/Geometry/PolygonTriangulation/TriangulateMonotone.hs @@ -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 diff --git a/src/Algorithms/Geometry/PolygonTriangulation/Types.hs b/src/Algorithms/Geometry/PolygonTriangulation/Types.hs --- a/src/Algorithms/Geometry/PolygonTriangulation/Types.hs +++ b/src/Algorithms/Geometry/PolygonTriangulation/Types.hs @@ -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 diff --git a/src/Algorithms/Geometry/RayShooting/Naive.hs b/src/Algorithms/Geometry/RayShooting/Naive.hs new file mode 100644 --- /dev/null +++ b/src/Algorithms/Geometry/RayShooting/Naive.hs @@ -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)) + + + +-------------------------------------------------------------------------------- diff --git a/src/Algorithms/Geometry/SSSP.hs b/src/Algorithms/Geometry/SSSP.hs --- a/src/Algorithms/Geometry/SSSP.hs +++ b/src/Algorithms/Geometry/SSSP.hs @@ -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) diff --git a/src/Algorithms/Geometry/VisibilityPolygon/Lee.hs b/src/Algorithms/Geometry/VisibilityPolygon/Lee.hs --- a/src/Algorithms/Geometry/VisibilityPolygon/Lee.hs +++ b/src/Algorithms/Geometry/VisibilityPolygon/Lee.hs @@ -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 diff --git a/src/Algorithms/Geometry/WSPD.hs b/src/Algorithms/Geometry/WSPD.hs --- a/src/Algorithms/Geometry/WSPD.hs +++ b/src/Algorithms/Geometry/WSPD.hs @@ -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..] diff --git a/src/Data/Geometry/Arrangement/Internal.hs b/src/Data/Geometry/Arrangement/Internal.hs --- a/src/Data/Geometry/Arrangement/Internal.hs +++ b/src/Data/Geometry/Arrangement/Internal.hs @@ -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 diff --git a/src/Data/Geometry/Ball.hs b/src/Data/Geometry/Ball.hs --- a/src/Data/Geometry/Ball.hs +++ b/src/Data/Geometry/Ball.hs @@ -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) diff --git a/src/Data/Geometry/BezierSpline.hs b/src/Data/Geometry/BezierSpline.hs --- a/src/Data/Geometry/BezierSpline.hs +++ b/src/Data/Geometry/BezierSpline.hs @@ -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 - - --------------------------------------------------------------------------------- - - - -------------------------------------------------------------------------------- diff --git a/src/Data/Geometry/Box.hs b/src/Data/Geometry/Box.hs --- a/src/Data/Geometry/Box.hs +++ b/src/Data/Geometry/Box.hs @@ -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 diff --git a/src/Data/Geometry/Box/Internal.hs b/src/Data/Geometry/Box/Internal.hs --- a/src/Data/Geometry/Box/Internal.hs +++ b/src/Data/Geometry/Box/Internal.hs @@ -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) diff --git a/src/Data/Geometry/Box/Sides.hs b/src/Data/Geometry/Box/Sides.hs --- a/src/Data/Geometry/Box/Sides.hs +++ b/src/Data/Geometry/Box/Sides.hs @@ -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 diff --git a/src/Data/Geometry/Duality.hs b/src/Data/Geometry/Duality.hs --- a/src/Data/Geometry/Duality.hs +++ b/src/Data/Geometry/Duality.hs @@ -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 diff --git a/src/Data/Geometry/HyperPlane.hs b/src/Data/Geometry/HyperPlane.hs --- a/src/Data/Geometry/HyperPlane.hs +++ b/src/Data/Geometry/HyperPlane.hs @@ -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 diff --git a/src/Data/Geometry/Interval.hs b/src/Data/Geometry/Interval.hs --- a/src/Data/Geometry/Interval.hs +++ b/src/Data/Geometry/Interval.hs @@ -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 diff --git a/src/Data/Geometry/KDTree.hs b/src/Data/Geometry/KDTree.hs --- a/src/Data/Geometry/KDTree.hs +++ b/src/Data/Geometry/KDTree.hs @@ -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 diff --git a/src/Data/Geometry/Line.hs b/src/Data/Geometry/Line.hs --- a/src/Data/Geometry/Line.hs +++ b/src/Data/Geometry/Line.hs @@ -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) diff --git a/src/Data/Geometry/Line/Internal.hs b/src/Data/Geometry/Line/Internal.hs --- a/src/Data/Geometry/Line/Internal.hs +++ b/src/Data/Geometry/Line/Internal.hs @@ -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 diff --git a/src/Data/Geometry/LineSegment.hs b/src/Data/Geometry/LineSegment.hs --- a/src/Data/Geometry/LineSegment.hs +++ b/src/Data/Geometry/LineSegment.hs @@ -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 diff --git a/src/Data/Geometry/LineSegment/Internal.hs b/src/Data/Geometry/LineSegment/Internal.hs --- a/src/Data/Geometry/LineSegment/Internal.hs +++ b/src/Data/Geometry/LineSegment/Internal.hs @@ -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) diff --git a/src/Data/Geometry/PlanarSubdivision.hs b/src/Data/Geometry/PlanarSubdivision.hs --- a/src/Data/Geometry/PlanarSubdivision.hs +++ b/src/Data/Geometry/PlanarSubdivision.hs @@ -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 diff --git a/src/Data/Geometry/PlanarSubdivision/Basic.hs b/src/Data/Geometry/PlanarSubdivision/Basic.hs --- a/src/Data/Geometry/PlanarSubdivision/Basic.hs +++ b/src/Data/Geometry/PlanarSubdivision/Basic.hs @@ -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 diff --git a/src/Data/Geometry/PlanarSubdivision/Dynamic.hs b/src/Data/Geometry/PlanarSubdivision/Dynamic.hs --- a/src/Data/Geometry/PlanarSubdivision/Dynamic.hs +++ b/src/Data/Geometry/PlanarSubdivision/Dynamic.hs @@ -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) diff --git a/src/Data/Geometry/PlanarSubdivision/Merge.hs b/src/Data/Geometry/PlanarSubdivision/Merge.hs --- a/src/Data/Geometry/PlanarSubdivision/Merge.hs +++ b/src/Data/Geometry/PlanarSubdivision/Merge.hs @@ -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] diff --git a/src/Data/Geometry/Point.hs b/src/Data/Geometry/Point.hs --- a/src/Data/Geometry/Point.hs +++ b/src/Data/Geometry/Point.hs @@ -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 diff --git a/src/Data/Geometry/Point/Class.hs b/src/Data/Geometry/Point/Class.hs --- a/src/Data/Geometry/Point/Class.hs +++ b/src/Data/Geometry/Point/Class.hs @@ -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 #-} diff --git a/src/Data/Geometry/Point/Internal.hs b/src/Data/Geometry/Point/Internal.hs --- a/src/Data/Geometry/Point/Internal.hs +++ b/src/Data/Geometry/Point/Internal.hs @@ -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 diff --git a/src/Data/Geometry/PointLocation/PersistentSweep.hs b/src/Data/Geometry/PointLocation/PersistentSweep.hs --- a/src/Data/Geometry/PointLocation/PersistentSweep.hs +++ b/src/Data/Geometry/PointLocation/PersistentSweep.hs @@ -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 -- diff --git a/src/Data/Geometry/PolyLine.hs b/src/Data/Geometry/PolyLine.hs --- a/src/Data/Geometry/PolyLine.hs +++ b/src/Data/Geometry/PolyLine.hs @@ -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. diff --git a/src/Data/Geometry/Polygon.hs b/src/Data/Geometry/Polygon.hs --- a/src/Data/Geometry/Polygon.hs +++ b/src/Data/Geometry/Polygon.hs @@ -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) diff --git a/src/Data/Geometry/Polygon/Core.hs b/src/Data/Geometry/Polygon/Core.hs --- a/src/Data/Geometry/Polygon/Core.hs +++ b/src/Data/Geometry/Polygon/Core.hs @@ -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." diff --git a/src/Data/Geometry/RangeTree.hs b/src/Data/Geometry/RangeTree.hs --- a/src/Data/Geometry/RangeTree.hs +++ b/src/Data/Geometry/RangeTree.hs @@ -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) diff --git a/src/Data/Geometry/RangeTree/Measure.hs b/src/Data/Geometry/RangeTree/Measure.hs --- a/src/Data/Geometry/RangeTree/Measure.hs +++ b/src/Data/Geometry/RangeTree/Measure.hs @@ -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 diff --git a/src/Data/Geometry/Slab.hs b/src/Data/Geometry/Slab.hs --- a/src/Data/Geometry/Slab.hs +++ b/src/Data/Geometry/Slab.hs @@ -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) = diff --git a/src/Data/Geometry/SubLine.hs b/src/Data/Geometry/SubLine.hs --- a/src/Data/Geometry/SubLine.hs +++ b/src/Data/Geometry/SubLine.hs @@ -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 :+ ()))) diff --git a/src/Data/Geometry/Transformation.hs b/src/Data/Geometry/Transformation.hs --- a/src/Data/Geometry/Transformation.hs +++ b/src/Data/Geometry/Transformation.hs @@ -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 diff --git a/src/Data/Geometry/Transformation/Internal.hs b/src/Data/Geometry/Transformation/Internal.hs --- a/src/Data/Geometry/Transformation/Internal.hs +++ b/src/Data/Geometry/Transformation/Internal.hs @@ -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 diff --git a/src/Data/Geometry/Triangle.hs b/src/Data/Geometry/Triangle.hs --- a/src/Data/Geometry/Triangle.hs +++ b/src/Data/Geometry/Triangle.hs @@ -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) diff --git a/src/Data/Geometry/Vector.hs b/src/Data/Geometry/Vector.hs --- a/src/Data/Geometry/Vector.hs +++ b/src/Data/Geometry/Vector.hs @@ -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 #-} diff --git a/src/Data/Geometry/Vector/VectorFamily.hs b/src/Data/Geometry/Vector/VectorFamily.hs --- a/src/Data/Geometry/Vector/VectorFamily.hs +++ b/src/Data/Geometry/Vector/VectorFamily.hs @@ -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) diff --git a/src/Data/Geometry/Vector/VectorFamilyPeano.hs b/src/Data/Geometry/Vector/VectorFamilyPeano.hs --- a/src/Data/Geometry/Vector/VectorFamilyPeano.hs +++ b/src/Data/Geometry/Vector/VectorFamilyPeano.hs @@ -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 diff --git a/src/Data/Geometry/Vector/VectorFixed.hs b/src/Data/Geometry/Vector/VectorFixed.hs --- a/src/Data/Geometry/Vector/VectorFixed.hs +++ b/src/Data/Geometry/Vector/VectorFixed.hs @@ -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) diff --git a/src/Data/Geometry/VerticalRayShooting/PersistentSweep.hs b/src/Data/Geometry/VerticalRayShooting/PersistentSweep.hs --- a/src/Data/Geometry/VerticalRayShooting/PersistentSweep.hs +++ b/src/Data/Geometry/VerticalRayShooting/PersistentSweep.hs @@ -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) diff --git a/src/Data/PlaneGraph/Core.hs b/src/Data/PlaneGraph/Core.hs --- a/src/Data/PlaneGraph/Core.hs +++ b/src/Data/PlaneGraph/Core.hs @@ -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 diff --git a/src/Data/PlaneGraph/IO.hs b/src/Data/PlaneGraph/IO.hs --- a/src/Data/PlaneGraph/IO.hs +++ b/src/Data/PlaneGraph/IO.hs @@ -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]