hgeometry-0.14: benchmark/Algorithms/Geometry/LineSegmentIntersection/BentleyOttmannNoExt.hs
{-# 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)