debugging and common halfspace intersection

hgeometry/hgeometry.cabal

@@ -259,6 +259,7 @@ library kernel
         HGeometry.Line.Class
         HGeometry.Line.PointAndVector
 
+        HGeometry.Line.General
         HGeometry.Line.LineEQ
         HGeometry.Line.NonVertical.Class
 
@@ -269,6 +270,7 @@ library kernel
         HGeometry.Duality
 
         HGeometry.HalfSpace
+        HGeometry.HalfSpace.Class
         HGeometry.HalfLine
 
         -- HGeometry.HalfPlane
@@ -378,7 +380,7 @@ library
 
         HGeometry.IntervalTree
 
-        -- HGeometry.HalfPlane.CommonIntersection
+        HGeometry.HalfPlane.CommonIntersection
 
         HGeometry.Polygon.Triangulation
         HGeometry.Polygon.Triangulation.TriangulateMonotone
@@ -667,6 +669,7 @@ test-suite with-ipe-hspec
       Polygon.Triangulation.TriangulateSpec
       LineSegment.Intersection.BentleyOttmannSpec
       PlaneGraph.RenderSpec
+      HalfPlane.CommonIntersectionSpec
 
   hs-source-dirs:   test-with-ipe/test
   build-depends:

hgeometry/kernel/src/HGeometry/HalfPlane.hs

@@ -19,6 +19,10 @@ data GHalfPlane line r = LeftOf  !r
 type HalfPlane r = GHalfPlane (LineEQ r) r
 type instance NumType (HalfPlane r) = r
 
+
+instance HalfPlane_ (GHalfPlane line r) r where
+  type
+
 -- type HalfPlane r = HalfSpace 2 r
 
 -- class

hgeometry/kernel/src/HGeometry/HalfSpace/Class.hs

@@ -12,8 +12,8 @@ import HGeometry.Sign
 
 -- | Types modelling halfspaces.
 class ( HyperPlane_ (BoundingHyperPlane halfSpace d r) d r
-      , Dimension halfSpace ~ d
-      , NumType halfSpace ~ r
+      , Dimension halfSpace ~ d, Dimension (BoundingHyperPlane halfSpace d r) ~ d
+      , NumType halfSpace ~ r,   NumType (BoundingHyperPlane halfSpace d r) ~ r
       )
        => HalfSpace_ halfSpace d r | halfSpace -> d,
                                      halfSpace -> r where

hgeometry/kernel/src/HGeometry/HyperPlane/Class.hs

