{-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE ScopedTypeVariables #-} -------------------------------------------------------------------------------- -- | -- Module : Data.Geometry.Arrangement.Internal -- Copyright : (C) Frank Staals -- License : see the LICENSE file -- Maintainer : Frank Staals -- -- Data type for representing an Arrangement of lines in \(\mathbb{R}^2\). -- -------------------------------------------------------------------------------- module Data.Geometry.Arrangement.Internal where import Control.Lens import qualified Data.CircularSeq as CSeq import Data.Ext import qualified Data.Foldable as F import Data.Geometry.Boundary import Data.Geometry.Box import Data.Geometry.Line import Data.Geometry.LineSegment import Data.Geometry.PlanarSubdivision import Data.Geometry.Point import Data.Geometry.Properties import qualified Data.List as List import Data.Maybe import Data.Ord (Down(..)) import Data.Proxy import Data.Sequence.Util import qualified Data.Vector as V import Data.Vinyl.CoRec -------------------------------------------------------------------------------- type ArrangementBoundary s e r = V.Vector (Point 2 r, VertexId' s, Maybe (Line 2 r :+ e)) -- | Data type representing a two dimensional planar arrangement data Arrangement s l v e f r = Arrangement { _inputLines :: V.Vector (Line 2 r :+ l) , _subdivision :: PlanarSubdivision s v e f r , _boundedArea :: Rectangle () r , _unboundedIntersections :: ArrangementBoundary s l r } deriving (Show,Eq) -- unboundedIntersections also stores the corners of the box. They are not -- associated with any line makeLenses ''Arrangement type instance NumType (Arrangement s l v e f r) = r type instance Dimension (Arrangement s l v e f r) = 2 -------------------------------------------------------------------------------- -- | 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 -- | 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 -- | 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 = Arrangement (V.fromList ls) subdiv rect (link parts' subdiv) where subdiv = fromConnectedSegments px segs & rawVertexData.traverse.dataVal .~ () (segs,parts') = computeSegsAndParts rect ls computeSegsAndParts :: (Ord r, Fractional r) => Rectangle () r -> [Line 2 r :+ l] -> ( [LineSegment 2 () r :+ Maybe l] , [(Point 2 r, Maybe (Line 2 r :+ l))] ) computeSegsAndParts rect ls = ( segs <> boundarySegs, parts') where segs = map (&extra %~ Just) . concatMap (\(l,ls') -> perLine rect l ls') $ makePairs ls boundarySegs = map (:+ Nothing) . toSegments . dupFirst $ map fst parts' dupFirst = \case [] -> [] xs@(x:_) -> xs ++ [x] parts' = unBoundedParts rect ls perLine :: forall r l. (Ord r, Fractional r) => Rectangle () r -> Line 2 r :+ l -> [Line 2 r :+ l] -> [LineSegment 2 () r :+ l] perLine b m ls = map (:+ m^.extra) . toSegments . rmDuplicates . List.sort $ vs <> vs' where rmDuplicates = map head . List.group vs = mapMaybe (m `intersectionPoint`) ls vs' = maybe [] (\(p,q) -> [p,q]) . asA (Proxy :: Proxy (Point 2 r, Point 2 r)) $ (m^.core) `intersect` (Boundary b) intersectionPoint :: (Ord r, Fractional r) => Line 2 r :+ l -> Line 2 r :+ l -> Maybe (Point 2 r) intersectionPoint (l :+ _) (m :+ _) = asA (Proxy :: Proxy (Point 2 r)) $ l `intersect` m toSegments :: Ord r => [Point 2 r] -> [LineSegment 2 () r] toSegments ps = let pts = map ext $ ps in zipWith ClosedLineSegment pts (tail pts) -- | Constructs a boundingbox containing all intersections -- -- running time: \(O(n^2)\), where \(n\) is the number of input lines makeBoundingBox :: (Ord r, Fractional r) => [Line 2 r :+ l] -> Rectangle () r makeBoundingBox = grow 1 . boundingBoxList' . intersections -- | Computes all intersections intersections :: (Ord r, Fractional r) => [Line 2 r :+ l] -> [Point 2 r] intersections = mapMaybe (uncurry intersectionPoint) . allPairs -- intersections :: forall p r. (Ord r, Fractional r) -- => [Line 2 r :+ p] -> Map.Map (Point 2 r) (NonEmpty (Line 2 r :+ p)) -- intersections = Map.map sortNub . collect -- . mapMaybe (\(l,m) -> (l, m,) <$> f l m) . allPairs -- where -- f (l :+ _) (m :+ _) = asA (Proxy :: Proxy (Point 2 r)) $ l `intersect` m -- collect :: Ord k => [(v,v,k)] -> Map.Map k (NonEmpty v) -- collect = foldr f mempty -- where -- f (l,m,p) = Map.insertWith (<>) p (NonEmpty.fromList [l,m]) -- sortNub :: Ord r => NonEmpty (Line 2 r :+ p) -> NonEmpty (Line 2 r :+ p) -- sortNub = fmap (NonEmpty.head) . groupLines -- groupLines :: Ord r => NonEmpty (Line 2 r :+ p) -- -> NonEmpty (NonEmpty (Line 2 r :+ p)) -- groupLines = NonEmpty.groupWith1 L2 . NonEmpty.sortWith L2 -- -- | Newtype wrapper that allows us to sort lines -- newtype L2 r p = L2 (Line 2 r :+ p) deriving (Show) -- instance Eq r => Eq (L2 r p) where -- (L2 (Line p u :+ _)) == (L2 (Line q v :+ _)) = (p,u) == (q,v) -- instance Ord r => Ord (L2 r p) where -- (L2 (Line p u :+ _)) `compare` (L2 (Line q v :+ _)) = p `compare` q <> u `compare` v -- -- | Collect the intersection points per line -- byLine :: Ord r -- => Map.Map (Point 2 r) (NonEmpty (Line 2 r :+ p)) -- -> Map.Map (L2 r p) (NonEmpty (Point 2 r)) -- byLine = foldr f mempty . flatten . Map.assocs -- where -- flatten = concatMap (\(p,ls) -> map (\l -> (L2 l,p)) $ NonEmpty.toList ls) -- f (l,p) = Map.insertWith (<>) l $ NonEmpty.fromList [p] -- | Computes the intersections with a particular side sideIntersections :: (Ord r, Fractional r) => [Line 2 r :+ l] -> LineSegment 2 q r -> [(Point 2 r, Line 2 r :+ l)] sideIntersections ls s = let l = supportingLine s :+ undefined in List.sortOn fst . filter (flip onSegment s . fst) . mapMaybe (\m -> (,m) <$> l `intersectionPoint` m) $ ls -- | Constructs the unbounded intersections. Reported in clockwise direction. unBoundedParts :: (Ord r, Fractional r) => Rectangle () r -> [Line 2 r :+ l] -> [(Point 2 r, Maybe (Line 2 r :+ l))] unBoundedParts rect ls = [tl] <> t <> [tr] <> reverse r <> [br] <> reverse b <> [bl] <> l where sideIntersections' = over (traverse._2) Just . sideIntersections ls (t,r,b,l) = map4 sideIntersections' $ sides rect (tl,tr,br,bl) = map4 ((,Nothing) . (^.core)) $ corners rect map4 :: (a -> b) -> (a,a,a,a) -> (b,b,b,b) map4 f (a,b',c,d) = (f a, f b', f c, f d) -- | Links the vertices of the outer boundary with those in the subdivision link :: Eq r => [(Point 2 r, a)] -> PlanarSubdivision s v (Maybe e) f r -> V.Vector (Point 2 r, VertexId' s, a) link vs ps = V.fromList . map (\((p,x),(_,y)) -> (p,y,x)) . F.toList . fromJust' $ alignWith (\(p,_) (q,_) -> p == q) (CSeq.fromList vs) vs' where vs' = CSeq.fromList . map (\v -> (ps^.locationOf v,v) ) . V.toList $ boundaryVertices (outerFaceId ps) ps fromJust' = fromMaybe (error "Data.Geometry.Arrangement.link: fromJust") -------------------------------------------------------------------------------- makePairs :: [a] -> [(a,[a])] makePairs = go where go [] = [] go (x:xs) = (x,xs) : map (\(y,ys) -> (y,x:ys)) (go xs) allPairs :: [a] -> [(a,a)] allPairs ys = go ys where go [] = [] go (x:xs) = map (x,) xs ++ go xs -- | Given a predicate that tests if two elements of a CSeq match, find a -- rotation of the seqs such at they match. -- -- Running time: \(O(n)\) alignWith :: (a -> b -> Bool) -> CSeq.CSeq a -> CSeq.CSeq b -> Maybe (CSeq.CSeq (a,b)) alignWith p xs ys = CSeq.zipL xs <$> CSeq.findRotateTo (p (CSeq.focus xs)) ys -------------------------------------------------------------------------------- -- | Given an Arrangement and a line in the arrangement, follow the line -- through he arrangement. -- traverseLine :: (Ord r, Fractional r) => Line 2 r -> Arrangement s l v (Maybe e) f r -> [Dart s] traverseLine l arr = let md = findStart l arr dup x = (x,x) in maybe [] (List.unfoldr (fmap dup . follow arr)) md -- | Find the starting point of the line the arrangement findStart :: (Ord r, Fractional r) => Line 2 r -> Arrangement s l v (Maybe e) f r -> Maybe (Dart s) findStart l arr = do (p,_) <- asA (Proxy :: Proxy (Point 2 r, Point 2 r)) $ l `intersect` (Boundary $ arr^.boundedArea) (_,v,_) <- findStartVertex p arr findStartDart (arr^.subdivision) v -- | Given a point on the boundary of the boundedArea box; find the vertex -- this point corresponds to. -- -- running time: \(O(\log n)\) -- -- basically; maps every point to a tuple of the point and the side the -- point occurs on. We then binary search to find the point we are looking -- for. findStartVertex :: (Ord r, Fractional r) => Point 2 r -> Arrangement s l v e f r -> Maybe (Point 2 r, VertexId' s, Maybe (Line 2 r :+ l)) findStartVertex p arr = do ss <- findSide p i <- binarySearchVec (pred' ss) (arr^.unboundedIntersections) pure $ arr^.unboundedIntersections.singular (ix i) where (t,r,b,l) = sides'' $ arr^.boundedArea sides'' = map4 (\(ClosedLineSegment a c) -> LineSegment (Closed a) (Open c)) . sides findSide q = fmap fst . List.find (onSegment q . snd) $ zip [1..] [t,r,b,l] pred' ss (q,_,_) = let Just j = findSide q x = before (ss,p) (j,q) in x == LT || x == EQ before (i,p') (j,q') = case i `compare` j of LT -> LT GT -> GT EQ | i == 2 || i == 3 -> Down p' `compare` Down q' | otherwise -> p' `compare` q' -- | Find the starting dart of the given vertex v. Reports a dart s.t. -- tailOf d = v -- -- running me: \(O(k)\) where \(k\) is the degree of the vertex findStartDart :: PlanarSubdivision s v (Maybe e) f r -> VertexId' s -> Maybe (Dart s) findStartDart ps v = V.find (\d -> isJust $ ps^.dataOf d) $ incidentEdges v ps -- the "real" dart is the one that has ata associated to it. -- | Given a dart d that incoming to v (headOf d == v), find the outgoing dart -- colinear with the incoming one. Again reports dart d' s.t. tailOf d' = v -- -- running time: \(O(k)\), where k is the degree of the vertex d points to follow :: (Ord r, Num r) => Arrangement s l v e f r -> Dart s -> Maybe (Dart s) follow arr d = V.find extends $ incidentEdges v ps where ps = arr^.subdivision v = headOf d ps (up,vp) = over both (^.location) $ endPointData d ps extends d' = let wp = ps^.locationOf (headOf d' ps) in d' /= twin d && ccw up vp wp == CoLinear -------------------------------------------------------------------------------- -- TODO: we can skip the findStart by just traversing from all boundary points -- computeFaceData :: (Arrangement s v e f r -> Dart s -> f') -- -> Arrangement s v e f r -> V.Vertex f' -- computeFaceData arr f = fmap fromJust . V.create $ do -- v <- MV.replicate (arr^.subdivision.to numFaces) Nothing -- mapM_ (computeFaceData' arr f v) $ arr^.inputLines -- pure v -- computeFaceData' arr f v l = mapM_ (assign ) traverseLine arr l --------------------------------------------------------------------------------