{-# LANGUAGE ScopedTypeVariables #-}
module Algorithms.Geometry.LineSegmentIntersection.BentleyOttmann where
import Algorithms.Geometry.LineSegmentIntersection.Types
import Control.Lens hiding (contains)
import Data.Ext
import qualified Data.Foldable as F
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 Data.OrdSeq (Compare)
import qualified Data.OrdSeq as SS
import qualified Data.Set as EQ
import Data.Vinyl
import Data.Vinyl.CoRec
intersections :: (Ord r, Fractional r)
=> [LineSegment 2 p r] -> Intersections p r
intersections ss = merge $ sweep pts mempty
where
pts = EQ.fromAscList . groupStarts . L.sort . concatMap asEventPts $ ss
interiorIntersections :: (Ord r, Fractional r)
=> [LineSegment 2 p r] -> Intersections p r
interiorIntersections = M.filter (not . isEndPointIntersection) . intersections
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)]
merge :: Ord r => [IntersectionPoint p r] -> Intersections p r
merge = foldr (\(IntersectionPoint p a) -> M.insertWith (<>) p a) M.empty
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
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 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
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
(Event p s) `compare` (Event q t) = case ordPoints p q of
EQ -> s `compare` t
x -> 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
startSegs :: Event p r -> [LineSegment 2 p r]
startSegs e = case eventType e of
Start ss -> NonEmpty.toList ss
_ -> []
ordAt :: (Fractional r, Ord r) => r -> Compare (LineSegment 2 p r)
ordAt y = comparing (xCoordAt y)
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
| otherwise = px + alpha * (qx - px)
where
alpha = (y - py) / (qy - py)
type EventQueue p r = EQ.Set (Event p r)
type StatusStructure p r = SS.OrdSeq (LineSegment 2 p r)
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
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'
starts' = filter (isClosedStart p) starts
toReport = case starts' ++ contains' of
(_:_:_) -> [IntersectionPoint p $ associated (starts' <> ends) contains]
_ -> []
ss' = before <> newSegs <> after
newSegs = toStatusStruct p $ starts ++ contains
eq' = foldr EQ.insert eq es
es | F.null newSegs = maybeToList $ app (findNewEvent p) sl sr
| otherwise = let s' = fst <$> SS.minView newSegs
s'' = fst <$> SS.maxView newSegs
in catMaybes [ app (findNewEvent p) sl s'
, app (findNewEvent p) s'' sr
]
sl = fst <$> SS.maxView before
sr = fst <$> SS.minView after
app f x y = do { x' <- x ; y' <- y ; f x' y'}
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 <> mid2, after)
where
(before, mid1, after') = SS.splitOn (xCoordAt $ p^.yCoord) (p^.xCoord) ss
(mid2, after) = SS.splitMonotonic (\s -> not $ p `onSegment` s) after'
toStatusStruct :: (Fractional r, Ord r)
=> Point 2 r -> [LineSegment 2 p r] -> StatusStructure p r
toStatusStruct p xs = ss <> hors
where
(hors',rest) = L.partition isHorizontal xs
ss = SS.fromListBy (ordAt $ maxY xs) rest
hors = SS.fromListBy (comparing rightEndpoint) hors'
isHorizontal s = s^.start.core.yCoord == s^.end.core.yCoord
maxY = maximum . filter (< p^.yCoord)
. concatMap (\s -> [s^.start.core.yCoord,s^.end.core.yCoord])
rightEndpoint :: Ord r => LineSegment 2 p r -> r
rightEndpoint s = (s^.start.core.xCoord) `max` (s^.end.core.xCoord)
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]
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 $ \NoIntersection -> Nothing)
:& (H $ \q -> if ordPoints q p == GT then Just (Event q Intersection)
else Nothing)
:& (H $ \_ -> Nothing)
:& RNil