From 6c478dcb3e2ce29ca0523a5c316b7d58d9d08012 Mon Sep 17 00:00:00 2001 From: Frank Staals Date: Wed, 24 Apr 2024 20:51:58 +0200 Subject: [PATCH] debugging and common halfspace intersection --- hgeometry/hgeometry.cabal | 5 +- hgeometry/kernel/src/HGeometry/HalfPlane.hs | 4 + .../kernel/src/HGeometry/HalfSpace/Class.hs | 4 +- .../kernel/src/HGeometry/HyperPlane/Class.hs | 16 +-- .../kernel/test/HGeometry/HalfSpaceSpec.hs | 7 +- .../kernel/test/HGeometry/HyperPlaneSpec.hs | 8 +- .../HGeometry/HalfPlane/CommonIntersection.hs | 117 ++++++++---------- todo.org | 23 +++- 8 files changed, 101 insertions(+), 83 deletions(-) diff --git a/hgeometry/hgeometry.cabal b/hgeometry/hgeometry.cabal index 5bdf68b01..8ec8c1bac 100644 --- a/hgeometry/hgeometry.cabal +++ b/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: diff --git a/hgeometry/kernel/src/HGeometry/HalfPlane.hs b/hgeometry/kernel/src/HGeometry/HalfPlane.hs index 6d89e8098..43ad7cff7 100644 --- a/hgeometry/kernel/src/HGeometry/HalfPlane.hs +++ b/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 diff --git a/hgeometry/kernel/src/HGeometry/HalfSpace/Class.hs b/hgeometry/kernel/src/HGeometry/HalfSpace/Class.hs index 2b3377826..bf429b994 100644 --- a/hgeometry/kernel/src/HGeometry/HalfSpace/Class.hs +++ b/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 diff --git a/hgeometry/kernel/src/HGeometry/HyperPlane/Class.hs b/hgeometry/kernel/src/HGeometry/HyperPlane/Class.hs index a4a6e6998..0596c2f88 100644 --- a/hgeometry/kernel/src/HGeometry/HyperPlane/Class.hs +++ b/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 #-} diff --git a/hgeometry/kernel/test/HGeometry/HalfSpaceSpec.hs b/hgeometry/kernel/test/HGeometry/HalfSpaceSpec.hs index 9cbcac4a6..55308ffa8 100644 --- a/hgeometry/kernel/test/HGeometry/HalfSpaceSpec.hs +++ b/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 diff --git a/hgeometry/kernel/test/HGeometry/HyperPlaneSpec.hs b/hgeometry/kernel/test/HGeometry/HyperPlaneSpec.hs index f6be1b0c2..e90515070 100644 --- a/hgeometry/kernel/test/HGeometry/HyperPlaneSpec.hs +++ b/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)" $ diff --git a/hgeometry/src/HGeometry/HalfPlane/CommonIntersection.hs b/hgeometry/src/HGeometry/HalfPlane/CommonIntersection.hs index effb83a4f..64f0ef84c 100644 --- a/hgeometry/src/HGeometry/HalfPlane/CommonIntersection.hs +++ b/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 diff --git a/todo.org b/todo.org index e711ea923..f506d2b2f 100644 --- a/todo.org +++ b/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