diff --git a/hgeometry/hgeometry.cabal b/hgeometry/hgeometry.cabal index ee1446307..ca190111d 100644 --- a/hgeometry/hgeometry.cabal +++ b/hgeometry/hgeometry.cabal @@ -403,6 +403,9 @@ library HGeometry.RangeTree HGeometry.RangeTree.Base + HGeometry.VerticalRayShooting + HGeometry.VerticalRayShooting.PersistentSweep + other-modules: HGeometry.Polygon.Simple.Implementation HGeometry.Polygon.Simple.InPolygon @@ -417,12 +420,12 @@ library HGeometry.LowerEnvelope.Type HGeometry.LowerEnvelope.Connected.Type HGeometry.LowerEnvelope.Connected.FromVertexForm + -- HGeometry.LowerEnvelope.Sample + -- HGeometry.LowerEnvelope.AtMostThree -- HGeometry.LowerEnvelope.Triangulate - -- HGeometry.LowerEnvelope.Sample - HGeometry.VoronoiDiagram.ViaLowerEnvelope @@ -629,6 +632,7 @@ test-suite hspec SegmentTreeSpec SegmentTree.R2Spec RangeTreeSpec + VerticalRayShootingSpec hs-source-dirs: test build-depends: diff --git a/hgeometry/kernel/src/HGeometry/Line/PointAndVector.hs b/hgeometry/kernel/src/HGeometry/Line/PointAndVector.hs index 35afc4dec..b7673769d 100644 --- a/hgeometry/kernel/src/HGeometry/Line/PointAndVector.hs +++ b/hgeometry/kernel/src/HGeometry/Line/PointAndVector.hs @@ -33,11 +33,12 @@ import qualified Data.Foldable as F import Data.Ord (comparing) import GHC.Generics (Generic) import GHC.TypeLits +import HGeometry.Ext import HGeometry.HyperPlane.Class import HGeometry.Intersection import HGeometry.Line.Class -import HGeometry.Line.LineEQ import HGeometry.Line.Intersection +import HGeometry.Line.LineEQ import HGeometry.Point -- import HGeometry.Point.EuclideanDistance -- import HGeometry.Point.Orientation.Degenerate @@ -243,6 +244,9 @@ instance ( Ord r class HasSupportingLine t where supportingLine :: t -> LinePV (Dimension t) (NumType t) +instance HasSupportingLine t => HasSupportingLine (t :+ extra) where + supportingLine = supportingLine . view core + instance HasSupportingLine (LinePV d r) where supportingLine = id diff --git a/hgeometry/src/HGeometry/LowerEnvelope/DivideAndConquer.hs b/hgeometry/src/HGeometry/LowerEnvelope/DivideAndConquer.hs index b15682c93..605e3b798 100644 --- a/hgeometry/src/HGeometry/LowerEnvelope/DivideAndConquer.hs +++ b/hgeometry/src/HGeometry/LowerEnvelope/DivideAndConquer.hs @@ -1,20 +1,23 @@ module HGeometry.LowerEnvelope.DivideAndConquer ( lowerEnvelope - , lowerEnvelopeWith ) where - --- import Control.Monad.Random -import HGeometry.LowerEnvelope.Type -import Witherable -import System.Random.Stateful -import Control.Monad.State.Class import Control.Lens +import qualified Data.List.NonEmpty as NonEmpty +import Data.List.NonEmpty (NonEmpty(..)) +import qualified Data.Map as Map import Data.Word +import HGeometry.HyperPlane.Class +import HGeometry.HyperPlane.NonVertical +import HGeometry.LowerEnvelope.AdjListForm +import HGeometry.LowerEnvelope.EpsApproximation import qualified HGeometry.LowerEnvelope.Naive as Naive import HGeometry.LowerEnvelope.Sample - --------------------------------------------------------------------------------- +import HGeometry.LowerEnvelope.Type +import HGeometry.LowerEnvelope.VertexForm +import HGeometry.Point +import HGeometry.Properties +import Witherable -------------------------------------------------------------------------------- @@ -29,45 +32,85 @@ eps = 1/8 -- TODO figure out what this had to be again. -- | divide and conquer algorithm -- --- expected running time: \(O(n \log n)\) -lowerEnvelope :: ( Foldable f, Ord r, Fractional r - , RandomGen gen, MonadState gen m - ) - => f (Plane r) -> m (LowerEnvelope [] Boxed.Vector r) -lowerEnvelope hs - | n <= nZero = pure $ Naive.lowerEnvelope hs - | otherwise = do ss <- sample p hs - (env, prisms) <- computePrisms hs ss - subEnvs <- mapM (over conflictList (SubEnv . lowerEnvelope)) prisms - merge env subEnvs +-- running time: \(O(n \log n)\) +lowerEnvelope :: ( Plane_ plane r + , Ord r, Fractional r, Foldable f, Functor f, Ord plane + , Show plane, Show r + ) => f plane -> LowerEnvelope plane +lowerEnvelope = fromVertexForm . lowerEnvelopeVertexForm + +-- FIXME: make sure not all planes are parallel first, otherwise the triangulatedEnvelope part is kind of weird. + +-- | Compute the vertices of the lower envelope +-- +-- +-- running time: \(O(n \log n)\) +lowerEnvelopeVertexForm :: forall f plane r. + ( Plane_ plane r + , Ord r, Fractional r, Foldable f, Ord plane + ) => f plane -> VertexForm plane +lowerEnvelopeVertexForm hs + | n <= nZero = Naive.lowerEnvelopeVertexForm hs + | otherwise = undefined where + r = undefined + s = undefined + + as = epsApproximation r hs + env = triangulatedLowerEnvelope as + + + conflictLists' = computeConflictLists env hs + + superCells = formSuperCells s env + + conflictLists = undefined --- combineConlictLists + + + + + -- do ss <- sample p hs + -- (env, prisms) <- computePrisms hs ss + -- subEnvs <- mapM (over conflictList lowerEnvelope) prisms + -- merge env subEnvs n = length hs - s = n ^^^ (1-eps) - p = probability s n + -- s = n ^^^ (1-eps) + -- p = probability s n -data SubEnv plane = SubEnv (LowerEnvelope [] Boxed.Vector (NumType plane)) --- | Run the divide and conquer algorithm with a given generator. --- --- expected running time: \(O(n \log n)\) -lowerEnvelopeWith :: ( Foldable f, Ord r, Fractional r - , RandomGen gen - ) - => gen -> f (Plane r) -> m (LowerEnvelope [] Boxed.Vector r) -lowerEnvelopeWith gen = runStateGen gen . lowerEnvelope +type TriangulatedLowerEnvelope plane = LowerEnvelope' plane --------------------------------------------------------------------------------- +triangulatedLowerEnvelope :: ( Plane_ plane r + , Ord r, Fractional r, Foldable f, Functor f, Ord plane + , Show plane, Show r + ) => f plane -> TriangulatedLowerEnvelope plane +triangulatedLowerEnvelope = undefined --- | represents the result after computing the lower envelope on the conflict lists -data SubEnv plane = SubEnv (LowerEnvelope [] Boxed.Vector (NumType plane)) --------------------------------------------------------------------------------- +type SuperCell plane = SuperCell' (NumType plane) plane +data SuperCell' r plane = SuperCell + -- { boundary :: SimplePolygon (Point 2 r) + -- , internalVertices :: [BoundedVertex plane] + -- } --- | Given the gloal envelope, and the envelopes of of the prisms, --- merge them into one big structure -merge :: LowerEnvelope [] Boxed.Vector r - -> f (Prism SubEnv r) - -> LowerEnvelope [] Boxed.Vector r -merge env _subEnvs = pure env -- TODO --------------------------------------------------------------------------------- + +formSuperCells :: Int -- ^ the number of triangles s in each super cell + -> TriangulatedLowerEnvelope plane + -> NonEmpty (SuperCell plane) +formSuperCells = undefined + + +-- newtype VerticesWithConflictLists + +newtype Vertex' plane = Vertex' (Point 3 (NumType plane)) + + +computeConflictLists :: TriangulatedLowerEnvelope plane + -> f plane + -> Map.Map (Vertex' plane) (f plane) +computeConflictLists = undefined + + +mergeConflictLists :: Map.Map (Vertex' plane) (f plane) +mergeConflictLists = undefined diff --git a/hgeometry/src/HGeometry/VerticalRayShooting.hs b/hgeometry/src/HGeometry/VerticalRayShooting.hs new file mode 100644 index 000000000..43e667e35 --- /dev/null +++ b/hgeometry/src/HGeometry/VerticalRayShooting.hs @@ -0,0 +1,12 @@ +-------------------------------------------------------------------------------- +-- | +-- Module : HGeometry.VerticalRayShooting +-- Copyright : (C) Frank Staals +-- License : see the LICENSE file +-- Maintainer : Frank Staals +-------------------------------------------------------------------------------- +module HGeometry.VerticalRayShooting + ( module HGeometry.VerticalRayShooting.PersistentSweep + ) where + +import HGeometry.VerticalRayShooting.PersistentSweep diff --git a/hgeometry/src/HGeometry/VerticalRayShooting/PersistentSweep.hs b/hgeometry/src/HGeometry/VerticalRayShooting/PersistentSweep.hs new file mode 100644 index 000000000..b5a52d5f8 --- /dev/null +++ b/hgeometry/src/HGeometry/VerticalRayShooting/PersistentSweep.hs @@ -0,0 +1,219 @@ +{-# Language TemplateHaskell #-} +-------------------------------------------------------------------------------- +-- | +-- Module : HGeometry.VerticalRayShooting.PersistentSweep +-- Copyright : (C) Frank Staals +-- License : see the LICENSE file +-- Maintainer : Frank Staals +-------------------------------------------------------------------------------- +module HGeometry.VerticalRayShooting.PersistentSweep + ( VerticalRayShootingStructure + , StatusStructure + -- , leftMost, sweepStruct + + -- * Building the Data Structure + , verticalRayShootingStructure + -- * Querying the Data Structure + , segmentAbove, segmentAboveOrOn + , findSlab + , lookupAbove, lookupAboveOrOn, searchInSlab + ) where + +import Control.Lens hiding (contains, below) +import Data.Foldable (toList) +import qualified Data.List as List +import Data.List.NonEmpty (NonEmpty(..)) +import qualified Data.List.NonEmpty as NonEmpty +import Data.Semigroup.Foldable +import qualified Data.Set as SS -- status struct +import qualified Data.Vector as V +import HGeometry.Algorithms.BinarySearch (binarySearchIn) +import HGeometry.Ext +import HGeometry.Line.PointAndVector +import HGeometry.LineSegment +import HGeometry.Point +import HGeometry.Properties +import qualified HGeometry.Set.Util as SS + +-------------------------------------------------------------------------------- + +-- | The vertical ray shooting data structure +type VerticalRayShootingStructure lineSegment = + VerticalRayShootingStructure' (NumType lineSegment) lineSegment + +-- | The implementatyion of the vertical ray shooting data structure +data VerticalRayShootingStructure' r lineSegment = + VerticalRayShootingStructure { _leftMost :: r + -- ^ x-coordinate of the leftmost vertex/endpoint + , _sweepStruct :: V.Vector (r :+ StatusStructure lineSegment) + -- ^ entry (r :+ s) means that "just" left of "r" the + -- status structure is 's', i.e up to 'r' + } deriving (Show,Eq) + +-- | The status structure +type StatusStructure lineSegment = SS.Set lineSegment + +makeLensesWith (lensRules&generateUpdateableOptics .~ False) ''VerticalRayShootingStructure' + +-------------------------------------------------------------------------------- +-- * Building the DS + +-- | Given a set of \(n\) interiorly pairwise disjoint *closed* segments, +-- compute a vertical ray shooting data structure. (i.e. the +-- endpoints of the segments may coincide). +-- +-- pre: no vertical segments +-- +-- running time: \(O(n\log n)\). +-- space: \(O(n\log n)\). +verticalRayShootingStructure :: ( LineSegment_ lineSegment point + , Point_ point 2 r + , Ord r, Fractional r, Foldable1 f) + => f lineSegment + -> VerticalRayShootingStructure lineSegment +verticalRayShootingStructure ss = VerticalRayShootingStructure (eventX e) (sweep' events) + where + events@(e :| _) = fmap combine + . NonEmpty.groupAllWith1 eventX + . foldMap1 (toEvents . orientLR) + $ ss + sweep' = V.fromList . toList . sweep + + toEvents seg = NonEmpty.fromList [ (seg^.start.xCoord) :+ Insert seg :| [] + , (seg^.end.xCoord) :+ Delete seg :| [] + ] + + +-- | Given a bunch of events happening at the same time, merge them into a single event +-- where we apply all actions. +combine :: NonEmpty (Event r lineSegment) -> Event r lineSegment +combine es@((x :+ _) :| _) = x :+ foldMap1 eventActions es + +---------------------------------------- + +data Action a = Insert a | Delete a deriving (Show,Eq) + +{- HLINT ignore "Avoid lambda using `infix`" -} +interpret :: Action a -> (a -> a -> Ordering) -> SS.Set a -> SS.Set a +interpret = \case + Insert s -> \cmp -> SS.insertBy cmp s + Delete s -> \cmp -> SS.deleteAllBy cmp s + + +type Event r lineSegment = r :+ NonEmpty (Action lineSegment) + +eventX :: Event r lineSegment -> r +eventX = view core + +eventActions :: Event r lineSegment -> NonEmpty (Action lineSegment) +eventActions = view extra + +---------------------------------------- + +-- | Runs the sweep building the data structure from left to right. +sweep :: ( LineSegment_ lineSegment point, Point_ point 2 r + , Ord r, Fractional r + ) + => NonEmpty (Event r lineSegment) -> NonEmpty (r :+ StatusStructure lineSegment) +sweep es = NonEmpty.fromList + . snd . List.mapAccumL h SS.empty + $ zip (toList es) (NonEmpty.tail es) + where + h ss evts = let x :+ ss' = handle ss evts in (ss',x :+ ss') + +-- | Given the current status structure (for left of the next event +-- 'l'), and the next two events (l,r); essentially defining the slab +-- between l and r, we construct the status structure for in the slab (l,r). +-- returns the right boundary and this status structure. +handle :: (Ord r, Fractional r, LineSegment_ lineSegment point, Point_ point 2 r) + => StatusStructure lineSegment + -> (Event r lineSegment, Event r lineSegment) + -> r :+ StatusStructure lineSegment +handle ss ( l :+ acts + , r :+ _) = let mid = (l+r)/2 + runActionAt x act = interpret act (ordAtX x) + in r :+ foldr (runActionAt mid) ss (orderActs acts) + -- run deletions first + +-- | 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 + Delete _ -> True + Insert _ -> False + ) acts + in NonEmpty.fromList $ ins <> dels + + +-------------------------------------------------------------------------------- +-- * Querying the DS + +-- | Find the segment vertically strictly above query point q, if it +-- exists. +-- +-- \(O(\log n)\) +segmentAbove :: ( LineSegment_ lineSegment point, Point_ point 2 r + , Point_ queryPoint 2 r + , Ord r, Num r, HasSupportingLine lineSegment + ) => queryPoint -> VerticalRayShootingStructure lineSegment + -> Maybe lineSegment +segmentAbove q ds = findSlab q ds >>= lookupAbove q + +-- | Find the segment vertically query point q, if it exists. +-- +-- \(O(\log n)\) +segmentAboveOrOn :: ( LineSegment_ lineSegment point, Point_ point 2 r + , Point_ queryPoint 2 r + , Ord r, Num r, HasSupportingLine lineSegment + ) => queryPoint -> VerticalRayShootingStructure lineSegment + -> Maybe lineSegment +segmentAboveOrOn q ds = findSlab q ds >>= lookupAboveOrOn q + + + +-- | Given a query point, find the (data structure of the) slab containing the query point +-- +-- \(O(\log n)\) +findSlab :: ( LineSegment_ lineSegment point, Point_ point 2 r + , Point_ queryPoint 2 r + , Ord r, Num r, HasSupportingLine lineSegment + ) + => queryPoint -> VerticalRayShootingStructure lineSegment + -> Maybe (StatusStructure lineSegment) +findSlab q ds | q^.xCoord < ds^.leftMost = Nothing + | otherwise = view extra + <$> binarySearchIn (q `leftOf `) (ds^.sweepStruct) + where + q' `leftOf` (r :+ _) = q'^.xCoord <= r + +-------------------------------------------------------------------------------- +-- * Querying in a single slab + +-- | Finds the segment containing or above the query point 'q' +-- +-- \(O(\log n)\) +lookupAboveOrOn :: ( LineSegment_ lineSegment point, Point_ point 2 r + , Point_ queryPoint 2 r + , Ord r, Num r, HasSupportingLine lineSegment + ) + => queryPoint -> StatusStructure lineSegment -> Maybe lineSegment +lookupAboveOrOn q = searchInSlab (not . (q `liesAbove`)) + +-- | Finds the first segment strictly above q +-- +-- \(O(\log n)\) +lookupAbove :: ( LineSegment_ lineSegment point, Point_ point 2 r + , Point_ queryPoint 2 r + , Ord r, Num r, HasSupportingLine lineSegment + ) + => queryPoint -> StatusStructure lineSegment -> Maybe lineSegment +lookupAbove q = searchInSlab (q `liesBelow`) + +-- | generic searching function +searchInSlab :: (LineSegment_ lineSegment point, Point_ point 2 r + , HasSupportingLine lineSegment, Num r) + => (LinePV 2 r -> Bool) + -> StatusStructure lineSegment -> Maybe lineSegment +searchInSlab p = binarySearchIn (p . supportingLine) + + +---------------------------------------------------------------------------------- diff --git a/hgeometry/test/VerticalRayShootingSpec.hs b/hgeometry/test/VerticalRayShootingSpec.hs new file mode 100644 index 000000000..652fecddb --- /dev/null +++ b/hgeometry/test/VerticalRayShootingSpec.hs @@ -0,0 +1,42 @@ +module VerticalRayShootingSpec where + +import HGeometry.VerticalRayShooting.PersistentSweep +import Control.Lens hiding (contains, below) +import HGeometry.Ext +import Data.List.NonEmpty (NonEmpty(..)) +import qualified Data.List.NonEmpty as NonEmpty +import HGeometry.Number.Real.Rational +import HGeometry.LineSegment +import HGeometry.Point +import Test.Hspec + +-------------------------------------------------------------------------------- + +type R = RealNumber 5 + +spec :: Spec +spec = describe "VerticalRayShooting Tests" $ do + it "manual queries on horizontal subidv" $ do + segmentAbove (Point2 5 0) test1 `shouldBe` + Just (ClosedLineSegment (Point2 0 2) (Point2 10 2) :+ 1) + segmentAbove (Point2 5 2) test1 `shouldBe` + Just (ClosedLineSegment (Point2 1 4) (Point2 12 4) :+ 2) + segmentAbove (Point2 5 1) test1 `shouldBe` + Just (ClosedLineSegment (Point2 0 2) (Point2 10 2) :+ 1) + segmentAbove (Point2 5 5) test1 `shouldBe` Nothing + segmentAbove (Point2 10 5) test1 `shouldBe` Nothing + segmentAbove (Point2 10 0) test1 `shouldBe` + Just (ClosedLineSegment (Point2 0 2) (Point2 10 2) :+ 1) + segmentAbove (Point2 10 2) test1 `shouldBe` + Just (ClosedLineSegment (Point2 1 4) (Point2 12 4) :+ 2) + segmentAbove (Point2 10 1) test1 `shouldBe` + Just (ClosedLineSegment (Point2 0 2) (Point2 10 2) :+ 1) + +test1 :: VerticalRayShootingStructure (ClosedLineSegment (Point 2 R) :+ Int) +test1 = verticalRayShootingStructure . NonEmpty.fromList $ zipWith (:+) + [ hor 2 0 10 + , hor 4 1 12 + , hor 2 10 14 + ] [1..] + where + hor y l r = ClosedLineSegment (Point2 l y) (Point2 r y)