From 53897d56a6bac02b74ca42a23d520f7eaa2a37de Mon Sep 17 00:00:00 2001 From: Frank Staals <991345+noinia@users.noreply.github.com> Date: Sun, 17 Nov 2024 11:53:29 +0100 Subject: [PATCH] Generate Arbitrary PlaneGraphs (#254) + adding a basic PLYWriter library * getting started on the random graph business * typo * using the random instances * working on generating random graphs * progress on makign it compile * generating random plane graphs * some cleaning up * remdering random plane graphs * creating a graph from connected segments * some cleaning up * trying to get the upper hull thing to work as well actually * uncommenting unused code for doctest * quick ply-writer * deleting the PlaneGraph drawing stuff * setting up some more tests * creating a demo to render a lower envelope * plywriter does something at least :) --- hgeometry-examples/draw/Main.hs | 32 ++++ hgeometry-examples/hgeometry-examples.cabal | 25 +++ hgeometry-examples/lowerEnv/Main.hs | 47 +++++ hgeometry/hgeometry.cabal | 19 +- hgeometry/ply/src/PLY/Writer.hs | 69 +++++++ .../HGeometry/PlaneGraph/Instances.hs | 177 ++++++++++++++++++ .../src/HGeometry/ConvexHull/R3/Naive/Dual.hs | 68 ++++--- .../Plane/LowerEnvelope/Connected.hs | 2 +- .../Plane/LowerEnvelope/Connected/Graph.hs | 81 ++++++-- .../LowerEnvelope/Connected/MonoidalMap.hs | 21 ++- .../LowerEnvelope/Connected/Separator.hs | 80 ++++---- .../Plane/LowerEnvelope/Connected/Type.hs | 40 ++++ .../Plane/LowerEnvelope/DivideAndConquer.hs | 57 +++--- .../Plane/LowerEnvelope/VertexForm.hs | 130 ------------- hgeometry/src/HGeometry/PlaneGraph.hs | 2 + hgeometry/src/HGeometry/PlaneGraph/Type.hs | 87 ++++++++- .../VoronoiDiagram/ViaLowerEnvelope.hs | 8 +- .../test/Plane/LowerEnvelopeSpec.hs | 8 +- hgeometry/test/ConvexHull/R3Spec.hs | 28 +++ hgeometry/test/Graphics/CameraSpec.hs | 1 - 20 files changed, 723 insertions(+), 259 deletions(-) create mode 100644 hgeometry-examples/draw/Main.hs create mode 100644 hgeometry-examples/lowerEnv/Main.hs create mode 100644 hgeometry/ply/src/PLY/Writer.hs create mode 100644 hgeometry/src-quickcheck/HGeometry/PlaneGraph/Instances.hs delete mode 100644 hgeometry/src/HGeometry/Plane/LowerEnvelope/VertexForm.hs create mode 100644 hgeometry/test/ConvexHull/R3Spec.hs diff --git a/hgeometry-examples/draw/Main.hs b/hgeometry-examples/draw/Main.hs new file mode 100644 index 000000000..e2748cb5a --- /dev/null +++ b/hgeometry-examples/draw/Main.hs @@ -0,0 +1,32 @@ +{-# LANGUAGE QuasiQuotes #-} +module Main(main) where + +import Control.Lens +import qualified Data.List.NonEmpty as NonEmpty +import HGeometry.Instances () +import HGeometry.Number.Real.Rational +import HGeometry.PlaneGraph +import HGeometry.PlaneGraph.Instances +import HGeometry.Point +import Ipe +import System.OsPath +import Test.QuickCheck + +-------------------------------------------------------------------------------- + +type R = RealNumber 5 + +renderGraph :: PlaneGraph QuickCheckWorld (Point 2 R) () () -> IpePage R +renderGraph gr = fromContent $ + concat [ [ iO $ defIO p | p <- gr^..vertices ] + , [ iO $ defIO seg | seg <- gr^..edgeSegments ] + ] + +main :: IO () +main = do + grs <- NonEmpty.fromList <$> sample' arbitrary + let outFp = [osp|foo.ipe|] + writeIpeFile outFp $ ipeFile (renderGraph <$> grs) + + -- (grs :: [PlaneGraph QuickCheckWorld (Point 2 R) () ()]) <- sample' arbitrary + -- mapM_ print grs diff --git a/hgeometry-examples/hgeometry-examples.cabal b/hgeometry-examples/hgeometry-examples.cabal index 361f0e6b8..0a76911cf 100644 --- a/hgeometry-examples/hgeometry-examples.cabal +++ b/hgeometry-examples/hgeometry-examples.cabal @@ -50,6 +50,7 @@ common setup , hgeometry:svg , hgeometry:miso , hgeometry:geojson + , hgeometry:quickcheck , hiraffe >= 0.1 , containers >= 0.6 , vector >= 0.13 @@ -76,6 +77,7 @@ common setup , file-io , filepath , geojson >= 4.1.1 + , QuickCheck >= 2.15.0 -- , dependent-map >= 0.4 -- , dependent-sum >= 0.7.1 @@ -269,3 +271,26 @@ executable hgeometry-skia Base GeometryStore GeometryStore.Helper + + +-------------------------------------------------------------------------------- +executable hgeometry-draw + import: setup, miso-setup + hs-source-dirs: draw + main-is: Main.hs + -- other-modules: + -- Options + + +-------------------------------------------------------------------------------- + +-- Renders a 3D model of the lower envelope of a set of planes +executable hgeometry-lowerEnv + import: setup, miso-setup + hs-source-dirs: lowerEnv + main-is: Main.hs + build-depends: + hgeometry:ply-writer + + -- other-modules: + -- Options diff --git a/hgeometry-examples/lowerEnv/Main.hs b/hgeometry-examples/lowerEnv/Main.hs new file mode 100644 index 000000000..234b58f84 --- /dev/null +++ b/hgeometry-examples/lowerEnv/Main.hs @@ -0,0 +1,47 @@ +{-# LANGUAGE QuasiQuotes #-} +module Main(main) where + +import Control.Lens +import qualified Data.Foldable as F +import qualified Data.List.NonEmpty as NonEmpty +import HGeometry.Ext +import HGeometry.Number.Real.Rational +import HGeometry.Plane.LowerEnvelope +import HGeometry.PlaneGraph +import HGeometry.PlaneGraph.Instances +import HGeometry.Point +import HGeometry.Triangle +import HGeometry.VoronoiDiagram.ViaLowerEnvelope +import Ipe +import PLY.Writer +import System.OsPath +import Test.QuickCheck + +-------------------------------------------------------------------------------- + +type R = RealNumber 5 + + +myPlanes = NonEmpty.fromList $ zipWith (\i p -> pointToPlane p :+ (i,p)) [0..] + [ Point2 16 80 + , Point2 64 48 + , Point2 208 128 + , Point2 176 48 + , Point2 96 112 + , Point2 128 80 + , Point2 48 144 + ] + +verticesOf = NonEmpty.fromList . foldMap F.toList . trianglesOf +trianglesOf _ = [ Triangle (origin :+ 0) (Point3 10 0 1 :+ 1) (Point3 0 10 2 :+ 2) ] + + + -- \case + -- ParallelStrips _ -> undefined + -- ConnectedEnvelope env -> undefined + +-- trianglesOf env = [] + +main :: IO () +main = renderOutputToFile [osp|myLowerEnv.ply|] (verticesOf $ lowerEnvelope myPlanes) + (trianglesOf $ lowerEnvelope myPlanes) diff --git a/hgeometry/hgeometry.cabal b/hgeometry/hgeometry.cabal index 90d84bfe0..9a8ff71db 100644 --- a/hgeometry/hgeometry.cabal +++ b/hgeometry/hgeometry.cabal @@ -50,6 +50,7 @@ common all-setup , hgeometry-combinatorial >= 1.0.0.0 && < 2 , hiraffe >= 0.1 && < 1 , containers >= 0.6 && < 1 + , nonempty-containers >= 0.3.4.5 && < 0.4 , dlist >= 0.7 && < 1.1 , bytestring >= 0.11 && < 1 , vector:vector >= 0.13 && < 1 @@ -391,7 +392,7 @@ library HGeometry.ConvexHull.Melkman -- HGeometry.ConvexHull.R3.Naive - -- HGeometry.ConvexHull.R3.Naive.Dual + HGeometry.ConvexHull.R3.Naive.Dual HGeometry.BezierSpline @@ -419,7 +420,7 @@ library HGeometry.VoronoiDiagram - + HGeometry.VoronoiDiagram.ViaLowerEnvelope HGeometry.SegmentTree HGeometry.SegmentTree.Base @@ -460,7 +461,6 @@ library -- HGeometry.Plane.LowerEnvelope.DivideAndConquer -- HGeometry.Plane.LowerEnvelope.EpsApproximation - HGeometry.VoronoiDiagram.ViaLowerEnvelope @@ -492,6 +492,7 @@ library ipe Ipe.IpeRender Ipe.IpeToIpe + -- HGeometry.PlaneGraph.Draw other-modules: @@ -570,6 +571,7 @@ library quickcheck exposed-modules: HGeometry.Instances HGeometry.Polygon.Instances + HGeometry.PlaneGraph.Instances other-modules: Paths_hgeometry autogen-modules: @@ -604,6 +606,16 @@ library vector-quickcheck build-depends: hgeometry:vector +library ply-writer + import: all-setup + visibility: public + hs-source-dirs: ply/src + exposed-modules: + PLY.Writer + build-depends: + hgeometry:point + , hgeometry:kernel + -------------------------------------------------------------------------------- -- * Test Suites -------------------------------------------------------------------------------- @@ -699,6 +711,7 @@ test-suite hspec Spec ConvexHull.ConvexHullSpec ConvexHull.MelkmanSpec + ConvexHull.R3Spec ClosestPair.ClosestPairSpec IntervalTreeSpec LowerEnvelope.RegionsSpec diff --git a/hgeometry/ply/src/PLY/Writer.hs b/hgeometry/ply/src/PLY/Writer.hs new file mode 100644 index 000000000..dcc2a6c10 --- /dev/null +++ b/hgeometry/ply/src/PLY/Writer.hs @@ -0,0 +1,69 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE QuasiQuotes #-} +-------------------------------------------------------------------------------- +-- | +-- Module : Ply.Writer +-- Copyright : (C) Frank Staals +-- License : see the LICENSE file +-- Maintainer : Frank Staals +-- +-- Helper module to write PLY files; so that we can export simple 3D scenes. +-- +-------------------------------------------------------------------------------- +module PLY.Writer + ( renderOutputToFile + , renderOutput + ) where + +import Control.Lens +import qualified Data.ByteString.Lazy.Char8 as Char8 +import qualified Data.Foldable as F +import Data.Foldable1 +import HGeometry.Ext +import HGeometry.Point +import HGeometry.Triangle +import qualified System.File.OsPath as File +import System.OsPath + +-------------------------------------------------------------------------------- + +-- | Writes the the points and triangles to a file in PLY format. +renderOutputToFile :: (Foldable1 f, Point_ point 3 r, Show r) + => OsPath -> f (point :+ Int) + -> [Triangle (point :+ Int)] + -> IO () +renderOutputToFile fp pts ts = File.writeFile fp $ renderOutput pts ts + +-- | Generates the content of the PLY file for the given non-empty list of points and the +-- list of triangles +-- assumes points are 0 indexed. +renderOutput :: (Foldable1 f, Point_ point 3 r, Show r) + => f (point :+ Int) -> [Triangle (point :+ Int)] + -> Char8.ByteString +renderOutput (F.toList -> pts) ts = + Char8.unlines $ hdr <> map renderPt pts <> map renderTri ts + where + hdr = ["ply" + , "format ascii 1.0" + ,"element vertex " <> (showT $ length pts) + ,"property float32 x" + ,"property float32 y" + ,"property float32 z" + ,"element face " <> (showT $ length ts) + ,"property list uchar int vertex_index" + ,"end_header" + ] + +-- | Writes a Point to ply format +renderPt :: (Point_ point 3 r, Show r) => (point :+ extra) -> Char8.ByteString +renderPt (p :+ _) = let Point3 x y z = over coordinates showT $ p^.asPoint + in Char8.unwords [x,y,z] + +-- | Writes a triangle to ply format +renderTri :: (Point_ point 3 r, Show r) + => Triangle (point :+ Int) -> Char8.ByteString +renderTri (Triangle p q r) = let i a = showT $ a^.extra in Char8.unwords ["3",i p, i q, i r] + +-- | Helper to output stuff to bytestrings +showT :: Show a => a -> Char8.ByteString +showT = Char8.pack . show diff --git a/hgeometry/src-quickcheck/HGeometry/PlaneGraph/Instances.hs b/hgeometry/src-quickcheck/HGeometry/PlaneGraph/Instances.hs new file mode 100644 index 000000000..36f7711a9 --- /dev/null +++ b/hgeometry/src-quickcheck/HGeometry/PlaneGraph/Instances.hs @@ -0,0 +1,177 @@ +-------------------------------------------------------------------------------- +-- | +-- Module : HGeometry.PlaneGraph.Instances +-- Copyright : (C) Frank Staals +-- License : see the LICENSE file +-- Maintainer : Frank Staals +-- +-- Arbitrary instance for a plane graph +-- +-------------------------------------------------------------------------------- +module HGeometry.PlaneGraph.Instances + ( arbitraryPlaneGraph + , QuickCheckWorld + ) where + +import Control.Lens +import Control.Monad.State +import Control.Subcategory.Functor +import Data.Coerce +import Data.Foldable1 +import qualified Data.List.NonEmpty as NonEmpty +import Data.List.NonEmpty(NonEmpty(..)) +import qualified Data.Map as Map +import qualified Data.Map.NonEmpty as NEMap +import Data.Ord(comparing) +import Data.Proxy +import qualified Data.Set as Set +import Data.Set(Set) +import Data.Tree +import qualified Data.Vector.NonEmpty as Vector +import HGeometry.Ext +import HGeometry.Foldable.Util +import HGeometry.HyperPlane.NonVertical +import HGeometry.Instances () +import HGeometry.Plane.LowerEnvelope.Connected(MinimizationDiagram(..)) +import HGeometry.Plane.LowerEnvelope.Connected.Graph +import HGeometry.PlaneGraph +import HGeometry.Point +import HGeometry.Properties +import HGeometry.VoronoiDiagram +import HGeometry.VoronoiDiagram.ViaLowerEnvelope +import Hiraffe.AdjacencyListRep.Map +import Hiraffe.BFS.Pure +import Hiraffe.PlanarGraph(planarGraph, vertexData) +import Hiraffe.PlanarGraph.Dart(Arc(..), Dart(..), Direction(..)) +import Prelude hiding (filter) +import Test.QuickCheck hiding (Positive, Negative) +import Test.QuickCheck.Instances () +import Witherable + + +import Debug.Trace +-------------------------------------------------------------------------------- + +data QuickCheckWorld + +instance ( Arbitrary r + , Ord r, Fractional r + , Show r + ) => Arbitrary (PlaneGraph QuickCheckWorld (Point 2 r) () ()) where + arbitrary = arbitraryPlaneGraph Proxy + +-- general strategy: +-- 1) generate a bunch of random points uniformly at random. +-- 2) Construct the Voronoi diagram +-- 3) turn it into a triangulated graph (on its bounded vertices) +-- 4) sample a subset of its edges +-- 5) return the largest connected component + + +-- | Generate an arbitrary Plane Graph +arbitraryPlaneGraph :: forall proxy s r. + ( Ord r, Fractional r, Arbitrary r + , Show r + ) + => proxy s -> Gen (PlaneGraph s (Point 2 r) () ()) +arbitraryPlaneGraph proxy = do + n <- scale (*2) arbitrary + (pts :: NonEmpty (Point 2 r)) <- genNDistinct (max 10 n) arbitrary + -- need at least a few vertices so that we generate at least a triangle in the planar graph + case voronoiDiagram pts of + AllColinear _ -> arbitraryPlaneGraph proxy -- retry + ConnectedVD vd -> do gr <- markWitherableEdges (toTriangulatedPlaneGraph' . asMD $ vd) + case traverseNeighbourOrder NEMap.nonEmptyMap $ largestComponent gr of + Nothing -> arbitraryPlaneGraph proxy -- retry + Just gr' -> pure $ toPlaneGraph proxy gr' + +-------------------------------------------------------------------------------- + + -- | Returns the largest connected component in the graph; i.e. it shrinks +-- the graph to contain only the vertices/edges in this connected component. +largestComponent :: (Ord i, Witherable f) + => GGraph f i v (Bool, e) -> GGraph f i v e +largestComponent gr = witherGraphTo tr gr + where + tr = maximumBy (comparing length) $ bff gr + + +-- | Convert a Vornooi diagram into a minimization diagram again +asMD :: ( Point_ point 2 r + , Num r, Ord r + ) => VoronoiDiagram' point -> MinimizationDiagram r (Plane r) +asMD = MinimizationDiagram . Map.mapKeys pointToPlane . coerce + + +-- | select a random subset of edges. I.e. it marks the edges we want to retain. +markWitherableEdges :: Ord i => GGraph f i v e -> Gen (GGraph f i v (Bool,e)) +markWitherableEdges gr = gr&edges %%~ \x -> (,x) <$> arbitrary' + where + arbitrary' = frequency [ (19, pure True) + , (1, pure False) + ] + +-- | Retain only the selected subset of the vertices, and the edges marked +witherGraphTo :: ( Foldable1 g, Witherable f, Ord i + ) => g i -> GGraph f i v (Bool, e) -> GGraph f i v e +witherGraphTo vs (Graph gr) = Graph $ fmap removeEdges m + where + -- retain only the vertices from vs + m = foldMap1 (\u -> NEMap.singleton u (gr NEMap.! u)) vs + -- remove edges to other components + removeEdges (VertexData x ns no) = let ns' = Map.foldMapWithKey p' ns + in VertexData x ns' (filter (`Map.member` ns') no) + -- test if we should retain the edge; i.e. if the edge is marked and the + -- other endpoint exists in our component + p' i (retain,x) | retain && NEMap.member i m = Map.singleton i x + | otherwise = Map.empty + + +-- | Given a connected plane graph in adjacency list format; convert it into an actual +-- PlaneGraph. +-- +-- \(O(n\log n)\) +toPlaneGraph :: (Ord r, Foldable1 f) + => proxy s + -> GGraph f (Point 2 r) v e -> PlaneGraph s (Point 2 r) () () +toPlaneGraph _ (Graph m) = PlaneGraph $ (planarGraph theDarts)&vertexData .~ vtxData + where + vtxData = Vector.fromNonEmptyN1 (length m) (NEMap.keys m) + + -- a non-empty list of vertices, with for each vertex the darts in order around the vertex + theDarts = evalState (sequence' theDarts') (0 :+ Map.empty) + theDarts' = toNonEmpty $ NEMap.mapWithKey toIncidentDarts m + -- turn the outgoing edges of u into darts + toIncidentDarts u (VertexData _ _ neighOrder) = + (\v -> (toDart u v, ())) <$> toNonEmpty neighOrder + -- create the dart corresponding to vertices u and v + + toDart u v | u <= v = flip Dart Positive <$> arc u v + | otherwise = flip Dart Negative <$> arc v u + + arc u v = gets (arcOf (u,v)) >>= \case + Just a -> pure a + Nothing -> do a <- nextArc + modify $ insertArc (u,v) a + pure a + + arcOf x = Map.lookup x . view extra + insertArc k v = over extra $ Map.insert k v + + nextArc = do i <- gets (view core) + modify $ over core (+1) + pure $ Arc i + + +sequence' :: Applicative m => NonEmpty (NonEmpty (m a, b)) -> m (NonEmpty (NonEmpty (a,b))) +sequence' = traverse $ traverse (\(fa,b) -> (,b) <$> fa) + +--------------------------------------------------------------------------------0 + +-- | pre: n >= 1 +genNDistinct :: Eq a => Int -> Gen a -> Gen (NonEmpty a) +genNDistinct n gen = NonEmpty.fromList <$> go [] n + where + go acc 0 = pure acc + go acc n' = do x <- gen `suchThat` (\x' -> all (/= x') acc) + go (x:acc) (n'-1) diff --git a/hgeometry/src/HGeometry/ConvexHull/R3/Naive/Dual.hs b/hgeometry/src/HGeometry/ConvexHull/R3/Naive/Dual.hs index 9402781e0..520378c2d 100644 --- a/hgeometry/src/HGeometry/ConvexHull/R3/Naive/Dual.hs +++ b/hgeometry/src/HGeometry/ConvexHull/R3/Naive/Dual.hs @@ -7,25 +7,29 @@ module HGeometry.ConvexHull.R3.Naive.Dual , facets ) where -import Control.Lens -import Data.Foldable (toList) -import Data.Foldable1 -import Data.List (find) -import Data.List.NonEmpty (NonEmpty(..)) -import Data.Maybe (isNothing) -import HGeometry.Combinatorial.Util -import HGeometry.Duality -import HGeometry.Ext -import HGeometry.HalfSpace -import HGeometry.HyperPlane -import HGeometry.HyperPlane.NonVertical -import HGeometry.Intersection (intersects) -import HGeometry.Plane.LowerEnvelope -import HGeometry.Plane.LowerEnvelope.Naive -import HGeometry.Point -import HGeometry.Properties -import HGeometry.Triangle -import HGeometry.Vector +import Control.Lens +import Data.Foldable (toList) +import Data.Foldable1 +import Data.List (find) +import Data.List.NonEmpty (NonEmpty(..)) +import qualified Data.List.NonEmpty as NonEmpty +import Data.Map (Map) +import qualified Data.Map as Map +import Data.Maybe (isNothing) +import HGeometry.Combinatorial.Util +import HGeometry.Duality +import HGeometry.Ext +import HGeometry.HalfSpace +import HGeometry.HyperPlane +import HGeometry.HyperPlane.NonVertical +import HGeometry.Intersection (intersects) +import HGeometry.Plane.LowerEnvelope +import HGeometry.Plane.LowerEnvelope.Connected.MonoidalMap (mapWithKeyMerge) +import HGeometry.Plane.LowerEnvelope.Naive +import HGeometry.Point +import HGeometry.Properties +import HGeometry.Triangle +import HGeometry.Vector -------------------------------------------------------------------------------- @@ -46,13 +50,25 @@ upperHull :: ( Point_ point 3 r upperHull = lowerEnvelope . fmap (\p -> dualHyperPlane p :+ p) -type Facet point = [point] +type Facet point = NonEmpty point -- | Outputs the facets of the upper hull. -facets :: UpperHull point -> [Facet point] +facets :: (Ord (NumType point)) => UpperHull point -> [Facet point] facets = \case - ParallelStrips _ -> error "facets: parallel strips; no bounded facets" - ConnectedEnvelope env -> env^..boundedVertices.traverse.to toFacet - where - toFacet :: BoundedVertexF f (plane :+ point) -> Facet point - toFacet v = (^.extra) <$> v^.definers.to toList + ParallelStrips _ -> [] -- error "facets: parallel strips; no bounded facets" + ConnectedEnvelope env -> toFacet <$> Map.elems theVertices + where + theVertices = mapWithKeyMerge (\h reg -> Map.fromList [ (v,NonEmpty.singleton h) + | v <- verticesOf reg ] + ) (asMap env) + verticesOf = \case + Bounded vertices -> vertices + Unbounded _ vertices _ -> NonEmpty.toList vertices + + toFacet hv = (^.extra) <$> hv + -- We want all vertices v of the lower envelope; let H_v denote all planes that + -- intersect in v. + -- + -- In the primal v is actually a plane; the facet is then defined by all points to the + -- planes in H_v. So, we want to compute H_v, and then take the convex hull of those + -- points (in the plane dual to v) diff --git a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected.hs b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected.hs index 066469f8d..6c43811b6 100644 --- a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected.hs +++ b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected.hs @@ -9,7 +9,7 @@ -- -------------------------------------------------------------------------------- module HGeometry.Plane.LowerEnvelope.Connected - ( MinimizationDiagram + ( MinimizationDiagram(..) , asMap , Region(..) , CircularList diff --git a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/Graph.hs b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/Graph.hs index 52dc05062..381ef6795 100644 --- a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/Graph.hs +++ b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/Graph.hs @@ -9,48 +9,74 @@ -- -------------------------------------------------------------------------------- module HGeometry.Plane.LowerEnvelope.Connected.Graph - ( PlaneGraph + ( PlaneGraph' , E(..) - , toPlaneGraph + , toTriangulatedPlaneGraph' + , toPlaneGraph' ) where import Data.List.NonEmpty (NonEmpty(..)) import Data.Map (Map) import qualified Data.Map as Map +import qualified Data.Map.NonEmpty as NEMap import Data.Semigroup (First(..)) import HGeometry.HyperPlane.Class import HGeometry.HyperPlane.NonVertical import HGeometry.Plane.LowerEnvelope.Connected.MonoidalMap import HGeometry.Plane.LowerEnvelope.Connected.Type +import HGeometry.PlaneGraph.Type (E(..)) import HGeometry.Point import HGeometry.Vector - +import Hiraffe.AdjacencyListRep.Map -------------------------------------------------------------------------------- -- | A Plane graph storing vertices of type v that are identified by keys of type k, and -- some ordered sequence of edges (which are ordered using e). -type PlaneGraph k v e = Map k (Map e k, v) +type PlaneGraph' k v e = GGraph (Map e) k v e + + + +-- | A Plane graph storing vertices of type v that are identified by keys of type k, and +-- some ordered sequence of edges (which are ordered using e). +type PlaneGraphMap k v e = Map k (Map e k, v) +-- | Produce a triangulated plane graph on the bounded vertices. every vertex is +-- represented by its point, it stores a list of its outgoing edges, and some data. +toTriangulatedPlaneGraph' :: (Plane_ plane r, Num r, Ord r) + => MinimizationDiagram r plane + -> PlaneGraph' (Point 2 r) (First r) (E r) +toTriangulatedPlaneGraph' = toPlaneGraph'' . toTriangulatedPlaneGraphMap -newtype E r = E (Vector 2 r) - deriving newtype (Show) +-- | Produce a plane graph on the bounded vertices. every vertex is represented by its +-- point, it stores a list of its outgoing edges, and some data. +toPlaneGraph' :: (Plane_ plane r, Num r, Ord r) + => MinimizationDiagram r plane -> PlaneGraph' (Point 2 r) (First r) (E r) +toPlaneGraph' = toPlaneGraph'' . toPlaneGraphMap + +-- | Helper to produce the graphs +toPlaneGraph'' :: (Ord r) + => PlaneGraphMap (Point 2 r) (First r) (E r) + -> PlaneGraph' (Point 2 r) (First r) (E r) +toPlaneGraph'' = Graph . NEMap.unsafeFromMap + . fmap (\(neighOrder, x) -> VertexData x (mkNeighMap neighOrder) neighOrder) + where + mkNeighMap = Map.foldMapWithKey (\e i -> Map.singleton i e) -instance (Ord r, Num r) => Eq (E r) where - a == b = a `compare` b == EQ -instance (Ord r, Num r) => Ord (E r) where - (E v) `compare` (E u) = ccwCmpAroundWith (Vector2 0 1) (origin :: Point 2 r) (Point v) (Point u) +-------------------------------------------------------------------------------- -- | Produce a triangulated plane graph on the bounded vertices. every vertex is -- represented by its point, it stores a list of its outgoing edges, and some data. -toPlaneGraph :: (Plane_ plane r, Num r, Ord r) - => MinimizationDiagram r plane -> PlaneGraph (Point 2 r) (First r) (E r) -toPlaneGraph = mapWithKeyMerge toTriangulatedGr . asMap +toTriangulatedPlaneGraphMap :: (Plane_ plane r, Num r, Ord r) + => MinimizationDiagram r plane + -> PlaneGraphMap (Point 2 r) (First r) (E r) +toTriangulatedPlaneGraphMap = mapWithKeyMerge toTriangulatedGr . asMap +-- | Helper function to construct a triangulated plane graph toTriangulatedGr :: (Plane_ plane r, Num r, Ord r) => plane -> Region r (Point 2 r) - -> PlaneGraph (Point 2 r) (First r) (E r) + -> PlaneGraphMap (Point 2 r) (First r) (E r) toTriangulatedGr h = Map.mapWithKey (\v adjs -> (adjs, First $ evalAt v h)) . \case Bounded vertices -> case vertices of (v0:v1:vs) -> triangulate v0 v1 vs @@ -69,3 +95,30 @@ toTriangulatedGr h = Map.mapWithKey (\v adjs -> (adjs, First $ evalAt v h)) . \c , (w, Map.fromList [ edge w u, edge w v]) ] edge u v = ((E $ v .-. u), v) + +-------------------------------------------------------------------------------- + +-- | Produce a plane graph on the bounded vertices. every vertex is represented by its +-- point, it stores a list of its outgoing edges (to other bounded vertices), and some +-- data. +toPlaneGraphMap :: (Plane_ plane r, Num r, Ord r) + => MinimizationDiagram r plane -> PlaneGraphMap (Point 2 r) (First r) (E r) +toPlaneGraphMap = mapWithKeyMerge toGr . asMap + + +-- | Helper function to construct a plane graph +toGr :: (Plane_ plane r, Num r, Ord r) + => plane -> Region r (Point 2 r) + -> PlaneGraphMap (Point 2 r) (First r) (E r) +toGr h = Map.mapWithKey (\v adjs -> (adjs, First $ evalAt v h)) . \case + Bounded vertices -> case vertices of + (_:v1:vs) -> Map.unionsWith (<>) $ zipWith mkEdge vertices (v1:vs) + _ -> error "toGR: absurd, <2 vertices" + Unbounded _ vertices _ -> case vertices of + _ :| [] -> Map.empty + u :| vs -> Map.unionsWith (<>) $ zipWith mkEdge (u:vs) vs + where + mkEdge u v = Map.fromList [ (u, (uncurry Map.singleton $ edge u v)) + , (v, (uncurry Map.singleton $ edge v u)) + ] + edge u v = ((E $ v .-. u), v) diff --git a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/MonoidalMap.hs b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/MonoidalMap.hs index 0e0118c23..a6fca3202 100644 --- a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/MonoidalMap.hs +++ b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/MonoidalMap.hs @@ -1,12 +1,17 @@ module HGeometry.Plane.LowerEnvelope.Connected.MonoidalMap - ( MonoidalMap, getMap + ( MonoidalMap(..) , unionsWithKey , mapWithKeyMerge + + , MonoidalNEMap(..) ) where import qualified Data.Foldable as F +import Data.Foldable1 import Data.Map (Map) import qualified Data.Map as Map +import Data.Map.NonEmpty (NEMap) +import qualified Data.Map.NonEmpty as NEMap -------------------------------------------------------------------------------- -- * Operations on Maps @@ -23,10 +28,22 @@ mapWithKeyMerge f = getMap . Map.foldMapWithKey (\k v -> MonoidalMap $ f k v) -- | A Map in which we combine conflicting elements by using their semigroup operation -- rather than picking the left value (as is done in the default Data.Map) newtype MonoidalMap k v = MonoidalMap { getMap :: Map k v } - deriving (Show) + deriving stock (Show) + deriving newtype (Functor,Foldable) instance (Ord k, Semigroup v) => Semigroup (MonoidalMap k v) where (MonoidalMap ma) <> (MonoidalMap mb) = MonoidalMap $ Map.unionWith (<>) ma mb instance (Ord k, Semigroup v) => Monoid (MonoidalMap k v) where mempty = MonoidalMap mempty + +-------------------------------------------------------------------------------- + +-- | A NonEmpty Map in which we combine conflicting elements by using their semigroup +-- operation rather than picking the left value (as is done in the default Data.Map) +newtype MonoidalNEMap k v = MonoidalNEMap { getNEMap :: NEMap k v } + deriving stock (Show) + deriving newtype (Functor,Foldable,Foldable1) + +instance (Ord k, Semigroup v) => Semigroup (MonoidalNEMap k v) where + (MonoidalNEMap ma) <> (MonoidalNEMap mb) = MonoidalNEMap $ NEMap.unionWith (<>) ma mb diff --git a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/Separator.hs b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/Separator.hs index 999f15d7c..148b4cbbf 100644 --- a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/Separator.hs +++ b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/Separator.hs @@ -13,6 +13,7 @@ module HGeometry.Plane.LowerEnvelope.Connected.Separator ( planarSeparator ) where +import Control.Lens ((^..),(^?), asIndex) import qualified Data.Foldable as F import Data.Kind (Type) import qualified Data.List as List @@ -20,47 +21,17 @@ import Data.List.NonEmpty (NonEmpty(..)) import qualified Data.List.NonEmpty as NonEmpty import Data.Map (Map) import qualified Data.Map as Map -import Data.Maybe (fromMaybe) +import Data.Maybe (fromMaybe, isJust) +import Data.Ord (comparing) import Data.Set (Set) import qualified Data.Set as Set import Data.Tree +import HGeometry.Foldable.Util () import HGeometry.Plane.LowerEnvelope.Connected.Graph import HGeometry.Vector - --------------------------------------------------------------------------------- --- * BFS - --- | Computes a breath first forest --- (O(n \log n)) -bff :: Ord k => PlaneGraph k v e -> [Tree k] -bff gr = go (Map.keysSet gr) - where - go remaining = case Set.minView remaining of - Nothing -> [] - Just (v, remaining') -> let (remaining'', t) = bfs gr v remaining' - in t : go remaining'' - --- | Turn the map into a tree. -toTree :: Ord k => Map k [k] -> k -> Tree k -toTree m = go - where - go s = Node s $ map go (fromMaybe [] $ Map.lookup s m) - --- | BFS from the given starting vertex, and the set of still-to-visit vertices. --- returns the remaining still-to-visit vertices and the tree. -bfs :: Ord k => PlaneGraph k v e -> k -> Set k -> (Set k, Tree k) -bfs gr s = fmap (flip toTree s) . bfs' [s] - where - bfs' lvl remaining = foldr visit (remaining, Map.empty) lvl - - visit v (remaining, m) = let chs = filter (flip Set.member remaining) $ neighs v - in (foldr Set.delete remaining chs, Map.insert v chs m) - neighs v = maybe [] (Map.elems . fst) $ Map.lookup v gr - --------------------------------------------------------------------------------- - - - +import Hiraffe.AdjacencyListRep.Map +import Hiraffe.BFS.Pure +import Hiraffe.Graph.Class -------------------------------------------------------------------------------- @@ -82,21 +53,30 @@ sqrt' = floor . sqrt . fromIntegral type Separator k = ([k],Vector 2 [k]) + +-- | Extracts the maximum from a non-empty list using the given ordering +viewMaximumBy :: (a -> a -> Ordering) -> NonEmpty a -> (a, [a]) +viewMaximumBy cmp (x0:|xs) = foldr (\x (m,rest) -> case x `cmp` m of + LT -> (m, x:rest) + EQ -> (m, x:rest) + GT -> (x, m:rest) + + ) (x0, []) xs + -- | Returns a pair (separator, Vector2 verticesSubGraphA verticesSubGraphB) -- so that -- -- 1) there are no edges connecting subGraph A and subgraph B, -- 2) the size of the separator is at most sqrt(n). -- 3) the vertex sets of A and B have weight at most 2/3 the total weight -planarSeparator :: Ord k => PlaneGraph k v e -> Separator k -planarSeparator gr = case trees of - [] -> ([],Vector2 [] []) - ((tr,m):rest) +planarSeparator :: Ord k => PlaneGraph' k v e -> Separator k +planarSeparator gr = case viewMaximumBy (comparing snd) trees of + ((tr,m),rest) | m <= twoThirds -> groupComponents | otherwise -> planarSeparator' tr m -- we should also add the remaining vertices where - trees = List.sortOn snd . map (\t -> (t, length t)) $ bff gr - n = sum $ map snd trees + trees = (\t -> (t, length t)) <$> bff gr + n = sum $ fmap snd trees half = n `div` 2 twoThirds = 2 * (n `div` 3) @@ -129,7 +109,7 @@ planarSeparator gr = case trees of -- | contracts the plane graph so that we get a spanning tree of diameter at most sqrt(n). -contract :: PlaneGraph k v e -> Tree k -> (PlaneGraph k v e, Tree k) +contract :: PlaneGraph' k v e -> Tree k -> (PlaneGraph' k v e, Tree k) contract = undefined trim _ _ tr = tr @@ -137,7 +117,7 @@ trim _ _ tr = tr -- | Given a spanning tree of the graph that has diameter r, compute -- a separator of size at most 2r+1 -planarSeparatorTree :: Ord k => PlaneGraph k v e -> Tree k -> Separator k +planarSeparatorTree :: Ord k => PlaneGraph' k v e -> Tree k -> Separator k planarSeparatorTree gr tr = (sep, foldMap F.toList <$> trees) -- FIXME: continue searching where @@ -283,14 +263,15 @@ type EndPoint a = (a, [Tree a], [Tree a]) -- | Split the leaf of the path splitLeaf :: Ord k - => PlaneGraph k v e + => PlaneGraph' k v e -> (k,k) -> SplitTree k (Tree k) -> SplitTree k (EndPoint k) splitLeaf gr (v',w') = fmap $ \(Node u chs) -> split u chs (if u == v' then w' else v') where split v chs w = case List.break (hasEdge v) chs of (before, _:after) -> (v, before, after) _ -> error "splitLeaf: absurd. edge not found!?" - hasEdge v w = maybe False (F.elem (root w) . fst) $ Map.lookup v gr + hasEdge v w = isJust $ gr^?edgeAt (v, root w) + -- maybe False (F.elem (root w) . fst) $ gr^?vertexAt v -- note: this is linear in the number of edges out of v -- | Turn the split tree into a separator, and the trees inside the cycle, and outside the @@ -373,8 +354,11 @@ annotate :: Tree k -> Tree (Weighted Int k) annotate (Node v chs) = let chs' = map annotate chs in Node (Weighted (1 + costs chs') v) chs' -graphEdges :: Ord k => PlaneGraph k v e -> Set (k,k) -graphEdges = Map.foldMapWithKey (\u (es,_) -> Set.fromList [ (u,v) | v <- Map.elems es, u <= v]) +graphEdges :: Ord k => PlaneGraph' k v e -> Set (k,k) +graphEdges gr = Set.fromList [ (u,v) | (u,v) <- gr^..edges.asIndex, u <= v ] + + + -- Map.foldMapWithKey (\u (es,_) -> Set.fromList [ (u,v) | v <- Map.elems es, u <= v]) treeEdges :: Ord k => Tree k -> Set (k,k) treeEdges (Node u chs) = Set.fromList [ orient (u,v) | Node v _ <- chs ] diff --git a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/Type.hs b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/Type.hs index 611a91994..c618042f9 100644 --- a/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/Type.hs +++ b/hgeometry/src/HGeometry/Plane/LowerEnvelope/Connected/Type.hs @@ -101,3 +101,43 @@ type CircularList a = [a] -------------------------------------------------------------------------------- + +-- data MDVertexIx = UnboundedVertexIx +-- | BoundedVertexIx {-# UNPACK #-}!Int +-- deriving (Show,Read,Eq,Ord) + +-- data MDVertex r plane = UnboundedVertex +-- | BoundedVertex { _location :: Point 2 r +-- , _definers :: Vector 3 plane +-- } +-- deriving (Show,Eq,Ord) + + +-- instance HasVertices' (MinimizationDiagram r plane) where +-- -- | invariant: vertexIx == UnboundedVertex <=> vertex = UnboundedVertex +-- type VertexIx (MinimizationDiagram r plane) = MDVertexIx +-- type Vertex (MinimizationDiagram r plane) = MDVertex r plane +-- vertexAt = \case +-- UnboundedVertexIx -> pure UnboundedVertex +-- BoundedVertexIx i -> + +-- instance HasVertices (MinimizationDiagram r plane) (MinimizationDiagram r plane) where +-- instance HasDarts' (MinimizationDiagram r plane) where +-- type DartIx (MinimizationDiagram r plane) = (MDVertexIx, MDVertexIx) +-- type Dart (MinimizationDiagram r plane) = (MDVertex r plane, MDVertex r plane) +-- dartAt = undefined + +-- instance HasDarts (MinimizationDiagram r plane) (MinimizationDiagram r plane) where + +-- instance DiGraph_ (MinimizationDiagram r plane) where +-- diGraphFromAdjacencyLists = undefined +-- endPoints = undefined +-- outNeighboursOf = undefined +-- twinDartOf = undefined + + +-- -- | Produce a triangulated plane graph on the bounded vertices. every vertex is +-- -- represented by its point, it stores a list of its outgoing edges, and some data. +-- toGraph :: (Plane_ plane r, Num r, Ord r) +-- => MinimizationDiagram r plane -> PlaneGraph (Point 2 r) (First r) (E r) +-- toGraph = mapWithKeyMerge toTriangulatedGr . asMap diff --git a/hgeometry/src/HGeometry/Plane/LowerEnvelope/DivideAndConquer.hs b/hgeometry/src/HGeometry/Plane/LowerEnvelope/DivideAndConquer.hs index f1afe09cf..fc5a7924a 100644 --- a/hgeometry/src/HGeometry/Plane/LowerEnvelope/DivideAndConquer.hs +++ b/hgeometry/src/HGeometry/Plane/LowerEnvelope/DivideAndConquer.hs @@ -2,36 +2,38 @@ module HGeometry.Plane.LowerEnvelope.DivideAndConquer ( lowerEnvelope ) where -import Control.Lens -import Data.Foldable1 -import qualified Data.List as List -import Data.List.NonEmpty (NonEmpty(..)) -import qualified Data.List.NonEmpty as NonEmpty -import qualified Data.Map as Map -import qualified Data.Set as Set -import Data.Word -import HGeometry.Ext -import HGeometry.HyperPlane.Class -import HGeometry.HyperPlane.NonVertical -import HGeometry.Intersection -import HGeometry.Line -import HGeometry.Line.General -import HGeometry.Line.LineEQ -import qualified HGeometry.Line.LowerEnvelope as LowerEnvelope -import HGeometry.Line.PointAndVector -import HGeometry.Plane.LowerEnvelope.AdjListForm -import HGeometry.Plane.LowerEnvelope.EpsApproximation -import qualified HGeometry.Plane.LowerEnvelope.Naive as Naive -import HGeometry.Plane.LowerEnvelope.Type -import HGeometry.Plane.LowerEnvelope.VertexForm -import HGeometry.Point -import HGeometry.Properties -import HGeometry.Vector -import Witherable +-- import Control.Lens +-- import Data.Foldable1 +-- import qualified Data.List as List +-- import Data.List.NonEmpty (NonEmpty(..)) +-- import qualified Data.List.NonEmpty as NonEmpty +-- import qualified Data.Map as Map +-- import qualified Data.Set as Set +-- import Data.Word +-- import HGeometry.Ext +-- import HGeometry.HyperPlane.Class +-- import HGeometry.HyperPlane.NonVertical +-- import HGeometry.Intersection +-- import HGeometry.Line +-- import HGeometry.Line.General +-- import HGeometry.Line.LineEQ +-- import qualified HGeometry.Line.LowerEnvelope as LowerEnvelope +-- import HGeometry.Line.PointAndVector +-- import HGeometry.Plane.LowerEnvelope.AdjListForm +-- import HGeometry.Plane.LowerEnvelope.EpsApproximation +-- import qualified HGeometry.Plane.LowerEnvelope.Naive as Naive +-- import HGeometry.Plane.LowerEnvelope.Type +-- import HGeometry.Plane.LowerEnvelope.VertexForm +-- import HGeometry.Point +-- import HGeometry.Properties +-- import HGeometry.Vector +-- import Witherable -------------------------------------------------------------------------------- +lowerEnvelope = undefined +{- -- | below this value we use the naive algorithm. nZero :: Int @@ -274,3 +276,6 @@ planarSeparator gr = case List.break (\lvl -> accumSize lvl < half) lvls of -- planarSeparators :: PlanarGraph_ planarGraph -- => planarGraph -> Tree + + +-} diff --git a/hgeometry/src/HGeometry/Plane/LowerEnvelope/VertexForm.hs b/hgeometry/src/HGeometry/Plane/LowerEnvelope/VertexForm.hs deleted file mode 100644 index 5298c137f..000000000 --- a/hgeometry/src/HGeometry/Plane/LowerEnvelope/VertexForm.hs +++ /dev/null @@ -1,130 +0,0 @@ -{-# LANGUAGE UndecidableInstances #-} --------------------------------------------------------------------------------- --- | --- Module : HGeometry.Plane.LowerEnvelope.VertexForm --- Copyright : (C) Frank Staals --- License : see the LICENSE file --- Maintainer : Frank Staals --- --- A Representation of the Lower envelope of planes in vertex form; --- i.e. storing only the vertices. --- --------------------------------------------------------------------------------- -module HGeometry.Plane.LowerEnvelope.VertexForm - ( VertexForm(VertexForm) - , hasVertices, vertices' - , singleton - , LEVertex, pattern LEVertex, Definers - , BoundedVertexF(Vertex) - , location, definers, location2 - - , intersectionLine - , intersectionPoint - ) where - -import Control.Lens -import qualified Data.Map as Map -import qualified Data.Set as Set -import HGeometry.Combinatorial.Util -import HGeometry.HyperPlane.NonVertical -import HGeometry.Intersection -import HGeometry.Line -import HGeometry.Line.General -import HGeometry.Plane.LowerEnvelope.Type -import HGeometry.Point -import HGeometry.Properties -import Hiraffe.Graph - --------------------------------------------------------------------------------- - --- | Vertices of a lower envelope in vertex form. -type LEVertex = BoundedVertexF (Const ()) - --- | Convenience Constructor for computing vertices of the lower --- envelope that are in Vertex form. -pattern LEVertex :: Point 3 (NumType plane) -> Definers plane -> LEVertex plane -pattern LEVertex v defs = Vertex v defs (Const ()) -{-# COMPLETE LEVertex #-} - --- | The definers of a vertex. -type Definers plane = Set.Set plane - --- | The lower envelope in vertex form -newtype VertexForm plane = - VertexForm (Map.Map (Point 3 (NumType plane)) - (Definers plane) - ) - --- | Iso to access the underlying Map. -_VertexFormMap :: Iso (VertexForm plane) (VertexForm plane') - (Map.Map (Point 3 (NumType plane)) (Definers plane)) - (Map.Map (Point 3 (NumType plane')) (Definers plane')) -_VertexFormMap = coerced - --- | Test if there are vertices -hasVertices :: VertexForm plane -> Bool -hasVertices = not . Map.null . view _VertexFormMap - -deriving instance ( Show plane, Show (NumType plane) - ) => Show (VertexForm plane) -deriving instance ( Eq plane - , Eq (NumType plane) - ) => Eq (VertexForm plane) - --- | Computes a lower envelope consisting of a single vertex. -singleton :: LEVertex plane -> VertexForm plane -singleton v = VertexForm $ Map.singleton (v^.location) (v^.definers) - -instance (Ord (NumType plane), Ord plane) => Semigroup (VertexForm plane) where - (VertexForm m) <> (VertexForm m') = VertexForm $ Map.unionWith (<>) m m' -instance (Ord (NumType plane), Ord plane) => Monoid (VertexForm plane) where - mempty = VertexForm mempty - -instance Ord (NumType plane) => HasVertices' (VertexForm plane) where - type Vertex (VertexForm plane) = Definers plane - type VertexIx (VertexForm plane) = Point 3 (NumType plane) - vertexAt i = _VertexFormMap . iix i - numVertices (VertexForm m) = Map.size m - --- instance Ord (NumType plane) => HasVertices (VertexForm plane) (VertexForm plane) where --- vertices = _VertexFormMap . itraversed - --- | Traversal (rather than a traversal1) of the vertices. -vertices' :: IndexedTraversal' (VertexIx (VertexForm plane)) - (VertexForm plane) (Vertex (VertexForm plane)) -vertices' = _VertexFormMap . itraversed - --------------------------------------------------------------------------------- - - --- | Given two planes, computes the line in which they intersect. -intersectionLine :: (Plane_ plane r, Fractional r, Eq r) - => plane -> plane -> Maybe (VerticalOrLineEQ r) -intersectionLine (Plane_ a1 b1 c1) (Plane_ a2 b2 c2) - | b1 /= b2 = Just $ NonVertical $ LineEQ ((a2 - a1) / diffB) ((c2 - c1) / diffB) - -- the two planes intersect in some normal line - | a1 /= a2 = Just $ VerticalLineThrough ((c2 -c1) / (a1 - a2)) - -- the planes intersect in a vertical line - | otherwise = Nothing - -- the planes don't intersect at all - where - diffB = b1 - b2 - --- | Computes there the three planes intersect -intersectionPoint :: ( Plane_ plane r, Ord r, Fractional r) - => Three plane -> Maybe (Point 3 r) -intersectionPoint (Three h1@(Plane_ a1 b1 c1) h2 h3) = - do l12 <- intersectionLine h1 h2 - l13 <- intersectionLine h1 h3 - case (l12,l13) of - (VerticalLineThrough _x12, VerticalLineThrough _x13) -> Nothing - -- if the xes are the same they would be the same plane even - (VerticalLineThrough x, NonVertical l) -> vertNonVertIntersect x l - (NonVertical l, VerticalLineThrough x) -> vertNonVertIntersect x l - (NonVertical l, NonVertical m) -> l `intersect` m >>= \case - Line_x_Line_Point (Point2 x y) -> Just $ Point3 x y (a1 * x + b1* y + c1) - Line_x_Line_Line _ -> Nothing - where - vertNonVertIntersect x l = let y = evalAt' x l - z = a1 * x + b1* y + c1 - in Just $ Point3 x y z diff --git a/hgeometry/src/HGeometry/PlaneGraph.hs b/hgeometry/src/HGeometry/PlaneGraph.hs index 9aa60d9c3..81b5fc0a0 100644 --- a/hgeometry/src/HGeometry/PlaneGraph.hs +++ b/hgeometry/src/HGeometry/PlaneGraph.hs @@ -13,6 +13,8 @@ module HGeometry.PlaneGraph , module HGeometry.PlaneGraph.Class , module Hiraffe.PlanarGraph.Class , PlaneGraph(..) + , fromAdjacencyRep + , fromConnectedSegments ) where import HGeometry.PlaneGraph.Class diff --git a/hgeometry/src/HGeometry/PlaneGraph/Type.hs b/hgeometry/src/HGeometry/PlaneGraph/Type.hs index 33e861c69..40b817032 100644 --- a/hgeometry/src/HGeometry/PlaneGraph/Type.hs +++ b/hgeometry/src/HGeometry/PlaneGraph/Type.hs @@ -1,7 +1,7 @@ {-# LANGUAGE UndecidableInstances #-} -------------------------------------------------------------------------------- -- | --- Module : HGeometryPlaneGraph.Type +-- Module : HGeometry.PlaneGraph.Type -- Copyright : (C) Frank Staals -- License : see the LICENSE file -- Maintainer : Frank Staals @@ -13,24 +13,41 @@ -------------------------------------------------------------------------------- module HGeometry.PlaneGraph.Type ( PlaneGraph(..) + , fromAdjacencyRep + , fromConnectedSegments -- , VertexData(VertexData), location + + , E(..) ) where import Control.Lens hiding (holes, holesOf, (.=)) import Data.Coerce import Data.Foldable1 +import Data.Foldable1.WithIndex +import qualified Data.List.NonEmpty as NonEmpty import qualified Data.Map as Map +import qualified Data.Map.NonEmpty as NEMap import qualified Data.Vector.NonEmpty as Vector import Data.YAML import GHC.Generics (Generic) import HGeometry.Box import HGeometry.Foldable.Sort (sortBy ) +import HGeometry.LineSegment +import HGeometry.Plane.LowerEnvelope.Connected.MonoidalMap import HGeometry.PlaneGraph.Class import HGeometry.Point import HGeometry.Properties import HGeometry.Transformation -import Hiraffe.PlanarGraph +import HGeometry.Vector +import Hiraffe.AdjacencyListRep.Map +import Hiraffe.Graph.Class +import Hiraffe.PlanarGraph ( PlanarGraph, World(..) + , DartId, VertexId, FaceId + ) import qualified Hiraffe.PlanarGraph as PG +import Hiraffe.PlanarGraph.Class +import qualified Hiraffe.PlanarGraph.Dart as Dart + -------------------------------------------------------------------------------- -- * The PlaneGraph type @@ -174,3 +191,69 @@ instance ( Point_ v 2 r -- boundingBox = boundingBoxList' . F.toList . fmap (^._2.location) . vertices + + +-------------------------------------------------------------------------------- + +-- | Constructs a connected plane graph +-- +-- pre: The segments form a single connected component +-- No two segments partially overlap. +-- +-- running time: \(O(n\log n)\) +fromConnectedSegments :: ( Foldable1 f, Ord r, Num r + , LineSegment_ lineSegment point + , Point_ point 2 r + ) + => f lineSegment + -> PlaneGraph s (NonEmpty.NonEmpty point) lineSegment () +fromConnectedSegments segs = PlaneGraph $ + (PG.planarGraph theDarts)&PG.vertexData .~ vtxData + where + -- to get the darts we simply convert the NEMap (_, NEMap _ (dart, seg)) into + -- a NonEmpty (NonEmpty (dart, seg)) + theDarts = toNonEmpty . snd <$> verts + vtxData = Vector.fromNonEmpty $ fst <$> verts + + -- Collects all edges per vertex + verts = toNonEmpty . ifoldMap1 f $ toNonEmpty segs + + -- Creates two vertices with one edge each ; combines them into a single Map + f i seg = let u = seg^.start + v = seg^.end + d = Dart.Dart (Dart.Arc i) Dart.Positive + in singleton (u^.asPoint) (vtx (d ,seg) u v) + <> singleton (v^.asPoint) (vtx (Dart.twin d,seg) v u) + + singleton k v = MonoidalNEMap $ NEMap.singleton k v + +-- | Helper type to represent the vertex data of a vertex. The NEMap +-- represents the edges ordered cyclically around the vertex +type VtxData v r e = (v, NEMap.NEMap (E r) e) + +-- | Creates the vertex data +vtx :: (Point_ point 2 r, Ord r, Num r) + => e -> point -> point -> VtxData (NonEmpty.NonEmpty point) r e +vtx e p q = (NonEmpty.singleton p, NEMap.singleton (E $ q .-. p) e) + +-------------------------------------------------------------------------------- + +-- | Given a connected plane graph in adjacency list format; convert it into an actual +-- PlaneGraph. +-- +-- \(O(n\log n)\) +fromAdjacencyRep :: (Point_ vertex 2 r, Ord i, Foldable1 f) + => proxy s -> GGraph f i vertex e -> PlaneGraph s vertex e () +fromAdjacencyRep proxy = PlaneGraph . PG.fromAdjacencyRep proxy + + +-------------------------------------------------------------------------------- + +-- | Helper type to sort vectors cyclically around the origine +newtype E r = E (Vector 2 r) + deriving newtype (Show) + +instance (Ord r, Num r) => Eq (E r) where + a == b = a `compare` b == EQ +instance (Ord r, Num r) => Ord (E r) where + (E v) `compare` (E u) = ccwCmpAroundWith (Vector2 0 1) (origin :: Point 2 r) (Point v) (Point u) diff --git a/hgeometry/src/HGeometry/VoronoiDiagram/ViaLowerEnvelope.hs b/hgeometry/src/HGeometry/VoronoiDiagram/ViaLowerEnvelope.hs index 81bec6937..2c8f54333 100644 --- a/hgeometry/src/HGeometry/VoronoiDiagram/ViaLowerEnvelope.hs +++ b/hgeometry/src/HGeometry/VoronoiDiagram/ViaLowerEnvelope.hs @@ -17,6 +17,7 @@ module HGeometry.VoronoiDiagram.ViaLowerEnvelope , voronoiDiagram , voronoiVertices -- , edgeGeometries + , pointToPlane ) where import Control.Lens @@ -104,9 +105,12 @@ voronoiDiagramWith :: ( Point_ point 2 r, Functor nonEmpty, Ord point voronoiDiagramWith lowerEnv pts = case lowerEnv . fmap (\p -> pointToPlane p :+ p) $ pts of ParallelStrips strips -> AllColinear $ fmap (^.extra) strips ConnectedEnvelope env -> ConnectedVD . VoronoiDiagram . cmap (^.extra) $ env + + +-- | lifts the point to a plane; so that the lower envelope corresponds to the VD +pointToPlane :: (Point_ point 2 r, Num r) => point -> Plane r +pointToPlane = flipZ . liftPointToPlane where - -- lifts the point to a plane; so that the lower envelope corresponds to the VD - pointToPlane = flipZ . liftPointToPlane flipZ = over (hyperPlaneCoefficients.traverse) negate diff --git a/hgeometry/test-with-ipe/test/Plane/LowerEnvelopeSpec.hs b/hgeometry/test-with-ipe/test/Plane/LowerEnvelopeSpec.hs index 5e6bebfa3..3e399ca7a 100644 --- a/hgeometry/test-with-ipe/test/Plane/LowerEnvelopeSpec.hs +++ b/hgeometry/test-with-ipe/test/Plane/LowerEnvelopeSpec.hs @@ -149,10 +149,10 @@ testIpe inFp outFp = do (addStyleSheet opacitiesStyle $ singlePageFromContent out) -theEdges :: PlaneGraph (Point 2 R) h (E R) -> IpeObject' Group R -theEdges = ipeGroup . Map.foldMapWithKey (\v (adjs, _) -> - foldMap (\w -> [ iO $ defIO (ClosedLineSegment v w) - ]) adjs) +-- theEdges :: PlaneGraph' (Point 2 R) h (E R) -> IpeObject' Group R +-- theEdges = ipeGroup . Map.foldMapWithKey (\v (adjs, _) -> +-- foldMap (\w -> [ iO $ defIO (ClosedLineSegment v w) +-- ]) adjs) -- -- build a triangulated graph from the points in the input file diff --git a/hgeometry/test/ConvexHull/R3Spec.hs b/hgeometry/test/ConvexHull/R3Spec.hs new file mode 100644 index 000000000..f31a2aa57 --- /dev/null +++ b/hgeometry/test/ConvexHull/R3Spec.hs @@ -0,0 +1,28 @@ +module ConvexHull.R3Spec (spec) where + +import Data.Foldable1 +import Data.List.NonEmpty (NonEmpty(..)) +import qualified Data.List.NonEmpty as NonEmpty +import qualified HGeometry.ConvexHull.R3.Naive.Dual as Dual +import HGeometry.Instances () +import HGeometry.Number.Real.Rational +import HGeometry.Point +import HGeometry.Triangle +import Test.Hspec +import Test.Hspec.QuickCheck +import Test.QuickCheck +import Test.QuickCheck.Instances () + + +import Debug.Trace +-------------------------------------------------------------------------------- + +type R = RealNumber 5 + + +spec = describe "3D convex hull through duality tests " $ do + fit "single triangle" $ do + let tri = Triangle origin (Point3 10 0 1) (Point3 0 10 2 :: Point 3 R) + Dual.facets (traceShowWith ("hull",) $ Dual.upperHull tri) `shouldBe` [] + -- FIXME: This is actually incorrect. I think it should be the thing below: + -- Dual.facets (traceShowWith ("hull",) $ Dual.upperHull tri) `shouldBe` [toNonEmpty tri] diff --git a/hgeometry/test/Graphics/CameraSpec.hs b/hgeometry/test/Graphics/CameraSpec.hs index b78377690..aab5fc776 100644 --- a/hgeometry/test/Graphics/CameraSpec.hs +++ b/hgeometry/test/Graphics/CameraSpec.hs @@ -1,7 +1,6 @@ module Graphics.CameraSpec where import Control.Lens -import HGeometry.Ext import HGeometry.Graphics.Camera import HGeometry.Point import HGeometry.Transformation