hgeometry-0.8.0.0: src/Data/Geometry/Arrangement/Internal.hs
{-# 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
--------------------------------------------------------------------------------