vertical ray shooting by persistent sweep :)

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:

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
 

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

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