diff --git a/hgeometry/src/HGeometry/LowerEnvelope/Connected/FromVertexForm.hs b/hgeometry/src/HGeometry/LowerEnvelope/Connected/FromVertexForm.hs index 78a16c3c1..7048700c6 100644 --- a/hgeometry/src/HGeometry/LowerEnvelope/Connected/FromVertexForm.hs +++ b/hgeometry/src/HGeometry/LowerEnvelope/Connected/FromVertexForm.hs @@ -235,8 +235,7 @@ sortAlongBoundary face = case mv0 of u' = projectPoint $ ivLoc u v' = projectPoint $ ivLoc v --- | given an edge (u,v) that has h to its left, and all remaining vertices of the face, --- sorted in CCW order around the face starting with *v*, compute all its edgesf +-- | given the vertices of the face in CCW order, compute all its edges of the face. faceToEdges :: forall plane r. (Plane_ plane r, Ord r, Fractional r, Ord plane , Show plane, Show r ) @@ -297,14 +296,14 @@ twoVertices (Vtx h ui up uDefs) (Vtx _ vi vp vDefs) = case extractEdgeDefs h up uDefs vp vDefs >>= withUVOrder of Nothing -> error "twoVertices. absurd, h and h' don't define an edge!?" Just (EdgeDefs hr hu hv, uBeforeV) -- -> traceShowWith ("twoVertices",u,v,uBeforeV,) $ - | uBeforeV -> NonEmpty.fromList [ (vi, Edge ui h hr) - , (ui, Edge unboundedVertexId h hu) - ] + | uBeforeV -> NonEmpty.fromList [ (vi, Edge ui h hr) + , (ui, Edge unboundedVertexId h hu) + ] -- the edge must be oriented from v to u so that h is on the left - | otherwise -> NonEmpty.fromList [ (ui, Edge vi h hr) - , (vi, Edge unboundedVertexId h hv) - ] - -- the edge must be oriented from u to v + | otherwise -> NonEmpty.fromList [ (ui, Edge vi h hr) + , (vi, Edge unboundedVertexId h hv) + ] + -- the edge must be oriented from u to v where -- determine if u lies before v in the order of u and v along the intersection line of -- h and hr. @@ -315,64 +314,6 @@ twoVertices (Vtx h ui up uDefs) (Vtx _ vi vp vDefs) = -- v if v does not lie in the left haflpalne ) <$> intersectionLine' h hr - --------------------------------------------------------------------------------- - --- | Vertices in of the lower envelope in adjacencylist form. -type BoundedVertex = BoundedVertexF Seq.Seq - - --- | Given a location of a vertex v, a pair of planes h1,h2 and the --- remaining defining planes of v, computes the outgoing half-line --- from v on which h1,h2 are the lowest (if such an halfline exists). -outgoingUnboundedEdge :: ( Plane_ plane r, Ord r, Fractional r - , Foldable1 f - ) - => Point 3 r -- ^ the location of the vertex v - -> Two plane -- ^ the pair of planes for which to compute - -- the halfine - -> f plane -- ^ the other planes intersecting at v - -> Maybe (HalfLine (Point 2 r) :+ EdgeDefiners plane) -outgoingUnboundedEdge v (Two h1 h2) h3s = - intersectionLineWithDefiners h1 h2 >>= toHalfLineFrom (projectPoint v) h3s - -- todo, if there are more planes, I guess we should check if the hl is not dominated by the other - -- planes either. - --- | Given : --- --- v : the projected location of the vertex --- hs : the remaining planes defining v (typically just one plane h3) --- l : the (projection of the) line l in which planes h1 and h2 intersect (containing v) --- --- we compute the half-line eminating from v in which h1 and h2 define --- an edge incident to v. -toHalfLineFrom :: (Plane_ plane r, Foldable1 f, Fractional r, Ord r) - => Point 2 r -- ^ vertex v - -> f plane -- ^ the remaining plane(s) hs - -> LinePV 2 r :+ EdgeDefiners plane -- ^ the line l - -> Maybe (HalfLine (Point 2 r) :+ EdgeDefiners plane) -toHalfLineFrom v hs ((LinePV _ w) :+ defs@(EdgeDefiners h1 h2)) = - validate w defs <|> validate ((-1) *^ w) (EdgeDefiners h2 h1) - -- We try both directions. Note that if we reverse the direction - -- of the line, what the left/right plane is changes. - -- - -- If neither direction works, then h1,h2 do not define a good - -- direction. This should happen only when there are more than 3 - -- planes intersecting in v, i.e. in degernate sitatuions - where - -- | test if direction d is a good direction, i.e. if towards direction d - -- h1 (and thus also h2) is actually lower than all remaining defining planes. - validate d defs' = let zVal = evalAt (v .+^ d) - in if all (\h3 -> zVal h1 < zVal h3) hs - then Just (HalfLine v d :+ defs') else Nothing - - - - - - - - -------------------------------------------------------------------------------- -- * Convenience functions @@ -394,7 +335,7 @@ maxBy cmp a b = case cmp a b of _ -> a --- | shift the vector by d +-- | shift the vector by d items to the right shiftR :: Int -> NonEmptyVector v -> NonEmptyVector v shiftR d v = let n = length v in NonEmptyV.generate1 n $ \i -> v NonEmptyV.! ((i+n+d) `mod` n) diff --git a/hgeometry/src/HGeometry/LowerEnvelope/Connected/Type.hs b/hgeometry/src/HGeometry/LowerEnvelope/Connected/Type.hs index 83096b447..3c8c85c9a 100644 --- a/hgeometry/src/HGeometry/LowerEnvelope/Connected/Type.hs +++ b/hgeometry/src/HGeometry/LowerEnvelope/Connected/Type.hs @@ -17,6 +17,7 @@ module HGeometry.LowerEnvelope.Connected.Type , singleton + , BoundedVertex , BoundedVertexF(Vertex) , location, definers, location2 @@ -28,7 +29,7 @@ module HGeometry.LowerEnvelope.Connected.Type , projectedEdgeGeometries, projectedEdgeGeometry - + , outgoingUnboundedEdge , edgeIntersectionLine , intersectionLine' , intersectionLineWithDefiners @@ -84,16 +85,16 @@ instance (Ord (NumType plane), Num (NumType plane)) => IsBoxable (LowerEnvelope' boundingBox env = boundingBox . NonEmpty.fromList $ env^..boundedVertices.traverse.location -- the fromList is safe since there is alwasy at least one vertex --- instance Functor LowerEnvelope' where --- instance Foldable LowerEnvelope' where --- instance Traversable LowerEnvelope' where - --- | Traversal of the planes in the lower envelope +-- | Traversal of the planes in the lower envelope. Since we need an Ord constraint we +-- can't make LowerEnvelope' an instance of Traversable. +-- +-- Be aware that this may destroy some of the invariants. So use this function with care. traverseLowerEnvelope :: ( Applicative f, NumType plane ~ NumType plane' , Ord plane' ) => (plane -> f plane') - -> LowerEnvelope' plane -> f (LowerEnvelope' plane') + -> LowerEnvelope' plane + -> f (LowerEnvelope' plane') traverseLowerEnvelope f (LowerEnvelope v0 vs) = LowerEnvelope <$> traverse f v0 <*> traverse (traverseBoundedV f) vs @@ -121,18 +122,6 @@ singleton v = LowerEnvelope v0 (Seq.singleton v') -------------------------------------------------------------------------------- - --------------------------------------------------------------------------------- - - --- | Helper type for edgedefs -data EdgeDefs plane = EdgeDefs { common :: plane - , uNeigh :: plane - , vNeigh :: plane - } deriving (Show,Eq) - --------------------------------------------------------------------------------- - -- | The unbounded vertex. newtype UnboundedVertex plane = UnboundedVertex { _incidentEdgesU :: Seq.Seq (LEEdge plane) } deriving (Show,Eq,Functor,Foldable,Traversable) @@ -270,8 +259,9 @@ edgeIntersectionLine e = case intersectionLine' (e^.leftPlane) (e^.rightPlane) o -------------------------------------------------------------------------------- --- | +-- | Types that have an IncidentEdges' field. class HasIncidentEdges t where + -- | Lens to access the incident edges field. incidentEdges' :: Lens' (t plane) (Seq.Seq (LEEdge plane)) instance HasIncidentEdges UnboundedVertex where @@ -450,17 +440,4 @@ projectedEdgeGeometry env (ui,vi) e = case env^?!vertexAt ui of -- the unbounded vertex into the bounded vertex that means that to construct the -- halfline we actually wish to flip the orientation of the line - - --- case --- -- traceShowWith ("unbEdge", e, "h=",e^.rightPlane, "h'=",e^.leftPlane, "line'",) $ --- intersectionLine' (e^.rightPlane) (e^.leftPlane) of --- Just l -> Left $ HalfLine (v^.location2) (l^.direction) --- -- edge e is oriented from the lowerId towards the higherId, --- -- so in particular *from* the unbounded vertex into the bounded vertex --- -- that means that to construct the halfline we actually wish to --- -- flip the left and right plane, so that the halfline is directed outwards. --- Nothing -> error "projectedEdgeGeometry: absurd, no intersection between the planes" - - -------------------------------------------------------------------------------- diff --git a/hgeometry/src/HGeometry/LowerEnvelope/Naive.hs b/hgeometry/src/HGeometry/LowerEnvelope/Naive.hs index 1ddc70368..960fbc1df 100644 --- a/hgeometry/src/HGeometry/LowerEnvelope/Naive.hs +++ b/hgeometry/src/HGeometry/LowerEnvelope/Naive.hs @@ -10,22 +10,14 @@ module HGeometry.LowerEnvelope.Naive -------------------------------------------------------------------------------- -import Control.Lens import Control.Monad (guard) -import qualified Data.Sequence as Seq import qualified Data.Set as Set import HGeometry.Combinatorial.Util -import HGeometry.Foldable.Sort import HGeometry.HyperPlane.Class import HGeometry.HyperPlane.NonVertical -import HGeometry.Intersection -import HGeometry.Line -import HGeometry.Line.LineEQ import HGeometry.LowerEnvelope.AdjListForm (LowerEnvelope, fromVertexForm) import HGeometry.LowerEnvelope.VertexForm import HGeometry.Point -import HGeometry.Properties -import Hiraffe.Graph -------------------------------------------------------------------------------- @@ -61,10 +53,7 @@ lowerEnvelope = fromVertexForm . lowerEnvelopeVertexForm lowerEnvelopeVertexForm :: ( Plane_ plane r , Ord r, Fractional r, Foldable f, Ord plane ) => f plane -> VertexForm plane -lowerEnvelopeVertexForm hs = foldMap (\t -> case asVertex hs t of - Nothing -> mempty - Just v -> singleton v - ) $ uniqueTriplets hs +lowerEnvelopeVertexForm hs = foldMap (maybe mempty singleton . asVertex hs) $ uniqueTriplets hs -- | Given all planes, and a triple, computes if the triple defines a -- vertex of the lower envelope, and if so returns it. @@ -75,8 +64,8 @@ asVertex hs t@(Three h1 h2 h3) = do v <- intersectionPoint t pure $ LEVertex v (Set.fromList [h1,h2,h3]) -- | test if v lies below (or on) all the planes in hs -belowAll :: (Plane_ plane r, Ord r, Num r, Foldable f) => Point 3 r -> f plane -> Bool -belowAll v hs = all (\h -> onSideTest v h /= GT) hs +belowAll :: (Plane_ plane r, Ord r, Num r, Foldable f) => Point 3 r -> f plane -> Bool +belowAll v = all (\h -> onSideTest v h /= GT) {-# INLINE belowAll #-} diff --git a/hgeometry/src/HGeometry/VoronoiDiagram/ViaLowerEnvelope.hs b/hgeometry/src/HGeometry/VoronoiDiagram/ViaLowerEnvelope.hs index 9b11020a1..9a546720c 100644 --- a/hgeometry/src/HGeometry/VoronoiDiagram/ViaLowerEnvelope.hs +++ b/hgeometry/src/HGeometry/VoronoiDiagram/ViaLowerEnvelope.hs @@ -20,7 +20,6 @@ module HGeometry.VoronoiDiagram.ViaLowerEnvelope import Control.Lens import Data.Default.Class -import qualified Data.Map as Map import HGeometry.Box import HGeometry.Duality import HGeometry.Ext @@ -31,7 +30,6 @@ import HGeometry.LowerEnvelope.Naive (lowerEnvelopeVertexForm) import HGeometry.LowerEnvelope.VertexForm (VertexForm) import HGeometry.Point import HGeometry.Properties -import HGeometry.Vector import Hiraffe.Graph -------------------------------------------------------------------------------- diff --git a/todo.org b/todo.org index 71205f27a..6aa857c9d 100644 --- a/todo.org +++ b/todo.org @@ -105,11 +105,20 @@ intersecting halflines with boxes seems to go wrong somehow. *** TODO plane graph * TODO 3d-lower-envelope/convex hull -** DONE naive +** TODO naive +*** TODO triangulated envelope +*** TODO handle degeneracies +*** TODO handle all colinear +*** TODO cyclic sorting of the edges + ** TODO define tests -*** TODO correctly render a lower envelope/vd with 1 vertex and 3 unbounded edges -*** TODO correctly render bounded edges of some larger set of points -*** TODO correctly render unbounded edges of some larger set of points +*** DONE correctly render a lower envelope/vd with 1 vertex and 3 unbounded edges +*** DONE correctly render bounded edges of some larger set of points +*** DONE correctly render unbounded edges of some larger set of points +*** TODO properyt test that every (bounded) face is convex + +** TODO some sort of benchmarking for the naive algorithm + ** TODO Set3 type to clean up and/or speed up the fromVertexForm code ?