@@ -120,7 +120,7 @@ class ( NumType hyperPlane ~ r
   {-# INLINE hyperPlaneEquation #-}
 
   -- | Get the normal vector of the hyperplane. The vector points into
-  -- the negative halfspace.
+  -- the positive halfspace.
   --
   -- >>> normalVector myVerticalLine
   -- Vector2 (-1.0) 0.0
@@ -130,7 +130,7 @@ class ( NumType hyperPlane ~ r
   default normalVector :: (KnownNat d, Num r, Eq r, 1 <= d)
                        => hyperPlane -> Vector d r
   normalVector h = let a = suffix $ hyperPlaneEquation h
-                   in if signum (a^.last) == 1 then negated a else a
+                   in if signum (a^.last) == 1 then a else negated a
   {-# INLINE normalVector #-}
   -- https://en.wikipedia.org/wiki/Normal_(geometry)#Hypersurfaces_in_n-dimensional_space
   -- states that: if the hyperplane is defined as the solution set of a single linear equation
@@ -159,9 +159,10 @@ class ( NumType hyperPlane ~ r
   q `onHyperPlane` h = (== 0) $ evalHyperPlaneEquation h q
   {-# INLINE onHyperPlane #-}
 
-  -- | Test if a point lies on a hyperplane. For non-vertical
-  -- hyperplanes, returns whether the point is *above* the hyperplane
-  -- or not with respect to the d^th dimension.
+  -- | Test if a point lies on a hyperplane. For non-vertical hyperplanes, returns whether
+  -- the point is *above* the hyperplane or not with respect to the d^th
+  -- dimension. (I.e. if (and only if) q lies above the hyperplane h, then q `onSideTest`
+  -- h return GT.)
   --
   -- For vertical hyperplanes (with respect to dimension d), we return 'LT' when the point
   -- is on the "left".
@@ -256,7 +257,8 @@ class HyperPlane_ hyperPlane d r
                          => Vector (d+1) r -> hyperPlane
 
 
-  -- | Construct a Hyperplane from a point and a normal.
+  -- | Construct a Hyperplane from a point and a normal. The normal points into the halfplane
+  -- for which the side-test is positive.
   --
   -- >>> myVerticalLine == fromPointAndNormal (Point2 5 30) (Vector2 (-1) 0)
   -- True
@@ -270,7 +272,7 @@ class HyperPlane_ hyperPlane d r
                              => point -> Vector d r -> hyperPlane
   fromPointAndNormal q n = hyperPlaneFromEquation $ cons a0 n
     where
-      a0 = negate $ (q^.vector) `dot` n
+      a0 = (q^.vector) `dot` n
   {-# INLINE fromPointAndNormal #-}
 
 

hgeometry/kernel/test/HGeometry/HalfSpaceSpec.hs

@@ -3,11 +3,11 @@ module HGeometry.HalfSpaceSpec
 
 import HGeometry.HalfSpace
 import HGeometry.HyperPlane
+import HGeometry.HyperPlane.NonVertical
 import HGeometry.Intersection
 import HGeometry.Kernel.Instances ()
 import HGeometry.Line
 import HGeometry.Point
--- import HGeometry.Vector hiding (head)
 import Test.Hspec
 import Test.Hspec.QuickCheck
 
@@ -44,6 +44,11 @@ spec = describe "halfspace Tests" $ do
              let n = normalVector h
                  p = pointOn h
              in (p .+^ n) `intersects` HalfSpace Positive h
+         prop "normal vector points into positive halfspace (NonVertical HyperPlane)" $
+           \(h :: NonVerticalHyperPlane 3 Rational) ->
+             let n = normalVector h
+                 p = pointOn h
+             in (p .+^ n) `intersects` HalfSpace Positive h
          prop "normal vector points into positive halfspace (LineEQ)" $
            \(h :: LineEQ Rational) ->
              let n = normalVector h

hgeometry/kernel/test/HGeometry/HyperPlaneSpec.hs

@@ -157,22 +157,22 @@ spec = describe "HyperPlane Tests" $ do
 
          prop "fromPointAnNormal and sideTest consistent for HyperPlane " $
            \(p :: Point 2 R) n ->
-             n /= Vector2 0 0 ==>
+             allOf components (>0) n ==>
              ((p .+^ n) `onSideTest` (fromPointAndNormal p n :: HyperPlane 2 R))
              `shouldBe` GT
          prop "fromPointAnNormal and sideTest consistent for NonVerticalHyperPlane " $
            \(p :: Point 2 R) n ->
-             n /= Vector2 0 0 ==>
+             allOf components (>0) n ==>
              ((p .+^ n) `onSideTest` (fromPointAndNormal p n :: NonVerticalHyperPlane 2 R))
              `shouldBe` GT
          prop "fromPointAnNormal and sideTest consistent for LineEQ" $
            \(p :: Point 2 R) n ->
-             n /= Vector2 0 0 ==>
+             allOf components (>0) n ==>
              ((p .+^ n) `onSideTest` (fromPointAndNormal p n :: LineEQ R))
              `shouldBe` GT
          prop "normalVector and fromPointAndNormal consistent (HyperPlane 2 R)" $
            \(p :: Point 2 R) n ->
-             n /= Vector2 0 0 ==>
+             allOf components (>0) n ==>
              normalVector (fromPointAndNormal p n :: HyperPlane 2 R) `shouldBe` n
 
          prop "nonVertical sidetest means above (NonVHyperplane 2)" $

hgeometry/src/HGeometry/HalfPlane/CommonIntersection.hs

@@ -9,27 +9,33 @@ module HGeometry.HalfPlane.CommonIntersection
 
 import           Control.Lens
 import           Data.Default.Class
+import qualified Data.Foldable as F
+import           Data.Foldable1
+import qualified Data.List as List
+import           Data.List.NonEmpty (NonEmpty(..))
 import qualified Data.List.NonEmpty as NonEmpty
 import           Data.Ord (comparing)
 import           Data.Sequence (Seq(..))
 import qualified Data.Sequence as Seq
-import qualified Data.Vector as V
+import qualified Data.Vector as Vector
 import           HGeometry.Duality
 import           HGeometry.Ext
 import           HGeometry.Foldable.Sort
+import           HGeometry.HalfSpace
 import           HGeometry.HyperPlane.Class
 import           HGeometry.HyperPlane.NonVertical
 import           HGeometry.Intersection
 import           HGeometry.Line
 import           HGeometry.Point
 import           HGeometry.Polygon.Convex
-
+import           HGeometry.Vector
 --------------------------------------------------------------------------------
 
 -- | Common intersection of a bunch of halfplanes
 data CommonIntersection halfPlane r =
     EmptyIntersection
   | Bounded (ConvexPolygon (Point 2 r :+ halfPlane))
+  | Slab halfPlane () -- TODO needs a boundingLLine
     -- ^ each vertex stores the interior halfplane of the CCW-edge it is incident to.
   | Unbounded (LowerChain halfPlane r)
   deriving (Show,Eq)
@@ -56,72 +62,40 @@ bimap' f g (LowerChain hs h) = LowerChain (fmap (bimap f g) hs) (f h)
 -- | Computes the common intersection of a \(n\) halfplanes.
 --
 -- running time: \(O(n\log n)\)
-commonIntersection     :: ( Foldable f, Functor f
+commonIntersection     :: ( Foldable1 f, Functor f
                           , HyperPlane_ halfPlane 2 r  -- this is not quite right yet
                           , Fractional r, Ord r
                           )
                        => f halfPlane -> CommonIntersection halfPlane r
-commonIntersection hs0 = undefined
-  where
-    hs = partitionHalfPlanes hs0
-    lb = fromMaybe . minimumByOf traverse (comparing (^.core)) $ lefts hs
-      -- common intersection of the left halfplanes
-    rb = fromMaybe . maximumByOf traverse (comparing (^.core)) $ rights hs
-      -- common intersection of the right halfplanes
-    fromMaybe = maybe EntirePlane BoundedBy
-
-    bb = lowerBoundary $ aboves hs
-    ub = upperBoundary $ belows hs
-
-
-
-
-
-
-
-data HalfPlanes verticals nonVerticals = HalfPlanes { lefts  :: !verticals
-                                                    , rights :: !verticals
-                                                    , aboves :: !nonVerticals
-                                                    , belows :: !nonVerticals
-                                                    }
-                                       deriving (Show,Read,Eq,Ord,Functor,Foldable,Traversable)
-instance (Semigroup verticals, Semigroup nonVerticals)
-         => Semigroup (HalfPlanes verticals nonVerticals) where
-  (HalfPlanes ls rs as bs) <> (HalfPlanes ls' rs' as' bs') =
-    HalfPlanes (ls <> ls') (rs <> rs') (as <> as') (bs <> bs')
-instance (Monoid verticals, Monoid nonVerticals)
-         => Monoid (HalfPlanes verticals nonVerticals) where
-  mempty = HalfPlanes mempty mempty mempty mempty
-
-
--- | Partition the halfplanes into left halfplanes (bounded from the right), right
--- halfplanes, aboves, and belows.
-partitionHalfPlanes :: (Foldable f
-                       )
-                    => f halfPlane -> HalfPlanes [r :+ halfPlane] [LineEQ r :+ halfPlane]
-partitionHalfPlanes = foldMap f mempty
-  where
-    f hl = undefined
+commonIntersection hs0 = case NonEmpty.nonEmpty <$> foldMap partitionHalfPlanes hs0 of
+  Vector2 Nothing       Nothing       -> error "commonIntersection: absurd, cannot happen"
+  Vector2 Nothing       (Just aboves) -> Unbounded $ lowerBoundary aboves
+  Vector2 (Just belows) Nothing       -> Unbounded $ upperBoundary belows
+  Vector2 (Just belows) (Just aboves) -> let bb = lowerBoundary aboves
+                                             ub = upperBoundary belows
+                                         in undefined
+
+-- | partitions the halfplanes into those with a negative sign (the first component of the
+-- result) and positive signs.
+partitionHalfPlanes :: (Foldable f, HalfSpace_ halfPlane d r
+                       ) => f halfPlane -> Vector 2 [halfPlane]
+partitionHalfPlanes = uncurry Vector2
+                    . List.partition (\hl -> hl^.halfSpaceSign == Negative) . F.toList
 
 --------------------------------------------------------------------------------
 
-data LowerBoundary chain = EntirePlane | BoundedBy chain deriving (Show,Eq,Functor)
-
 -- | Given the bounding lines of a bunch of halfplanes that are all bounded from below,
 -- computes their common intersection.
 --
---
 -- running time: O(n\log n)
-lowerBoundary :: ( NonVerticalHyperPlane_ boundingLine 2 r
-                 , Foldable f, Fractional r, Ord r
+lowerBoundary :: ( HalfPlane_ halfSpace r
+                 , Foldable1 f, Fractional r, Ord r
                  )
-              => f boundingLine -> LowerBoundary (LowerChain boundingLine r)
-lowerBoundary = initialize . dropParallel . sortOnCheap @V.Vector dualPoint
+              => f halfSpace -> LowerChain (BoundingHyperPlane halfSpace 2 r) r
+lowerBoundary = initialize . dropParallel . sortOnCheap @Vector.Vector increasingSlope
                 -- we sort lexicographically on increasing slope and decreasing intercept
   where
-    initialize = \case
-      []   -> EntirePlane
-      h:hs -> BoundedBy $ go (LowerChain mempty h) hs
+    initialize (h :| hs) = go (LowerChain mempty h) hs
 
     -- we go through the halfplanes by increasing slope. That means the newest bounding
     -- line is the steepest, and therefore is guaranteed to appear as the rightmost
@@ -142,6 +116,15 @@ lowerBoundary = initialize . dropParallel . sortOnCheap @V.Vector dualPoint
       where
         toLineEQ = MkLineEQ . NonVerticalHyperPlane . view hyperPlaneCoefficients
 
+-- | computes the slope of the bounding line of the halfplane; vertical lines are
+-- represented first (in left to right order) ,and then lines are ordered
+-- lexicographically on increasing slope. (Hence the slightly weird return type).
+increasingSlope ::  ( HalfPlane_ halfSpace r, Fractional r, Ord r
+                    ) => halfSpace -> Either r (Vector 2 r)
+increasingSlope hl = undefined -- case hl^.boundingHyperPlane
+  -- to dualPoint.vector of
+  --                      Vector2
+
 
 -- | Drop the edges whose left-endpoint lies below h'
 dropFrom                        :: (HyperPlane_ boundingLine 2 r, Ord r, Num r)
@@ -159,15 +142,17 @@ dropFrom (LowerChain hs0 h0) h' = go h0 hs0
 
 --------------------------------------------------------------------------------
 
--- | Given the bounding lines of a bunch of halfplanes that are all bounded from above,
--- computes their common intersection.
---
--- running time: O(n\log n)
-upperBoundary :: ( NonVerticalHyperPlane_ boundingLine 2 r
-                 , Fractional r, Ord r, Foldable f, Functor f
-                 )
-              => f boundingLine -> LowerBoundary (LowerChain boundingLine r)
-upperBoundary = fmap (bimap' flipY (over yCoord negate)) . lowerBoundary . fmap flipY
-                -- by flipping the plane, and using the existing lowerBoundary machinery
-  where
-    flipY = over (hyperPlaneCoefficients.traverse) negate
+upperBoundary = undefined
+
+-- -- | Given the bounding lines of a bunch of halfplanes that are all bounded from above,
+-- -- computes their common intersection.
+-- --
+-- -- running time: O(n\log n)
+-- upperBoundary :: ( NonVerticalHyperPlane_ boundingLine 2 r
+--                  , Fractional r, Ord r, Foldable f, Functor f
+--                  )
+--               => f boundingLine -> LowerBoundary (LowerChain boundingLine r)
+-- upperBoundary = fmap (bimap' flipY (over yCoord negate)) . lowerBoundary . fmap flipY
+--                 -- by flipping the plane, and using the existing lowerBoundary machinery
+--   where
+--     flipY = over (hyperPlaneCoefficients.traverse) negate

todo.org

@@ -29,6 +29,22 @@ this can't be right:
 
 we should really use the in plane point p
 
+* TODO normal vector business
+
+want:
+
+- normalVector points into the positive halfspace
+- normalVector (fromPointAnNormal p n) == n
+
+
+- onSideTest p
+
+
+
+
+
+
+
 * DONE write some tests for testing line x line intersections of the various forms
 
 basically test that converting to a diff line type results in
@@ -92,9 +108,12 @@ intersecting halflines with boxes seems to go wrong somehow.
 
 * TODO polygon triangulation
 ** DONE triangulate monotone
-** TODO triangiulate non-monotone
+** DONE triangiulate non-monotone
 *** DONE split into non-monotone parts
-*** TODO graph representation
+*** DONE graph representation
+** TODO triangulate world demo/benchmark
+** TODO polygons with holes
+*** TODO represent polygons with holes
 * DONE polyline simplification
 ** DONE imai iri
 ** DONE DP