diff --git a/benchmark/Algorithms/Geometry/PolygonTriangulation/Bench.hs b/benchmark/Algorithms/Geometry/PolygonTriangulation/Bench.hs --- a/benchmark/Algorithms/Geometry/PolygonTriangulation/Bench.hs +++ b/benchmark/Algorithms/Geometry/PolygonTriangulation/Bench.hs @@ -9,7 +9,7 @@ import Data.Ext import Test.Tasty.Bench import qualified Data.Foldable as F -import Data.Geometry.Ipe +import Ipe import Data.Geometry.LineSegment import Data.Geometry.Polygon import Data.Geometry.PlanarSubdivision diff --git a/changelog b/changelog --- a/changelog +++ b/changelog @@ -2,6 +2,24 @@ * Changelog +** 0.13 + +- Implementation of Logaritmic Method, wich allows us to transform a + static data structure into an insertion only data structure +- Moved 'intersects' from the HasIntersectionWith class into a new + class IsIntersectableWith. This allows separate (weaker) constraints + for checking *if* geometries intersect rather than computing exact + intersections. +- New BezierSpline features. +- "Zoom to fit" transformation. +- Many fixes related to PlaneGraph/PlanarSubdivison; i.e. bugs in + which order the vertices/darts where reported when traversing a + face. The polygon representing the outer boundary now is some area + inside a bounding polygon. +- Fixed a bug in the DelaunayTriangulation. +- Preliminary implementations for updating planar subdivisions + (e.g. subdividing edges). + ** 0.12 - New website: https://hgeometry.org/ diff --git a/changelog.org b/changelog.org --- a/changelog.org +++ b/changelog.org @@ -2,6 +2,24 @@ * Changelog +** 0.13 + +- Implementation of Logaritmic Method, wich allows us to transform a + static data structure into an insertion only data structure +- Moved 'intersects' from the HasIntersectionWith class into a new + class IsIntersectableWith. This allows separate (weaker) constraints + for checking *if* geometries intersect rather than computing exact + intersections. +- New BezierSpline features. +- "Zoom to fit" transformation. +- Many fixes related to PlaneGraph/PlanarSubdivison; i.e. bugs in + which order the vertices/darts where reported when traversing a + face. The polygon representing the outer boundary now is some area + inside a bounding polygon. +- Fixed a bug in the DelaunayTriangulation. +- Preliminary implementations for updating planar subdivisions + (e.g. subdividing edges). + ** 0.12 - New website: https://hgeometry.org/ diff --git a/docs/Data/Geometry/PlanarSubdivision/mySubdiv.jpg b/docs/Data/Geometry/PlanarSubdivision/mySubdiv.jpg new file mode 100644 Binary files /dev/null and b/docs/Data/Geometry/PlanarSubdivision/mySubdiv.jpg differ diff --git a/docs/Data/PlaneGraph/planegraph.png b/docs/Data/PlaneGraph/planegraph.png new file mode 100644 Binary files /dev/null and b/docs/Data/PlaneGraph/planegraph.png differ diff --git a/hgeometry.cabal b/hgeometry.cabal --- a/hgeometry.cabal +++ b/hgeometry.cabal @@ -1,5 +1,6 @@ +cabal-version: 2.4 name: hgeometry -version: 0.12.0.4 +version: 0.13 synopsis: Geometric Algorithms, Data structures, and Data types. description: HGeometry provides some basic geometry types, and geometric algorithms and @@ -8,7 +9,7 @@ asymptotic running time guarantees. Note that HGeometry is still highly experimental, don't be surprised to find bugs. homepage: https://fstaals.net/software/hgeometry -license: BSD3 +license: BSD-3-Clause license-file: LICENSE author: Frank Staals maintainer: frank@fstaals.net @@ -23,10 +24,8 @@ changelog changelog.org -Extra-doc-files: docs/Data/PlaneGraph/small.png - -- docs/**/*.png - -cabal-version: 2.0 +Extra-doc-files: docs/**/*.png + docs/**/*.jpg source-repository head type: git location: https://github.com/noinia/hgeometry @@ -99,6 +98,8 @@ Data.Geometry.PlanarSubdivision.Raw Data.Geometry.PlanarSubdivision.Basic Data.Geometry.PlanarSubdivision.Merge + Data.Geometry.PlanarSubdivision.Dynamic + Data.Geometry.PlanarSubdivision.TreeRep Data.Geometry.Arrangement Data.Geometry.Arrangement.Internal @@ -156,7 +157,7 @@ Algorithms.Geometry.Diameter.ConvexHull -- Algorithms.Geometry.Sweep - + Algorithms.Geometry.PolygonTriangulation Algorithms.Geometry.PolygonTriangulation.Types Algorithms.Geometry.PolygonTriangulation.Triangulate Algorithms.Geometry.PolygonTriangulation.MakeMonotone @@ -197,6 +198,7 @@ other-modules: Data.Geometry.Matrix.Internal + Data.Geometry.Transformation.Internal -- * Implementation Internals of Polygons Data.Geometry.Polygon.Core @@ -226,7 +228,7 @@ -- other-extensions: build-depends: base >= 4.11 && < 5 - , hgeometry-combinatorial >= 0.12.0.3 + , hgeometry-combinatorial >= 0.13 , bifunctors >= 4.1 , bytestring >= 0.10 diff --git a/src/Algorithms/Geometry/DelaunayTriangulation/Types.hs b/src/Algorithms/Geometry/DelaunayTriangulation/Types.hs --- a/src/Algorithms/Geometry/DelaunayTriangulation/Types.hs +++ b/src/Algorithms/Geometry/DelaunayTriangulation/Types.hs @@ -86,11 +86,11 @@ -- showDT :: (Show p, Show r) => Triangulation p r -> IO () -- showDT = mapM_ print . edgesAsPoints -{- HLINT ignore edgesAsPoints -} + -- | List add edges as point pairs. edgesAsPoints :: Triangulation p r -> [(Point 2 r :+ p, Point 2 r :+ p)] edgesAsPoints t = let pts = _positions t - in map (\(u,v) -> (pts V.! u, pts V.! v)) . edgesAsVertices $ t + in map (bimap (pts V.!) (pts V.!)) . edgesAsVertices $ t -- | List add edges as VertexID pairs. edgesAsVertices :: Triangulation p r -> [(VertexID,VertexID)] @@ -127,5 +127,5 @@ toPlaneGraph _ tr = PG.PlaneGraph $ g&PPG.vertexData .~ vtxData where g = PPG.fromAdjacencyLists . V.toList . V.imap f $ tr^.neighbours - f i adj = (VertexId i, VertexId <$> adj) + f i adj = (VertexId i, C.leftElements $ VertexId <$> adj) -- report in CCW order vtxData = (\(loc :+ p) -> VertexData loc p) <$> tr^.positions diff --git a/src/Algorithms/Geometry/PolyLineSimplification/ImaiIri.hs b/src/Algorithms/Geometry/PolyLineSimplification/ImaiIri.hs --- a/src/Algorithms/Geometry/PolyLineSimplification/ImaiIri.hs +++ b/src/Algorithms/Geometry/PolyLineSimplification/ImaiIri.hs @@ -25,8 +25,8 @@ import qualified Data.Vector as V import Witherable -import Data.RealNumber.Rational -type R = RealNumber 5 +-- import Data.RealNumber.Rational +-- type R = RealNumber 5 -------------------------------------------------------------------------------- @@ -122,17 +122,17 @@ -------------------------------------------------------------------------------- -tr :: Tree Int -tr = Node 0 [Node 1 [], Node 2 [Node 3 [], Node 2 [], Node 4 [Node 5 []]]] +-- tr :: Tree Int +-- tr = Node 0 [Node 1 [], Node 2 [Node 3 [], Node 2 [], Node 4 [Node 5 []]]] -poly :: PolyLine 2 Int R -poly = case fromPoints [origin :+ 0, Point2 1 1 :+ 1, Point2 2 2 :+ 2, Point2 3 3 :+ 3] of - Just p -> p +-- poly :: PolyLine 2 Int R +-- poly = case fromPoints [origin :+ 0, Point2 1 1 :+ 1, Point2 2 2 :+ 2, Point2 3 3 :+ 3] of +-- Just p -> p -test = Seq.fromList [0..5] +-- test = Seq.fromList [0..5] -myTree :: Tree Int -myTree = Node {rootLabel = 0, subForest = [Node {rootLabel = 1, subForest = []} - ,Node {rootLabel = 2, subForest = []} - ,Node {rootLabel = 3, subForest = []}] - } +-- myTree :: Tree Int +-- myTree = Node {rootLabel = 0, subForest = [Node {rootLabel = 1, subForest = []} +-- ,Node {rootLabel = 2, subForest = []} +-- ,Node {rootLabel = 3, subForest = []}] +-- } diff --git a/src/Algorithms/Geometry/PolygonTriangulation.hs b/src/Algorithms/Geometry/PolygonTriangulation.hs new file mode 100644 --- /dev/null +++ b/src/Algorithms/Geometry/PolygonTriangulation.hs @@ -0,0 +1,14 @@ +-------------------------------------------------------------------------------- +-- | +-- Module : Algorithms.Geometry.PolygonTriangulation +-- Copyright : (C) Frank Staals +-- License : see the LICENSE file +-- Maintainer : Frank Staals +-------------------------------------------------------------------------------- +module Algorithms.Geometry.PolygonTriangulation + ( triangulate + , triangulate' + , computeDiagonals + ) where + +import Algorithms.Geometry.PolygonTriangulation.Triangulate diff --git a/src/Algorithms/Geometry/PolygonTriangulation/Triangulate.hs b/src/Algorithms/Geometry/PolygonTriangulation/Triangulate.hs --- a/src/Algorithms/Geometry/PolygonTriangulation/Triangulate.hs +++ b/src/Algorithms/Geometry/PolygonTriangulation/Triangulate.hs @@ -71,8 +71,6 @@ . lefts . map (^._2.core) . filter (\mp -> mp^._2.extra == Inside) -- triangulate only the insides -- . filter (\f -> f^._1 /= outerFaceId') - . F.toList . rawFacePolygons $ monotoneP - - -- -- we alredy know we get the polgyons in *clockwise* order, so skip the - -- -- check if it is counter clockwise - -- toCounterClockWiseOrder'' = reverseOuterBoundary + . F.toList . internalFacePolygons $ monotoneP + -- FIXME: we should already get all polygons in CCW order, so no + -- need for the toClockwiseOrder' call diff --git a/src/Data/Geometry.hs b/src/Data/Geometry.hs --- a/src/Data/Geometry.hs +++ b/src/Data/Geometry.hs @@ -26,7 +26,7 @@ import Data.Geometry.LineSegment import Data.Geometry.Point import Data.Geometry.PolyLine hiding (fromPoints) -import Data.Geometry.Polygon hiding (fromPoints, maximumBy, minimumBy) +import Data.Geometry.Polygon hiding (fromPoints) import Data.Geometry.Properties import Data.Geometry.Transformation -- import Linear.Affine hiding (Point, Vector, origin) diff --git a/src/Data/Geometry/Ball.hs b/src/Data/Geometry/Ball.hs --- a/src/Data/Geometry/Ball.hs +++ b/src/Data/Geometry/Ball.hs @@ -31,6 +31,7 @@ import Linear.Matrix import Linear.V3 (V3(..)) + -------------------------------------------------------------------------------- -- * A d-dimensional ball @@ -192,16 +193,36 @@ newtype Touching p = Touching p deriving (Show,Eq,Ord,Functor,F.Foldable,T.Traversable) -- | No intersection, one touching point, or two points -type instance IntersectionOf (Line 2 r) (Circle p r) = [ NoIntersection - , Touching (Point 2 r) - , (Point 2 r, Point 2 r) - ] +type instance IntersectionOf (Line d r) (Sphere d p r) = [ NoIntersection + , Touching (Point d r) + , (Point d r, Point d r) + ] +instance {-# OVERLAPPABLE #-} (Ord r, Fractional r, Arity d) + => Line d r `HasIntersectionWith` Sphere d q r where + l `intersects` (Sphere (c :+ _) r) = let closest = pointClosestTo c l + in squaredEuclideanDist c closest <= r +instance {-# OVERLAPPING #-} (Ord r, Num r) => Line 2 r `HasIntersectionWith` Circle p r where + (Line p' v) `intersects` (Circle (c :+ _) r) = discr >= 0 + where + (Vector2 vx vy) = v + -- (px, py) is the vector/point after translating the circle s.t. it is centered at the + -- origin + (Vector2 px py) = p' .-. c + + -- let q lambda be the intersection point. We solve the following equation + -- solving the equation (q_x)^2 + (q_y)^2 = r^2 then yields the equation + -- L^2(vx^2 + vy^2) + L2(px*vx + py*vy) + px^2 + py^2 = 0 + -- where L = \lambda + aa = vx^2 + vy^2 + bb = 2 * (px * vx + py * vy) + cc = px^2 + py^2 - r^2 + discr = bb^2 - 4*aa*cc + instance (Ord r, Floating r) => Line 2 r `IsIntersectableWith` Circle p r where nonEmptyIntersection = defaultNonEmptyIntersection - (Line p' v) `intersect` (Circle (c :+ _) r) = case discr `compare` 0 of LT -> coRec NoIntersection EQ -> coRec . Touching $ q' (lambda (+)) @@ -233,12 +254,19 @@ -- | A line segment may not intersect a circle, touch it, or intersect it -- properly in one or two points. -type instance IntersectionOf (LineSegment 2 p r) (Circle q r) = [ NoIntersection - , Touching (Point 2 r) - , Point 2 r - , (Point 2 r, Point 2 r) - ] +type instance IntersectionOf (LineSegment d p r) (Sphere d q r) = [ NoIntersection + , Touching (Point d r) + , Point d r + , (Point d r, Point d r) + ] +instance (Ord r, Fractional r, Arity d) + => LineSegment d p r `HasIntersectionWith` Sphere d q r where + seg `intersects` (Sphere (c :+ _) r) = let closest = pointClosestTo c (supportingLine seg) + in case squaredEuclideanDist c closest `compare` r of + LT -> True + EQ -> closest `intersects` seg + GT -> False instance (Ord r, Floating r) => LineSegment 2 p r `IsIntersectableWith` Circle q r where diff --git a/src/Data/Geometry/BezierSpline.hs b/src/Data/Geometry/BezierSpline.hs --- a/src/Data/Geometry/BezierSpline.hs +++ b/src/Data/Geometry/BezierSpline.hs @@ -8,40 +8,60 @@ -- Maintainer : Frank Staals -------------------------------------------------------------------------------- module Data.Geometry.BezierSpline( - BezierSpline (BezierSpline) + BezierSpline (BezierSpline, Bezier2, Bezier3) , controlPoints , fromPointSeq + , endPoints + , Data.Geometry.BezierSpline.reverse , evaluate , split + , splitMany + , splitMonotone + , splitByPoints + , extension + , extend + , growTo + , merge , subBezier , tangent , approximate , parameterOf , snap - - , pattern Bezier2, pattern Bezier3 - + , intersectB , colinear - , lineApproximate , quadToCubic ) where -import Control.Lens hiding (Empty) -import qualified Data.Foldable as F +import Algorithms.Geometry.ConvexHull.GrahamScan +import Algorithms.Geometry.SmallestEnclosingBall.RIC +import Algorithms.Geometry.SmallestEnclosingBall.Types +import Control.Lens hiding (Empty) +import Data.Ext +import qualified Data.Foldable as F +import Data.Geometry.Ball +import Data.Geometry.Box.Internal import Data.Geometry.Line +import Data.Geometry.LineSegment hiding (endPoints) import Data.Geometry.Point +import Data.Geometry.PolyLine (PolyLine(..)) +import Data.Geometry.Polygon +import Data.Geometry.Polygon.Convex hiding (merge) import Data.Geometry.Properties import Data.Geometry.Transformation -import Data.Geometry.Vector -import Data.LSeq (LSeq) -import qualified Data.LSeq as LSeq -import Data.Sequence (Seq (..)) -import qualified Data.Sequence as Seq -import Data.Traversable (fmapDefault, foldMapDefault) +import Data.Geometry.Vector hiding (init) +import Data.LSeq (LSeq) +import qualified Data.LSeq as LSeq +import Data.List (sort) +import qualified Data.List.NonEmpty as NonEmpty +import Data.Sequence (Seq(..)) +import qualified Data.Sequence as Seq +import Data.Traversable (fmapDefault,foldMapDefault) import GHC.TypeNats -import qualified Test.QuickCheck as QC +import qualified Test.QuickCheck as QC +-- import Debug.Trace + -------------------------------------------------------------------------------- -- | Datatype representing a Bezier curve of degree \(n\) in \(d\)-dimensional space. @@ -49,7 +69,8 @@ -- makeLenses ''BezierSpline -- | Bezier control points. With n degrees, there are n+1 control points. -controlPoints :: Iso (BezierSpline n1 d1 r1) (BezierSpline n2 d2 r2) (LSeq (1+n1) (Point d1 r1)) (LSeq (1+n2) (Point d2 r2)) +controlPoints :: Iso (BezierSpline n1 d1 r1) (BezierSpline n2 d2 r2) + (LSeq (1+n1) (Point d1 r1)) (LSeq (1+n2) (Point d2 r2)) controlPoints = iso _controlPoints BezierSpline -- | Quadratic Bezier Spline @@ -66,19 +87,32 @@ Bezier3 p q r s = fromPointSeq . Seq.fromList $ [p,q,r,s] {-# COMPLETE Bezier3 #-} +-- | Constructs the Bezier Spline from a given sequence of points. +fromPointSeq :: Seq (Point d r) -> BezierSpline n d r +fromPointSeq = BezierSpline . LSeq.promise . LSeq.fromSeq + + deriving instance (Arity d, Eq r) => Eq (BezierSpline n d r) type instance Dimension (BezierSpline n d r) = d type instance NumType (BezierSpline n d r) = r - instance (Arity n, Arity d, QC.Arbitrary r) => QC.Arbitrary (BezierSpline n d r) where arbitrary = fromPointSeq . Seq.fromList <$> QC.vector (fromIntegral . (1+) . natVal $ C @n) --- | Constructs the Bezier Spline from a given sequence of points. -fromPointSeq :: Seq (Point d r) -> BezierSpline n d r -fromPointSeq = BezierSpline . LSeq.promise . LSeq.fromSeq +{- +instance (Arity n, Arity d, QC.Arbitrary r, Ord r) => QC.Arbitrary (BezierSpline n d r) where + arbitrary = fromPointSeq . Seq.fromList <$> allDifferent (fromIntegral . (1+) . natVal $ C @n) +-- | Generates a set of unique items. +allDifferent :: (Ord a, QC.Arbitrary a) => Int -> QC.Gen [a] +allDifferent n = take n . Set.toList . go maxattempts mempty <$> QC.infiniteList + where + maxattempts = 100 + go 0 s _ = s -- too many attempts + go t s (x:xs) | Set.size s == n = s + | otherwise = go (t-1) (Set.insert x s) xs +-} instance (Arity d, Show r) => Show (BezierSpline n d r) where show (BezierSpline ps) = @@ -100,38 +134,79 @@ instance PointFunctor (BezierSpline n d) where pmap f = over controlPoints (fmap f) +-------------------------------------------------------------------------------- + +-- | Convert a quadratic bezier to a cubic bezier. +quadToCubic :: Fractional r => BezierSpline 2 2 r -> BezierSpline 3 2 r +quadToCubic (Bezier2 a (Point b) c) = + Bezier3 a (Point $ (1/3)*^ (toVec a ^+^ 2*^b)) (Point $ (1/3)*^ (2*^ b ^+^ toVec c)) c + +-------------------------------------------------------------------------------- + +-- | Reverse a BezierSpline +reverse :: (Arity d, Ord r, Num r) => BezierSpline n d r -> BezierSpline n d r +reverse = controlPoints %~ LSeq.reverse + + -- | Evaluate a BezierSpline curve at time t in [0, 1] -- -- pre: \(t \in [0,1]\) -evaluate :: (Arity d, Ord r, Num r) => BezierSpline n d r -> r -> Point d r +evaluate :: (Arity d, Eq r, Num r) => BezierSpline n d r -> r -> Point d r +evaluate b 0 = fst $ endPoints b +evaluate b 1 = snd $ endPoints b evaluate b t = evaluate' (b^.controlPoints.to LSeq.toSeq) where evaluate' = \case (p :<| Empty) -> p pts@(_ :<| tl) -> let (ini :|> _) = pts in evaluate' $ Seq.zipWith blend ini tl _ -> error "evaluate: absurd" - blend p q = p .+^ t *^ (q .-. p) --- | Tangent to the bezier spline at the starting point. +-- | Extract a tangent vector from the first to the second control point. tangent :: (Arity d, Num r, 1 <= n) => BezierSpline n d r -> Vector d r -tangent b = b^?!controlPoints.ix 1 .-. b^?!controlPoints.ix 0 +tangent b = b^?!controlPoints.ix 1 .-. b^?!controlPoints.ix 0 --- | Restrict a Bezier curve to th,e piece between parameters t < u in [0, 1]. +-- | Return the endpoints of the Bezier spline. +endPoints :: BezierSpline n d r -> (Point d r, Point d r) +endPoints b = let (p LSeq.:<| _) = b^.controlPoints + (_ LSeq.:|> q) = b^.controlPoints + in (p,q) + + + + +-- | Restrict a Bezier curve to the piece between parameters t < u in [0, 1]. subBezier :: (KnownNat n, Arity d, Ord r, Num r) => r -> r -> BezierSpline n d r -> BezierSpline n d r subBezier t u = fst . split u . snd . split t + +-- | Compute the convex hull of the control polygon of a 2-dimensional Bezier curve. +-- Should also work in any dimension, but convex hull is not yet implemented. +convexHullB :: (Ord r, Fractional r) => BezierSpline n 2 r -> ConvexPolygon () r +convexHullB = convexHull . NonEmpty.fromList . fmap ext . F.toList . _controlPoints + +-------------------------------------------------------------------------------- + -- | Split a Bezier curve at time t in [0, 1] into two pieces. -split :: forall n d r. (KnownNat n, Arity d, Ord r, Num r) - => r -> BezierSpline n d r -> (BezierSpline n d r, BezierSpline n d r) -split t b | t < 0 || t > 1 = error "Split parameter out of bounds." - | otherwise = let n = fromIntegral $ natVal (C @n) - ps = collect t $ b^.controlPoints - in ( fromPointSeq . Seq.take (n + 1) $ ps - , fromPointSeq . Seq.drop n $ ps - ) +split :: forall n d r. (KnownNat n, Arity d, Ord r, Num r) + => r -> BezierSpline n d r -> (BezierSpline n d r, BezierSpline n d r) +split t b | t < 0 = error "split: t < 0" -- ++ show t ++ " < 0" + | t > 1 = error "split: t > 1" -- ++ show t ++ " > 1" + | otherwise = splitRaw t b + +-- | Split without parameter check. If t outside [0,1], will actually extend the curve +-- rather than split it. +splitRaw :: forall n d r. (KnownNat n, Arity d, Ord r, Num r) + => r -> BezierSpline n d r -> (BezierSpline n d r, BezierSpline n d r) +splitRaw t b = let n = fromIntegral $ natVal (C @n) + ps = collect t $ b^.controlPoints + in ( fromPointSeq . Seq.take (n + 1) $ ps + , fromPointSeq . Seq.drop (n + 0) $ ps + ) + +-- | implementation of splitRaw collect :: (Arity d, Ord r, Num r) => r -> LSeq n (Point d r) -> Seq (Point d r) collect t = go . LSeq.toSeq where @@ -142,48 +217,168 @@ blend p q = p .+^ t *^ (q .-. p) --- {- +-- | Split a Bezier curve into many pieces. +-- Todo: filter out duplicate parameter values! +splitMany :: forall n d r. (KnownNat n, Arity d, Ord r, Fractional r) + => [r] -> BezierSpline n d r -> [BezierSpline n d r] +splitMany = splitManySorted . sort . map (restrict "splitMany" 0 1) --- -- | Merge to Bezier pieces. Assumes they can be merged into a single piece of the same degree --- -- (as would e.g. be the case for the result of a 'split' operation). --- -- Does not test whether this is the case! --- merge :: (Arity d, Ord r, Num r) => (Bezier d r, Bezier d r) -> Bezier d r + where splitManySorted [] b' = [b'] + splitManySorted (t : ts) b' = let (a,c) = split t b' + in a : splitManySorted (map (rescale t) ts) c + rescale :: r -> r -> r + rescale 1 _ = 1 + rescale t u = (u - t) / (1 - t) --- -} --- | Approximate Bezier curve by Polyline with given resolution. -approximate :: forall n d r. (KnownNat n, Arity d, Ord r, Fractional r) - => r -> BezierSpline n d r -> [Point d r] -approximate eps b - | squaredEuclideanDist p q < eps^2 = [p,q] - | otherwise = let (b1, b2) = split 0.5 b - in approximate eps b1 ++ tail (approximate eps b2) - where - p = b^.controlPoints.to LSeq.head - q = b^.controlPoints.to LSeq.last +-- | Cut a Bezier curve into $x_i$-monotone pieces. +-- Can only be solved exactly for degree 4 or smaller. +-- Only gives rational result for degree 2 or smaller. +-- Currentlly implemented for degree 3. +splitMonotone :: (Arity d, Ord r, Enum r, Floating r) => Int -> BezierSpline 3 d r -> [BezierSpline 3 d r] +splitMonotone i b = splitMany (locallyExtremalParameters i b) b --- | Given a point on (or close to) a Bezier curve, return the corresponding parameter value. --- (For points far away from the curve, the function will return the parameter value of --- an approximate locally closest point to the input point.) -parameterOf :: (Arity d, Ord r, Fractional r) => BezierSpline n d r -> Point d r -> r -parameterOf b p = binarySearch (qdA p . evaluate b) treshold (1 - treshold) - where treshold = 0.0001 +{- +type family RealTypeConstraint (n :: Nat) (r :: *) :: Constraint where + RealTypeConstraint 1 r = (Fractional r) + RealTypeConstraint 2 r = (Fractional r) + RealTypeConstraint 3 r = (Floating r) + RealTypeConstraint 4 r = (Floating r) + RealTypeConstraint 5 r = (Floating r) + RealTypeConstraint n r = TypeError "" +-} -binarySearch :: (Ord r, Fractional r) => (r -> r) -> r -> r -> r -binarySearch f l r | abs (f l - f r) < treshold = m - | derivative f m > 0 = binarySearch f l m - | otherwise = binarySearch f m r - where m = (l + r) / 2 - treshold = 0.0001 +-- | Report all parameter values at which the derivative of the $i$th coordinate is 0. +locallyExtremalParameters :: (Arity d, Ord r, Enum r, Floating r) + => Int -> BezierSpline 3 d r -> [r] +locallyExtremalParameters i curve = + let [x1, x2, x3, x4] = map (view $ unsafeCoord i) $ F.toList $ _controlPoints curve + a = 3 * x4 - 9 * x3 + 9 * x2 - 3 * x1 + b = 6 * x1 - 12 * x2 + 6 * x3 + c = 3 * x2 - 3 * x1 + in filter (\j -> 0 <= j && j <= 1) $ solveQuadraticEquation a b c -derivative :: Fractional r => (r -> r) -> r -> r -derivative f x = (f (x + delta) - f x) / delta - where delta = 0.00001 --- | Snap a point close to a Bezier curve to the curve. -snap :: (Arity d, Ord r, Fractional r) => BezierSpline n d r -> Point d r -> Point d r -snap b = evaluate b . parameterOf b +-- | Subdivide a curve based on a sequence of points. +-- Assumes these points are all supposed to lie on the curve, and +-- snaps endpoints of pieces to these points. +-- (higher dimensions would work, but depends on convex hull) +splitByPoints :: (KnownNat n, Ord r, RealFrac r) + => r -> [Point 2 r] -> BezierSpline n 2 r -> [BezierSpline n 2 r] +splitByPoints treshold points curve = + let a = fst $ endPoints curve + b = snd $ endPoints curve + intern = filter (\p -> p /= a && p /= b) points + times = map (parameterOf treshold curve) intern + tipos = sort $ zip times intern + pieces = splitMany (map fst tipos) curve + stapts = a : map snd tipos + endpts = map snd tipos ++ [b] + in zipWith3 snapEndpoints stapts endpts pieces +-------------------------------------------------------------------------------- + +-- | Extend a Bezier curve to a parameter value t outside the interval [0,1]. +-- For t < 0, returns a Bezier representation of the section of the underlying curve +-- from parameter value t until paramater value 0. For t > 1, the same from 1 to t. +-- +-- pre: t outside [0,1] +extension :: forall n d r. (KnownNat n, Arity d, Ord r, Num r) + => r -> BezierSpline n d r -> BezierSpline n d r +extension t b | t > 0 && t < 1 = error "extension: 0 < t < 1" -- ++ show t ++ " < 1" + | t <= 0 = fst $ splitRaw t b + | otherwise {- t >= 1-} = snd $ splitRaw t b + +-- | Extend a Bezier curve to a parameter value t outside the interval [0,1]. +-- For t < 0, returns a Bezier representation of the section of the underlying curve +-- from parameter value t until paramater value 1. For t > 1, the same from 0 to t. +-- +-- pre: t outside [0,1] +extend :: forall n d r. (KnownNat n, Arity d, Ord r, Num r) + => r -> BezierSpline n d r -> BezierSpline n d r +extend t b | t > 0 && t < 1 = error "extend: 0 < t < 1" -- ++ show t ++ " < 1" + | t <= 0 = snd $ splitRaw t b + | otherwise {- t >= 1 -} = fst $ splitRaw t b + + +-- | Extend a Bezier curve to a point not on the curve, but on / close +-- to the extended underlying curve. +growTo :: (KnownNat n, Arity d, Ord r, Fractional r) + => r -> Point d r -> BezierSpline n d r -> BezierSpline n d r +growTo treshold p b = + let t = extendedParameterOf treshold b p + r | t < 0 = extend t b + | t > 1 = extend t b + | otherwise = b + in r + +{- + +-- | Tries to fit a degree n Bezier curve through a list of points, with error parameter eps. +-- Either returns an appropriate curve, or fails. +fit :: r -> [Point 2 r] -> Maybe (Bezier n d r) +fit eps pts + +-} + + +-------------------------------------------------------------------------------- + +-- | Merge two Bezier pieces. Assumes they can be merged into a single piece of the same degree +-- (as would e.g. be the case for the result of a 'split' operation). +-- Does not test whether this is the case! +merge :: (KnownNat n, Arity d, Ord r, Fractional r) + => r -> BezierSpline n d r -> BezierSpline n d r -> BezierSpline n d r +merge treshold b1 b2 = let (p1, q1) = endPoints b1 + (p2, q2) = endPoints b2 + result | q1 /= p2 = error "merge: something is wrong, maybe need to flip one of the curves?" + | otherwise = snapEndpoints p1 q2 $ growTo treshold p1 b2 + in result + +-- need distance function between polyBeziers... + + +-------------------------------------------------------------------------------- + + +-- | Approximate Bezier curve by Polyline with given resolution. That +-- is, every point on the approximation will have distance at most res +-- to the Bezier curve. +approximate :: (KnownNat n, Arity d, Ord r, Fractional r) + => r -> BezierSpline n d r -> PolyLine d () r +approximate res = PolyLine . fmap ext . approximate' res + +-- | implementation of approximate; returns the polyline as an LSeq +approximate' :: (KnownNat n, Arity d, Ord r, Fractional r) + => r -> BezierSpline n d r -> LSeq 2 (Point d r) +approximate' res = LSeq.promise . LSeq.fromSeq . go + where + go b | flat res b = let (p,q) = endPoints b in Seq.fromList [p,q] + | otherwise = let (b1, b2) = split 0.5 b in go b1 <> Seq.drop 1 (go b2) + +-- | Test whether a Bezier curve can be approximated by a single line segment, +-- given the resolution parameter. +flat :: (KnownNat n, Arity d, Ord r, Fractional r) => r -> BezierSpline n d r -> Bool +flat r b = let p = fst $ endPoints b + q = snd $ endPoints b + s = ClosedLineSegment (p :+ ()) (q :+ ()) + e t = sqDistanceToSeg (evaluate b t) s < r ^ 2 + in qdA p q < r ^ 2 || all e [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9] + +-- seems this is now covered by approximate +-- +-- +-- -- | Approximate curve as line segments where no point on the curve is further away +-- -- from the nearest line segment than the given tolerance. +-- lineApproximate :: (Ord r, Fractional r) => r -> BezierSpline 3 2 r -> [Point 2 r] +-- lineApproximate eps bezier +-- | colinear eps bezier = +-- [ bezier^.controlPoints.to LSeq.head +-- , bezier^.controlPoints.to LSeq.last ] +-- | otherwise = +-- let (b1, b2) = split 0.5 bezier +-- in lineApproximate eps b1 ++ tail (lineApproximate eps b2) + -- If both control points are on the same side of the straight line from the start and end -- points then the curve is guaranteed to be within 3/4 of the distance from the straight line -- to the furthest control point. @@ -206,18 +401,251 @@ | sameSide = 9/16 * maxDist | otherwise = 16/81 * maxDist --- | Approximate curve as line segments where no point on the curve is further away --- from the nearest line segment than the given tolerance. -lineApproximate :: (Ord r, Fractional r) => r -> BezierSpline 3 2 r -> [Point 2 r] -lineApproximate eps bezier - | colinear eps bezier = - [ bezier^.controlPoints.to LSeq.head - , bezier^.controlPoints.to LSeq.last ] - | otherwise = - let (b1, b2) = split 0.5 bezier - in lineApproximate eps b1 ++ tail (lineApproximate eps b2) +-------------------------------------------------------------------------------- --- | Convert a quadratic bezier to a cubic bezier. -quadToCubic :: (Fractional r) => BezierSpline 2 2 r -> BezierSpline 3 2 r -quadToCubic (Bezier2 a b c) = - Bezier3 a ((1/3)*^(Point (toVec a ^+^ 2*^toVec b))) ((1/3)*^(Point (2*^ toVec b ^+^ toVec c))) c +-- general d depends on convex hull +-- parameterOf :: (Arity d, Ord r, Fractional r) => BezierSpline n d r -> Point d r -> r +-- +-- | Given a point on (or within distance treshold to) a Bezier curve, return the parameter value +-- of some point on the curve within distance treshold from p. +-- For points farther than treshold from the curve, the function will attempt to return the +-- parameter value of an approximate locally closest point to the input point, but no guarantees. +parameterOf :: (KnownNat n, Ord r, RealFrac r) => r -> BezierSpline n 2 r -> Point 2 r -> r +parameterOf treshold b p | closeEnough treshold p $ fst $ endPoints b = 0 + | closeEnough treshold p $ snd $ endPoints b = 1 + | otherwise = parameterInterior treshold b p + +-- parameterInterior is slow, look into algebraic solution? + +-- general d depends on convex hull +parameterInterior :: (KnownNat n, Ord r, RealFrac r) => r -> BezierSpline n 2 r -> Point 2 r -> r +parameterInterior treshold b p | sqrad (F.toList $ view controlPoints b) < (0.5 * treshold)^2 = 0.5 + | otherwise = + let (b1, b2) = split 0.5 b + recurse1 = 0.5 * parameterInterior treshold b1 p + recurse2 = 0.5 + 0.5 * parameterInterior treshold b2 p + chb1 = _simplePolygon $ convexHullB b1 + chb2 = _simplePolygon $ convexHullB b2 + in1 = sqDistanceToPolygon p chb1 < treshold^2 + in2 = sqDistanceToPolygon p chb2 < treshold^2 + result | in1 && in2 = betterFit b p recurse1 recurse2 + | in2 && not in2 = recurse1 + | not in2 && in2 = recurse2 + | sqDistanceToPolygon p chb1 < sqDistanceToPolygon p chb2 = recurse1 + | otherwise = recurse2 + in result + +-- | Given a point on (or close to) the extension of a Bezier curve, return the corresponding +-- parameter value, which might also be smaller than 0 or larger than 1. +-- (For points far away from the curve, the function will return the parameter value of +-- an approximate locally closest point to the input point.) +-- +-- This implementation is not robust: might return a locally closest point on the curve +-- even though the point lies on another part of the curve. For points on the actual +-- curve, use parameterOf instead. +extendedParameterOf :: (Arity d, KnownNat n, Ord r, Fractional r) + => r -> BezierSpline n d r -> Point d r -> r +extendedParameterOf treshold b p | p == fst (endPoints b) = 0 + | p == snd (endPoints b) = 1 + | otherwise = binarySearch treshold (qdA p . evaluate b) (-100) 100 + +---------------------------------------- +-- * Stuff to implement parameterOf and extendedParameterOf + +betterFit :: (KnownNat n, Arity d, Ord r, Fractional r) + => BezierSpline n d r -> Point d r -> r -> r -> r +betterFit b p t u = + let q = evaluate b t + r = evaluate b u + in if qdA q p < qdA r p then t else u + +sqDistanceToPolygon :: (Ord r, Fractional r) => Point 2 r -> SimplePolygon p r -> r +sqDistanceToPolygon point poly | insidePolygon point poly = 0 + | otherwise = minimum $ map (sqDistanceToSeg point) $ listEdges poly + + +-------------------------------------------------------------------------------- + + + + +-------------------------------------------------------------------------------- + +-- | Given two Bezier curves, list all intersection points. +-- Not exact, since for degree >= 3 there is no closed form. +-- (In principle, this algorithm works in any dimension +-- but this requires convexHull, area/volume, and intersect.) +intersectB :: (KnownNat n, Ord r, RealFrac r) => r -> BezierSpline n 2 r -> BezierSpline n 2 r -> [Point 2 r] +intersectB treshold a b + | a == b = [fst $ endPoints b, snd $ endPoints b] -- should really return the whole curve + | otherwise = let [a1, _a2, _a3, a4] = F.toList $ _controlPoints a + [b1, _b2, _b3, b4] = F.toList $ _controlPoints b + in intersectPointsPoints treshold [a1, a4] [b1, b4] + ++ intersectPointsInterior treshold [a1, a4] b + ++ intersectPointsInterior treshold [b1, b4] a + ++ intersectInteriorInterior treshold [a1, a4, b1, b4] a b + + +closeEnough :: (Arity d, Ord r, Fractional r) => r -> Point d r -> Point d r -> Bool +closeEnough treshold p q = qdA p q < treshold ^ 2 + +intersectPointsPoints :: (Ord r, Fractional r) => r -> [Point 2 r] -> [Point 2 r] -> [Point 2 r] +intersectPointsPoints treshold ps = filter (\q -> any (closeEnough treshold q) ps) + +intersectPointsInterior :: (KnownNat n, Ord r, RealFrac r) => r -> [Point 2 r] -> BezierSpline n 2 r -> [Point 2 r] +intersectPointsInterior treshold ps b = + let [b1, _b2, _b3, b4] = F.toList $ _controlPoints b + nearc p = closeEnough treshold (snap treshold b p) p + near1 = closeEnough treshold b1 + near4 = closeEnough treshold b4 + in filter (\p -> nearc p && not (near1 p) && not (near4 p)) ps + + +intersectInteriorInterior :: (KnownNat n, Ord r, RealFrac r) => r -> [Point 2 r] -> BezierSpline n 2 r -> BezierSpline n 2 r -> [Point 2 r] +intersectInteriorInterior treshold forbidden a b = + let cha = _simplePolygon $ convexHullB a + chb = _simplePolygon $ convexHullB b + (a1, a2) = split 0.5 a + (b1, b2) = split 0.5 b + points = F.toList (view controlPoints a) + ++ F.toList (view controlPoints b) + approx = average points + done | not (cha `intersectsP` chb) = True + | sqrad points < treshold^2 = True + | otherwise = False + result | not (cha `intersectsP` chb) = [] + | any (closeEnough treshold approx) forbidden = [] + | otherwise = [approx] + recurse = intersectInteriorInterior treshold forbidden a1 b1 + ++ intersectInteriorInterior treshold forbidden a1 b2 + ++ intersectInteriorInterior treshold forbidden a2 b1 + ++ intersectInteriorInterior treshold forbidden a2 b2 + in if done then result else recurse + +sqrad :: (Ord r, RealFrac r) => [Point 2 r] -> r +sqrad points | length points < 2 = error "sqrad: not enough points" +sqrad points | otherwise = + let rationalPoints :: [Point 2 Rational] -- smallestEnclosingDisk fails on Floats + rationalPoints = map (traverse %~ realToFrac) points + (a : b : cs) = map (:+ ()) rationalPoints + diskResult = smallestEnclosingDisk' a b cs + in realToFrac $ view squaredRadius $ view enclosingDisk $ diskResult + +average :: (Functor t, Foldable t, Arity d, Fractional r) => t (Point d r) -> Point d r +average ps = origin .+^ foldr1 (^+^) (fmap toVec ps) ^/ realToFrac (length ps) + +{- +type instance IntersectionOf (BezierSpline n 2 r) (BezierSpline n 2 r) = [ NoIntersection + , [Point 2 r] + , BezierSpline n 2 r + ] + + +instance (KnownNat n, Ord r, Fractional r) => (BezierSpline n 2 r) `IsIntersectableWith` (BezierSpline n 2 r) where + nonEmptyIntersection = defaultNonEmptyIntersection + a `intersect` b = a `intersectB` b +-} + + +-- function to test whether two convex polygons intersect +-- for speed, first test bounding boxes +-- maybe would be faster to directly compare bounding boxes of points, rather than +-- call convex hull first? +intersectsP :: (Ord r, Fractional r) => SimplePolygon p r -> SimplePolygon p r -> Bool +intersectsP p q | not $ boundingBox p `intersects` boundingBox q = False + | otherwise = or [a `intersects` b | a <- p & listEdges, b <- q & listEdges] + || (any (flip insidePolygon p) $ map _core $ F.toList $ polygonVertices q) + || (any (flip insidePolygon q) $ map _core $ F.toList $ polygonVertices p) + -- first test bounding box? + + +{- + +instance (Arity d, Floating r) => IsBoxable (BezierSpline 3 d r) where + boundingBox b = foldr1 (<>) $ map (\i -> boundingBox (extremal True i b) <> boundingBox (extremal False i b)) [1 .. d] + +-- | Find extremal points on curve in the $i$th dimension. +extremal :: Floating r => Bool -> Int -> BezierSpline 3 d r -> Point d r +extremal pos i b = + let [p1, _, _, p4] = F.toList $ view controlPoints b + ps = map evaluate $ locallyExtremalParameters i b + candidates = [p1, p4] ++ ps + result | pos = maximumBy (unsafeCoord i . snd) candidates + | not pos = minimumBy (unsafeCoord i . snd) candidates + in result + +-} + + +-------------------------------------------------------------------------------- + +snapEndpoints :: (KnownNat n, Arity d, Ord r, Fractional r) + => Point d r -> Point d r -> BezierSpline n d r -> BezierSpline n d r +snapEndpoints p q curve = + let points = F.toList $ _controlPoints curve + middle = tail . init $ points + new = [p] ++ middle ++ [q] + in fromPointSeq $ Seq.fromList new + + +-- | Snap a point close to a Bezier curve to the curve. +snap :: (KnownNat n, Ord r, RealFrac r) => r -> BezierSpline n 2 r -> Point 2 r -> Point 2 r +snap treshold b = evaluate b . parameterOf treshold b + +-------------------------------------------------------------------------------- +-- * Helper functions + +-- | Solve equation of the form ax^2 + bx + c = 0. +-- If there are multiple solutions, report in ascending order. +-- Attempt at a somewhat robust implementation. +solveQuadraticEquation :: (Ord r, Enum r, Floating r) => r -> r -> r -> [r] +solveQuadraticEquation 0 0 0 = [0..] -- error "infinite solutions" +solveQuadraticEquation _ 0 0 = [0] +solveQuadraticEquation 0 _ 0 = [0] +solveQuadraticEquation 0 0 _ = [] +solveQuadraticEquation a b 0 = sort [0, -b / a] +solveQuadraticEquation a 0 c | (-c / a) < 0 = [] + | (-c / a) == 0 = [0] + | (-c / a) > 0 = [sqrt (-c / a)] +solveQuadraticEquation 0 b c = [-c / b] +solveQuadraticEquation a b c | almostzero a || almostzero (a / b) || almostzero (a / c) = solveQuadraticEquation 0 b c +solveQuadraticEquation a b c = + let d = b^2 - 4 * a * c + result | d == 0 = [-b / (2 * a)] + | d > 0 = [(-b - sqrt d) / (2 * a), (-b + sqrt d) / (2 * a)] + | otherwise = [] + in result + -- trace ("soving equation " ++ show a ++ "x^2 + " ++ show b ++ "x + " ++ show c ++ " = 0") $ result + +-- | Test whether a floating point number is close enough to zero, taking rounding errors into account. +almostzero :: (Floating r, Ord r) => r -> Bool +almostzero x = abs x < epsilon + +-- | Treshold for rounding errors in almostzero test. +-- TODO: Should be different depending on the type. +epsilon :: Floating r => r +epsilon = 0.0001 + + + +-- | This function tests whether a value lies within bounds of a given interval. +-- If not, graciously continues with value snapped to interval. +-- This should never happen, but apparently it sometimes does? +restrict :: (Ord r) => String -> r -> r -> r -> r +restrict f l r x | l > r = error $ f <> ": restrict [l,r] is not an interval" --error $ f ++ ": restrict: [" ++ show l ++ ", " ++ show r ++ "] is not an interval" + -- | x < l = trace (f ++ ": restricting " ++ show x ++ " to [" ++ show l ++ ", " ++ show r ++ "]") l + -- | x > r = trace (f ++ ": restricting " ++ show x ++ " to [" ++ show l ++ ", " ++ show r ++ "]") r + | otherwise = x + + +binarySearch :: (Ord r, Fractional r) + => r -> (r -> r) -> r -> r -> r +binarySearch treshold f l r + | abs (f l - f r) < treshold = restrict "binarySearch" l r m + | derivative f m > 0 = restrict "binarySearch" l r $ binarySearch treshold f l m + | otherwise = restrict "binarySearch" l r $ binarySearch treshold f m r + where m = (l + r) / 2 + +derivative :: Fractional r => (r -> r) -> r -> r +derivative f x = (f (x + delta) - f x) / delta + where delta = 0.0000001 diff --git a/src/Data/Geometry/Box/Corners.hs b/src/Data/Geometry/Box/Corners.hs --- a/src/Data/Geometry/Box/Corners.hs +++ b/src/Data/Geometry/Box/Corners.hs @@ -23,7 +23,9 @@ -------------------------------------------------------------------------------- --- | A Quadrant data type +-- | A data type rperesenting the corners of a box. the order of the +-- Corners is 'northWest, northEast, southEast, southWest', i.e. in +-- clockwise order starting from the topleft. data Corners a = Corners { _northWest :: !a , _northEast :: !a , _southEast :: !a diff --git a/src/Data/Geometry/Box/Internal.hs b/src/Data/Geometry/Box/Internal.hs --- a/src/Data/Geometry/Box/Internal.hs +++ b/src/Data/Geometry/Box/Internal.hs @@ -22,7 +22,7 @@ import qualified Data.Foldable as F import Data.Geometry.Point import Data.Geometry.Properties -import Data.Geometry.Transformation +import Data.Geometry.Transformation.Internal import Data.Geometry.Vector import qualified Data.Geometry.Vector as V import qualified Data.List.NonEmpty as NE @@ -106,6 +106,8 @@ type instance IntersectionOf (Box d p r) (Box d q r) = '[ NoIntersection, Box d () r] +instance (Ord r, Arity d) => Box d p r `HasIntersectionWith` Box d q r + instance (Ord r, Arity d) => Box d p r `IsIntersectableWith` Box d q r where nonEmptyIntersection = defaultNonEmptyIntersection @@ -144,6 +146,9 @@ type instance IntersectionOf (Point d r) (Box d p r) = '[ NoIntersection, Point d r] + +instance (Arity d, Ord r) => Point d r `HasIntersectionWith` Box d p r where + intersects = inBox instance (Arity d, Ord r) => Point d r `IsIntersectableWith` Box d p r where nonEmptyIntersection = defaultNonEmptyIntersection diff --git a/src/Data/Geometry/HalfLine.hs b/src/Data/Geometry/HalfLine.hs --- a/src/Data/Geometry/HalfLine.hs +++ b/src/Data/Geometry/HalfLine.hs @@ -120,6 +120,7 @@ , Point d r ] +instance (Ord r, Fractional r) => HalfLine 2 r `HasIntersectionWith` Line 2 r instance (Ord r, Fractional r) => HalfLine 2 r `IsIntersectableWith` Line 2 r where nonEmptyIntersection = defaultNonEmptyIntersection @@ -130,6 +131,7 @@ :& RNil +instance (Ord r, Fractional r) => HalfLine 2 r `HasIntersectionWith` HalfLine 2 r instance (Ord r, Fractional r) => HalfLine 2 r `IsIntersectableWith` HalfLine 2 r where nonEmptyIntersection = defaultNonEmptyIntersection @@ -155,6 +157,8 @@ ) :& RNil +instance (Ord r, Fractional r) => LineSegment 2 () r `HasIntersectionWith` HalfLine 2 r + instance (Ord r, Fractional r) => LineSegment 2 () r `IsIntersectableWith` HalfLine 2 r where nonEmptyIntersection = defaultNonEmptyIntersection @@ -173,9 +177,11 @@ :& RNil +instance (Ord r, Fractional r, Arity d) => Point d r `HasIntersectionWith` HalfLine d r where + intersects = onHalfLine + instance (Ord r, Fractional r, Arity d) => Point d r `IsIntersectableWith` HalfLine d r where nonEmptyIntersection = defaultNonEmptyIntersection - intersects = onHalfLine p `intersect` hl | p `intersects` hl = coRec p | otherwise = coRec NoIntersection @@ -188,6 +194,8 @@ , Point 2 r , LineSegment 2 () r ] +instance (Ord r, Fractional r) + => HalfLine 2 r `HasIntersectionWith` Boundary (Rectangle p r) instance (Ord r, Fractional r) => HalfLine 2 r `IsIntersectableWith` Boundary (Rectangle p r) where @@ -207,6 +215,8 @@ (True,False) -> coRec $ ClosedLineSegment (ext o) p (True,True) -> coRec s) :& RNil +instance (Ord r, Fractional r) + => HalfLine 2 r `HasIntersectionWith` Rectangle p r instance (Ord r, Fractional r) => HalfLine 2 r `IsIntersectableWith` Rectangle p r where diff --git a/src/Data/Geometry/HalfSpace.hs b/src/Data/Geometry/HalfSpace.hs --- a/src/Data/Geometry/HalfSpace.hs +++ b/src/Data/Geometry/HalfSpace.hs @@ -98,19 +98,19 @@ type instance IntersectionOf (Point d r) (HalfSpace d r) = [NoIntersection, Point d r] -instance (Num r, Ord r, Arity d) => Point d r `IsIntersectableWith` HalfSpace d r where - nonEmptyIntersection = defaultNonEmptyIntersection - +instance (Num r, Ord r, Arity d) => Point d r `HasIntersectionWith` HalfSpace d r where q `intersects` h = q `inHalfSpace` h /= Outside +instance (Num r, Ord r, Arity d) => Point d r `IsIntersectableWith` HalfSpace d r where + nonEmptyIntersection = defaultNonEmptyIntersection q `intersect` h | q `intersects` h = coRec q | otherwise = coRec NoIntersection - type instance IntersectionOf (Line d r) (HalfSpace d r) = [NoIntersection, HalfLine d r, Line d r] +instance (Fractional r, Ord r) => Line 2 r `HasIntersectionWith` HalfSpace 2 r instance (Fractional r, Ord r) => Line 2 r `IsIntersectableWith` HalfSpace 2 r where nonEmptyIntersection = defaultNonEmptyIntersection diff --git a/src/Data/Geometry/HyperPlane.hs b/src/Data/Geometry/HyperPlane.hs --- a/src/Data/Geometry/HyperPlane.hs +++ b/src/Data/Geometry/HyperPlane.hs @@ -48,9 +48,11 @@ type instance IntersectionOf (Point d r) (HyperPlane d r) = [NoIntersection, Point d r] +instance (Num r, Eq r, Arity d) => Point d r `HasIntersectionWith` HyperPlane d r where + q `intersects` (HyperPlane p n) = n `dot` (q .-. p) == 0 + instance (Num r, Eq r, Arity d) => Point d r `IsIntersectableWith` HyperPlane d r where nonEmptyIntersection = defaultNonEmptyIntersection - q `intersects` (HyperPlane p n) = n `dot` (q .-. p) == 0 q `intersect` h | q `intersects` h = coRec q | otherwise = coRec NoIntersection @@ -103,6 +105,8 @@ GT -> Above type instance IntersectionOf (Line 3 r) (Plane r) = [NoIntersection, Point 3 r, Line 3 r] + +instance (Eq r, Fractional r) => Line 3 r `HasIntersectionWith` Plane r instance (Eq r, Fractional r) => Line 3 r `IsIntersectableWith` Plane r where nonEmptyIntersection = defaultNonEmptyIntersection diff --git a/src/Data/Geometry/Interval.hs b/src/Data/Geometry/Interval.hs --- a/src/Data/Geometry/Interval.hs +++ b/src/Data/Geometry/Interval.hs @@ -129,6 +129,7 @@ type instance IntersectionOf (Interval a r) (Interval a r) = [NoIntersection, Interval a r] +instance Ord r => Interval a r `HasIntersectionWith` Interval a r instance Ord r => Interval a r `IsIntersectableWith` Interval a r where nonEmptyIntersection = defaultNonEmptyIntersection diff --git a/src/Data/Geometry/Line.hs b/src/Data/Geometry/Line.hs --- a/src/Data/Geometry/Line.hs +++ b/src/Data/Geometry/Line.hs @@ -48,25 +48,27 @@ type instance IntersectionOf (Point d r) (Line d r) = [NoIntersection, Point d r] -instance (Eq r, Fractional r, Arity d) => Point d r `IsIntersectableWith` Line d r where - nonEmptyIntersection = defaultNonEmptyIntersection +instance (Eq r, Fractional r, Arity d) => Point d r `HasIntersectionWith` Line d r where intersects = onLine +instance {-# OVERLAPPING #-} (Ord r, Num r) => Point 2 r `HasIntersectionWith` Line 2 r where + intersects = onLine2 + +instance (Eq r, Fractional r, Arity d) => Point d r `IsIntersectableWith` Line d r where + nonEmptyIntersection = defaultNonEmptyIntersection p `intersect` l | p `intersects` l = coRec p | otherwise = coRec NoIntersection -instance {-# OVERLAPPING #-} (Ord r, Num r) - => Point 2 r `IsIntersectableWith` Line 2 r where +instance {-# OVERLAPPING #-} (Ord r, Num r) => Point 2 r `IsIntersectableWith` Line 2 r where nonEmptyIntersection = defaultNonEmptyIntersection - intersects = onLine2 p `intersect` l | p `intersects` l = coRec p | otherwise = coRec NoIntersection - type instance IntersectionOf (Line 2 r) (Boundary (Rectangle p r)) = [ NoIntersection, Point 2 r, (Point 2 r, Point 2 r) , LineSegment 2 () r] - instance (Ord r, Fractional r) + => Line 2 r `HasIntersectionWith` Boundary (Rectangle p r) +instance (Ord r, Fractional r) => Line 2 r `IsIntersectableWith` Boundary (Rectangle p r) where nonEmptyIntersection = defaultNonEmptyIntersection @@ -102,7 +104,8 @@ type instance IntersectionOf (Line 2 r) (Rectangle p r) = [ NoIntersection, Point 2 r, LineSegment 2 () r] - +instance (Ord r, Fractional r) + => Line 2 r `HasIntersectionWith` Rectangle p r instance (Ord r, Fractional r) => Line 2 r `IsIntersectableWith` Rectangle p r where nonEmptyIntersection = defaultNonEmptyIntersection diff --git a/src/Data/Geometry/Line/Internal.hs b/src/Data/Geometry/Line/Internal.hs --- a/src/Data/Geometry/Line/Internal.hs +++ b/src/Data/Geometry/Line/Internal.hs @@ -171,6 +171,8 @@ , Line 2 r ] +instance (Eq r, Fractional r) => Line 2 r `HasIntersectionWith` Line 2 r + instance (Eq r, Fractional r) => Line 2 r `IsIntersectableWith` Line 2 r where @@ -234,6 +236,21 @@ :& (H $ \(Point2 _ b) -> Just (vy / vx,b)) :& (H $ \_ -> Nothing) -- l is a vertical line (through x=0) :& RNil + +-- -- | get values a,b,c s.t. the input line is described by ax + by + c = 0 +-- toLinearFunction' :: Line 2 r -> (r,r,r) +-- toLinearFunction' () + +-- | Given a point p and a line l, computes the point q on l closest to p. +pointClosestTo :: (Fractional r, Arity d) => Point d r -> Line d r -> Point d r +pointClosestTo p (Line a m) = a .+^ (t0 *^ m) + where + -- see https://monkeyproofsolutions.nl/wordpress/how-to-calculate-the-shortest-distance-between-a-point-and-a-line/ + t0 = numerator / divisor + numerator = (p .-. a) `dot` m + divisor = m `dot` m + + -- | Result of a side test diff --git a/src/Data/Geometry/LineSegment/Internal.hs b/src/Data/Geometry/LineSegment/Internal.hs --- a/src/Data/Geometry/LineSegment/Internal.hs +++ b/src/Data/Geometry/LineSegment/Internal.hs @@ -43,7 +43,7 @@ import Data.Geometry.Point import Data.Geometry.Properties import Data.Geometry.SubLine -import Data.Geometry.Transformation +import Data.Geometry.Transformation.Internal import Data.Geometry.Vector import Data.Ord (comparing) import Data.Vinyl @@ -234,16 +234,23 @@ instance {-# OVERLAPPING #-} (Ord r, Num r) + => Point 2 r `HasIntersectionWith` LineSegment 2 p r where + intersects = onSegment2 + +instance {-# OVERLAPPING #-} (Ord r, Num r) => Point 2 r `IsIntersectableWith` LineSegment 2 p r where nonEmptyIntersection = defaultNonEmptyIntersection - intersects = onSegment2 p `intersect` seg | p `intersects` seg = coRec p | otherwise = coRec NoIntersection + instance {-# OVERLAPPABLE #-} (Ord r, Fractional r, Arity d) + => Point d r `HasIntersectionWith` LineSegment d p r where + intersects = onSegment + +instance {-# OVERLAPPABLE #-} (Ord r, Fractional r, Arity d) => Point d r `IsIntersectableWith` LineSegment d p r where nonEmptyIntersection = defaultNonEmptyIntersection - intersects = onSegment p `intersect` seg | p `intersects` seg = coRec p | otherwise = coRec NoIntersection @@ -265,7 +272,8 @@ -- work in higher dimensions that might allow us to drop the -- Fractional constraint - +instance (Ord r, Fractional r) => + LineSegment 2 p r `HasIntersectionWith` LineSegment 2 p r instance (Ord r, Fractional r) => LineSegment 2 p r `IsIntersectableWith` LineSegment 2 p r where @@ -276,6 +284,9 @@ :& H coRec :& H (coRec . subLineToSegment) :& RNil + +instance (Ord r, Fractional r) => + LineSegment 2 p r `HasIntersectionWith` Line 2 r where instance (Ord r, Fractional r) => LineSegment 2 p r `IsIntersectableWith` Line 2 r where diff --git a/src/Data/Geometry/PlanarSubdivision.hs b/src/Data/Geometry/PlanarSubdivision.hs --- a/src/Data/Geometry/PlanarSubdivision.hs +++ b/src/Data/Geometry/PlanarSubdivision.hs @@ -22,6 +22,7 @@ import qualified Data.List.NonEmpty as NonEmpty import Data.Geometry.PlanarSubdivision.Basic import Data.Geometry.PlanarSubdivision.Merge +import Data.Geometry.PlanarSubdivision.TreeRep import Data.Geometry.Polygon import Data.Proxy @@ -156,3 +157,26 @@ -- mySubDiv = fromSimplePolygons (Id Test) -- 0 -- (NonEmpty.fromList [simplePg' :+ 1, trianglePG :+ 2]) + + + + +-- type R = Int +-- data MyWorld + +-- mySubDiv :: PlanarSubdivision MyWorld Int (Int,Int) String R +-- mySubDiv = undefined + +-- faceData xs = FaceData (Seq.fromList xs) + + + +-- fromTreeRep :: TreeRep v e f r -> PlanarSubdivision s v e f r +-- fromTreeRep (PlanarSD of' (InnerSD ajs fs)) = undefined + + +-- fromInnerRep :: forall s v e f r. (Ord r, Fractional r) +-- => InnerRep v e f r -> PlanarSubdivision s v e () r +-- fromInnerRep f (InnerSD ajs fs) = fromConnectedSegments (Proxy @s) segs +-- where +-- segs = adjs diff --git a/src/Data/Geometry/PlanarSubdivision/Basic.hs b/src/Data/Geometry/PlanarSubdivision/Basic.hs --- a/src/Data/Geometry/PlanarSubdivision/Basic.hs +++ b/src/Data/Geometry/PlanarSubdivision/Basic.hs @@ -1,6 +1,4 @@ {-# LANGUAGE TemplateHaskell #-} -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE PartialTypeSignatures #-} {-# LANGUAGE ScopedTypeVariables #-} -------------------------------------------------------------------------------- -- | @@ -35,14 +33,16 @@ , components, component , vertices', vertices , edges', edges - , faces', faces, internalFaces + , faces', internalFaces', faces, internalFaces , darts' - -- , traverseVertices, traverseDarts, traverseFaces + , traverseVertices, traverseDarts, traverseFaces + , mapVertices, mapDarts, mapFaces , headOf, tailOf, twin, endPoints , incidentEdges, incomingEdges, outgoingEdges - , nextIncidentEdge + , nextIncidentEdge, prevIncidentEdge + , nextIncidentEdgeFrom, prevIncidentEdgeFrom , neighboursOf , leftFace, rightFace @@ -58,8 +58,10 @@ , faceDataOf , edgeSegment, edgeSegments - , rawFacePolygon, rawFaceBoundary - , rawFacePolygons + , faceBoundary + , internalFacePolygon, internalFacePolygons + , outerFacePolygon, outerFacePolygon' + , facePolygons , VertexId(..), FaceId(..), Dart, World(..) @@ -68,9 +70,14 @@ , dataVal , dartMapping, Raw(..) + + , asLocalD, asLocalV, asLocalF + , Incident (incidences) + , common, commonVertices, commonDarts, commonFaces ) where import Control.Lens hiding (holes, holesOf, (.=)) +import Data.Bifunctor (first, second) import Data.Coerce import Data.Ext import qualified Data.Foldable as F @@ -92,6 +99,7 @@ , HasDataOf(..) ) import qualified Data.Sequence as Seq +import qualified Data.Set as Set import qualified Data.Vector as V import qualified Data.Vector.Mutable as MV import GHC.Generics (Generic) @@ -136,6 +144,7 @@ boundingBox = boundingBoxList' . V.toList . _components +-- | Lens to access a particular component of the planar subdivision. component :: ComponentId s -> Lens' (PlanarSubdivision s v e f r) (Component s r) component ci = components.singular (ix $ unCI ci) @@ -333,6 +342,10 @@ faces' ps = let n = numFaces ps in V.fromList $ map (FaceId . VertexId) [0..n-1] +-- | \( O(n) \). Vector of all primal faces. +internalFaces' :: PlanarSubdivision s v e f r -> V.Vector (FaceId' s) +internalFaces' = V.tail . faces' + -- | \( O(n) \). Vector of all primal faces with associated data. faces :: PlanarSubdivision s v e f r -> V.Vector (FaceId' s, FaceData (Dart s) f) faces ps = (\fi -> (fi,ps^.faceDataOf fi)) <$> faces' ps @@ -416,8 +429,8 @@ in (\d -> g^.dataOf d) <$> ds --- | Given a dart d that points into some vertex v, report the next --- dart e in the cyclic order around v. +-- | Given a dart d that points into some vertex v, report the next dart in the +-- cyclic (counterclockwise) order around v. -- -- running time: \(O(1)\) nextIncidentEdge :: Dart s -> PlanarSubdivision s v e f r -> Dart s @@ -425,6 +438,40 @@ d'' = PG.nextIncidentEdge d' g in g^.dataOf d'' +-- | Given a dart d that points into some vertex v, report the +-- previous dart in the cyclic (counterclockwise) order around v. +-- +-- running time: \(O(1)\) +-- +-- >>> prevIncidentEdge (dart 1 "+1") smallG +-- Dart (Arc 3) +1 +prevIncidentEdge :: Dart s -> PlanarSubdivision s v e f r -> Dart s +prevIncidentEdge d ps = let (_,d',g) = asLocalD d ps + d'' = PG.prevIncidentEdge d' g + in g^.dataOf d'' + +-- | Given a dart d that points away from some vertex v, report the +-- next dart in the cyclic (counterclockwise) order around v. +-- +-- +-- running time: \(O(1)\) +-- +nextIncidentEdgeFrom :: Dart s -> PlanarSubdivision s v e f r -> Dart s +nextIncidentEdgeFrom d ps = let (_,d',g) = asLocalD d ps + d'' = PG.nextIncidentEdgeFrom d' g + in g^.dataOf d'' + +-- | Given a dart d that points into away from vertex v, report the previous dart in the +-- cyclic (counterclockwise) order around v. +-- +-- running time: \(O(1)\) +-- +prevIncidentEdgeFrom :: Dart s -> PlanarSubdivision s v e f r -> Dart s +prevIncidentEdgeFrom d ps = let (_,d',g) = asLocalD d ps + d'' = PG.prevIncidentEdgeFrom d' g + in g^.dataOf d'' + + -- | All edges incident to vertex v in incoming direction -- (i.e. pointing into v) in counterclockwise order around v. -- @@ -467,10 +514,11 @@ fi = PG.rightFace d' g in g^.dataOf fi --- | The darts on the outer boundary of the face, for internal faces --- the darts are in clockwise order. For the outer face the darts are --- in counterclockwise order, and the darts from various components are in no particular order. --- +-- | The darts on the outer boundary of this face. The darts are +-- reported in order along the face. This means that for internal +-- faces the darts are reported in *clockwise* order along the +-- boundary, whereas for the outer face the darts are reported in +-- counter clockwise order. -- -- running time: \(O(k)\), where \(k\) is the output size. outerBoundaryDarts :: FaceId' s -> PlanarSubdivision s v e f r -> V.Vector (Dart s) @@ -478,6 +526,7 @@ where single (_,f',g) = (\d -> g^.dataOf d) <$> PG.boundary f' g + -- | Get the local face and component from a given face. asLocalF :: FaceId' s -> PlanarSubdivision s v e f r -> NonEmpty (ComponentId s, FaceId' (Wrap s), Component s r) @@ -487,8 +536,9 @@ where toLocalF d = let (ci,d',c) = asLocalD d ps in (ci,PG.leftFace d' c,c) --- | The vertices of the outer boundary of the face, for internal faces in --- clockwise order, for the outer face in counter clockwise order. +-- | The vertices of the outer boundary of the face, for internal +-- faces in clockwise order, for the outer face in counter clockwise +-- order. -- -- -- running time: \(O(k)\), where \(k\) is the output size. @@ -526,7 +576,10 @@ asLocalV (VertexId v) ps = let (Raw ci v' _) = ps^?!rawVertexData.ix v in (ci,v',ps^.component ci) --- | Note that using the setting part of this lens may be very expensive!! +-- | Lens to access the vertex data +-- +-- Note that using the setting part of this lens may be very +-- expensive!! (O(n)) vertexDataOf :: VertexId' s -> Lens' (PlanarSubdivision s v e f r ) (VertexData r v) vertexDataOf (VertexId vi) = lens get' set'' @@ -538,10 +591,17 @@ in ps&rawVertexData.ix vi.dataVal .~ (x^.vData) &component ci.PG.vertexDataOf wvdi.location .~ (x^.location) + +-- | Get the location of a vertex in the planar subdivision. +-- +-- Note that the setting part of this lens may be very expensive! +-- Moreover, use with care (as this may destroy planarity etc.) locationOf :: VertexId' s -> Lens' (PlanarSubdivision s v e f r ) (Point 2 r) locationOf v = vertexDataOf v.location +-- | Lens to get the face data of a particular face. Note that the +-- setting part of this lens may be very expensive! (O(n)) faceDataOf :: FaceId' s -> Lens' (PlanarSubdivision s v e f r) (FaceData (Dart s) f) faceDataOf fi = lens getF setF @@ -562,33 +622,54 @@ type DataOf (PlanarSubdivision s v e f r) (FaceId' s) = f dataOf f = faceDataOf f.fData +-- | Traverse the vertices of the planar subdivision +traverseVertices :: Applicative g + => (VertexId' s -> v -> g v') + -> PlanarSubdivision s v e f r -> g (PlanarSubdivision s v' e f r) +traverseVertices h = traverseOf rawVertexData (traverseWith VertexId h) --- -- | Traverse the vertices --- -- --- traverseVertices :: Applicative m --- => (VertexId' s -> v -> m v') --- -> PlanarSubdivision s v e f r --- -> m (PlanarSubdivision s v' e f r) --- traverseVertices f = itraverseOf (vertexData.itraversed) (\i -> f (VertexId i)) +-- | Traverse the darts of the Planar subdivision +traverseDarts :: Applicative g + => (Dart s -> e -> g e') + -> PlanarSubdivision s v e f r -> g (PlanarSubdivision s v e' f r) +traverseDarts h = traverseOf rawDartData (traverseWith toEnum h) --- -- | Traverses the darts --- -- --- traverseDarts :: Applicative m --- => (Dart s -> e -> m e') --- -> PlanarSubdivision s v e f r --- -> m (PlaneGraph s v e' f r) --- traverseDarts f = traverseOf (dart) (PG.traverseDarts f) +-- | Traverse the faces of the planar subdivision. +traverseFaces :: Applicative g + => (FaceId' s -> f -> g f') + -> PlanarSubdivision s v e f r -> g (PlanarSubdivision s v e f' r) +traverseFaces h = traverseOf rawFaceData (traverseFaces' h) + where + traverseFaces' h' = itraverse (\i -> traverse (h' (FaceId . VertexId $ i))) --- -- | Traverses the faces --- -- --- traverseFaces :: Applicative m --- => (FaceId' s -> f -> m f') --- -> PlaneGraph s v e f r --- -> m (PlaneGraph s v e f' r) --- traverseFaces f = traverseOf graph (PG.traverseFaces f) +-- | Helper function to implement traver(vertertices|darts|faces) +traverseWith :: Applicative g + => (Int -> w s) + -> (w s -> v -> g v') + -> V.Vector (Raw ci i v) + -> g (V.Vector (Raw ci i v')) +traverseWith mkIdx h = itraverse (\i -> traverse (h $ mkIdx i)) +-------------------------------------------------------------------------------- +-- | Map with index over all faces +mapFaces :: (FaceId' s -> t -> f') + -> PlanarSubdivision s v e t r -> PlanarSubdivision s v e f' r +mapFaces h = runIdentity . traverseFaces (\i x -> Identity $ h i x) + +-- | Map with index over all vertices +mapVertices :: (VertexId' s -> t -> v') + -> PlanarSubdivision s t e f r -> PlanarSubdivision s v' e f r +mapVertices h = runIdentity . traverseVertices (\i x -> Identity $ h i x) + +-- | Map with index over all darts +mapDarts :: (Dart s -> t -> e') + -> PlanarSubdivision s v t f r -> PlanarSubdivision s v e' f r +mapDarts h = runIdentity . traverseDarts (\i x -> Identity $ h i x) + +-------------------------------------------------------------------------------- + -- | Getter for the data at the endpoints of a dart -- -- running time: \(O(1)\) @@ -621,7 +702,9 @@ edgeSegments ps = (\d -> (d,edgeSegment d ps)) <$> edges' ps --- | Given a dart and the subdivision constructs the line segment representing it +-- | Given a dart and the subdivision constructs the line segment +-- representing it. The segment \(\overline{uv})\) is has \(u\) as its +-- tail and \(v\) as its head. -- -- \(O(1)\) edgeSegment :: Dart s -> PlanarSubdivision s v e f r -> LineSegment 2 v r :+ e @@ -629,45 +712,106 @@ in ClosedLineSegment p q :+ ps^.dataOf d --- | Generates the darts incident to a face, starting with the given dart. +-- | Given a dart d, generates the darts on (the current component of) +-- the boundary of the the face that is to the right of the given +-- dart. The darts are reported in order along the face. This means +-- that for -- +-- - (the outer boundary of an) internal faces the darts are reported +-- in *clockwise* order along the boundary, +-- - the "inner" boundary of a face, i.e. the boundary of ahole, the +-- darts are reported in *counter clockwise* order. -- +-- Note that this latter case means that in the darts of a a component +-- of the outer face are reported in counter clockwise order. +-- -- \(O(k)\), where \(k\) is the number of darts reported boundary' :: Dart s -> PlanarSubdivision s v e f r -> V.Vector (Dart s) boundary' d ps = let (_,d',g) = asLocalD d ps in (\d'' -> g^.dataOf d'') <$> PG.boundary' d' g - --- | Constructs the outer boundary of the face +-- | The outerboundary of the face as a simple polygon. For internal +-- faces the polygon that is reported has its vertices stored in CCW +-- order (as expected). -- --- \(O(k)\), where \(k\) is the complexity of the outer boundary of the face -rawFaceBoundary :: FaceId' s -> PlanarSubdivision s v e f r -> SimplePolygon v r :+ f -rawFaceBoundary i ps = unsafeFromPoints pts :+ (ps^.dataOf i) +-- pre: FaceId refers to an internal face. +-- +-- \(O(k)\), where \(k\) is the complexity of the outer boundary of +-- the face +faceBoundary :: FaceId' s -> PlanarSubdivision s v e f r -> SimplePolygon v r :+ f +faceBoundary i ps = unsafeFromPoints (reverse pts) :+ (ps^.dataOf i) where d = V.head $ outerBoundaryDarts i ps pts = (\d' -> PG.vtxDataToExt $ ps^.vertexDataOf (headOf d' ps)) <$> V.toList (boundary' d ps) - + -- for internal faces boundary' produces the boundary darts in + -- clockwise order. Hence, we reverse the sequence of points we + -- obtain to get the points/vertices in CCW order, so that we can + -- construct a simplepolygon out of them. --- | Constructs the boundary of the given face +-- | Constructs the boundary of the given face. -- -- \(O(k)\), where \(k\) is the complexity of the face -rawFacePolygon :: FaceId' s -> PlanarSubdivision s v e f r - -> SomePolygon v r :+ f -rawFacePolygon i ps = case F.toList $ holesOf i ps of +internalFacePolygon :: FaceId' s -> PlanarSubdivision s v e f r + -> SomePolygon v r :+ f +internalFacePolygon i ps = case F.toList $ holesOf i ps of [] -> Left res :+ x hs -> Right (MultiPolygon res $ map toHole hs) :+ x where - res :+ x = rawFaceBoundary i ps - toHole d = rawFaceBoundary (leftFace d ps) ps ^. core + res :+ x = faceBoundary i ps + toHole d = faceBoundary (leftFace d ps) ps ^. core +-- TODO: Verify that holes are in the right orientation. --- | Lists all *internal* faces of the planar subdivision. -rawFacePolygons :: PlanarSubdivision s v e f r - -> V.Vector (FaceId' s, SomePolygon v r :+ f) -rawFacePolygons ps = fmap (\(i,_) -> (i,rawFacePolygon i ps)) . internalFaces $ ps +-- | Returns a sufficiently large, rectangular, polygon that contains +-- the entire planar subdivision. Each component corresponds to a hole +-- in this polygon. +outerFacePolygon :: (Num r, Ord r) + => PlanarSubdivision s v e f r -> MultiPolygon (Maybe v) r :+ f +outerFacePolygon ps = outerFacePolygon' outer ps & core %~ first (either (const Nothing) Just) + where + outer = rectToPolygon . grow 1 . boundingBox $ ps + rectToPolygon = unsafeFromPoints . reverse . F.toList . corners +-- | Given a sufficiently large outer boundary, draw the outerface as +-- a polygon with a hole. +outerFacePolygon' :: SimplePolygon v' r + -> PlanarSubdivision s v e f r -> MultiPolygon (Either v' v) r :+ f +outerFacePolygon' outer ps = MultiPolygon (first Left outer) holePgs :+ ps^.dataOf i + where + i = outerFaceId ps + holePgs = map getBoundary . F.toList $ holesOf i ps + -- get the bondary of a hole. Note that for holes, the function + -- 'boundary' promisses to report the darts, and therefore the + -- vertices in CCW order. Hence, we can directly construct a SimplePolygon out of it. + getBoundary d = unsafeFromPoints . fmap (second Right) $ faceBoundary' (twin d) + faceBoundary' d = (\d' -> PG.vtxDataToExt $ ps^.vertexDataOf (headOf d' ps)) + <$> V.toList (boundary' d ps) +-- | Procuces a polygon for each *internal* face of the planar +-- subdivision. +internalFacePolygons :: PlanarSubdivision s v e f r + -> V.Vector (FaceId' s, SomePolygon v r :+ f) +internalFacePolygons ps = fmap (\(i,_) -> (i,internalFacePolygon i ps)) . internalFaces $ ps + + +-- | Procuces a polygon for each face of the planar subdivision. +facePolygons :: (Num r, Ord r) + => PlanarSubdivision s v e f r + -> V.Vector (FaceId' s, SomePolygon (Maybe v) r :+ f) +facePolygons ps = V.cons (outerFaceId ps, first Right $ outerFacePolygon ps) ifs + where + ifs = wrapJust <$> internalFacePolygons ps + g :: Bifunctor g => g a b -> g (Maybe a) b + g = first Just + + wrapJust :: (FaceId' s, SomePolygon v r :+ f) + -> (FaceId' s, SomePolygon (Maybe v) r :+ f) + wrapJust (i,(spg :+ f)) = (i,bimap g g spg :+ f) + + + +-- | Mapping between the internal and extenral darts dartMapping :: PlanarSubdivision s v e f r -> V.Vector (Dart (Wrap s), Dart s) dartMapping ps = ps^.component (ComponentId 0).PG.dartData @@ -683,3 +827,63 @@ -- $ trianglePG -- trianglePG = fromPoints . map ext $ [origin, Point2 10 0, Point2 10 10] + + + + + + + + + + + + + + +-------------------------------------------------------------------------------- + + +-- | A class for describing which features (vertex, edge, face) of a planar subdivision +-- can be incident to each other. +class Incident s a b where + incidences :: PlanarSubdivision s v e f r -> a -> [b] + +instance Incident s (VertexId' s) (Dart s) where + incidences psd i = V.toList (incidentEdges i psd) ++ map twin (V.toList $ incidentEdges i psd) + +instance Incident s (VertexId' s) (FaceId' s) where + incidences psd i = map ((flip leftFace) psd) $ V.toList $ incidentEdges i psd + +instance Incident s (Dart s) (VertexId' s) where + incidences psd i = [headOf i psd, tailOf i psd] + +instance Incident s (Dart s) (FaceId' s) where + incidences psd i = [leftFace i psd, rightFace i psd] + +instance Incident s (FaceId' s) (VertexId' s) where + incidences psd i = V.toList $ boundaryVertices i psd + +instance Incident s (FaceId' s) (Dart s) where + incidences psd i = V.toList (outerBoundaryDarts i psd) ++ map twin (V.toList $ outerBoundaryDarts i psd) + +-- | Given two features (vertex, edge, or face) of a subdivision, +-- report all features of a given type that are incident to both. +common :: (Incident s a c, Incident s b c, Ord c) => PlanarSubdivision s v e f r -> a -> b -> [c] +common psd a b = Set.toList $ Set.intersection (Set.fromList $ incidences psd a) (Set.fromList $ incidences psd b) + +-- | Given two features (edge or face) of a subdivision, report all +-- vertices that are incident to both. +commonVertices :: (Incident s a (VertexId' s), Incident s b (VertexId' s)) => PlanarSubdivision s v e f r -> a -> b -> [VertexId' s] +commonVertices = common + +-- | Given two features (vertex or face) of a subdivision, report all +-- edges that are incident to both. Returns both darts of each +-- qualifying edge. +commonDarts :: (Incident s a (Dart s), Incident s b (Dart s)) => PlanarSubdivision s v e f r -> a -> b -> [Dart s] +commonDarts = common + +-- | Given two features (vertex or edge) of a subdivision, report all +-- faces that are incident to both. +commonFaces :: (Incident s a (FaceId' s), Incident s b (FaceId' s)) => PlanarSubdivision s v e f r -> a -> b -> [FaceId' s] +commonFaces = common diff --git a/src/Data/Geometry/PlanarSubdivision/Dynamic.hs b/src/Data/Geometry/PlanarSubdivision/Dynamic.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/PlanarSubdivision/Dynamic.hs @@ -0,0 +1,530 @@ +module Data.Geometry.PlanarSubdivision.Dynamic + ( splitEdge, unSplitEdge + , sproutIntoFace + , splitFace + ) where + +import Control.Lens + +import Data.Vector (Vector, toList, (//), empty) +import qualified Data.Vector as V +import Data.List (sort, sortOn, findIndex) + +import Data.Functor.Identity +import Data.Ext +import Data.Geometry hiding (Vector, head, imap) +import Data.Geometry.PlanarSubdivision +import Data.Geometry.PlanarSubdivision.Raw + +import Data.PlanarGraph (Dart (Dart), Arc (Arc), VertexId (VertexId), FaceId (FaceId), Direction (Positive, Negative)) +import Data.PlaneGraph (PlaneGraph) +import qualified Data.PlaneGraph as PG +import Data.PlaneGraph.AdjRep hiding (id, vData, faces) +import qualified Data.PlaneGraph.AdjRep as AR (id, vData, fData, faces, Face (..)) + +import Data.List.NonEmpty (NonEmpty(..)) +import qualified Data.List.NonEmpty as NonEmpty + +import Debug.Trace + +import Data.Geometry.PlanarSubdivision.Basic + +tracingOn = False + +tr :: Show a => String -> a -> a +tr s a | tracingOn = trace ("\9608 " ++ s ++ ": " ++ show a) a + | otherwise = a + + +-- TO DO: +-- ADD EDGE JOINING TWO COMPONENTS +-- CREATE NEW COMPONENT (SINGLE VERTEX) +-- DELETIONS + + +-- | Splits a given edge of a planar subdivision by inserting a new vertex on the edges. +-- Increases #vertices and #edges by 1. +splitEdge + :: (Show v, Show e, Show f, Show r) + => VertexId' s + -> VertexId' s + -> Point 2 r + -> v + -> (e -> (e, e)) + -> PlanarSubdivision s v e f r + -> PlanarSubdivision s v e f r + +splitEdge a b p v f d = + let (_, la, _) = asLocalV a d + (_, lb, _) = asLocalV b d + v' = (freeVertexId d, v) + fd = freeDart d + f' (Dart i Positive, e) = ((Dart i Positive, fst $ f e), (fd, snd $ f e)) + f' (Dart i Negative, e) = ((twin fd, fst $ f e), (Dart i Negative, snd $ f e)) + in tr "splitEdge" $ d & components' %~ fmap (splitEdgeInPlaneGraph la lb p v' f') + +-- | Sprouts a new edge from a given vertex into the interior of a given (incident) face. +-- Increases #vertices and #edges by 1. +sproutIntoFace + :: (Show v, Show e, Show f, Show r) + => VertexId' s + -> FaceId' s + -> Point 2 r + -> v + -> (e, e) + -> PlanarSubdivision s v e f r + -> PlanarSubdivision s v e f r + +sproutIntoFace a f p v (e1, e2) d = + let [ea] = tr "[ea]" $ filter (\e -> headOf e d == a && leftFace e d == f) $ commonDarts d a f + (_, la, _) = asLocalV a d + (_, lc, _) = asLocalV (tailOf ea d) d + v' = (freeVertexId d, v) + fd = freeDart d + e1' = (fd, e1) + e2' = (twin fd, e2) + in tr "sproutIntoFace" $ d & components' %~ fmap (sproutIntoFaceInPlaneGraph la lc p v' (e1', e2')) + +-- | Inserts a new edge between two given vertices, adjacent to a common face. +-- Increases #edges and #faces by 1. +splitFace + :: (Show v, Show e, Show f, Show r) + => VertexId' s + -> VertexId' s + -> (e, e) + -> (f -> (f, f)) + -> PlanarSubdivision s v e f r + -> PlanarSubdivision s v e f r + +splitFace a b e g d = + let (ca, _, _) = asLocalV a d + (cb, _, _) = asLocalV b d + in if ca == cb then splitFaceSameComponent a b e g d + else splitFaceDifferentComponents a b e g d + +splitFaceSameComponent a b e g d = + let fs = commonFaces d a b + f | length fs == 1 = tr "f(a)" $ headTrace "splitFaceSameComponent f" fs + | otherwise = tr "f(b)" $ headTrace "splitFaceSameComponent f" $ filter (not . isOuterFace) fs + [ea] = tr "[ea]" $ filter (\e -> headOf e d == a && leftFace e d == f) $ commonDarts d a f + [eb] = tr "[eb]" $ filter (\e -> headOf e d == b && leftFace e d == f) $ commonDarts d b f + (_, la, _) = asLocalV a d + (_, lb, _) = asLocalV b d + (_, lc, _) = asLocalV (tailOf ea d) d + (_, ld, _) = asLocalV (tailOf eb d) d + (_, lf, _) :| [] = asLocalF f d + fd = freeDart d + e' = ((fd, fst e), (twin fd, snd e)) + tf = freeFaceId d + g' (ef, x) = ((ef, fst $ g x), (tf, snd $ g x)) + in tr "splitFaceSameComponent" $ d & components' %~ fmap (splitFaceInPlaneGraph (tr "la" la) (tr "lb" lb) (tr "lc" lc) (tr "ld" ld) (tr "lf" lf) e' g') + +splitFaceDifferentComponents = undefined + + +-- | Splits a given edge of a planar subdivision by inserting a new vertex on the edges. +-- Increases #vertices and #edges by 1. +unSplitEdge + :: (Show v, Show e, Show f, Show r) + => VertexId' s + -> ((e, e) -> e) + -> PlanarSubdivision s v e f r + -> PlanarSubdivision s v e f r + +unSplitEdge b f d = + let [a, c] = tr "[a, c]" $ toList $ neighboursOf b d + (_, la, _) = asLocalV a d + (_, lb, _) = asLocalV b d + (_, lc, _) = asLocalV c d + [dab] = filter (\e -> tailOf e d == a) $ commonDarts d a b + [dcb] = filter (\e -> tailOf e d == c) $ commonDarts d b c + f' ((di, ei), (dj, ej)) | di == dab = ( dab, f (ei, ej)) + | di == dcb = (twin dab, f (ei, ej)) + | otherwise = error "you shouldn't call f' on any other dart" + -- no longer used: vertex id b and dart id dcb + in tr "unSplitEdge" $ d & components' %~ fmap (unSplitEdgeInPlaneGraph la lb lc f') +-- globally, need to restore VertexId and DartIds ??? + + + + + +-- nodig: + +freeVertexId :: PlanarSubdivision s v e f r -> VertexId' s +freeDart :: PlanarSubdivision s v e f r -> Dart s +freeFaceId :: PlanarSubdivision s v e f r -> FaceId' s + +freeVertexId = VertexId . numVertices +freeDart = flip Dart Positive . Arc . numEdges +freeFaceId = FaceId . VertexId . numFaces + +components' :: (Show v, Show e, Show f, Show r) => Lens' (PlanarSubdivision s v e f r) (Vector (Component' s v e f r)) +type Component' s v e f r = PlaneGraph (Wrap s) (VertexId' s, v) (Dart s, e) (FaceId' s, f) r +components' = lens getComponents' setComponents' + +getComponents' :: PlanarSubdivision s v e f r -> Vector (Component' s v e f r) +getComponents' p = fmap (addExtraData p) $ p ^. components + +addExtraData :: PlanarSubdivision s v e f r -> Component s r -> Component' s v e f r +addExtraData p c = c & PG.vertexData . traverse %~ (\i -> (i, p ^. dataOf i)) + & PG.rawDartData . traverse %~ (\i -> (i, p ^. dataOf i)) + & PG.faceData . traverse %~ (\i -> (i, p ^. dataOf i)) + +setComponents' :: (Show v, Show e, Show f, Show r) => PlanarSubdivision s v e f r -> Vector (Component' s v e f r) -> PlanarSubdivision s v e f r +setComponents' p cs = p & components .~ fmap remExtraData cs + & rawVertexData .~ (tr "rawVertexData" . vectorise $ getRawVertexData cs) + & rawDartData .~ (tr "rawDartData" . vectorise $ getRawEdgeData cs) + & rawFaceData .~ (tr "rawFaceData" . vectorise $ getRawFaceData cs) + +getRawVertexData :: Vector (Component' s v e f r) + -> [(VertexId' s, Raw s (VertexId' (Wrap s)) v)] +getRawVertexData = concat . imap (\ci g -> map (\(li, VertexData _ (gi, v)) -> (gi, Raw (toEnum ci) li v)) $ toList $ PG.vertices g) . toList + +--getEdgeData :: Vector (Component' s v e f r) -> [(Dart s, (Dart s, e))] +--getEdgeData = map (\(a, b) -> (a, (a, b))) . concatMap (toList . (^. PG.rawDartData)) . toList + +getRawEdgeData :: Vector (Component' s v e f r) + -> [(Dart s, Raw s (Dart (Wrap s)) e)] +getRawEdgeData = concat . imap (\ci g -> map (\(li, (gi, e)) -> (gi, Raw (toEnum ci) li e)) $ toList $ PG.darts g) . toList + + +--getFaceData :: Vector (Component' s v e f r) -> [(FaceId' s, f)] +--getFaceData = concatMap (toList . (^. PG.faceData)) . toList + + +-- data RawFace s f +-- _faceIdx :: !(Maybe (ComponentId s, FaceId' (Wrap s))) +-- _faceDataVal :: !(FaceData (Dart s) f) + +-- | Something in this implementation is not right. It makes asLocalF produce an error. +getRawFaceData :: Vector (Component' s v e f r) + -> [(FaceId' s, RawFace s f)] +getRawFaceData = concat . imap (\ci g -> map (bla ci) $ toList $ PG.faces g) . toList + where + bla ci (li, (gi, f)) | isOuterFace gi = (gi, RawFace Nothing (FaceData Empty f)) + | otherwise = (gi, RawFace (Just (toEnum ci, li)) (FaceData Empty f)) +-- holes are always empty! (where to get them from?) + +isOuterFace :: FaceId' s -> Bool +isOuterFace i = fromEnum i == 0 + +remExtraData :: Component' s v e f r -> Component s r +remExtraData c = c & PG.vertexData . traverse %~ fst + & PG.rawDartData . traverse %~ fst + & PG.faceData . traverse %~ fst + + +vectorise :: (Enum i, Show i) => [(i, a)] -> Vector a +vectorise vs = V.replicate (length vs) undefined // map (\(i, a) -> (fromEnum i, a)) vs + + + + +------------------ +-- PLANE GRAPHS -- +------------------ + + +-- INSERTIONS -- + + +splitEdgeInPlaneGraph + :: (Show v, Show e, Show f, Show r) + => VertexId' s + -> VertexId' s + -> Point 2 r + -> v + -> (e -> (e, e)) + -> PlaneGraph s v e f r + -> PlaneGraph s v e f r +-- LET OP! TEST OF a EN b WEL VOORKOMEN! +splitEdgeInPlaneGraph a b p v f + = tr "splitEdgeInPlaneGraph" + . PG.fromAdjRep undefined + . splitEdgeInAdjRep (fromEnum a) (fromEnum b) p v f + . PG.toAdjRep + +sproutIntoFaceInPlaneGraph + :: (Show v, Show e, Show f, Show r) + => VertexId' s + -> VertexId' s + -> Point 2 r + -> v + -> (e, e) + -> PlaneGraph s v e f r + -> PlaneGraph s v e f r +sproutIntoFaceInPlaneGraph a c p v e g = + let ai = fromEnum a + ci = fromEnum c + in tr "splitEdgeInPlaneGraph" + $ PG.fromAdjRep undefined + $ sproutInAdjRep ai ci p v e + $ PG.toAdjRep g + + +-- PG.toAdjRep :: PlaneGraph s v e f r -> Gr (Vtx v e r) (Face f) +-- PG.fromAdjRep :: proxy s -> Gr (Vtx v e r) (Face f) -> PlaneGraph s v e f r + + +splitFaceInPlaneGraph + :: (Show v, Show e, Show f, Show r) + => VertexId' s -- index van vertex a + -> VertexId' s -- index van vertex b + -> VertexId' s -- index van vertex c + -> VertexId' s -- index van vertex d + -> FaceId' s -- index van te splitsen face + -> (e, e) -- extra data voor nieuwe edge ab + -> (f -> (f, f)) -- functie om face data in twee stukken te knippen + -> PlaneGraph s v e f r -- input graaf + -> PlaneGraph s v e f r -- output graaf + +splitFaceInPlaneGraph a b c d f e h g = + let ai = fromEnum a + bi = fromEnum b + ci = fromEnum c + di = fromEnum d + fi = fromEnum $ tr "fi" $ traceShow (g ^. dataOf f) $ PG.tailOf (PG.boundaryDart f g) g + fj = fromEnum $ tr "fj" $ PG.headOf (PG.boundaryDart f g) g + -- ^ boundaryDart seems not working either + in tr "splitFaceInPlaneGraph" + $ PG.fromAdjRep undefined + $ splitFaceInAdjRep ai bi ci di fi fj e h + $ PG.toAdjRep g + + +-- DELETIONS -- + + +unSplitEdgeInPlaneGraph + :: (Show v, Show e, Show f, Show r) + => VertexId' s + -> VertexId' s + -> VertexId' s + -> ((e, e) -> e) + -> PlaneGraph s v e f r + -> PlaneGraph s v e f r + +unSplitEdgeInPlaneGraph a b c f + = tr "unSplitEdgeInPlaneGraph" + . PG.fromAdjRep undefined + . unSplitEdgeInAdjRep (fromEnum a) (fromEnum b) (fromEnum c) f + . PG.toAdjRep + + +------------- +-- ADJREPS -- +------------- + +-- Gr +-- adjacencies :: [v] +-- faces :: [f] + +-- Vtx +-- id :: Int +-- loc :: Point 2 r +-- adj :: [(Int, e)] +-- vData :: v + +-- Face +-- incidentEdge :: (Int, Int) +-- fData :: f + +--deriving instance (Show v, Show f) => Show (Gr v f) +--deriving instance (Show v, Show e, Show r) => Show (Vtx v e r) +--deriving instance Show f => Show (Face f) + + +-- instance {-# OVERLAPS #-} Show (VertexId s Primal) where show i = 'v' : show (fromEnum i) +-- instance {-# OVERLAPS #-} Show (FaceId s Primal) where show i = 'f' : show (fromEnum i) +-- instance {-# OVERLAPS #-} Show (Dart s, v) where +-- show (Dart (Arc s) Positive, _) = 'd' : show (fromEnum s) ++ "+" +-- show (Dart (Arc s) Negative, _) = 'd' : show (fromEnum s) ++ "-" + +-- instance Show f => Show (Face f) where show f = (show $ AR.fData f) ++ "~>" ++ (show $ incidentEdge f) +-- instance (Show e, Show r) => Show (Vtx v e r) where show v = (show $ AR.id v) ++ "~>" ++ (show $ adj v) +-- instance (Show v, Show f) => Show (Gr v f) where show g = "Gr " ++ (show $ adjacencies g) ++ " " ++ (show $ AR.faces g) + +-- ik heb: +splitEdgeInAdjRep + :: (Show v, Show e, Show f, Show r) + => Int -- index van vertex a + -> Int -- index van vertex b + -> Point 2 r -- locatie voor nieuwe vertex c + -> v -- extra data voor vertex c + -> (e -> (e, e)) -- functie om edge data in twee stukken te knippen + -> Gr (Vtx v e r) (Face f) -- input graaf + -> Gr (Vtx v e r) (Face f) -- output graaf + +splitEdgeInAdjRep a b p v f g = + let n = length $ adjacencies g + -- first find vertices a and b + oa = headTrace "splitEdgeInAdjRep oa" $ filter ((== a) . AR.id) $ adjacencies g + ob = headTrace "splitEdgeInAdjRep ob" $ filter ((== b) . AR.id) $ adjacencies g + os = filter ((lift (&&) (/= a) (/= b)) . AR.id) $ adjacencies g + -- find edge data + e1 = snd $ headTrace "splitEdgeInAdjRep e1" $ filter ((== b) . fst) $ adj oa + e2 = snd $ headTrace "splitEdgeInAdjRep e2" $ filter ((== a) . fst) $ adj ob + -- create new adjacencies to c in a and b + na = oa {adj = replace ((== b) . fst) (const (n, fst $ f e1)) $ adj oa} + nb = ob {adj = replace ((== a) . fst) (const (n, fst $ f e2)) $ adj ob} + -- create new vertex c + nc = Vtx {AR.id = n, loc = p, adj = [(a, snd $ f e2), (b, snd $ f e1)], AR.vData = v} + -- update faces (only if incidentEdge happens to point to ab) + nf = replace ((== (a, b)) . incidentEdge) (\f -> f {incidentEdge = (a, n)}) + $ replace ((== (b, a)) . incidentEdge) (\f -> f {incidentEdge = (b, n)}) + $ AR.faces g + in tr "splitEdgeInAdjRep" $ (tr "original" g) {adjacencies = sortOn AR.id $ na : nb : nc : os, AR.faces = nf} + + +sproutInAdjRep + :: (Show v, Show e, Show f, Show r) + => Int -- index van vertex a + -> Int -- index van vertex c (andere kant van edge a) + -> Point 2 r -- locatie voor nieuwe vertex c + -> v -- extra data voor vertex c + -> (e, e) -- extra data voor nieuwe edge + -> Gr (Vtx v e r) (Face f) -- input graaf + -> Gr (Vtx v e r) (Face f) -- output graaf + +sproutInAdjRep a c p v e g = + let n = length $ adjacencies g + -- first find vertex a + oa = tr "oa" $ headTrace "sproutInAdjRep oa" $ filter ((== a) . AR.id) $ adjacencies g + os = tr "os" $ filter ((/= a) . AR.id) $ adjacencies g + -- need to find index of c + fj (Just x) = x + fj Nothing = error "splitFaceInAdjRep got Nothing" + ci = tr "ci" $ fj $ findIndex ((== c) . fst) $ adj oa + -- create new adjacency to new vertex z in a + na = tr "na" $ oa {adj = take ci (adj oa) ++ (n, fst e) : drop ci (adj oa)} + -- create new vertex z + nz = Vtx {AR.id = n, loc = p, adj = [(a, snd e)], AR.vData = v} + in tr "splitFaceInAdjRep" $ (tr "original" g) {adjacencies = sortOn AR.id $ na : nz : os} + +splitFaceInAdjRep + :: (Show v, Show e, Show f, Show r) + => Int -- index van vertex a + -> Int -- index van vertex b + -> Int -- index van vertex c (andere kant van edge a) + -> Int -- index van vertex d (andere kant van edge b) + -> Int -- index van face edge start + -> Int -- index van face edge eind + -> (e, e) -- extra data voor nieuwe edge ab + -> (f -> (f, f)) -- functie om face data in twee stukken te knippen + -> Gr (Vtx v e r) (Face f) -- input graaf + -> Gr (Vtx v e r) (Face f) -- output graaf + +-- is it easier to split a vertex than a face? + +splitFaceInAdjRep a b c d u v e f g = + let + -- first find vertices a and b + oa = tr "oa" $ headTrace "splitFaceInAdjRep oa" $ filter ((== a) . AR.id) $ adjacencies g + ob = tr "ob" $ headTrace "splitFaceInAdjRep ob" $ filter ((== b) . AR.id) $ adjacencies g + os = tr "os" $ filter ((lift (&&) (/= a) (/= b)) . AR.id) $ adjacencies g + -- insert new adjacency between a and b + fj (Just x) = x + fj Nothing = error "splitFaceInAdjRep got Nothing" + -- need to find indices c and d! + ci = tr "ci" $ fj $ findIndex ((== c) . fst) $ adj oa + di = tr "di" $ fj $ findIndex ((== d) . fst) $ adj ob + -- insert new adjacencies to each other in a and b + na = tr "na" $ oa {adj = take ci (adj oa) ++ (b, fst e) : drop ci (adj oa)} + nb = tr "nb" $ ob {adj = take di (adj ob) ++ (a, snd e) : drop di (adj ob)} + -- find the face that is incident to both a and b + i = tr "i" $ fj $ findIndex ((== (u, v)) . incidentEdge) $ AR.faces g + fd = tr "fd" $ AR.fData $ AR.faces g !! i + ef = tr "ef" $ take i (AR.faces g) ++ drop (i + 1) (AR.faces g) + f1 = tr "f1" $ AR.Face {incidentEdge = (a, b), AR.fData = fst $ f fd} + f2 = tr "f2" $ AR.Face {incidentEdge = (b, a), AR.fData = snd $ f fd} + in tr "splitFaceInAdjRep" $ (tr "original" g) {adjacencies = sortOn AR.id $ na : nb : os, AR.faces = ef ++ [f1, f2]} + + + + + +unSplitEdgeInAdjRep + :: (Show v, Show e, Show f, Show r) + => Int -- index van vertex a + -> Int -- index van vertex b (te verwijderen) + -> Int -- index van vertex c + -> ((e, e) -> e) -- functie om edge data te mergen + -> Gr (Vtx v e r) (Face f) -- input graaf + -> Gr (Vtx v e r) (Face f) -- output graaf + +unSplitEdgeInAdjRep a b c f g = + let n = length $ adjacencies g + -- first find vertices a, b and c + oa = head $ filter ((== a) . AR.id) $ adjacencies g + ob = head $ filter ((== b) . AR.id) $ adjacencies g + oc = head $ filter ((== c) . AR.id) $ adjacencies g + os = filter ((\i -> i /= a && i /= b && i /= c) . AR.id) $ adjacencies g + -- find edge data + eab = snd $ head $ filter ((== b) . fst) $ adj oa + eba = snd $ head $ filter ((== a) . fst) $ adj ob + ebc = snd $ head $ filter ((== c) . fst) $ adj ob + ecb = snd $ head $ filter ((== b) . fst) $ adj oc + -- create new adjacencies between a and c + na = oa {adj = replace ((== b) . fst) (const (c, f (eab, ebc))) $ adj oa} + nc = oc {adj = replace ((== b) . fst) (const (a, f (ecb, eba))) $ adj oc} + nv = sortOn AR.id $ na : nc : os + -- update faces (only if incidentEdge happens to point to ab or bc) + nf = replace ((== (a, b)) . incidentEdge) (\f -> f {incidentEdge = (a, c)}) + $ replace ((== (b, a)) . incidentEdge) (\f -> f {incidentEdge = (c, a)}) + $ replace ((== (b, c)) . incidentEdge) (\f -> f {incidentEdge = (a, c)}) + $ replace ((== (c, b)) . incidentEdge) (\f -> f {incidentEdge = (c, a)}) + $ AR.faces g + -- restore consecutive numbering: replace vertex n-1 by b + ng = replaceIndex (n - 1) b $ (tr "original" g) {adjacencies = nv, AR.faces = nf} + in tr "unSplitEdgeInAdjRep" $ ng + +-- Gr +-- adjacencies :: [v] +-- faces :: [f] + +-- Vtx +-- id :: Int +-- loc :: Point 2 r +-- adj :: [(Int, e)] +-- vData :: v + +-- Face +-- incidentEdge :: (Int, Int) +-- fData :: f + +replaceIndex :: Int -> Int -> Gr (Vtx v e r) (Face f) -> Gr (Vtx v e r) (Face f) +replaceIndex i j g = g { adjacencies = map (replaceIndexAdjacency i j) $ adjacencies g + , AR.faces = map (replaceIndexFace i j) $ AR.faces g + } + +replaceIndexAdjacency :: Int -> Int -> Vtx v e r -> Vtx v e r +replaceIndexAdjacency i j v = v { AR.id = if AR.id v == i then j else AR.id v + , adj = replace ((== i) . fst) (set _1 j) $ adj v + } + +replaceIndexFace :: Int -> Int -> Face f -> Face f +replaceIndexFace i j f | fst (incidentEdge f) == i = f {incidentEdge = incidentEdge f & set _1 j} + | snd (incidentEdge f) == i = f {incidentEdge = incidentEdge f & set _2 j} + | otherwise = f + + +------------- +-- HELPERS -- +------------- + +replace :: (a -> Bool) -> (a -> a) -> [a] -> [a] +replace f g = map $ replace' f g + +replace' :: (a -> Bool) -> (a -> a) -> a -> a +replace' f g x | f x = g x + | otherwise = x + +lift :: (a -> b -> c) -> (d -> a) -> (d -> b) -> d -> c +lift f g h x = f (g x) (h x) + + + +headTrace :: String -> [a] -> a +headTrace s xs | null xs = error $ s ++ ": head of empty list" + | otherwise = head xs diff --git a/src/Data/Geometry/PlanarSubdivision/Raw.hs b/src/Data/Geometry/PlanarSubdivision/Raw.hs --- a/src/Data/Geometry/PlanarSubdivision/Raw.hs +++ b/src/Data/Geometry/PlanarSubdivision/Raw.hs @@ -39,6 +39,14 @@ instance (ToJSON ia, ToJSON a) => ToJSON (Raw s ia a) where toEncoding = genericToEncoding defaultOptions +instance FunctorWithIndex i (Raw ci i) where + imap f (Raw ci i x) = Raw ci i (f i x) +instance FoldableWithIndex i (Raw ci i) where + ifoldMap f (Raw _ i x) = f i x +instance TraversableWithIndex i (Raw ci i) where + itraverse f (Raw ci i x) = Raw ci i <$> f i x + + -- | get the dataVal of a Raw dataVal :: Lens (Raw s ia a) (Raw s ia b) a b dataVal = lens (\(Raw _ _ x) -> x) (\(Raw c i _) y -> Raw c i y) diff --git a/src/Data/Geometry/PlanarSubdivision/TreeRep.hs b/src/Data/Geometry/PlanarSubdivision/TreeRep.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/PlanarSubdivision/TreeRep.hs @@ -0,0 +1,110 @@ +-------------------------------------------------------------------------------- +-- | +-- Module : Data.Geometry.PlanarSubdivision.TreeRep +-- Copyright : (C) Frank Staals +-- License : see the LICENSE file +-- Maintainer : Frank Staals +-- +-- Data types that help encode/decode a planegraph as a JSON/YAML file. +-- +-------------------------------------------------------------------------------- +module Data.Geometry.PlanarSubdivision.TreeRep( PlanarSD(..) + , Vtx(..) + , myTreeRep + ) where + +-- FIXME; uncomment myTreeRep + +import Data.Aeson +import Data.PlaneGraph.AdjRep (Vtx(..)) +import GHC.Generics (Generic) + +import Data.Geometry.Point +import Data.RealNumber.Rational + +-------------------------------------------------------------------------------- + + + +-- | Specify the planar subdivison as a tree of components +data PlanarSD v e f r = PlanarSD + { outerFace :: f -- ^ outer face + , inner :: InnerSD v e f r + } deriving (Show,Eq,Functor,Generic) + +instance (ToJSON r, ToJSON v, ToJSON e, ToJSON f) => ToJSON (PlanarSD v e f r) where + toEncoding = genericToEncoding defaultOptions +instance (FromJSON r, FromJSON v, FromJSON e, FromJSON f) => FromJSON (PlanarSD v e f r) + + +data InnerSD v e f r = InnerSD + { adjs :: [Vtx v e r] -- ^ list of vertices and edges in the + -- components incident to the outer + -- face + , faces :: [(f, [InnerSD v e f r])] -- ^ for each internal + -- face in the component described by adjs its data, + -- and possible holes + } deriving (Show,Eq,Functor,Generic) + +instance (ToJSON r, ToJSON v, ToJSON e, ToJSON f) => ToJSON (InnerSD v e r f) where + toEncoding = genericToEncoding defaultOptions +instance (FromJSON r, FromJSON v, FromJSON e, FromJSON f) => FromJSON (InnerSD v e r f) + + + +-------------------------------------------------------------------------------- + +-- | This represents the following Planar subdivision. Note that the +-- graph is undirected, the arrows are just to indicate what the +-- Positive direction of the darts is. +-- +-- ![mySubDiv](docs/Data/Geometry/PlanarSubdivision/mySubDiv.jpg) +myTreeRep :: PlanarSD Int () String (RealNumber 3) +myTreeRep = PlanarSD "f_infty" (InnerSD ads fs) + where + fs = [ ("f_1", []) + , ("f_2", [f5, f6]) + , ("f_3", []) + , ("f_4", [f7]) + ] + + f5 = InnerSD [ vtx 16 (Point2 3 8) [e 17, e 18] + , vtx 17 (Point2 0 7) [e 16, e 18] + , vtx 18 (Point2 (-1) 4) [e 16, e 17] + ] [("f_5",[])] + + f6 = InnerSD [ vtx 15 (Point2 3 3) [e 14, e 13] + , vtx 13 (Point2 6 4) [e 14, e 15] + , vtx 14 (Point2 3 6) [e 13, e 15] + ] [("f_6",[])] + + f7 = InnerSD [ vtx 19 (Point2 0 9) [e 20, e 23] + , vtx 20 (Point2 0 4) [e 19, e 21] + , vtx 21 (Point2 15 2) [e 20, e 22] + , vtx 22 (Point2 17 5) [e 21, e 23] + , vtx 23 (Point2 15 8) [e 19, e 22] + ] [("f_7",[f8])] + + f8 = InnerSD [ vtx 24 (Point2 14 6) [e 25, e 26] + , vtx 25 (Point2 13 8) [e 24, e 26] + , vtx 26 (Point2 12 5) [e 24, e 25] + ] [("f_8",[])] + + ads = [ vtx 0 (Point2 0 0) [e 1, e 4] + , vtx 1 (Point2 10 2) [e 0, e 5] + , vtx 2 (Point2 9 9) [e 1, e 7, e 3] + , vtx 3 (Point2 0 10) [e 2, e 4] + , vtx 4 (Point2 (-4) 5) [e 0, e 3] + , vtx 5 (Point2 15 3) [e 1, e 6] + , vtx 6 (Point2 20 6) [e 5, e 7] + , vtx 7 (Point2 10 14) [e 2, e 6, e 8] + , vtx 8 (Point2 4 13) [e 7, e 3] + , vtx 9 (Point2 4 (-4)) [e 10, e 11] + , vtx 10 (Point2 8 (-4)) [e 11, e 9] + , vtx 11 (Point2 11 (-2)) [e 10, e 12] + , vtx 12 (Point2 7 (-1)) [e 9, e 11] + ] + + e i = (i,()) + + vtx i p as = Vtx i p as i diff --git a/src/Data/Geometry/Polygon.hs b/src/Data/Geometry/Polygon.hs --- a/src/Data/Geometry/Polygon.hs +++ b/src/Data/Geometry/Polygon.hs @@ -1,3 +1,4 @@ +{-# OPTIONS_GHC -fno-warn-orphans #-} -------------------------------------------------------------------------------- -- | -- Module : Data.Geometry.Polygon @@ -121,9 +122,11 @@ type instance IntersectionOf (Point 2 r) (Polygon t p r) = [NoIntersection, Point 2 r] +instance (Fractional r, Ord r) => Point 2 r `HasIntersectionWith` Polygon t p r where + q `intersects` pg = q `inPolygon` pg /= Outside + instance (Fractional r, Ord r) => Point 2 r `IsIntersectableWith` Polygon t p r where nonEmptyIntersection = defaultNonEmptyIntersection - q `intersects` pg = q `inPolygon` pg /= Outside q `intersect` pg | q `intersects` pg = coRec q | otherwise = coRec NoIntersection diff --git a/src/Data/Geometry/Polygon/Bezier.hs b/src/Data/Geometry/Polygon/Bezier.hs --- a/src/Data/Geometry/Polygon/Bezier.hs +++ b/src/Data/Geometry/Polygon/Bezier.hs @@ -9,10 +9,13 @@ import Control.Lens import Data.Ext -import Data.Geometry.BezierSpline (BezierSpline, lineApproximate, pattern Bezier3) +import Data.Geometry.BezierSpline (BezierSpline, pattern Bezier3) +import qualified Data.Geometry.BezierSpline as Bezier import Data.Geometry.Point +import Data.Geometry.PolyLine(points) import Data.Geometry.Polygon import qualified Data.Vector.Circular as CV +import qualified Data.Foldable as F data PathJoin r = JoinLine @@ -46,8 +49,8 @@ where f :: (Point 2 r :+ PathJoin r, Point 2 r :+ PathJoin r) -> CV.CircularVector (Point 2 r :+ ()) f (a :+ JoinLine, _) = CV.singleton (ext a) - f (a :+ JoinCurve b c, d :+ _) = - CV.unsafeFromList $ map ext $ init (lineApproximate eps (Bezier3 a b c d)) + f (a :+ JoinCurve b c, d :+ _) = let poly = Bezier.approximate eps (Bezier3 a b c d) + in CV.unsafeFromList . init . F.toList $ poly^.points approximateSome :: (Ord r, Fractional r) => r -> SomePolygon (PathJoin r) r -> SomePolygon () r approximateSome eps (Left p) = Left $ approximate eps p diff --git a/src/Data/Geometry/QuadTree/Cell.hs b/src/Data/Geometry/QuadTree/Cell.hs --- a/src/Data/Geometry/QuadTree/Cell.hs +++ b/src/Data/Geometry/QuadTree/Cell.hs @@ -48,6 +48,9 @@ type instance IntersectionOf (Point 2 r) (Cell r) = '[ NoIntersection, Point 2 r] +instance (Ord r, Fractional r) => Point 2 r `HasIntersectionWith` Cell r where + p `intersects` c = p `intersects` toBox c + instance (Ord r, Fractional r) => Point 2 r `IsIntersectableWith` Cell r where nonEmptyIntersection = defaultNonEmptyIntersection p `intersect` c = p `intersect` toBox c diff --git a/src/Data/Geometry/Slab.hs b/src/Data/Geometry/Slab.hs --- a/src/Data/Geometry/Slab.hs +++ b/src/Data/Geometry/Slab.hs @@ -63,6 +63,8 @@ '[Rectangle (a,a) r] +instance Ord r => Slab o a r `HasIntersectionWith` Slab o a r + instance Ord r => Slab o a r `IsIntersectableWith` Slab o a r where nonEmptyIntersection = defaultNonEmptyIntersection @@ -71,6 +73,9 @@ :& H (\i'' -> coRec (Slab i'' :: Slab o a r)) :& RNil +instance Slab Horizontal a r `HasIntersectionWith` Slab Vertical a r where + _ `intersects` _ = True + instance Slab Horizontal a r `IsIntersectableWith` Slab Vertical a r where nonEmptyIntersection _ _ _ = True @@ -105,6 +110,9 @@ [NoIntersection, Line 2 r, LineSegment 2 a r] instance (Fractional r, Ord r, HasBoundingLines o) => + Line 2 r `HasIntersectionWith` Slab o a r + +instance (Fractional r, Ord r, HasBoundingLines o) => Line 2 r `IsIntersectableWith` Slab o a r where nonEmptyIntersection = defaultNonEmptyIntersection @@ -132,6 +140,9 @@ [NoIntersection, SubLine 2 () s r] instance (Fractional r, Ord r, HasBoundingLines o) => + SubLine 2 a r r `HasIntersectionWith` Slab o a r + +instance (Fractional r, Ord r, HasBoundingLines o) => SubLine 2 a r r `IsIntersectableWith` Slab o a r where nonEmptyIntersection = defaultNonEmptyIntersection @@ -151,6 +162,9 @@ type instance IntersectionOf (LineSegment 2 p r) (Slab o a r) = [NoIntersection, LineSegment 2 () r] + +instance (Fractional r, Ord r, HasBoundingLines o) => + LineSegment 2 a r `HasIntersectionWith` Slab o a r instance (Fractional r, Ord r, HasBoundingLines o) => LineSegment 2 a r `IsIntersectableWith` Slab o a r where diff --git a/src/Data/Geometry/SubLine.hs b/src/Data/Geometry/SubLine.hs --- a/src/Data/Geometry/SubLine.hs +++ b/src/Data/Geometry/SubLine.hs @@ -141,6 +141,8 @@ , Point 2 r , SubLine 2 p s r ] +instance (Ord r, Fractional r) => + SubLine 2 p r r `HasIntersectionWith` SubLine 2 p r r {- HLINT ignore "Redundant bracket" -} instance (Ord r, Fractional r) => @@ -164,6 +166,9 @@ s'' = asProperInterval . first (^.extra) $ s'&start.core .~ toOffset' (s'^.start.extra.core) l &end.core .~ toOffset' (s'^.end.extra.core) l + +instance (Ord r, Fractional r) => + SubLine 2 p (UnBounded r) r `HasIntersectionWith` SubLine 2 p (UnBounded r) r instance (Ord r, Fractional r) => SubLine 2 p (UnBounded r) r `IsIntersectableWith` SubLine 2 p (UnBounded r) r where diff --git a/src/Data/Geometry/Transformation.hs b/src/Data/Geometry/Transformation.hs --- a/src/Data/Geometry/Transformation.hs +++ b/src/Data/Geometry/Transformation.hs @@ -1,5 +1,3 @@ -{-# LANGUAGE Unsafe #-} -{-# LANGUAGE UndecidableInstances #-} -------------------------------------------------------------------------------- -- | -- Module : Data.Geometry.Transformation @@ -7,186 +5,55 @@ -- License : see the LICENSE file -- Maintainer : Frank Staals -------------------------------------------------------------------------------- -module Data.Geometry.Transformation where - -import Control.Lens (iso,set,Iso,imap) -import Data.Geometry.Matrix -import Data.Geometry.Matrix.Internal (mkRow) -import Data.Geometry.Point -import Data.Geometry.Properties -import Data.Geometry.Vector -import qualified Data.Geometry.Vector as V -import Data.Proxy -import GHC.TypeLits - -{- $setup ->>> import Data.Geometry.LineSegment ->>> import Data.Ext --} - --------------------------------------------------------------------------------- --- * Transformations - --- | A type representing a Transformation for d dimensional objects -newtype Transformation d r = Transformation { _transformationMatrix :: Matrix (d + 1) (d + 1) r } - --- | Transformations and Matrices are isomorphic. -transformationMatrix :: Iso (Transformation d r) (Transformation d s) - (Matrix (d + 1) (d + 1) r) (Matrix (d + 1) (d + 1) s) -transformationMatrix = iso _transformationMatrix Transformation +module Data.Geometry.Transformation + ( Transformation(Transformation) + , transformationMatrix + , (|.|), identity, inverseOf -deriving instance (Show r, Arity (d + 1)) => Show (Transformation d r) -deriving instance (Eq r, Arity (d + 1)) => Eq (Transformation d r) -deriving instance (Ord r, Arity (d + 1)) => Ord (Transformation d r) -deriving instance Arity (d + 1) => Functor (Transformation d) -deriving instance Arity (d + 1) => Foldable (Transformation d) -deriving instance Arity (d + 1) => Traversable (Transformation d) + , IsTransformable(..) + , transformAllBy + , transformPointFunctor -type instance NumType (Transformation d r) = r + , translation, scaling, uniformScaling --- | Compose transformations (right to left) -(|.|) :: (Num r, Arity (d + 1)) => Transformation d r -> Transformation d r -> Transformation d r -(Transformation f) |.| (Transformation g) = Transformation $ f `multM` g + , translateBy, scaleBy, scaleUniformlyBy + , rotateTo --- if it exists? + , skewX, rotation, reflection, reflectionV, reflectionH --- | Compute the inverse transformation --- --- >>> inverseOf $ translation (Vector2 (10.0) (5.0)) --- Transformation {_transformationMatrix = Matrix (Vector3 (Vector3 1.0 0.0 (-10.0)) (Vector3 0.0 1.0 (-5.0)) (Vector3 0.0 0.0 1.0))} -inverseOf :: (Fractional r, Invertible (d + 1) r) - => Transformation d r -> Transformation d r -inverseOf = Transformation . inverse' . _transformationMatrix + , fitToBox + , fitToBoxTransform + ) where +import Control.Lens +import Data.Ext +import Data.Geometry.Box (Rectangle, IsBoxable) +import qualified Data.Geometry.Box as Box +import Data.Geometry.Properties +import Data.Geometry.Point +import Data.Geometry.Transformation.Internal +import Data.Geometry.Vector -------------------------------------------------------------------------------- --- * Transformable geometry objects --- | A class representing types that can be transformed using a transformation -class IsTransformable g where - transformBy :: Transformation (Dimension g) (NumType g) -> g -> g - --- | Apply a transformation to a collection of objects. --- --- >>> transformAllBy (uniformScaling 2) [Point1 1, Point1 2, Point1 3] --- [Point1 2.0,Point1 4.0,Point1 6.0] -transformAllBy :: (Functor c, IsTransformable g) - => Transformation (Dimension g) (NumType g) -> c g -> c g -transformAllBy t = fmap (transformBy t) - --- | Apply transformation to a PointFunctor, ie something that contains --- points. Polygons, triangles, line segments, etc, are all PointFunctors. --- --- >>> transformPointFunctor (uniformScaling 2) $ OpenLineSegment (Point1 1 :+ ()) (Point1 2 :+ ()) --- OpenLineSegment (Point1 2.0 :+ ()) (Point1 4.0 :+ ()) -transformPointFunctor :: ( PointFunctor g, Fractional r, d ~ Dimension (g r) - , Arity d, Arity (d + 1) - ) => Transformation d r -> g r -> g r -transformPointFunctor t = pmap (transformBy t) - -instance (Fractional r, Arity d, Arity (d + 1)) - => IsTransformable (Point d r) where - transformBy t = Point . transformBy t . toVec - -instance (Fractional r, Arity d, Arity (d + 1)) - => IsTransformable (Vector d r) where - transformBy (Transformation m) v = f $ m `mult` snoc v 1 - where - f u = (/ V.last u) <$> V.init u - - --------------------------------------------------------------------------------- --- * Common transformations - --- | Create translation transformation from a vector. --- --- >>> transformBy (translation $ Vector2 1 2) $ Point2 2 3 --- Point2 3.0 5.0 -translation :: (Num r, Arity d, Arity (d + 1)) - => Vector d r -> Transformation d r -translation v = Transformation . Matrix $ imap transRow (snoc v 1) - --- | Create scaling transformation from a vector. --- --- >>> transformBy (scaling $ Vector2 2 (-1)) $ Point2 2 3 --- Point2 4.0 (-3.0) -scaling :: (Num r, Arity d, Arity (d + 1)) - => Vector d r -> Transformation d r -scaling v = Transformation . Matrix $ imap mkRow (snoc v 1) - --- | Create scaling transformation from a scalar that is applied --- to all dimensions. --- --- >>> transformBy (uniformScaling 5) $ Point2 2 3 --- Point2 10.0 15.0 --- >>> uniformScaling 5 == scaling (Vector2 5 5) --- True --- >>> uniformScaling 5 == scaling (Vector3 5 5 5) --- True -uniformScaling :: (Num r, Arity d, Arity (d + 1)) => r -> Transformation d r -uniformScaling = scaling . pure - - --------------------------------------------------------------------------------- --- * Functions that execute transformations - --- | Translate a given point. --- --- >>> translateBy (Vector2 1 2) $ Point2 2 3 --- Point2 3.0 5.0 -translateBy :: ( IsTransformable g, Num (NumType g) - , Arity (Dimension g), Arity (Dimension g + 1) - ) => Vector (Dimension g) (NumType g) -> g -> g -translateBy = transformBy . translation - --- | Scale a given point. --- --- >>> scaleBy (Vector2 2 (-1)) $ Point2 2 3 --- Point2 4.0 (-3.0) -scaleBy :: ( IsTransformable g, Num (NumType g) - , Arity (Dimension g), Arity (Dimension g + 1) - ) => Vector (Dimension g) (NumType g) -> g -> g -scaleBy = transformBy . scaling - - --- | Scale a given point uniformly in all dimensions. --- --- >>> scaleUniformlyBy 5 $ Point2 2 3 --- Point2 10.0 15.0 -scaleUniformlyBy :: ( IsTransformable g, Num (NumType g) - , Arity (Dimension g), Arity (Dimension g + 1) - ) => NumType g -> g -> g -scaleUniformlyBy = transformBy . uniformScaling - - --- | Row in a translation matrix --- transRow :: forall n r. ( Arity n, Arity (n- 1), ((n - 1) + 1) ~ n --- , Num r) => Int -> r -> Vector n r --- transRow i x = set (V.element (Proxy :: Proxy (n-1))) x $ mkRow i 1 - -transRow :: forall n r. (Arity n, Arity (n + 1), Num r) - => Int -> r -> Vector (n + 1) r -transRow i x = set (V.element (Proxy :: Proxy n)) x $ mkRow i 1 - --------------------------------------------------------------------------------- --- * 3D Rotations - --- | Given three new unit-length basis vectors (u,v,w) that map to (x,y,z), --- construct the appropriate rotation that does this. --- --- -rotateTo :: Num r => Vector 3 (Vector 3 r) -> Transformation 3 r -rotateTo (Vector3 u v w) = Transformation . Matrix $ Vector4 (snoc u 0) - (snoc v 0) - (snoc w 0) - (Vector4 0 0 0 1) - --------------------------------------------------------------------------------- --- * 2D Transformations +-- | Given a rectangle r and a geometry g with its boundingbox, +-- transform the g to fit r. +fitToBox :: forall g r q. + ( IsTransformable g, IsBoxable g, NumType g ~ r, Dimension g ~ 2 + , Ord r, Fractional r + ) => Rectangle q r -> g -> g +fitToBox r g = transformBy (fitToBoxTransform r g) g --- | Skew transformation that keeps the y-coordinates fixed and shifts --- the x coordinates. -skewX :: Num r => r -> Transformation 2 r -skewX lambda = Transformation . Matrix $ Vector3 (Vector3 1 lambda 0) - (Vector3 0 1 0) - (Vector3 0 0 1) +-- | Given a rectangle r and a geometry g with its boundingbox, +-- compute a transformation can fit g to r. +fitToBoxTransform :: forall g r q. ( IsTransformable g, IsBoxable g + , NumType g ~ r, Dimension g ~ 2 + , Ord r, Fractional r + ) => Rectangle q r -> g -> Transformation 2 r +fitToBoxTransform r g = translation v2 |.| uniformScaling lam |.| translation v1 + where + b = Box.boundingBox g + v1 :: Vector 2 r + v1 = negate <$> b^.to Box.minPoint.core.vector + v2 = r^.to Box.minPoint.core.vector + lam = minimum $ (/) <$> Box.size r <*> Box.size b diff --git a/src/Data/Geometry/Transformation/Internal.hs b/src/Data/Geometry/Transformation/Internal.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Transformation/Internal.hs @@ -0,0 +1,219 @@ +{-# LANGUAGE UndecidableInstances #-} +-------------------------------------------------------------------------------- +-- | +-- Module : Data.Geometry.Transformation.Internal +-- Copyright : (C) Frank Staals +-- License : see the LICENSE file +-- Maintainer : Frank Staals +-------------------------------------------------------------------------------- +module Data.Geometry.Transformation.Internal where + +import Control.Lens (iso,set,Iso,imap) +import Data.Geometry.Matrix +import Data.Geometry.Matrix.Internal (mkRow) +import Data.Geometry.Point +import Data.Geometry.Properties +import Data.Geometry.Vector +import qualified Data.Geometry.Vector as V +import Data.Proxy +import GHC.TypeLits + +{- $setup +>>> import Data.Geometry.LineSegment +>>> import Data.Ext +-} + +-------------------------------------------------------------------------------- +-- * Transformations + +-- | A type representing a Transformation for d dimensional objects +newtype Transformation d r = Transformation { _transformationMatrix :: Matrix (d + 1) (d + 1) r } + +-- | Transformations and Matrices are isomorphic. +transformationMatrix :: Iso (Transformation d r) (Transformation d s) + (Matrix (d + 1) (d + 1) r) (Matrix (d + 1) (d + 1) s) +transformationMatrix = iso _transformationMatrix Transformation + +deriving instance (Show r, Arity (d + 1)) => Show (Transformation d r) +deriving instance (Eq r, Arity (d + 1)) => Eq (Transformation d r) +deriving instance (Ord r, Arity (d + 1)) => Ord (Transformation d r) +deriving instance Arity (d + 1) => Functor (Transformation d) +deriving instance Arity (d + 1) => Foldable (Transformation d) +deriving instance Arity (d + 1) => Traversable (Transformation d) + +type instance NumType (Transformation d r) = r + +-- | Compose transformations (right to left) +(|.|) :: (Num r, Arity (d + 1)) => Transformation d r -> Transformation d r -> Transformation d r +(Transformation f) |.| (Transformation g) = Transformation $ f `multM` g + +-- | Identity transformation; i.e. the transformation which does not change anything. +identity :: (Num r, Arity (d + 1)) => Transformation d r +identity = Transformation identityMatrix + +-- if it exists? + +-- | Compute the inverse transformation +-- +-- >>> inverseOf $ translation (Vector2 (10.0) (5.0)) +-- Transformation {_transformationMatrix = Matrix (Vector3 (Vector3 1.0 0.0 (-10.0)) (Vector3 0.0 1.0 (-5.0)) (Vector3 0.0 0.0 1.0))} +inverseOf :: (Fractional r, Invertible (d + 1) r) + => Transformation d r -> Transformation d r +inverseOf = Transformation . inverse' . _transformationMatrix + +-------------------------------------------------------------------------------- +-- * Transformable geometry objects + +-- | A class representing types that can be transformed using a transformation +class IsTransformable g where + transformBy :: Transformation (Dimension g) (NumType g) -> g -> g + +-- | Apply a transformation to a collection of objects. +-- +-- >>> transformAllBy (uniformScaling 2) [Point1 1, Point1 2, Point1 3] +-- [Point1 2.0,Point1 4.0,Point1 6.0] +transformAllBy :: (Functor c, IsTransformable g) + => Transformation (Dimension g) (NumType g) -> c g -> c g +transformAllBy t = fmap (transformBy t) + +-- | Apply transformation to a PointFunctor, ie something that contains +-- points. Polygons, triangles, line segments, etc, are all PointFunctors. +-- +-- >>> transformPointFunctor (uniformScaling 2) $ OpenLineSegment (Point1 1 :+ ()) (Point1 2 :+ ()) +-- OpenLineSegment (Point1 2.0 :+ ()) (Point1 4.0 :+ ()) +transformPointFunctor :: ( PointFunctor g, Fractional r, d ~ Dimension (g r) + , Arity d, Arity (d + 1) + ) => Transformation d r -> g r -> g r +transformPointFunctor t = pmap (transformBy t) + +instance (Fractional r, Arity d, Arity (d + 1)) + => IsTransformable (Point d r) where + transformBy t = Point . transformBy t . toVec + +instance (Fractional r, Arity d, Arity (d + 1)) + => IsTransformable (Vector d r) where + transformBy (Transformation m) v = f $ m `mult` snoc v 1 + where + f u = (/ V.last u) <$> V.init u + + +-------------------------------------------------------------------------------- +-- * Common transformations + +-- | Create translation transformation from a vector. +-- +-- >>> transformBy (translation $ Vector2 1 2) $ Point2 2 3 +-- Point2 3.0 5.0 +translation :: (Num r, Arity d, Arity (d + 1)) + => Vector d r -> Transformation d r +translation v = Transformation . Matrix $ imap transRow (snoc v 1) + +-- | Create scaling transformation from a vector. +-- +-- >>> transformBy (scaling $ Vector2 2 (-1)) $ Point2 2 3 +-- Point2 4.0 (-3.0) +scaling :: (Num r, Arity d, Arity (d + 1)) + => Vector d r -> Transformation d r +scaling v = Transformation . Matrix $ imap mkRow (snoc v 1) + +-- | Create scaling transformation from a scalar that is applied +-- to all dimensions. +-- +-- >>> transformBy (uniformScaling 5) $ Point2 2 3 +-- Point2 10.0 15.0 +-- >>> uniformScaling 5 == scaling (Vector2 5 5) +-- True +-- >>> uniformScaling 5 == scaling (Vector3 5 5 5) +-- True +uniformScaling :: (Num r, Arity d, Arity (d + 1)) => r -> Transformation d r +uniformScaling = scaling . pure + + +-------------------------------------------------------------------------------- +-- * Functions that execute transformations + +-- | Translate a given point. +-- +-- >>> translateBy (Vector2 1 2) $ Point2 2 3 +-- Point2 3.0 5.0 +translateBy :: ( IsTransformable g, Num (NumType g) + , Arity (Dimension g), Arity (Dimension g + 1) + ) => Vector (Dimension g) (NumType g) -> g -> g +translateBy = transformBy . translation + +-- | Scale a given point. +-- +-- >>> scaleBy (Vector2 2 (-1)) $ Point2 2 3 +-- Point2 4.0 (-3.0) +scaleBy :: ( IsTransformable g, Num (NumType g) + , Arity (Dimension g), Arity (Dimension g + 1) + ) => Vector (Dimension g) (NumType g) -> g -> g +scaleBy = transformBy . scaling + + +-- | Scale a given point uniformly in all dimensions. +-- +-- >>> scaleUniformlyBy 5 $ Point2 2 3 +-- Point2 10.0 15.0 +scaleUniformlyBy :: ( IsTransformable g, Num (NumType g) + , Arity (Dimension g), Arity (Dimension g + 1) + ) => NumType g -> g -> g +scaleUniformlyBy = transformBy . uniformScaling + + +-- | Row in a translation matrix +-- transRow :: forall n r. ( Arity n, Arity (n- 1), ((n - 1) + 1) ~ n +-- , Num r) => Int -> r -> Vector n r +-- transRow i x = set (V.element (Proxy :: Proxy (n-1))) x $ mkRow i 1 + +transRow :: forall n r. (Arity n, Arity (n + 1), Num r) + => Int -> r -> Vector (n + 1) r +transRow i x = set (V.element (Proxy :: Proxy n)) x $ mkRow i 1 + +-------------------------------------------------------------------------------- +-- * 3D Rotations + +-- | Given three new unit-length basis vectors (u,v,w) that map to (x,y,z), +-- construct the appropriate rotation that does this. +-- +-- +rotateTo :: Num r => Vector 3 (Vector 3 r) -> Transformation 3 r +rotateTo (Vector3 u v w) = Transformation . Matrix $ Vector4 (snoc u 0) + (snoc v 0) + (snoc w 0) + (Vector4 0 0 0 1) + +-------------------------------------------------------------------------------- +-- * 2D Transformations + +-- | Skew transformation that keeps the y-coordinates fixed and shifts +-- the x coordinates. +skewX :: Num r => r -> Transformation 2 r +skewX lambda = Transformation . Matrix $ Vector3 (Vector3 1 lambda 0) + (Vector3 0 1 0) + (Vector3 0 0 1) + +-- | Create a matrix that corresponds to a rotation by 'a' radians counter-clockwise +-- around the origin. +rotation :: Floating r => r -> Transformation 2 r +rotation a = Transformation . Matrix $ Vector3 (Vector3 (cos a) (- sin a) 0) + (Vector3 (sin a) ( cos a) 0) + (Vector3 0 0 1) + +-- | Create a matrix that corresponds to a reflection in a line through the origin +-- which makes an angle of 'a' radians with the positive 'x'-asis, in counter-clockwise +-- orientation. +reflection :: Floating r => r -> Transformation 2 r +reflection a = rotation a |.| reflectionV |.| rotation (-a) + +-- | Vertical reflection +reflectionV :: Num r => Transformation 2 r +reflectionV = Transformation . Matrix $ Vector3 (Vector3 1 0 0) + (Vector3 0 (-1) 0) + (Vector3 0 0 1) + +-- | Horizontal reflection +reflectionH :: Num r => Transformation 2 r +reflectionH = Transformation . Matrix $ Vector3 (Vector3 (-1) 0 0) + (Vector3 0 1 0) + (Vector3 0 0 1) diff --git a/src/Data/Geometry/Triangle.hs b/src/Data/Geometry/Triangle.hs --- a/src/Data/Geometry/Triangle.hs +++ b/src/Data/Geometry/Triangle.hs @@ -182,6 +182,8 @@ type instance IntersectionOf (Line 2 r) (Triangle 2 p r) = [ NoIntersection, Point 2 r, LineSegment 2 () r ] +instance (Fractional r, Ord r) => Line 2 r `HasIntersectionWith` Triangle 2 p r + instance (Fractional r, Ord r) => Line 2 r `IsIntersectableWith` Triangle 2 p r where nonEmptyIntersection = defaultNonEmptyIntersection @@ -206,6 +208,8 @@ type instance IntersectionOf (Line 3 r) (Triangle 3 p r) = [ NoIntersection, Point 3 r, LineSegment 3 () r ] + +instance (Fractional r, Ord r) => Line 3 r `HasIntersectionWith` Triangle 3 p r {- HLINT ignore "Use const" -} instance (Fractional r, Ord r) => Line 3 r `IsIntersectableWith` Triangle 3 p r where diff --git a/src/Data/PlaneGraph.hs b/src/Data/PlaneGraph.hs --- a/src/Data/PlaneGraph.hs +++ b/src/Data/PlaneGraph.hs @@ -11,7 +11,8 @@ -- embedding. -- -------------------------------------------------------------------------------- -module Data.PlaneGraph( PlaneGraph(PlaneGraph), graph +module Data.PlaneGraph( -- $setup + PlaneGraph(PlaneGraph), graph , PlanarGraph , VertexData(VertexData), vData, location, vtxDataToExt , fromSimplePolygon, fromConnectedSegments @@ -31,7 +32,7 @@ , incidentEdges, incomingEdges, outgoingEdges , neighboursOf , nextIncidentEdge, prevIncidentEdge - + , nextIncidentEdgeFrom, prevIncidentEdgeFrom , leftFace, rightFace , nextEdge, prevEdge @@ -44,15 +45,93 @@ , vertexData, faceData, dartData, rawDartData , edgeSegment, edgeSegments - , rawFacePolygon, rawFaceBoundary - , rawFacePolygons + , faceBoundary, internalFacePolygon + , outerFacePolygon, outerFacePolygon' + , facePolygons, facePolygons' , VertexId(..), FaceId(..), Dart, World(..), VertexId', FaceId' - , withEdgeDistances , writePlaneGraph, readPlaneGraph ) where import Data.PlaneGraph.IO import Data.PlaneGraph.Core + + +-------------------------------------------------------------------------------- + +-- $setup +-- >>> import Data.Proxy +-- >>> import Data.PlaneGraph.AdjRep(Gr(Gr),Face(Face),Vtx(Vtx)) +-- >>> import Data.PlaneGraph.IO(fromAdjRep) +-- >>> import Data.PlanarGraph.Dart(Dart(..),Arc(..)) +-- >>> :{ +-- let dart i s = Dart (Arc i) (read s) +-- small :: Gr (Vtx Int String Int) (Face String) +-- small = Gr [ Vtx 0 (Point2 0 0) [ (2,"0->2") +-- , (1,"0->1") +-- , (3,"0->3") +-- ] 0 +-- , Vtx 1 (Point2 2 2) [ (0,"1->0") +-- , (2,"1->2") +-- , (3,"1->3") +-- ] 1 +-- , Vtx 2 (Point2 2 0) [ (0,"2->0") +-- , (1,"2->1") +-- ] 2 +-- , Vtx 3 (Point2 (-1) 4) [ (0,"3->0") +-- , (1,"3->1") +-- ] 3 +-- ] +-- [ Face (2,1) "OuterFace" +-- , Face (0,1) "A" +-- , Face (1,0) "B" +-- ] +-- smallG = fromAdjRep (Proxy :: Proxy ()) small +-- :} +-- +-- +-- This represents the following graph. Note that the graph is undirected, the +-- arrows are just to indicate what the Positive direction of the darts is. +-- +-- ![myGraph](docs/Data/PlaneGraph/small.png) +-- +-- +-- Here is also a slightly larger example graph: +-- ![myGraph](docs/Data/PlaneGraph/planegraph.png) +-- +-- >>> import Data.RealNumber.Rational +-- >>> data MyWorld +-- >>> :{ +-- let myPlaneGraph :: PlaneGraph MyWorld Int () String (RealNumber 5) +-- myPlaneGraph = fromAdjRep (Proxy @MyWorld) myPlaneGraphAdjrep +-- myPlaneGraphAdjrep :: Gr (Vtx Int () (RealNumber 5)) (Face String) +-- myPlaneGraphAdjrep = Gr [ vtx 0 (Point2 0 0 ) [e 9, e 5, e 1, e 2] +-- , vtx 1 (Point2 4 4 ) [e 0, e 5, e 12] +-- , vtx 2 (Point2 3 7 ) [e 0, e 3] +-- , vtx 3 (Point2 0 5 ) [e 4, e 2] +-- , vtx 4 (Point2 3 8 ) [e 3, e 13] +-- , vtx 5 (Point2 8 1 ) [e 0, e 6, e 8, e 1] +-- , vtx 6 (Point2 6 (-1)) [e 5, e 9] +-- , vtx 7 (Point2 9 (-1)) [e 8, e 11] +-- , vtx 8 (Point2 12 1 ) [e 7, e 12, e 5] +-- , vtx 9 (Point2 8 (-5)) [e 0, e 10, e 6] +-- , vtx 10 (Point2 12 (-3)) [e 9, e 11] +-- , vtx 11 (Point2 14 (-1)) [e 10, e 7] +-- , vtx 12 (Point2 10 4 ) [e 1, e 8, e 13, e 14] +-- , vtx 13 (Point2 9 6 ) [e 4, e 14, e 12] +-- , vtx 14 (Point2 8 5 ) [e 13, e 12] +-- ] +-- [ Face (0,9) "OuterFace" +-- , Face (0,5) "A" +-- , Face (0,1) "B" +-- , Face (0,2) "C" +-- , Face (14,13) "D" +-- , Face (1,12) "E" +-- , Face (5,8) "F" +-- ] +-- where +-- e i = (i,()) +-- vtx i p es = Vtx i p es i +-- :} diff --git a/src/Data/PlaneGraph/AdjRep.hs b/src/Data/PlaneGraph/AdjRep.hs --- a/src/Data/PlaneGraph/AdjRep.hs +++ b/src/Data/PlaneGraph/AdjRep.hs @@ -27,7 +27,7 @@ -- edge. Adjacencies are given in -- arbitrary order , vData :: v - } deriving (Generic, Functor) + } deriving (Generic, Show, Eq, Functor) instance (ToJSON r, ToJSON v, ToJSON e) => ToJSON (Vtx v e r) where toEncoding = genericToEncoding defaultOptions diff --git a/src/Data/PlaneGraph/Core.hs b/src/Data/PlaneGraph/Core.hs --- a/src/Data/PlaneGraph/Core.hs +++ b/src/Data/PlaneGraph/Core.hs @@ -25,7 +25,7 @@ , vertices', vertices , edges', edges - , faces', faces, internalFaces, faces'' + , faces', internalFaces', faces, internalFaces, faces'' , darts', darts , traverseVertices, traverseDarts, traverseFaces @@ -34,6 +34,7 @@ , incidentEdges, incomingEdges, outgoingEdges , neighboursOf , nextIncidentEdge, prevIncidentEdge + , nextIncidentEdgeFrom, prevIncidentEdgeFrom , leftFace, rightFace @@ -47,41 +48,43 @@ , vertexData, faceData, dartData, rawDartData , edgeSegment, edgeSegments - , rawFacePolygon, rawFaceBoundary - , rawFacePolygons + , faceBoundary, internalFacePolygon + , outerFacePolygon, outerFacePolygon' + , facePolygons, facePolygons', internalFacePolygons , VertexId(..), FaceId(..), Dart, World(..), VertexId', FaceId' - , withEdgeDistances -- , writePlaneGraph, readPlaneGraph ) where -import Control.Lens hiding (holes, holesOf, (.=)) +import Control.Lens hiding (holes, holesOf, (.=)) import Data.Aeson +import Data.Bifunctor (first) import Data.Ext -import qualified Data.Foldable as F -import Data.Function (on) +import qualified Data.Foldable as F +import Data.Function (on) import Data.Geometry.Box import Data.Geometry.Interval -import Data.Geometry.Line (cmpSlope, supportingLine) +import Data.Geometry.Line (cmpSlope, supportingLine) import Data.Geometry.LineSegment hiding (endPoints) import Data.Geometry.Point import Data.Geometry.Polygon import Data.Geometry.Properties -import qualified Data.List.NonEmpty as NonEmpty -import qualified Data.Map as M -import Data.Ord (comparing) +import qualified Data.List.NonEmpty as NonEmpty +import qualified Data.Map as M +import Data.Ord (comparing) import Data.PlanarGraph (Arc (..), Dart (..), Direction (..), FaceId (..), FaceId', HasDataOf (..), PlanarGraph, VertexId (..), VertexId', World (..), dual, planarGraph, twin) -import qualified Data.PlanarGraph as PG +import qualified Data.PlanarGraph as PG import Data.Util -import qualified Data.Vector as V -import Data.Vector.Circular (CircularVector) -import GHC.Generics (Generic) +import qualified Data.Vector as V +import Data.Vector.Circular (CircularVector) +import GHC.Generics (Generic) + -------------------------------------------------------------------------------- -- $setup @@ -119,7 +122,45 @@ -- arrows are just to indicate what the Positive direction of the darts is. -- -- ![myGraph](docs/Data/PlaneGraph/small.png) - +-- +-- +-- Here is also a slightly larger example graph: +-- ![myGraph](docs/Data/PlaneGraph/planegraph.png) +-- +-- >>> import Data.RealNumber.Rational +-- >>> data MyWorld +-- >>> :{ +-- let myPlaneGraph :: PlaneGraph MyWorld Int () String (RealNumber 5) +-- myPlaneGraph = fromAdjRep (Proxy @MyWorld) myPlaneGraphAdjrep +-- myPlaneGraphAdjrep :: Gr (Vtx Int () (RealNumber 5)) (Face String) +-- myPlaneGraphAdjrep = Gr [ vtx 0 (Point2 0 0 ) [e 9, e 5, e 1, e 2] +-- , vtx 1 (Point2 4 4 ) [e 0, e 5, e 12] +-- , vtx 2 (Point2 3 7 ) [e 0, e 3] +-- , vtx 3 (Point2 0 5 ) [e 4, e 2] +-- , vtx 4 (Point2 3 8 ) [e 3, e 13] +-- , vtx 5 (Point2 8 1 ) [e 0, e 6, e 8, e 1] +-- , vtx 6 (Point2 6 (-1)) [e 5, e 9] +-- , vtx 7 (Point2 9 (-1)) [e 8, e 11] +-- , vtx 8 (Point2 12 1 ) [e 7, e 12, e 5] +-- , vtx 9 (Point2 8 (-5)) [e 0, e 10, e 6] +-- , vtx 10 (Point2 12 (-3)) [e 9, e 11] +-- , vtx 11 (Point2 14 (-1)) [e 10, e 7] +-- , vtx 12 (Point2 10 4 ) [e 1, e 8, e 13, e 14] +-- , vtx 13 (Point2 9 6 ) [e 4, e 14, e 12] +-- , vtx 14 (Point2 8 5 ) [e 13, e 12] +-- ] +-- [ Face (0,9) "OuterFace" +-- , Face (0,5) "A" +-- , Face (0,1) "B" +-- , Face (0,2) "C" +-- , Face (14,13) "D" +-- , Face (1,12) "E" +-- , Face (5,8) "F" +-- ] +-- where +-- e i = (i,()) +-- vtx i p es = Vtx i p es i +-- :} -------------------------------------------------------------------------------- -- * Vertex Data @@ -131,6 +172,7 @@ ,Functor,Foldable,Traversable) makeLenses ''VertexData +-- | Convert to an Ext vtxDataToExt :: VertexData r v -> Point 2 r :+ v vtxDataToExt (VertexData p v) = p :+ v @@ -235,6 +277,8 @@ -- -- >>> numVertices smallG -- 4 +-- >>> numVertices myPlaneGraph +-- 15 numVertices :: PlaneGraph s v e f r -> Int numVertices = PG.numVertices . _graph @@ -242,6 +286,7 @@ -- -- >>> numDarts smallG -- 10 +-- numDarts :: PlaneGraph s v e f r -> Int numDarts = PG.numDarts . _graph @@ -256,6 +301,8 @@ -- -- >>> numFaces smallG -- 3 +-- >>> numFaces myPlaneGraph +-- 7 numFaces :: PlaneGraph s v e f r -> Int numFaces = PG.numFaces . _graph @@ -281,9 +328,12 @@ darts' = PG.darts' . _graph -- | Get all darts together with their data +-- +-- darts :: PlaneGraph s v e f r -> V.Vector (Dart s, e) darts = PG.darts . _graph + -- | Enumerate all edges. We report only the Positive darts edges' :: PlaneGraph s v e f r -> V.Vector (Dart s) edges' = PG.edges' . _graph @@ -327,16 +377,30 @@ faces' :: PlaneGraph s v e f r -> V.Vector (FaceId' s) faces' = PG.faces' . _graph + +-- | face Ids of all internal faces in the plane graph +-- +-- running time: \(O(n)\) +internalFaces' :: (Ord r, Fractional r) => PlaneGraph s v e f r -> V.Vector (FaceId' s) +internalFaces' g = let i = outerFaceId g in V.filter (/= i) $ faces' g + -- | All faces with their face data. -- -- >>> mapM_ print $ faces smallG -- (FaceId 0,"OuterFace") -- (FaceId 1,"A") -- (FaceId 2,"B") +-- >>> mapM_ print $ faces myPlaneGraph +-- (FaceId 0,"OuterFace") +-- (FaceId 1,"A") +-- (FaceId 2,"B") +-- (FaceId 3,"C") +-- (FaceId 4,"E") +-- (FaceId 5,"F") +-- (FaceId 6,"D") faces :: PlaneGraph s v e f r -> V.Vector (FaceId' s, f) faces = PG.faces . _graph - -- | Reports the outerface and all internal faces separately. -- running time: \(O(n)\) faces'' :: (Ord r, Fractional r) @@ -384,6 +448,11 @@ -- -- >>> incidentEdges (VertexId 1) smallG -- [Dart (Arc 1) -1,Dart (Arc 4) +1,Dart (Arc 3) +1] +-- >>> mapM_ print $ incidentEdges (VertexId 5) myPlaneGraph +-- Dart (Arc 1) -1 +-- Dart (Arc 7) +1 +-- Dart (Arc 10) +1 +-- Dart (Arc 4) -1 incidentEdges :: VertexId' s -> PlaneGraph s v e f r -> V.Vector (Dart s) incidentEdges v = PG.incidentEdges v . _graph @@ -417,31 +486,69 @@ -- -- >>> neighboursOf (VertexId 1) smallG -- [VertexId 0,VertexId 2,VertexId 3] +-- >>> neighboursOf (VertexId 5) myPlaneGraph +-- [VertexId 0,VertexId 6,VertexId 8,VertexId 1] neighboursOf :: VertexId' s -> PlaneGraph s v e f r -> V.Vector (VertexId' s) neighboursOf v = PG.neighboursOf v . _graph -- | Given a dart d that points into some vertex v, report the next dart in the --- cyclic order around v in clockwise direction. +-- cyclic (counterclockwise) order around v. -- -- running time: \(O(1)\) -- -- >>> nextIncidentEdge (dart 1 "+1") smallG --- Dart (Arc 2) +1 +-- Dart (Arc 4) +1 +-- >>> nextIncidentEdge (dart 1 "+1") myPlaneGraph +-- Dart (Arc 7) +1 +-- >>> nextIncidentEdge (dart 17 "-1") myPlaneGraph +-- Dart (Arc 15) -1 nextIncidentEdge :: Dart s -> PlaneGraph s v e f r -> Dart s nextIncidentEdge d = PG.nextIncidentEdge d . _graph --- | Given a dart d that points into some vertex v, report the next dart in the --- cyclic order around v (in clockwise order) +-- | Given a dart d that points into some vertex v, report the previous dart in the +-- cyclic (counterclockwise) order around v. -- -- running time: \(O(1)\) -- -- >>> prevIncidentEdge (dart 1 "+1") smallG --- Dart (Arc 0) +1 +-- Dart (Arc 3) +1 +-- >>> prevIncidentEdge (dart 1 "+1") myPlaneGraph +-- Dart (Arc 4) -1 +-- >>> prevIncidentEdge (dart 7 "-1") myPlaneGraph +-- Dart (Arc 1) -1 prevIncidentEdge :: Dart s -> PlaneGraph s v e f r -> Dart s prevIncidentEdge d = PG.prevIncidentEdge d . _graph +-- | Given a dart d that points away from some vertex v, report the +-- next dart in the cyclic (counterclockwise) order around v. +-- +-- +-- running time: \(O(1)\) +-- +-- >>> nextIncidentEdgeFrom (dart 1 "+1") smallG +-- Dart (Arc 2) +1 +-- >>> nextIncidentEdgeFrom (dart 1 "+1") myPlaneGraph +-- Dart (Arc 2) +1 +-- >>> nextIncidentEdgeFrom (dart 4 "+1") myPlaneGraph +-- Dart (Arc 15) +1 +nextIncidentEdgeFrom :: Dart s -> PlaneGraph s v e f r -> Dart s +nextIncidentEdgeFrom d = PG.nextIncidentEdgeFrom d . _graph + +-- | Given a dart d that points into away from vertex v, report the previous dart in the +-- cyclic (counterclockwise) order around v. +-- +-- running time: \(O(1)\) +-- +-- >>> prevIncidentEdgeFrom (dart 1 "+1") smallG +-- Dart (Arc 0) +1 +-- >>> prevIncidentEdgeFrom (dart 4 "+1") myPlaneGraph +-- Dart (Arc 2) -1 +prevIncidentEdgeFrom :: Dart s -> PlaneGraph s v e f r -> Dart s +prevIncidentEdgeFrom d = PG.prevIncidentEdgeFrom d . _graph + + -- | The face to the left of the dart -- -- running time: \(O(1)\). @@ -494,20 +601,34 @@ prevEdge d = PG.prevEdge d . _graph --- | The darts bounding this face, for internal faces in clockwise order, for --- the outer face in counter clockwise order. --- +-- | The darts bounding this face. The darts are reported in order +-- along the face. This means that for internal faces the darts are +-- reported in *clockwise* order along the boundary, whereas for the +-- outer face the darts are reported in counter clockwise order. -- -- running time: \(O(k)\), where \(k\) is the output size. -- +-- >>> boundary (FaceId $ VertexId 2) smallG -- around face B +-- [Dart (Arc 2) +1,Dart (Arc 3) -1,Dart (Arc 1) -1] +-- >>> boundary (FaceId $ VertexId 0) smallG -- around outer face +-- [Dart (Arc 0) +1,Dart (Arc 4) -1,Dart (Arc 3) +1,Dart (Arc 2) -1] -- boundary :: FaceId' s -> PlaneGraph s v e f r -> V.Vector (Dart s) boundary f = PG.boundary f . _graph --- | Generates the darts incident to a face, starting with the given dart. +-- | Given a dart d, generates the darts bounding the face that is to +-- the right of the given dart. The darts are reported in order along +-- the face. This means that for internal faces the darts are reported +-- in *clockwise* order along the boundary, whereas for the outer face +-- the darts are reported in counter clockwise order. -- +-- running time: \(O(k)\), where \(k\) is the number of darts reported -- --- \(O(k)\), where \(k\) is the number of darts reported +-- >>> boundary' (dart 2 "+1") smallG -- around face B +-- [Dart (Arc 2) +1,Dart (Arc 3) -1,Dart (Arc 1) -1] +-- >>> boundary' (dart 0 "+1") smallG -- around outer face +-- [Dart (Arc 0) +1,Dart (Arc 4) -1,Dart (Arc 3) +1,Dart (Arc 2) -1] +-- boundary' :: Dart s -> PlaneGraph s v e f r -> V.Vector (Dart s) boundary' d = PG.boundary' d . _graph @@ -519,8 +640,24 @@ -- | The vertices bounding this face, for internal faces in clockwise order, for -- the outer face in counter clockwise order. -- --- -- running time: \(O(k)\), where \(k\) is the output size. +-- +-- >>> boundaryVertices (FaceId $ VertexId 2) smallG -- around B +-- [VertexId 0,VertexId 3,VertexId 1] +-- >>> boundaryVertices (FaceId $ VertexId 0) smallG -- around outerface +-- [VertexId 0,VertexId 2,VertexId 1,VertexId 3] +-- >>> mapM_ print $ boundaryVertices (FaceId $ VertexId 0) myPlaneGraph +-- VertexId 0 +-- VertexId 9 +-- VertexId 10 +-- VertexId 11 +-- VertexId 7 +-- VertexId 8 +-- VertexId 12 +-- VertexId 13 +-- VertexId 4 +-- VertexId 3 +-- VertexId 2 boundaryVertices :: FaceId' s -> PlaneGraph s v e f r -> V.Vector (VertexId' s) boundaryVertices f = PG.boundaryVertices f . _graph @@ -529,9 +666,18 @@ -------------------------------------------------------------------------------- -- * Access data + +-- | Lens to access the vertex data +-- +-- Note that using the setting part of this lens may be very +-- expensive!! (O(n)) vertexDataOf :: VertexId' s -> Lens' (PlaneGraph s v e f r ) (VertexData r v) vertexDataOf v = graph.PG.dataOf v +-- | Get the location of a vertex in the plane graph +-- +-- Note that the setting part of this lens may be very expensive! +-- Moreover, use with care (as this may destroy planarity etc.) locationOf :: VertexId' s -> Lens' (PlaneGraph s v e f r ) (Point 2 r) locationOf v = vertexDataOf v.location @@ -652,6 +798,27 @@ -- (Dart (Arc 2) +1,ClosedLineSegment (Point2 0 0 :+ 0) (Point2 (-1) 4 :+ 3) :+ "0->3") -- (Dart (Arc 4) +1,ClosedLineSegment (Point2 2 2 :+ 1) (Point2 2 0 :+ 2) :+ "1->2") -- (Dart (Arc 3) +1,ClosedLineSegment (Point2 2 2 :+ 1) (Point2 (-1) 4 :+ 3) :+ "1->3") +-- >>> mapM_ print $ edgeSegments myPlaneGraph +-- (Dart (Arc 0) +1,ClosedLineSegment (Point2 0 0 :+ 0) (Point2 8 (-5) :+ 9) :+ ()) +-- (Dart (Arc 1) +1,ClosedLineSegment (Point2 0 0 :+ 0) (Point2 8 1 :+ 5) :+ ()) +-- (Dart (Arc 2) +1,ClosedLineSegment (Point2 0 0 :+ 0) (Point2 4 4 :+ 1) :+ ()) +-- (Dart (Arc 3) +1,ClosedLineSegment (Point2 0 0 :+ 0) (Point2 3 7 :+ 2) :+ ()) +-- (Dart (Arc 4) +1,ClosedLineSegment (Point2 4 4 :+ 1) (Point2 8 1 :+ 5) :+ ()) +-- (Dart (Arc 15) +1,ClosedLineSegment (Point2 4 4 :+ 1) (Point2 10 4 :+ 12) :+ ()) +-- (Dart (Arc 5) +1,ClosedLineSegment (Point2 3 7 :+ 2) (Point2 0 5 :+ 3) :+ ()) +-- (Dart (Arc 6) +1,ClosedLineSegment (Point2 0 5 :+ 3) (Point2 3 8 :+ 4) :+ ()) +-- (Dart (Arc 18) +1,ClosedLineSegment (Point2 3 8 :+ 4) (Point2 9 6 :+ 13) :+ ()) +-- (Dart (Arc 7) +1,ClosedLineSegment (Point2 8 1 :+ 5) (Point2 6 (-1) :+ 6) :+ ()) +-- (Dart (Arc 10) +1,ClosedLineSegment (Point2 8 1 :+ 5) (Point2 12 1 :+ 8) :+ ()) +-- (Dart (Arc 12) +1,ClosedLineSegment (Point2 6 (-1) :+ 6) (Point2 8 (-5) :+ 9) :+ ()) +-- (Dart (Arc 8) +1,ClosedLineSegment (Point2 9 (-1) :+ 7) (Point2 12 1 :+ 8) :+ ()) +-- (Dart (Arc 14) +1,ClosedLineSegment (Point2 9 (-1) :+ 7) (Point2 14 (-1) :+ 11) :+ ()) +-- (Dart (Arc 9) +1,ClosedLineSegment (Point2 12 1 :+ 8) (Point2 10 4 :+ 12) :+ ()) +-- (Dart (Arc 11) +1,ClosedLineSegment (Point2 8 (-5) :+ 9) (Point2 12 (-3) :+ 10) :+ ()) +-- (Dart (Arc 13) +1,ClosedLineSegment (Point2 12 (-3) :+ 10) (Point2 14 (-1) :+ 11) :+ ()) +-- (Dart (Arc 16) +1,ClosedLineSegment (Point2 10 4 :+ 12) (Point2 9 6 :+ 13) :+ ()) +-- (Dart (Arc 17) +1,ClosedLineSegment (Point2 10 4 :+ 12) (Point2 8 5 :+ 14) :+ ()) +-- (Dart (Arc 19) +1,ClosedLineSegment (Point2 9 6 :+ 13) (Point2 8 5 :+ 14) :+ ()) edgeSegments :: PlaneGraph s v e f r -> V.Vector (Dart s, LineSegment 2 v r :+ e) edgeSegments ps = fmap withSegment . edges $ ps where @@ -671,31 +838,87 @@ seg = let (p,q) = bimap vtxDataToExt vtxDataToExt $ ps^.endPointsOf d in ClosedLineSegment p q --- | The polygon describing the face + +-- | The boundary of the face as a simple polygon. For internal faces +-- the polygon that is reported has its vertices stored in CCW order +-- (as expected). -- --- runningtime: \(O(k)\), where \(k\) is the size of the face. +-- pre: FaceId refers to an internal face. -- +-- For the other face this prodcuces a polygon in CW order (this may +-- lead to unexpected results.) -- -rawFaceBoundary :: FaceId' s -> PlaneGraph s v e f r - -> SimplePolygon v r :+ f -rawFaceBoundary i ps = pg :+ (ps^.dataOf i) +-- runningtime: \(O(k)\), where \(k\) is the size of the face. +faceBoundary :: FaceId' s -> PlaneGraph s v e f r -> SimplePolygon v r :+ f +faceBoundary i ps = pg :+ (ps^.dataOf i) where - pg = unsafeFromPoints . F.toList . fmap (\j -> ps^.graph.dataOf j.to vtxDataToExt) + pg = unsafeFromVector . V.reverse . fmap (\j -> ps^.graph.dataOf j.to vtxDataToExt) . boundaryVertices i $ ps + -- polygons are stored in CCW order, the boundaryVertices of + -- internal faces are reported in CW order we reverse them. --- | Alias for rawFace Boundary +-------------------------------------------------------------------------------- + +-- | The boundary of the face as a simple polygon. For internal faces +-- the polygon that is reported has its vertices stored in CCW order +-- (as expected). -- +-- pre: FaceId refers to an internal face. +-- +-- For the other face this prodcuces a polygon in CW order (this may +-- lead to unexpected results.) +-- -- runningtime: \(O(k)\), where \(k\) is the size of the face. -rawFacePolygon :: FaceId' s -> PlaneGraph s v e f r -> SimplePolygon v r :+ f -rawFacePolygon = rawFaceBoundary +internalFacePolygon :: FaceId' s -> PlaneGraph s v e f r -> SimplePolygon v r :+ f +internalFacePolygon = faceBoundary --- | Lists all faces of the plane graph. -rawFacePolygons :: PlaneGraph s v e f r - -> V.Vector (FaceId' s, SimplePolygon v r :+ f) -rawFacePolygons ps = fmap (\i -> (i,rawFacePolygon i ps)) . faces' $ ps +-- | Given the outerFaceId and the graph, construct a sufficiently +-- large rectangular multipolygon ith a hole containing the boundary +-- of the outer face. +outerFacePolygon :: (Num r, Ord r) + => FaceId' s -> PlaneGraph s v e f r -> MultiPolygon (Maybe v) r :+ f +outerFacePolygon i pg = + outerFacePolygon' i outer pg & core %~ first (either (const Nothing) Just) + where + outer = rectToPolygon . grow 1 . boundingBox $ pg + rectToPolygon = unsafeFromPoints . reverse . F.toList . corners +-- | Given the outerface id, and a sufficiently large outer boundary, +-- draw the outerface as a polygon with a hole. +outerFacePolygon' :: FaceId' s -> SimplePolygon v' r + -> PlaneGraph s v e f r -> MultiPolygon (Either v' v) r :+ f +outerFacePolygon' i outer pg = MultiPolygon (first Left outer) [hole] :+ pg^.dataOf i + where + hole = reverseOuterBoundary . first Right . view core $ faceBoundary i pg + -- if we call faceBoundary on the outerface we get a polygon in + -- the wrong orientation. So reverse it. + -------------------------------------------------------------------------------- +-- | Given the outerFace Id, construct polygons for all faces. We +-- construct a polygon with a hole for the outer face. +-- +facePolygons :: (Num r, Ord r) => FaceId' s -> PlaneGraph s v e f r + -> ( (FaceId' s, MultiPolygon (Maybe v) r :+ f) + , V.Vector (FaceId' s, SimplePolygon v r :+ f) + ) +facePolygons i ps = ((i, outerFacePolygon i ps), facePolygons' i ps) + +-- | Given the outerFace Id, lists all internal faces of the plane +-- graph with their boundaries. +facePolygons' :: FaceId' s -> PlaneGraph s v e f r + -> V.Vector (FaceId' s, SimplePolygon v r :+ f) +facePolygons' i ps = fmap (\j -> (j,internalFacePolygon j ps)) . V.filter (/= i) . faces' $ ps + + +-- | lists all internal faces of the plane graph with their +-- boundaries. +internalFacePolygons :: (Ord r, Fractional r) + => PlaneGraph s v e f r -> V.Vector (FaceId' s, SimplePolygon v r :+ f) +internalFacePolygons pg = facePolygons' (outerFaceId pg) pg + +-------------------------------------------------------------------------------- + -- | Labels the edges of a plane graph with their distances, as specified by -- the distance function. withEdgeDistances :: (Point 2 r -> Point 2 r -> a) @@ -703,3 +926,7 @@ withEdgeDistances f g = g&graph.PG.dartData %~ fmap (\(d,x) -> (d,len d :+ x)) where len d = uncurry f . over both (^.location) $ endPointData d g + + + +-------------------------------------------------------------------------------- diff --git a/src/Data/PlaneGraph/IO.hs b/src/Data/PlaneGraph/IO.hs --- a/src/Data/PlaneGraph/IO.hs +++ b/src/Data/PlaneGraph/IO.hs @@ -21,13 +21,18 @@ import qualified Data.PlanarGraph.AdjRep as PGA import qualified Data.PlanarGraph.IO as PGIO import Data.PlaneGraph.Core -import Data.PlaneGraph.AdjRep (Face,Vtx(Vtx),Gr(Gr)) +import Data.PlaneGraph.AdjRep import Data.Proxy import qualified Data.Vector as V import qualified Data.Vector.Mutable as MV import Data.Yaml (ParseException) import Data.Yaml.Util + +import Data.RealNumber.Rational +-- import Data.PlanarGraph.Dart +-- import Data.PlaneGraph.AdjRep + -------------------------------------------------------------------------------- -- $setup @@ -68,10 +73,10 @@ -- * Reading and Writing the Plane Graph -- | Reads a plane graph from a bytestring -readPlaneGraph :: (FromJSON v, FromJSON e, FromJSON f, FromJSON r) - => proxy s -> B.ByteString - -> Either ParseException (PlaneGraph s v e f r) -readPlaneGraph _ = decodeYaml +readPlaneGraph :: forall s v e f r. (FromJSON v, FromJSON e, FromJSON f, FromJSON r) + => B.ByteString + -> Either ParseException (PlaneGraph s v e f r) +readPlaneGraph = decodeYaml -- | Writes a plane graph to a bytestring writePlaneGraph :: (ToJSON v, ToJSON e, ToJSON f, ToJSON r) @@ -132,3 +137,71 @@ -- hence, no need to pick a secondary comparison -------------------------------------------------------------------------------- + +-- smallG = fromAdjRep (Proxy :: Proxy ()) small +-- where +-- small :: Gr (Vtx Int String Int) (Face String) +-- small = Gr [ Vtx 0 (Point2 0 0) [ (2,"0->2") +-- , (1,"0->1") +-- , (3,"0->3") +-- ] 0 +-- , Vtx 1 (Point2 2 2) [ (0,"1->0") +-- , (2,"1->2") +-- , (3,"1->3") +-- ] 1 +-- , Vtx 2 (Point2 2 0) [ (0,"2->0") +-- , (1,"2->1") +-- ] 2 +-- , Vtx 3 (Point2 (-1) 4) [ (0,"3->0") +-- , (1,"3->1") +-- ] 3 +-- ] +-- [ Face (2,1) "OuterFace" +-- , Face (0,1) "A" +-- , Face (1,0) "B" +-- ] + +-- dart i s = Dart (Arc i) (read s) + +data MyWorld + +-- ![myGraph](docs/Data/PlaneGraph/planegraph.png) +myPlaneGraph :: PlaneGraph MyWorld Int () String (RealNumber 5) +myPlaneGraph = fromAdjRep (Proxy @MyWorld) myPlaneGraphAdjrep + +myPlaneGraphAdjrep :: Gr (Vtx Int () (RealNumber 5)) (Face String) +myPlaneGraphAdjrep = Gr [ vtx 0 (Point2 0 0 ) [e 9, e 5, e 1, e 2] + , vtx 1 (Point2 4 4 ) [e 0, e 5, e 12] + , vtx 2 (Point2 3 7 ) [e 0, e 3] + , vtx 3 (Point2 0 5 ) [e 4, e 2] + , vtx 4 (Point2 3 8 ) [e 3, e 13] + , vtx 5 (Point2 8 1 ) [e 0, e 6, e 8, e 1] + , vtx 6 (Point2 6 (-1)) [e 5, e 9] + , vtx 7 (Point2 9 (-1)) [e 8, e 11] + , vtx 8 (Point2 12 1 ) [e 7, e 12, e 5] + , vtx 9 (Point2 8 (-5)) [e 0, e 10, e 6] + , vtx 10 (Point2 12 (-3)) [e 9, e 11] + , vtx 11 (Point2 14 (-1)) [e 10, e 7] + , vtx 12 (Point2 10 4 ) [e 1, e 8, e 13, e 14] + , vtx 13 (Point2 9 6 ) [e 4, e 14, e 12] + , vtx 14 (Point2 8 5 ) [e 13, e 12] + ] + [ Face (0,9) "OuterFace" + , Face (0,5) "A" + , Face (0,1) "B" + , Face (0,2) "C" + , Face (14,13) "D" + , Face (1,12) "E" + , Face (5,8) "F" + ] + where + e i = (i,()) + vtx i p es = Vtx i p es i + + + + +-- myPlaneGraph' :: IO (PlaneGraph MyWorld () () () (RealNumber 5)) +-- myPlaneGraph' = let err x = error $ show x +-- in either err id . readPlaneGraph +-- <$> B.readFile "docs/Data/PlaneGraph/myPlaneGraph.yaml"