diff --git a/changelog.org b/changelog.org --- a/changelog.org +++ b/changelog.org @@ -2,6 +2,40 @@ * Changelog +** 0.11 + +- Removed Functor instance from Triangle and replaced it with Bifunctor/Bifoldable/Bitraversable +- Testing if a point lies above/below a line is now in a typeclass, + moreover there now is also an instance of this typeclass for + planes. Hence, we can test if a point in R^3 lies above or below a + plane. +- Bugfixes in the incomingEdges and outgoingEdges functions in + Planar/Plane graphs and Planar subdivisions +- Added separate data types for Sides and Corners of Rectangles. +- More functionality for working with Halfspaces +- Fixed a bug in computing the intersection of overlapping + linesegments +- PolyLine.fromPoints now returns a Maybe PolyLine rather than a + Polyine. Use fromPointsUnsafe for the old behavior. +- Interval now no longer exports its constructor. Use the provided + patterns instead. +- Added an OpenLineSegment pattern/constructor +- The corners and sides functions in Box now return specific types + representing those rather than four tuples. +- Added a BezierSpline module and data type (Thanks to Maarten). +- Added a QuadTree implementation. It can be built from a set of + points, and to represent the zeroset of some function. +- Added a Naive implementation of Convex hull in R^3. Note however + that it works only for points in general position. In particular, no + four points should be coplanar. +- Added a Data.Geometry.Directions module that defines cardinal and + InterCardinal directions. +- Added an Ellipse type (mostly so that hgeometry-ipe can read + ellipses) +- Added FunctorWithIndex, FoldableWithIndex, and TraversableWithIndex + instances for Vector, and removed specifically exporting imap; we + can now just use those functions from the Lens package. + ** 0.10 - renamed the smallest enclosing ball to RIC diff --git a/doctests.hs b/doctests.hs --- a/doctests.hs +++ b/doctests.hs @@ -63,6 +63,11 @@ , "Data.Geometry.Polygon" , "Data.Geometry.Ball" , "Data.Geometry.Box" + , "Data.Geometry.HyperPlane" -- , "Algorithms.Geometry.HiddenSurfaceRemoval.HiddenSurfaceRemoval" + , "Algorithms.Geometry.ConvexHull.Naive" + , "Algorithms.Geometry.ConvexHull.JarvisMarch" + + , "Algorithms.Geometry.SoS.Orientation" ] diff --git a/hgeometry.cabal b/hgeometry.cabal --- a/hgeometry.cabal +++ b/hgeometry.cabal @@ -1,5 +1,5 @@ name: hgeometry -version: 0.10.0.0 +version: 0.11.0.0 synopsis: Geometric Algorithms, Data structures, and Data types. description: HGeometry provides some basic geometry types, and geometric algorithms and @@ -56,12 +56,17 @@ ghc-options: -O2 -Wall -fno-warn-unticked-promoted-constructors -fno-warn-type-defaults exposed-modules: + -- * Primitives; Simulating General Position + Algorithms.Geometry.SoS + Algorithms.Geometry.SoS.Symbolic + -- * Generic Geometry Data.Geometry Data.Geometry.Properties Data.Geometry.Transformation Data.Geometry.Boundary Data.Geometry.Duality + Data.Geometry.Directions -- * Basic Geometry Types Data.Geometry.Vector @@ -69,10 +74,13 @@ Data.Geometry.Vector.VectorFamily Data.Geometry.Vector.VectorFamilyPeano + Data.Geometry.Matrix + -- Data.Geometry.Vector.Vinyl Data.Geometry.Interval Data.Geometry.Interval.Util Data.Geometry.Point + Data.Geometry.Line Data.Geometry.Line.Internal Data.Geometry.LineSegment @@ -86,10 +94,17 @@ Data.Geometry.Slab Data.Geometry.Box Data.Geometry.Box.Internal + Data.Geometry.Box.Sides + Data.Geometry.Box.Corners + Data.Geometry.Ball + Data.Geometry.Ellipse + Data.Geometry.Polygon Data.Geometry.Polygon.Convex + Data.Geometry.BezierSpline + -- * Geometric Data Structures Data.Geometry.IntervalTree Data.Geometry.SegmentTree @@ -111,13 +126,21 @@ Data.Geometry.PrioritySearchTree + Data.Geometry.QuadTree + Data.Geometry.QuadTree.Cell + Data.Geometry.QuadTree.Quadrants + Data.Geometry.QuadTree.Split + Data.Geometry.QuadTree.Tree + + -- * Algorithms -- * Geometric Algorithms Algorithms.Geometry.ConvexHull.GrahamScan Algorithms.Geometry.ConvexHull.DivideAndConquer Algorithms.Geometry.ConvexHull.QuickHull - -- Algorithms.Geometry.ConvexHull.JarvisMarch + Algorithms.Geometry.ConvexHull.JarvisMarch + Algorithms.Geometry.ConvexHull.Naive Algorithms.Geometry.LowerEnvelope.DualCH @@ -174,14 +197,32 @@ Graphics.Render other-modules: + Data.Geometry.Matrix.Internal + -- * Implementation Internals of Polygons Data.Geometry.Polygon.Core Data.Geometry.Polygon.Extremes + + Data.Geometry.Point.Internal + Data.Geometry.Point.Orientation + Data.Geometry.Point.Quadrants + Data.Geometry.Point.Orientation.Degenerate + Data.Geometry.Point.Class + + Algorithms.Geometry.SoS.Expr + Algorithms.Geometry.SoS.AsPoint + Algorithms.Geometry.SoS.Internal + Algorithms.Geometry.SoS.Orientation + Algorithms.Geometry.SoS.Determinant + Algorithms.Geometry.SoS.Sign + + + -- other-extensions: build-depends: base >= 4.11 && < 5 - , hgeometry-combinatorial >= 0.10.0.0 + , hgeometry-combinatorial >= 0.11.0.0 , bifunctors >= 4.1 , bytestring >= 0.10 @@ -198,10 +239,13 @@ , deepseq >= 1.1 , fingertree >= 0.1 , MonadRandom >= 0.5 + , random >= 1.1 , QuickCheck >= 2.5 , quickcheck-instances >= 0.3 , reflection >= 2.1 , primitive >= 0.6.3.0 + , hashable >= 1.2 + -- , singleton-typelits >= 0.1.0.0 -- , ghc-typelits-natnormalise >= 0.6 @@ -220,7 +264,7 @@ , hspec, QuickCheck, quickcheck-instances - hs-source-dirs: src test + hs-source-dirs: src default-language: Haskell2010 diff --git a/src/Algorithms/Geometry/ConvexHull/DivideAndConquer.hs b/src/Algorithms/Geometry/ConvexHull/DivideAndConquer.hs --- a/src/Algorithms/Geometry/ConvexHull/DivideAndConquer.hs +++ b/src/Algorithms/Geometry/ConvexHull/DivideAndConquer.hs @@ -16,7 +16,6 @@ import Algorithms.DivideAndConquer import Control.Arrow ((&&&)) -import Control.Lens ((^.), to) import Data.Ext import Data.Geometry.Point import Data.Geometry.Polygon diff --git a/src/Algorithms/Geometry/ConvexHull/GrahamScan.hs b/src/Algorithms/Geometry/ConvexHull/GrahamScan.hs --- a/src/Algorithms/Geometry/ConvexHull/GrahamScan.hs +++ b/src/Algorithms/Geometry/ConvexHull/GrahamScan.hs @@ -1,6 +1,8 @@ module Algorithms.Geometry.ConvexHull.GrahamScan( convexHull - , upperHull - , lowerHull + , upperHull, upperHull' + , lowerHull, lowerHull' + + , upperHullFromSorted, upperHullFromSorted' ) where import Control.Lens ((^.)) @@ -23,13 +25,54 @@ in ConvexPolygon . fromPoints . reverse $ lh ++ uh -- | Computes the upper hull. The upper hull is given from left to right. +-- +-- Specifically. A pair of points defines an edge of the upper hull +-- iff all other points are strictly to the right of its supporting +-- line. +-- +-- Note that this definition implies that the segment may be +-- vertical. Use 'upperHull'' if such an edge should not be reported. +-- +-- running time: \(O(n\log n)\) upperHull :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p) upperHull = NonEmpty.reverse . hull id --- | Computes the upper hull. The upper hull is given from left to right +-- | Computes the upper hull, making sure that there are no vertical segments. +-- +-- The upper hull is given from left to right +-- +upperHull' :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p) +upperHull' = NonEmpty.reverse . dropVertical . hull id + +-- | Helper function to remove vertical segments from the hull. +-- +-- Tests if the first two points are on a vertical line, if so removes +-- the first point. +dropVertical :: Eq r => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p) +dropVertical = \case + h@(_ :| []) -> h + h@(p :| (q : rest)) | p^.core.xCoord == q^.core.xCoord -> q :| rest + | otherwise -> h + + +-- | Computes the upper hull. The upper hull is given from left to right. +-- +-- Specifically. A pair of points defines an edge of the lower hull +-- iff all other points are strictly to the left of its supporting +-- line. +-- +-- Note that this definition implies that the segment may be +-- vertical. Use 'lowerHull'' if such an edge should not be reported. +-- +-- running time: \(O(n\log n)\) lowerHull :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p) lowerHull = hull reverse +-- | Computes the lower hull, making sure there are no vertical +-- segments. (Note that the only such segment could be the first +-- segment). +lowerHull' :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p) +lowerHull' = dropVertical . hull reverse -- | Helper function so that that can compute both the upper or the lower hull, depending -- on the function f @@ -43,6 +86,32 @@ incXdecY :: Ord r => (Point 2 r) :+ p -> (Point 2 r) :+ q -> Ordering incXdecY (Point2 px py :+ _) (Point2 qx qy :+ _) = compare px qx <> compare qy py + + +-- | Given a sequence of points that is sorted on increasing +-- x-coordinate and decreasing y-coordinate, computes the upper +-- hull, in *right to left order*. +-- +-- Specifically. A pair of points defines an edge of the upper hull +-- iff all other points are strictly to the right of its supporting +-- line. +-- +-- +-- Note that In constrast to the 'upperHull' function, the result is +-- returned *from right to left* !!! +-- +-- running time: \(O(n)\). +upperHullFromSorted :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p) +upperHullFromSorted = \case + h@(_ :| []) -> h + pts -> hull' $ NonEmpty.toList pts + +-- | Computes the upper hull from a sorted input. Removes the last vertical segment. +-- +-- +-- running time: \(O(n)\). +upperHullFromSorted' :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p) +upperHullFromSorted' = dropVertical . upperHullFromSorted -- | Precondition: The list of input points is sorted diff --git a/src/Algorithms/Geometry/ConvexHull/JarvisMarch.hs b/src/Algorithms/Geometry/ConvexHull/JarvisMarch.hs new file mode 100644 --- /dev/null +++ b/src/Algorithms/Geometry/ConvexHull/JarvisMarch.hs @@ -0,0 +1,145 @@ +module Algorithms.Geometry.ConvexHull.JarvisMarch( + convexHull + + , upperHull, upperHull' + , lowerHull, lowerHull' + , steepestCcwFrom, steepestCwFrom + ) where + +import Control.Lens ((^.)) +import Data.Bifunctor +import Data.Either (either) +import Data.Ext +import Data.Foldable +import Data.Geometry.Point +import Data.Geometry.Polygon +import Data.Geometry.Polygon.Convex (ConvexPolygon(..)) +import Data.Geometry.Vector +import qualified Data.List as List +import Data.List.NonEmpty (NonEmpty(..), (<|)) +import qualified Data.List.NonEmpty as NonEmpty +import Data.Ord (comparing, Down(..)) +import Data.Semigroup.Foldable + +-------------------------------------------------------------------------------- + +-- | Compute the convexhull using JarvisMarch. The resulting polygon +-- is given in clockwise order. +-- +-- running time: \(O(nh)\), where \(n\) is the number of input points +-- and \(h\) is the complexity of the hull. +convexHull :: (Ord r, Num r) + => NonEmpty (Point 2 r :+ p) -> ConvexPolygon p r +convexHull (p :| []) = ConvexPolygon . fromPoints $ [p] +convexHull pts = ConvexPolygon . fromPoints $ uh <> reverse lh + where + lh = case NonEmpty.nonEmpty (NonEmpty.init $ lowerHull pts) of + Nothing -> [] + Just (_:|lh') -> lh' + uh = toList $ upperHull pts + + -- note that fromList is afe since ps contains at least two elements + -- where + -- SP p@(c :+ _) pts = minViewBy incXdecY ps + -- takeWhile' pf (x :| xs) = x : takeWhile pf xs + +upperHull :: (Num r, Ord r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p) +upperHull pts = repeatedly cmp steepestCwFrom s rest + where + (s:_ :+ rest) = extractMinimaBy cmp (NonEmpty.toList pts) + cmp = comparing (\(Point2 x y :+ _) -> (x, Down y)) + -- start from the topmost point that has minimum x-coord + -- also use cmp as the comparator, so that we also select the last + -- vertical segment. + +-- | Upepr hull from left to right, without any vertical segments. +upperHull' :: (Num r, Ord r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p) +upperHull' pts = pruneVertical $ repeatedly cmp steepestCwFrom s rest + where + (s:_ :+ rest) = extractMinimaBy cmp0 (NonEmpty.toList pts) + cmp0 = comparing (\(Point2 x y :+ _) -> (x, Down y)) + -- start from the topmost point that has minimum x-coord + cmp = comparing (^.core) + -- for the rest select them in normal + -- lexicographic order, this causes the last + -- vertical segment to be ignored. + +-- | Computes the lower hull, from left to right. Includes vertical +-- segments at the start. +-- +-- running time: \(O(nh)\), where \(h\) is the complexity of the hull. +lowerHull :: (Num r, Ord r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p) +lowerHull pts = pruneVertical $ repeatedly cmp steepestCcwFrom s rest + where + (s:_ :+ rest) = extractMinimaBy cmp0 (NonEmpty.toList pts) + cmp0 = comparing (\(Point2 x y :+ _) -> (x, Down y)) + -- start from the topmost point that has minimum x-coord + cmp = comparing (^.core) + -- for the rest of the comparions use the normal + -- lexicographic comparing order. + +-- | Jarvis March to compute the lower hull, without any vertical segments. +-- +-- +-- running time: \(O(nh)\), where \(h\) is the complexity of the hull. +lowerHull' :: (Num r, Ord r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p) +lowerHull' pts = pruneVertical $ repeatedly cmp steepestCcwFrom s rest + where + (s:_ :+ rest) = extractMinimaBy cmp (NonEmpty.toList pts) + cmp = comparing (^.core) + + +-- | Find the next point in counter clockwise order, i.e. the point +-- with minimum slope w.r.t. the given point. +steepestCcwFrom :: (Ord r, Num r) + => (Point 2 r :+ a) -> NonEmpty (Point 2 r :+ b) -> Point 2 r :+ b +steepestCcwFrom p = List.minimumBy (ccwCmpAroundWith (Vector2 0 (-1)) p) + +-- | Find the next point in clockwise order, i.e. the point +-- with maximum slope w.r.t. the given point. +steepestCwFrom :: (Ord r, Num r) + => (Point 2 r :+ a) -> NonEmpty (Point 2 r :+ b) -> Point 2 r :+ b +steepestCwFrom p = List.minimumBy (cwCmpAroundWith (Vector2 0 1) p) + +repeatedly :: (a -> a -> Ordering) -> (a -> NonEmpty a -> a) -> a -> [a] -> NonEmpty a +repeatedly cmp f = go + where + go m xs' = case NonEmpty.nonEmpty xs' of + Nothing -> m :| [] + Just xs -> let p = f m xs + in m <| go p (NonEmpty.filter (\x -> p `cmp` x == LT) xs) + + +-- | Removes the topmost vertical points, if they exist +pruneVertical :: Eq r => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p) +pruneVertical = either id id . foldr1With f (\q -> Left $ q:|[]) + where + f p = \case + Left (q:|qs) | p^.core.xCoord == q^.core.xCoord -> Left (p :| qs) + | otherwise -> Right (p :| q:qs) + Right pts -> Right (p <| pts) + +-- | Foldr, but start by applying some function on the rightmost +-- element to get the starting value. +foldr1With :: Foldable1 f => (a -> b -> b) -> (a -> b) -> f a -> b +foldr1With f b = go . toNonEmpty + where + go (x :| xs) = case NonEmpty.nonEmpty xs of + Nothing -> b x + Just xs' -> x `f` (go xs') + +-- | extracts all minima from the list. The result consists of the +-- list of minima, and all remaining points. Both lists are returned +-- in the order in which they occur in the input. +-- +-- >>> extractMinimaBy compare [1,2,3,0,1,2,3,0,1,2,0,2] +-- [0,0,0] :+ [2,3,1,2,3,1,2,1,2] +extractMinimaBy :: (a -> a -> Ordering) -> [a] -> [a] :+ [a] +extractMinimaBy cmp = \case + [] -> [] :+ [] + (x:xs) -> first NonEmpty.toList $ foldr (\y (mins@(m:|_) :+ rest) -> + case m `cmp` y of + LT -> mins :+ y:rest + EQ -> (y NonEmpty.<| mins) :+ rest + GT -> (y:|[]) :+ NonEmpty.toList mins <> rest + ) ((x:|[]) :+ []) xs diff --git a/src/Algorithms/Geometry/ConvexHull/Naive.hs b/src/Algorithms/Geometry/ConvexHull/Naive.hs new file mode 100644 --- /dev/null +++ b/src/Algorithms/Geometry/ConvexHull/Naive.hs @@ -0,0 +1,93 @@ +module Algorithms.Geometry.ConvexHull.Naive( ConvexHull + , lowerHull', lowerHullAll + + , isValidTriangle, upperHalfSpaceOf + ) where + +import Control.Lens +import Data.Ext +import Data.Foldable (toList) +import Data.Geometry.HalfSpace +import Data.Geometry.HyperPlane +import Data.Geometry.Line +import Data.Geometry.Point +import Data.Geometry.Triangle +import Data.Geometry.Vector +import Data.Intersection(intersects) +import qualified Data.List as List +import Data.List.NonEmpty (NonEmpty(..)) +import qualified Data.List.NonEmpty as NonEmpty +import Data.Maybe (listToMaybe, isNothing) +import Data.Util +-------------------------------------------------------------------------------- + +type ConvexHull d p r = [Triangle 3 p r] + +-- | Computes the lower hull without its vertical triangles. +-- +-- pre: The points are in general position. In particular, no four +-- points should be coplanar. +-- +-- running time: \(O(n^4)\) +lowerHull' :: forall r p. (Ord r, Fractional r, Show r) + => NonEmpty (Point 3 r :+ p) -> ConvexHull 3 p r +lowerHull' = filter (not . isVertical) . lowerHullAll + where + isVertical (Triangle p q r) = + ccw' (p&core %~ projectPoint) (q&core %~ projectPoint) (r&core %~ projectPoint) == CoLinear + +-- | Generates a set of triangles to be used to construct a complete +-- convex hull. In particular, it may contain vertical triangles. +-- +-- pre: The points are in general position. In particular, no four +-- points should be coplanar. +-- +-- running time: \(O(n^4)\) +lowerHullAll :: forall r p. (Ord r, Fractional r, Show r) + => NonEmpty (Point 3 r :+ p) -> ConvexHull 3 p r +lowerHullAll (toList -> pts) = let mkT (Three p q r) = Triangle p q r in + [ t | t <- mkT <$> uniqueTriplets pts, isNothing (isValidTriangle t pts) ] + + + +killOverlapping :: (Ord r, Fractional r) => [Triangle 3 p r] -> [Triangle 3 p r] +killOverlapping = foldr keepIfNotOverlaps [] + where + keepIfNotOverlaps t ts | any (t `overlaps`) ts = ts + | otherwise = t:ts + + +t1@(Triangle p q r) `overlaps` t2@(Triangle a b c) = upperHalfSpaceOf t1 == upperHalfSpaceOf t2 + && False + + + +-- | Tests if this is a valid triangle for the lower envelope. That +-- is, if all point lie above the plane through these points. Returns +-- a Maybe; if the result is a Nothing the triangle is valid, if not +-- it returns a counter example. +-- +-- >>> let t = (Triangle (ext origin) (ext $ Point3 1 0 0) (ext $ Point3 0 1 0)) +-- >>> isValidTriangle t [ext $ Point3 5 5 0] +-- Nothing +-- >>> let t = (Triangle (ext origin) (ext $ Point3 1 0 0) (ext $ Point3 0 1 0)) +-- >>> isValidTriangle t [ext $ Point3 5 5 (-10)] +-- Just (Point3 [5,5,-10] :+ ()) +isValidTriangle :: (Num r, Ord r) + => Triangle 3 p r -> [Point 3 r :+ q] -> Maybe (Point 3 r :+ q) +isValidTriangle t = listToMaybe . filter (\a -> not $ (a^.core) `intersects` h) + where + h = upperHalfSpaceOf t + + +-- | Computes the halfspace above the triangle. +-- +-- >>> upperHalfSpaceOf (Triangle (ext $ origin) (ext $ Point3 10 0 0) (ext $ Point3 0 10 0)) +-- HalfSpace {_boundingPlane = HyperPlane {_inPlane = Point3 [0,0,0], _normalVec = Vector3 [0,0,100]}} +upperHalfSpaceOf :: (Ord r, Num r) => Triangle 3 p r -> HalfSpace 3 r +upperHalfSpaceOf (Triangle p q r) = HalfSpace h + where + h' = from3Points (p^.core) (q^.core) (r^.core) + c = p&core.zCoord -~ 1 + h = if (c^.core) `liesBelow` h' then h' else h'&normalVec %~ ((-1) *^) + a `liesBelow` plane = (a `onSideUpDown` plane) == Below diff --git a/src/Algorithms/Geometry/DelaunayTriangulation/DivideAndConquer.hs b/src/Algorithms/Geometry/DelaunayTriangulation/DivideAndConquer.hs --- a/src/Algorithms/Geometry/DelaunayTriangulation/DivideAndConquer.hs +++ b/src/Algorithms/Geometry/DelaunayTriangulation/DivideAndConquer.hs @@ -8,21 +8,22 @@ import Control.Monad.State import Data.BinaryTree import qualified Data.CircularList as CL -import qualified Data.CircularSeq as CS import qualified Data.CircularList.Util as CU +import qualified Data.CircularSeq as CS import Data.Ext import qualified Data.Foldable as F import Data.Function (on) import Data.Geometry hiding (rotateTo) import Data.Geometry.Ball (disk, insideBall) import Data.Geometry.Polygon -import qualified Data.Geometry.Polygon.Convex as Convex import Data.Geometry.Polygon.Convex (ConvexPolygon(..), simplePolygon) +import qualified Data.Geometry.Polygon.Convex as Convex import qualified Data.IntMap.Strict as IM import qualified Data.List as L import qualified Data.List.NonEmpty as NonEmpty import qualified Data.Map as M import Data.Maybe (fromJust, fromMaybe) +import Data.Measured.Size import qualified Data.Vector as V ------------------------------------------------------------------------------- diff --git a/src/Algorithms/Geometry/LineSegmentIntersection/BentleyOttmann.hs b/src/Algorithms/Geometry/LineSegmentIntersection/BentleyOttmann.hs --- a/src/Algorithms/Geometry/LineSegmentIntersection/BentleyOttmann.hs +++ b/src/Algorithms/Geometry/LineSegmentIntersection/BentleyOttmann.hs @@ -26,7 +26,6 @@ import qualified Data.Map as M import Data.Maybe import Data.Ord (Down(..), comparing) -import Data.OrdSeq (Compare) import qualified Data.Set as SS -- status struct import qualified Data.Set.Util as SS -- status struct import qualified Data.Set as EQ -- event queue diff --git a/src/Algorithms/Geometry/LineSegmentIntersection/Types.hs b/src/Algorithms/Geometry/LineSegmentIntersection/Types.hs --- a/src/Algorithms/Geometry/LineSegmentIntersection/Types.hs +++ b/src/Algorithms/Geometry/LineSegmentIntersection/Types.hs @@ -16,6 +16,8 @@ -------------------------------------------------------------------------------- +type Compare a = a -> a -> Ordering + -- get the endpoints of a line segment endPoints' :: (HasEnd s, HasStart s) => s -> (StartCore s, EndCore s) endPoints' s = (s^.start.core,s^.end.core) diff --git a/src/Algorithms/Geometry/LinearProgramming/LP2DRIC.hs b/src/Algorithms/Geometry/LinearProgramming/LP2DRIC.hs --- a/src/Algorithms/Geometry/LinearProgramming/LP2DRIC.hs +++ b/src/Algorithms/Geometry/LinearProgramming/LP2DRIC.hs @@ -86,8 +86,8 @@ , _current :: !(Point d r) } -deriving instance (Arity d, Show r) => Show (LPState d r) -deriving instance (Arity d, Eq r) => Eq (LPState d r) +deriving instance (Arity d, Show r) => Show (LPState d r) +deriving instance (Arity d, Eq r, Fractional r) => Eq (LPState d r) obj :: Lens' (LPState d r) (Vector d r) obj = lens _obj (\(LPState _ s p) o -> LPState o s p) diff --git a/src/Algorithms/Geometry/LinearProgramming/Types.hs b/src/Algorithms/Geometry/LinearProgramming/Types.hs --- a/src/Algorithms/Geometry/LinearProgramming/Types.hs +++ b/src/Algorithms/Geometry/LinearProgramming/Types.hs @@ -26,14 +26,14 @@ | UnBounded (HalfLine d r) makePrisms ''LPSolution -deriving instance (Arity d, Show r) => Show (LPSolution d r) -deriving instance (Arity d, Eq r) => Eq (LPSolution d r) +deriving instance (Arity d, Show r) => Show (LPSolution d r) +deriving instance (Arity d, Eq r, Fractional r) => Eq (LPSolution d r) data LinearProgram d r = LinearProgram { _objective :: !(Vector d r) , _constraints :: [HalfSpace d r] } makeLenses ''LinearProgram -deriving instance Arity d => Functor (LinearProgram d) -deriving instance (Arity d, Show r) => Show (LinearProgram d r) -deriving instance (Arity d, Eq r) => Eq (LinearProgram d r) +deriving instance Arity d => Functor (LinearProgram d) +deriving instance (Arity d, Show r) => Show (LinearProgram d r) +deriving instance (Arity d, Fractional r, Eq r) => Eq (LinearProgram d r) diff --git a/src/Algorithms/Geometry/PolyLineSimplification/DouglasPeucker.hs b/src/Algorithms/Geometry/PolyLineSimplification/DouglasPeucker.hs --- a/src/Algorithms/Geometry/PolyLineSimplification/DouglasPeucker.hs +++ b/src/Algorithms/Geometry/PolyLineSimplification/DouglasPeucker.hs @@ -22,15 +22,15 @@ douglasPeucker :: (Ord r, Fractional r, Arity d) => r -> PolyLine d p r -> PolyLine d p r douglasPeucker eps pl - | dst <= (eps*eps) = fromPoints [a,b] + | dst <= (eps*eps) = fromPointsUnsafe [a,b] -- at least two points, so we are fine. | otherwise = douglasPeucker eps pref `merge` douglasPeucker eps subf where - pts = pl^.points - a = LSeq.head pts - b = LSeq.last pts - (i,dst) = maxDist pts (ClosedLineSegment a b) + pts = pl^.points + a = LSeq.head pts + b = LSeq.last pts + (i,dst) = maxDist pts (ClosedLineSegment a b) - (pref,subf) = split i pl + (pref,subf) = split i pl -------------------------------------------------------------------------------- -- * Internal functions diff --git a/src/Algorithms/Geometry/SmallestEnclosingBall/Naive.hs b/src/Algorithms/Geometry/SmallestEnclosingBall/Naive.hs --- a/src/Algorithms/Geometry/SmallestEnclosingBall/Naive.hs +++ b/src/Algorithms/Geometry/SmallestEnclosingBall/Naive.hs @@ -22,8 +22,8 @@ import Data.List (minimumBy) import Data.Function (on) import Data.Maybe (fromMaybe) -import Data.Util(STR(..),SP(..), uniquePairs, uniqueTriplets) - +import Data.Util(uniquePairs, uniqueTriplets) +import qualified Data.Util as Util -------------------------------------------------------------------------------- -- | Horrible O(n^4) implementation that simply tries all disks, checks if they @@ -38,11 +38,11 @@ pairs :: Fractional r => [Point 2 r :+ p] -> [DiskResult p r] pairs pts = [ DiskResult (fromDiameter (a^.core) (b^.core)) (Two a b) - | SP a b <- uniquePairs pts] + | Util.Two a b <- uniquePairs pts] triplets :: (Ord r, Fractional r) => [Point 2 r :+ p] -> [DiskResult p r] triplets pts = [DiskResult (disk' a b c) (Three a b c) - | STR a b c <- uniqueTriplets pts] + | Util.Three a b c <- uniqueTriplets pts] disk' :: (Ord r, Fractional r) => Point 2 r :+ p -> Point 2 r :+ p -> Point 2 r :+ p -> Disk () r diff --git a/src/Algorithms/Geometry/SmallestEnclosingBall/RIC.hs b/src/Algorithms/Geometry/SmallestEnclosingBall/RIC.hs --- a/src/Algorithms/Geometry/SmallestEnclosingBall/RIC.hs +++ b/src/Algorithms/Geometry/SmallestEnclosingBall/RIC.hs @@ -31,6 +31,7 @@ import Data.Ord (comparing) import System.Random.Shuffle (shuffle) +import Data.RealNumber.Rational import Debug.Trace -------------------------------------------------------------------------------- @@ -149,16 +150,32 @@ -------------------------------------------------------------------------------- -test :: Maybe (DiskResult () Rational) -test = smallestEnclosingDiskWithPoints p q myPts - where - p = ext $ Point2 0 (-6) - q = ext $ Point2 0 6 +-- test :: Maybe (DiskResult () Rational) +-- test = smallestEnclosingDiskWithPoints p q myPts +-- where +-- p = ext $ Point2 0 (-6) +-- q = ext $ Point2 0 6 -myPts = map ext [Point2 5 1, Point2 3 3, Point2 (-2) 2, Point2 (-4) 5] +-- myPts = map ext [Point2 5 1, Point2 3 3, Point2 (-2) 2, Point2 (-4) 5] -disk'' r = (r:+) <$> disk (p^.core) (q^.core) (r^.core) - where - p = ext $ Point2 0 (-6) - q = ext $ Point2 0 6 +-- disk'' r = (r:+) <$> disk (p^.core) (q^.core) (r^.core) +-- where +-- p = ext $ Point2 0 (-6) +-- q = ext $ Point2 0 6 + + +-- maartenBug :: DiskResult () Double +-- maartenBug = let (p:q:rest) = maartenBug' +-- in smallestEnclosingDisk' p q rest + +-- maartenBug' :: [Point 2 Double :+ ()] +-- maartenBug' = [ Point2 (7.2784424e-3) (249.23) :+ () +-- , Point2 (-5.188493 ) (249.23) :+ () +-- , Point2 (-10.382694 ) (249.23) :+ () +-- , Point2 (-15.575621 ) (249.23) :+ () +-- , Point2 (0.0 ) (249.23) :+ () +-- , Point2 (0.0 ) (239.9031) :+ () +-- , Point2 (0.0 ) (230.37791) :+ () +-- , Point2 (0.0 ) (220.67882) :+ () +-- ] diff --git a/src/Algorithms/Geometry/SoS.hs b/src/Algorithms/Geometry/SoS.hs new file mode 100644 --- /dev/null +++ b/src/Algorithms/Geometry/SoS.hs @@ -0,0 +1,241 @@ +-------------------------------------------------------------------------------- +-- | +-- Module : Algorithms.Geometry.SoS +-- Copyright : (C) Frank Staals +-- License : see the LICENSE file +-- Maintainer : Frank Staals +-- +-- Implementation of +-- Simulation of Simplicity: A Technique to Cope with Degenerate Cases in Geometric Algorithms +-- +-- By +-- Herbert Edelsbrunner and Ernst Peter Mucke +-- +-------------------------------------------------------------------------------- +module Algorithms.Geometry.SoS + ( module Algorithms.Geometry.SoS.Sign + , module Algorithms.Geometry.SoS.Orientation + , module Algorithms.Geometry.SoS.Determinant + ) where + +-- import Algorithms.Geometry.SoS.Internal +import Algorithms.Geometry.SoS.Orientation +import Algorithms.Geometry.SoS.Determinant +import Algorithms.Geometry.SoS.Sign +import Control.CanAquire +import Control.Lens +import Data.Ext +import Data.Geometry.Point.Internal +import Data.Geometry.Properties +import Data.Geometry.Vector + +-------------------------------------------------------------------------------- + +-------------------------------------------------------------------------------- + + + +-- sideTest' :: ( SoS p, Dimension p ~ 2, r ~ NumType p +-- , Eq r, Num r +-- ) => [p] -> Sign +-- sideTest' (q:p1:p2:_) = sideTest q (Vector2 p1 p2) + + + + + + +-------------------------------------------------------------------------------- + + +---------------------------------------- +-------------------------------------------------------------------------------- + +-------------------------------------------------------------------------------- + + + + + +-- instance (i `CanAquire` Point d r, Arity d) => P i d r `CanAquire` Point d (R i) where +-- aquire (P i) = Point $ pure () + + + + +-------------------------------------------------------------------------------- + + + +-------------------------------------------------------------------------------- + + +-------------------------------------------------------------------------------- + + +-- -- TODO: Remove this one +-- instance HasIndex (Point d r :+ Int) where +-- indexOf = view extra + + +-- test1 :: Sign +-- test1 = sideTest (Point1 1 :+ 0 :: Point 1 Int :+ Int) (Vector1 $ Point1 5 :+ 1) + +-- test2 :: Sign +-- test2 = sideTest (Point1 5 :+ 0 :: Point 1 Int :+ Int) (Vector1 $ Point1 5 :+ 1) + + +-- test3 :: Sign +-- test3 = sideTest (Point2 (-1) 5 :+ 0 :: Point 2 Int :+ Int) (Vector2 (Point2 0 0 :+ 1) +-- (Point2 0 10 :+ 2) +-- ) + + +-- pattern Point1 x = Point (Vector1 x) + + +-- testV :: Sign +-- testV = simulateSimplicity sideTest' [ Point2 (-1) 5 +-- , Point2 0 0 +-- , Point2 0 10 +-- ] + + + + + +-------------------------------------------------------------------------------- + + + + + + + + +-- cmpSignificance :: Ord k => Bag k -> Bag k -> Ordering +-- cmpSignificance (Bag e1) (Bag e2) = e1e2 `compare` e2e1 +-- where +-- e1e2 = fmap fst . Map.lookupMax $ e1 `Map.difference` e2 +-- e2e1 = fmap fst . Map.lookupMax $ e2 `Map.difference` e1 + + + +-- -- | Represents a Sum of terms, i.e. a value that has the form: +-- -- +-- -- \[ +-- -- \sum c \Pi_{(i,j)} \varepsilon(i,j) +-- -- \] +-- newtype Symbolic i j r = Symbolic [Term i j r] deriving (Show,Eq,Functor) + +-- instance (Ord i, Ord j, Num r) => Num (Symbolic i j r) where +-- (Symbolic ts) + (Symbolic ts') = Symbolic (ts `addTerms` ts') +-- negate = fmap negate +-- (Symbolic ts) * (Symbolic ts') = Symbolic $ multiplyTerms ts ts' +-- fromInteger x = constant (fromInteger x) +-- -- abs x | signum x == -1 = (-1)*x +-- -- | oterwise = x + +-- -- signum = undefined + + + + + + + + + + +-- -- | Adds two lists of terms +-- addTerms :: forall i j r. (Ord i, Ord j, Num r) +-- => [Term i j r] -> [Term i j r] -> [Term i j r] +-- addTerms ts ts' = (\(eps,c) -> Term c eps) <$> Map.toList m +-- where +-- m :: Map.Map (EpsFold i j) r +-- m = Map.fromListWith (+) [ (eps,c) | (Term c eps) <- ts <> ts' ] + +-- multiplyTerms :: forall i j r. (Ord i, Ord j, Num r) +-- => [Term i j r] -> [Term i j r] -> [Term i j r] +-- multiplyTerms ts ts' = (\(eps,c) -> Term c eps) <$> Map.toList m +-- where +-- m :: Map.Map (EpsFold i j) r +-- m = Map.fromListWith (+) [ (es <> es',c*d) | (Term c es) <- ts, (Term d es') <- ts' ] + + + + +-- orderedTerms :: (Ord i, Ord j) => Symbolic i j r -> [Term i j r] +-- orderedTerms (Symbolic ts) = List.sortBy (\(Term _ e1) (Term _ e2) -> cmpSignificance e1 e2) ts + + + + + + + + + + + + + + + + + + -- zipWith (\j x -> Term x $ singleton (i,j)) [0..] . toList + + + + + + +-- orderTerms :: (Ord i, Ord j) => Symbolic i j r -> Symbolic i j r +-- orderTerms (Symbolic ts) = Symbolic $ List.sortBy cmpSignificance ts + + + +-- fromPoint' :: Foldable f => i -> f r -> Symbolic i Int r +-- fromPoint' i = Symbolic . zipWith (\j x -> Term x [(i,j)]) [0..] . toList + + + +-- testZ :: Symbolic Int Int Int +-- testZ = (5 + 6) * + + + + + + -- case sign i of + -- (-1) -> Negative $ fromInteger i + -- 0 -> Zero + -- _ -> Positive $ fromInteger i + -- negate = \case + -- Negative c -> Positive c + -- Positive c -> Negative c + + +-- newtype N = N String deriving (Show,Eq) + + +-- instance Num N where +-- (N x) + (N y) = N $ x <> "+" <> y +-- (N x) * (N y) = N $ x <> y +-- negate (N x) = N ("negate(" <> x <> ")") +-- fromInteger = N . show + + +-- n :: (Ord i, Ord j) => String -> i -> j -> Symbolic i j N +-- n x i j = Symbolic [Term (N x) mempty, Term 1 (singleton (i,j))] + + + + + +-- testM3 = det33 $ V3 (fromPoint' [N "px", N "py"] <> 1) +-- (fromPoint' [N "px", N "py"] <> 1) +-- (fromPoint' [N "px", N "py"] <> 1) +-- -- (V3 (N "qx") (N "qy") 1) +-- -- (V3 (N "rx") (N "ry") 1) diff --git a/src/Algorithms/Geometry/SoS/AsPoint.hs b/src/Algorithms/Geometry/SoS/AsPoint.hs new file mode 100644 --- /dev/null +++ b/src/Algorithms/Geometry/SoS/AsPoint.hs @@ -0,0 +1,29 @@ +module Algorithms.Geometry.SoS.AsPoint where + +import Control.CanAquire +import Control.Lens +import Data.Ext +import Data.Geometry.Point.Class +import Data.Geometry.Point.Internal +import Data.Geometry.Properties +import Data.Geometry.Vector + +-------------------------------------------------------------------------------- +-- | a P is a 'read only' point in d dimensions +newtype P i d r = P i deriving (Eq, Show) + +-- | Indxec type that can disambiguate points +newtype SoSIndex i = SoSIndex i deriving (Show,Eq,Ord) + +instance HasIndex (P i d r) i where + indexOf (P i) = i + +instance Int `CanAquire` (Point d r) => (P Int d r) `CanAquire` (Point d r) where + aquire (P i) = aquire i + +type instance NumType (P i d r) = r +type instance Dimension (P i d r) = d + +asPointWithIndex :: (Arity d, i `CanAquire` Point d r) + => P i d r -> Point d r :+ SoSIndex i +asPointWithIndex (P i) = aquire i :+ (SoSIndex i) diff --git a/src/Algorithms/Geometry/SoS/Determinant.hs b/src/Algorithms/Geometry/SoS/Determinant.hs new file mode 100644 --- /dev/null +++ b/src/Algorithms/Geometry/SoS/Determinant.hs @@ -0,0 +1,13 @@ +module Algorithms.Geometry.SoS.Determinant where + +import Algorithms.Geometry.SoS.Sign +import Algorithms.Geometry.SoS.Symbolic +import Data.Geometry.Matrix + + +-- | pre: computes the sign of the determinant +signDet :: (HasDeterminant d, Ord i, Num r, Ord r) => Matrix d d (Symbolic i r) -> Sign +signDet m = case det m `compare` 0 of + LT -> Negative + GT -> Positive + EQ -> error "signDet: determinant is zero! this should not happen!" diff --git a/src/Algorithms/Geometry/SoS/Expr.hs b/src/Algorithms/Geometry/SoS/Expr.hs new file mode 100644 --- /dev/null +++ b/src/Algorithms/Geometry/SoS/Expr.hs @@ -0,0 +1,77 @@ +{-# LANGUAGE TemplateHaskell #-} +module Algorithms.Geometry.SoS.Expr where + +import Control.Lens +import qualified Data.List as List +import Data.List.NonEmpty (NonEmpty(..),nonEmpty) + +-------------------------------------------------------------------------------- + +data Expr v r = Constant r + | Negate (Expr v r) + | Sum [Expr v r] + | Prod [Expr v r] + | Var v + deriving (Show,Eq) +makePrisms ''Expr + + +foldExpr :: (r -> b) -> (b -> b) -> ([b] -> b) -> ([b] -> b) -> (v -> b) -> Expr v r -> b +foldExpr con' neg' sum' prod' var' = go + where + go = \case + Constant c -> con' c + Negate e -> neg' $ go e + Sum es -> sum' $ map go es + Prod es -> prod' $ map go es + Var v -> var' v + +-- | Test if the expression has any variables. +hasVariables :: Expr v r -> Bool +hasVariables = foldExpr (const False) + id + or + or + (const True) + +instance (Num r) => Num (Expr i r) where + fromInteger = Constant . fromInteger + negate = \case + Negate e -> e + e -> Negate e + + (Sum es) + (Sum es') = Sum $ es <> es' + (Sum es) + e = Sum (e:es) + e + (Sum es) = Sum (e:es) + e + e' = Sum [e,e'] + + (Prod es) * (Prod es') = Prod $ es <> es' + (Prod es) * e = Prod (e:es) + e * (Prod es) = Prod (e:es) + e * e' = Prod [e,e'] + + +simplify :: (Num r, Eq r) => Expr v r -> Expr v r +simplify = \case + Prod es -> case filter (isn't $ _Constant.only 1) es of + [] -> Constant 1 + es' -> Prod $ map simplify es' + Sum es -> case filter (isn't $ _Constant.only 0) es of + [] -> Constant 0 + es' -> Sum $ map simplify es' + Negate e -> Negate $ simplify e + e -> e + +prettyP :: (Show r, Show v) => Expr v r -> String +prettyP = \case + Constant c -> show c + Negate e -> "(-1)*(" <> prettyP e <> ")" + Prod es -> mconcat [ "(" + , List.intercalate ")*(" (prettyP <$> es) + , ")" + ] + Sum es -> mconcat [ "(" + , List.intercalate ") + (" (prettyP <$> es) + , ")" + ] + Var v -> show v diff --git a/src/Algorithms/Geometry/SoS/Internal.hs b/src/Algorithms/Geometry/SoS/Internal.hs new file mode 100644 --- /dev/null +++ b/src/Algorithms/Geometry/SoS/Internal.hs @@ -0,0 +1,28 @@ +module Algorithms.Geometry.SoS.Internal where + +import Algorithms.Geometry.SoS.AsPoint +import Algorithms.Geometry.SoS.Orientation +import Control.CanAquire +import Data.Geometry.Point.Internal + +-------------------------------------------------------------------------------- + +-- simulateSimplicity :: forall t d r b. (Traversable t, SoSD d) +-- => (forall p. ( AsPoint p, HasIndex p +-- , d ~ Dimension p, r ~ NumType p +-- ) => t p -> b) +-- -> t (Point d r) -> b +-- simulateSimplicity = simulateSimplicity' + + +-- | The actual implementation of SoS +simulateSimplicity' :: forall t d r b. (Traversable t, SoS d) + => (forall i. ( CanAquire i (Point d r) + , SoS d + ) => t (P i d r) -> b) + -> t (Point d r) -> b +simulateSimplicity' alg = runAcquire alg' + where + alg' :: forall i. CanAquire i (Point d r) => t i -> b + alg' = alg . fmap (P @i @d @r) + -- ideally the fmap would just be a coerce, but GHC does not want to do that. diff --git a/src/Algorithms/Geometry/SoS/Orientation.hs b/src/Algorithms/Geometry/SoS/Orientation.hs new file mode 100644 --- /dev/null +++ b/src/Algorithms/Geometry/SoS/Orientation.hs @@ -0,0 +1,83 @@ +module Algorithms.Geometry.SoS.Orientation( SoS + + , sideTest + , sideTest' + + , toSymbolic + ) where + +import Algorithms.Geometry.SoS.Determinant +import Algorithms.Geometry.SoS.Sign +import Algorithms.Geometry.SoS.Symbolic +import Control.Lens hiding (snoc,cons) +import Data.Ext +import Data.Geometry.Matrix +import Data.Geometry.Point +import Data.Geometry.Vector +import GHC.TypeNats + +-------------------------------------------------------------------------------- + + + +-- | A dimension d has support for SoS when we can: compute a +-- dterminant of a d+1 by d+1 dimensional matrix. +type SoS d = (Arity d, HasDeterminant (d+1)) + +-- | Given a query point q, and a vector of d points defining a +-- hyperplane test if q lies above or below the hyperplane. Each point +-- is assumed to have an unique index of type i that can be used to +-- disambiguate it in case of degeneracies. +-- +-- some 1D examples: +-- +-- >>> sideTest (Point1 0 :+ 0) (Vector1 $ Point1 2 :+ 1) +-- Negative +-- >>> sideTest (Point1 10 :+ 0) (Vector1 $ Point1 2 :+ 1) +-- Positive +-- >>> sideTest (Point1 2 :+ 0) (Vector1 $ Point1 2 :+ 1) +-- Positive +-- >>> sideTest (Point1 2 :+ 3) (Vector1 $ Point1 2 :+ 1) +-- Negative +-- +-- some 2D examples: +-- +-- >>> sideTest (Point2 1 2 :+ 0) $ Vector2 (Point2 0 0 :+ 1) (Point2 2 2 :+ 3) +-- Positive +-- >>> sideTest (Point2 1 (-2) :+ 0) $ Vector2 (Point2 0 0 :+ 1) (Point2 2 2 :+ 3) +-- Negative +-- >>> sideTest (Point2 1 1 :+ 0) $ Vector2 (Point2 0 0 :+ 1) (Point2 2 2 :+ 3) +-- Positive +-- >>> sideTest (Point2 1 1 :+ 10) $ Vector2 (Point2 0 0 :+ 1) (Point2 2 2 :+ 3) +-- Negative +-- >>> sideTest (Point2 1 1 :+ 10) $ Vector2 (Point2 0 0 :+ 3) (Point2 2 2 :+ 1) +-- Negative +sideTest :: (SoS d, Num r, Ord r, Ord i) + => Point d r :+ i -> Vector d (Point d r :+ i) -> Sign +sideTest q ps = sideTest'' . fmap toSymbolic $ cons q ps + +-- | Given an input point, transform its number type to include +-- symbolic $\varepsilon$ expressions so that we can use SoS. +toSymbolic :: (Ord i, Arity d) => Point d r :+ i -> Point d (Symbolic (i,Int) r) +toSymbolic (p :+ i) = p&vector' %~ imap (\j x -> symbolic x (i,j)) + +-- | Given a point q and a vector of d points defining a hyperplane, +-- test on which side of the hyperplane q lies. +-- +-- TODO: Specify what the sign means +sideTest' :: (Num r, Ord r, Ord i, HasDeterminant (d+1), Arity d, Arity (d+1)) + => Point d (Symbolic i r) -> Vector d (Point d (Symbolic i r)) -> Sign +sideTest' q ps = sideTest'' $ cons q ps + +-- | Given a vector of points, tests if the point encoded in the first +-- row is above/below the hyperplane defined by the remaining points +-- (rows). +sideTest'' :: (Num r, Ord r, Ord i, HasDeterminant (d+1), Arity d, Arity (d+1)) + => Vector (d+1) (Point d (Symbolic i r)) -> Sign +sideTest'' = signDet . Matrix . fmap mkLambdaRow + +-- | Given a point produces the vector/row corresponding to this point +-- in a homogeneous matrix represetnation. I.e. we add a 1 as an +-- additonal column at the end. +mkLambdaRow :: (Num r, Arity d, Arity (d+1)) => Point d r -> Vector (d+1) r +mkLambdaRow = flip snoc 1 . view vector' diff --git a/src/Algorithms/Geometry/SoS/Sign.hs b/src/Algorithms/Geometry/SoS/Sign.hs new file mode 100644 --- /dev/null +++ b/src/Algorithms/Geometry/SoS/Sign.hs @@ -0,0 +1,30 @@ +module Algorithms.Geometry.SoS.Sign where + +import qualified Data.List as List +import Data.Maybe + +-------------------------------------------------------------------------------- + +-- | The sign of an expression +data Sign = Negative | Positive deriving (Show,Eq,Ord,Enum,Bounded) + +flipSign :: Sign -> Sign +flipSign = \case + Negative -> Positive + Positive -> Negative + +-------------------------------------------------------------------------------- + +-- | Given the terms, in decreasing order of significance, computes the sign +-- +-- i.e. expects a list of terms, we base the sign on the sign of the first non-zero term. +-- +-- pre: the list contains at least one such a term. +signFromTerms :: (Num r, Eq r) => [r] -> Sign +signFromTerms = List.head . mapMaybe signum' + where + signum' x = case signum x of + -1 -> Just Negative + 0 -> Nothing + 1 -> Just Positive + _ -> error "signum': absurd" diff --git a/src/Algorithms/Geometry/SoS/Symbolic.hs b/src/Algorithms/Geometry/SoS/Symbolic.hs new file mode 100644 --- /dev/null +++ b/src/Algorithms/Geometry/SoS/Symbolic.hs @@ -0,0 +1,352 @@ +module Algorithms.Geometry.SoS.Symbolic( + EpsFold + , eps, mkEpsFold + , hasNoPertubation + , factors + , suitableBase + + , Term(..), term, constantFactor + + , Symbolic + , constant, symbolic, perturb + + , toTerms + , signOf + ) where + +import Algorithms.Geometry.SoS.Sign (Sign(..)) +import Control.Lens +import Data.Foldable (toList) +import qualified Data.List as List +import qualified Data.Map as Map +import qualified Data.Map.Merge.Strict as Map +import Data.Maybe (isNothing) +import Data.Word +import Test.QuickCheck (Arbitrary(..), listOf, suchThat) +import Test.QuickCheck.Instances () + +-------------------------------------------------------------------------------- +-- * EpsFolds + +{- +Let \(\mathcal{I}\) be a bag with indices, let \(c\) be an upper +bound on the number of times a single item may occur in +\(\mathcal{I}\), and let \(\varepsilon\) be a function mapping indices +to real numbers that satisfies: + +1. \(0 < \varepsilon(j) < 1\), for all \(1 \leq j\), +2. \(\prod_{0 \leq i \leq j} \varepsilon(i)^c > \varepsilon(k)\), for all \(1 \leq j < k\) + +Note that such a function exists: + +\begin{lemma} + \label{lem:condition_2} + Let \(\delta \in (0,1)\) and \(d \geq c+1\). The function + \(\varepsilon(i) = \delta^{d^i}\) satisfies condition 2. +\end{lemma} + +\begin{proof} + By transitivity it suffices to argue this for \(k=j+1\): + + \begin{align*} + &\qquad \prod_{0 \leq i \leq j} \varepsilon(i)^c > \varepsilon(j+1) \\ + \equiv &\qquad \prod_{0 \leq i \leq j} (\delta^{d^i})^c > \delta^{d^{j+1}}\\ + \equiv &\qquad \prod_{0 \leq i \leq j} \delta^{cd^i} > \delta^{d^{j+1}} \\ + \equiv &\qquad \delta^{\sum_{0 \leq i \leq j} cd^i} > \delta^{d^{j+1}} & + \text{using + } + \delta \in (0,1)\\ + \equiv &\qquad \sum_{0 \leq i \leq j} cd^i < d^{j+1} \\ + \equiv &\qquad c\sum_{0 \leq i \leq j} d^i < d^{j+1} \\ + \end{align*} + + We prove this by induction. + + For the base case \(j=0\): we have \(0 < 1\), which is trivially true. + + For the step case we have the induction hypothesis + \(c\sum_{0 \leq i \leq j} d^i < d^{j+1}\), and we have to prove that + \(c\sum_{0 \leq i \leq j+1} d^i < d^{j+2}\): + + \begin{align*} + c\sum_{0 \leq i \leq j+1} d^i + &= cd^{j+1} + c\sum_{0 \leq i \leq j} d^i \\ + &< cd^{j+1} + d^{j+1} & \text{using IH} \\ + &= (c+1)d^{j+1} & \text{using that } c+1 \leq d \\ + &\leq dd^{j+1} \\ + &=d^{j+2} + \end{align*} + This completes the proof. +\end{proof} + + + + + + +An EpsFold now represents the term + +\[ \prod_{i \in \mathcal{I}} \varepsilon(i) \] + +for some bag \(\mathcal{I}\). + + +Let \(\mathcal{J}\) be some sub-bag of \(\mathcal{I}\). Note that +condition 2 implies that: + +\(\prod_{i \in \mathcal{J}} \varepsilon(i) > \varepsilon(k)\), for all \(1 \leq j < k\) + +This means that when comparing two EpsFolds, say \(e_1\) and \(e_2\), +representing bags \(\mathcal{I}_1\) and \(\mathcal{I}_2\), +respectively. It suffices to compare the largest index +\(j \in \mathcal{I}_1\setminus\mathcal{I}_2\) with the largest index +\(k \in \mathcal{I}_2\setminus\mathcal{I}_1\). We have that + +\(e_1 > e_2\) if and only if \(j < k\). +-} +newtype EpsFold i = Pi (Bag i) deriving (Semigroup,Monoid) + +-- | Gets the factors +factors :: EpsFold i -> Bag i +factors (Pi is) = is + +-- | Creates the term \(\varepsilon(i)\) +eps :: i -> EpsFold i +eps = Pi . singleton + +mkEpsFold :: Ord i => [i] -> EpsFold i +mkEpsFold = Pi . foldMap singleton + + + +-- | computes a base 'd' that can be used as: +-- +-- \( \varepsilon(i) = \varepsilon^{d^i} \) +suitableBase :: EpsFold i -> Int +suitableBase = max 2 . (1+) . maxMultiplicity . factors + +instance Show i => Show (EpsFold i) where + showsPrec d (Pi b) = showParen (d > app_prec) $ + showString "Pi " . showsPrec d (toList b) + where + app_prec = 10 + + +instance Ord i => Eq (EpsFold i) where + e1 == e2 = (e1 `compare` e2) == EQ + +instance Ord i => Ord (EpsFold i) where + (Pi e1) `compare` (Pi e2) = k `compare` j -- note that k and j are flipped here + where + j = maximum' $ e1 `difference` e2 + k = maximum' $ e2 `difference` e1 + -- note: If the terms are all the same, the difference of the bags is empty + -- and thus both e1e2 and e2e1 are Nothing and thus equal. + + -- otherwise, let j be the largest term that is in e1 but not in e2. + -- If e2 does not have any terms at all (Nothing) it will be bigger than e1 + -- + -- if e2 does have a term, let k be the largest one, then the + -- biggest of those terms is the pair whose indices comes first. + +instance (Arbitrary i, Ord i) => Arbitrary (EpsFold i) where + arbitrary = (mkEpsFold . take 4) <$> listOf arbitrary + + +-- | Test if the epsfold has no pertubation at all (i.e. if it is \(\Pi_{\emptyset}\) +hasNoPertubation :: EpsFold i -> Bool +hasNoPertubation (Pi b) = null b + + +-------------------------------------------------------------------------------- +-- * Terms + +-- | A term 'Term c es' represents a term: +-- +-- \[ c \Pi_{i \in es} \varepsilon(i) +-- \] +-- +-- for a constant c and an arbitrarily small value \(\varepsilon\), +-- parameterized by i. +data Term i r = Term r (EpsFold i) deriving (Eq,Functor) + +-- | Lens to access the constant 'c' in the term. +constantFactor :: Lens' (Term i r) r +constantFactor = lens (\(Term c _) -> c) (\(Term _ es) c -> Term c es) + + +instance (Show i, Show r) => Show (Term i r) where + showsPrec d (Term c es) = showParen (d > up_prec) $ + showsPrec (up_prec + 1) c + . showString " * " + . showsPrec (up_prec + 1) es + where + up_prec = 5 + + +-- | Creates a singleton term +term :: r -> i -> Term i r +term r i = Term r $ eps i + +instance (Ord i, Ord r, Num r) => Ord (Term i r) where + (Term c e1) `compare` (Term d e2) = case (hasNoPertubation e1, hasNoPertubation e2) of + (True,True) -> c `compare` d + _ -> case (signum c, signum d) of + (-1,-1) -> e2 `compare` e1 + (0,0) -> e1 `compare` e2 + (1,1) -> e1 `compare` e2 + (-1,_) -> LT + (_,-1) -> GT + _ -> error "SoS: Term.ord absurd" + -- If both the eps folds are zero, and thus we just have constants + -- then we should compare the individual terms. + + -- if *one* of the two has an eps term, then we can choose eps to be + -- arbitrarily small, i.e. small enough so that that terms is + -- actually smaller than the other term. this is reflected since + -- findMax will then return a Noting, which is smaller than anything + -- else + + -- if both terms have epsilon terms, we first look at the sign. If + -- they have non-negative signs we compare the eps-folds as in the + -- paper. (Lemma 3.3). If both are negative, that reverses the + -- ordering. If the signs are different then we can base the + -- ordering on that. + +instance (Arbitrary r, Arbitrary (EpsFold i), Ord i) => Arbitrary (Term i r) where + arbitrary = Term <$> arbitrary <*> arbitrary + +-------------------------------------------------------------------------------- +-- * Symbolic + +-- | Represents a Sum of terms, i.e. a value that has the form: +-- +-- \[ +-- \sum c \Pi_i \varepsilon(i) +-- \] +-- +-- The terms are represented in order of decreasing significance. +-- +-- The main idea in this type is that, if symbolic values contains +-- \(\varepsilon(i)\) terms we can always order them. That is, two +-- Symbolic terms will be equal only if: +-- +-- - they contain *only* a constant term (that is equal) +-- - they contain the exact same \(\varepsilon\)-fold. +-- +newtype Symbolic i r = Sum (Map.Map (EpsFold i) r) deriving (Functor) + +-- | Produces a list of terms, in decreasing order of significance +toTerms :: Symbolic i r -> [Term i r] +toTerms (Sum m) = map (\(i,c) -> Term c i) . Map.toDescList $ m + +-- | Computing the Sign of an expression. (Nothing represents zero) +signOf :: (Num r, Eq r) => Symbolic i r -> Maybe Sign +signOf e = case List.dropWhile (== 0) . map (\(Term c _) -> signum c) $ toTerms e of + [] -> Nothing + (-1:_) -> Just Negative + _ -> Just Positive + +instance (Ord i, Eq r, Num r) => Eq (Symbolic i r) where + e1 == e2 = isNothing $ signOf (e1 - e2) + +instance (Ord i, Ord r, Num r) => Ord (Symbolic i r) where + e1 `compare` e2 = case signOf (e1 - e2) of + Nothing -> EQ + Just Negative -> LT + Just Positive -> GT + +instance (Ord i, Num r, Eq r) => Num (Symbolic i r) where + (Sum e1) + (Sum e2) = Sum $ Map.merge Map.preserveMissing -- insert things only in e1 + Map.preserveMissing -- insert things only in e2 + combine + e1 e2 + where + -- if things are in both e1 and e2, we add the constant terms. If they are non-zero + -- we use this value in the map. Otherwise we drop it. + combine = Map.zipWithMaybeMatched + (\_ c d -> let x = c + d in if x /= 0 then Just x else Nothing) + -- Symbolic $ Map.unionWith (+) ts ts' + + negate = fmap negate + + (Sum ts) * (Sum ts') = Sum $ Map.fromListWith (+) [ (es <> es',c*d) + | (es, c) <- Map.toList ts + , (es',d) <- Map.toList ts' + , c*d /= 0 + ] + + fromInteger x = constant (fromInteger x) + + signum s = case signOf s of + Nothing -> 0 + Just Negative -> (-1) + Just Positive -> 1 + + abs x | signum x == -1 = (-1)*x + | otherwise = x + + +instance (Show i, Show r) => Show (Symbolic i r) where + showsPrec d s = showParen (d > app_prec) $ + showString "Sum " . showsPrec d (toTerms s) + where + app_prec = 10 + +instance (Arbitrary r, Ord i, Arbitrary (EpsFold i)) => Arbitrary (Symbolic i r) where + arbitrary = Sum <$> arbitrary + +---------------------------------------- + +-- | Creates a constant symbolic value +constant :: Ord i => r -> Symbolic i r +constant c = Sum $ Map.singleton mempty c + +-- | Creates a symbolic vlaue with a single indexed term. If you just need a constant (i.e. non-indexed), use 'constant' +symbolic :: Ord i => r -> i -> Symbolic i r +symbolic r i = Sum $ Map.singleton (eps i) r + +-- | given the value c and the index i, creates the perturbed value +-- \(c + \varepsilon(i)\) +perturb :: (Num r, Ord i) => r -> i -> Symbolic i r +perturb c i = Sum $ Map.fromAscList [ (mempty,c) , (eps i,1) ] + + +-------------------------------------------------------------------------------- + +-- | The word specifiies how many *duplicates* there are. I.e. If the +-- Bag maps k to i, then k has multiplicity i+1. +newtype Bag a = Bag (Map.Map a Int) deriving (Show,Eq,Ord,Arbitrary) + +singleton :: k -> Bag k +singleton x = Bag $ Map.singleton x 0 + + +instance Foldable Bag where + -- ^ Takes multiplicity into account. + foldMap f (Bag m) = + Map.foldMapWithKey (\k d -> foldMap f (List.replicate (fromIntegral d+1) k)) m + null (Bag m) = Map.null m + +instance Ord k => Semigroup (Bag k) where + (Bag m) <> (Bag m') = Bag $ Map.unionWith (\d d' -> d + d' + 1) m m' + +instance Ord k => Monoid (Bag k) where + mempty = Bag $ Map.empty + +-- | Computes the difference of the two maps +difference :: Ord a => Bag a -> Bag a -> Bag a +difference (Bag m1) (Bag m2) = Bag $ Map.differenceWith updateCount m1 m2 + where + updateCount i j = let d = i - j -- note that we should actually compare (i+1) and (j+1) + in if d <= 0 then Nothing -- we have no copies left + else Just $ d - 1 + + +maximum' :: Bag b -> Maybe b +maximum' (Bag m) = fmap fst . Map.lookupMax $ m + + +-- | maximum multiplicity of an element in the bag +maxMultiplicity :: Bag a -> Int +maxMultiplicity (Bag m) = maximum . (0:) . map (1+) . Map.elems $ m diff --git a/src/Algorithms/Geometry/WellSeparatedPairDecomposition/Types.hs b/src/Algorithms/Geometry/WellSeparatedPairDecomposition/Types.hs --- a/src/Algorithms/Geometry/WellSeparatedPairDecomposition/Types.hs +++ b/src/Algorithms/Geometry/WellSeparatedPairDecomposition/Types.hs @@ -19,6 +19,7 @@ import Data.Geometry.Point import Data.Geometry.Vector import qualified Data.LSeq as LSeq +import Data.Measured.Class import qualified Data.Sequence as S import qualified Data.Traversable as Tr diff --git a/src/Algorithms/Geometry/WellSeparatedPairDecomposition/WSPD.hs b/src/Algorithms/Geometry/WellSeparatedPairDecomposition/WSPD.hs --- a/src/Algorithms/Geometry/WellSeparatedPairDecomposition/WSPD.hs +++ b/src/Algorithms/Geometry/WellSeparatedPairDecomposition/WSPD.hs @@ -50,7 +50,7 @@ => NonEmpty.NonEmpty (Point d r :+ p) -> SplitTree d p r () fairSplitTree pts = foldUp node' Leaf $ fairSplitTree' n pts' where - pts' = GV.imap sortOn . pure . g $ pts + pts' = imap sortOn . pure . g $ pts n = length $ pts'^.GV.element (C :: C 0) sortOn' i = NonEmpty.sortWith (^.core.unsafeCoord i) @@ -354,7 +354,7 @@ -- -- pre: points are sorted according to their dimension extends :: Arity d => GV.Vector d (PointSeq d p r) -> GV.Vector d (Range r) -extends = GV.imap (\i pts -> +extends = imap (\i pts -> ClosedRange ((LSeq.head pts)^.core.unsafeCoord (i + 1)) ((LSeq.last pts)^.core.unsafeCoord (i + 1))) diff --git a/src/Data/Geometry/Arrangement/Internal.hs b/src/Data/Geometry/Arrangement/Internal.hs --- a/src/Data/Geometry/Arrangement/Internal.hs +++ b/src/Data/Geometry/Arrangement/Internal.hs @@ -11,6 +11,7 @@ -------------------------------------------------------------------------------- module Data.Geometry.Arrangement.Internal where +import Algorithms.BinarySearch import Control.Lens import qualified Data.CircularSeq as CSeq import Data.Ext @@ -25,7 +26,6 @@ import qualified Data.List as List import Data.Maybe import Data.Ord (Down(..)) -import Data.Sequence.Util import qualified Data.Vector as V import Data.Vinyl.CoRec @@ -192,13 +192,10 @@ unBoundedParts rect ls = [tl] <> t <> [tr] <> reverse r <> [br] <> reverse b <> [bl] <> l where sideIntersections' = over (traverse._2) Just . sideIntersections ls - (t,r,b,l) = map4 sideIntersections' $ sides rect - (tl,tr,br,bl) = map4 ((,Nothing) . (^.core)) $ corners rect + Sides t r b l = fmap sideIntersections' $ sides rect + Corners tl tr br bl = fmap ((,Nothing) . (^.core)) $ corners rect -map4 :: (a -> b) -> (a,a,a,a) -> (b,b,b,b) -map4 f (a,b',c,d) = (f a, f b', f c, f d) - -- | Links the vertices of the outer boundary with those in the subdivision link :: Eq r => [(Point 2 r, a)] -> PlanarSubdivision s v (Maybe e) f r -> V.Vector (Point 2 r, VertexId' s, a) @@ -270,8 +267,8 @@ i <- binarySearchVec (pred' ss) (arr^.unboundedIntersections) pure $ arr^.unboundedIntersections.singular (ix i) where - (t,r,b,l) = sides'' $ arr^.boundedArea - sides'' = map4 (\(ClosedLineSegment a c) -> LineSegment (Closed a) (Open c)) . sides + Sides t r b l = sides'' $ arr^.boundedArea + sides'' = fmap (\(ClosedLineSegment a c) -> LineSegment (Closed a) (Open c)) . sides findSide q = fmap fst . List.find (onSegment q . snd) $ zip [1..] [t,r,b,l] 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 @@ -123,21 +123,27 @@ pattern Sphere c r = Boundary (Ball c r) {-# COMPLETE Sphere #-} - - +-- | +_BallSphere :: Iso (Disk p r) (Disk p s) (Circle p r) (Circle p s) +_BallSphere = _Boundary -------------------------------------------------------------------------------- -- * Disks and Circles, aka 2-dimensional Balls and Spheres type Disk p r = Ball 2 p r +-- | Given the center and the squared radius, constructs a disk pattern Disk :: Point 2 r :+ p -> r -> Disk p r pattern Disk c r = Ball c r {-# COMPLETE Disk #-} - type Circle p r = Sphere 2 p r +-- | Iso for converting between Disks and Circles, i.e. forgetting the boundary +_DiskCircle :: Iso (Disk p r) (Disk p s) (Circle p r) (Circle p s) +_DiskCircle = _BallSphere + +-- | Given the center and the squared radius, constructs a circle pattern Circle :: Point 2 r :+ p -> r -> Circle p r pattern Circle c r = Sphere c r {-# COMPLETE Circle #-} diff --git a/src/Data/Geometry/BezierSpline.hs b/src/Data/Geometry/BezierSpline.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/BezierSpline.hs @@ -0,0 +1,169 @@ +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE TemplateHaskell #-} +module Data.Geometry.BezierSpline( + BezierSpline (BezierSpline) + , controlPoints + , fromPointSeq + + , evaluate + , split + , subBezier + , tangent + , approximate + , parameterOf + , snap + + , pattern Bezier2, pattern Bezier3 + ) where + +import Control.Lens hiding (Empty) +import qualified Data.Foldable as F +import Data.Geometry.Point +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 GHC.TypeNats +import qualified Test.QuickCheck as QC + +-------------------------------------------------------------------------------- + +-- | Datatype representing a Bezier curve of degree \(n\) in \(d\)-dimensional space. +newtype BezierSpline n d r = BezierSpline { _controlPoints :: LSeq (1+n) (Point d r) } +makeLenses ''BezierSpline + +-- | Quadratic Bezier Spline +pattern Bezier2 :: Point d r -> Point d r -> Point d r -> BezierSpline 2 d r +pattern Bezier2 p q r <- ((F.toList . LSeq.take 3 . _controlPoints) -> [p,q,r]) + where + Bezier2 p q r = fromPointSeq . Seq.fromList $ [p,q,r] +{-# COMPLETE Bezier2 #-} + +-- | Cubic Bezier Spline +pattern Bezier3 :: Point d r -> Point d r -> Point d r -> Point d r -> BezierSpline 3 d r +pattern Bezier3 p q r s <- ((F.toList . LSeq.take 4 . _controlPoints) -> [p,q,r,s]) + where + Bezier3 p q r s = fromPointSeq . Seq.fromList $ [p,q,r,s] +{-# COMPLETE Bezier3 #-} + +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 d, Show r) => Show (BezierSpline n d r) where + show (BezierSpline ps) = + mconcat [ "BezierSpline", show $ length ps - 1, " ", show (F.toList ps) ] + +instance Arity d => Functor (BezierSpline n d) where + fmap = fmapDefault + +instance Arity d => Foldable (BezierSpline n d) where + foldMap = foldMapDefault + +instance Arity d => Traversable (BezierSpline n d) where + traverse f (BezierSpline ps) = BezierSpline <$> traverse (traverse f) ps + +instance (Fractional r, Arity d, Arity (d + 1), Arity n) + => IsTransformable (BezierSpline n d r) where + transformBy = transformPointFunctor + +instance PointFunctor (BezierSpline n d) where + pmap f = over controlPoints (fmap f) + +-- | 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 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 :: (Arity d, Num r, 1 <= n) => BezierSpline n d r -> Vector d r +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]. +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 + +-- | 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 + ) + +collect :: (Arity d, Ord r, Num r) => r -> LSeq n (Point d r) -> Seq (Point d r) +collect t = go . LSeq.toSeq + where + go = \case + ps@(_ :<| Empty) -> ps + ps@(p :<| tl) -> let (ini :|> q) = ps in (p :<| go (Seq.zipWith blend ini tl)) :|> q + _ -> error "collect: absurd" + + blend p q = p .+^ t *^ (q .-. p) + +-- {- + +-- -- | 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 + +-- -} + +-- | 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 + +-- | 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 + +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 + +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 diff --git a/src/Data/Geometry/Boundary.hs b/src/Data/Geometry/Boundary.hs --- a/src/Data/Geometry/Boundary.hs +++ b/src/Data/Geometry/Boundary.hs @@ -1,7 +1,8 @@ module Data.Geometry.Boundary where -import Data.Geometry.Properties -import Data.Geometry.Transformation +import Control.Lens (iso,Iso) +import Data.Geometry.Properties +import Data.Geometry.Transformation -------------------------------------------------------------------------------- @@ -12,6 +13,10 @@ type instance NumType (Boundary g) = NumType g type instance Dimension (Boundary g) = Dimension g + +-- | Iso for converting between things with a boundary and without its boundary +_Boundary :: Iso g h (Boundary g) (Boundary h) +_Boundary = iso Boundary (\(Boundary b) -> b) -- | Result of a query that asks if something is Inside a g, *on* the boundary diff --git a/src/Data/Geometry/Box.hs b/src/Data/Geometry/Box.hs --- a/src/Data/Geometry/Box.hs +++ b/src/Data/Geometry/Box.hs @@ -1,5 +1,3 @@ -{-# LANGUAGE TemplateHaskell #-} -{-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE DeriveAnyClass #-} {-# OPTIONS_GHC -fno-warn-orphans #-} @@ -14,56 +12,16 @@ -- -------------------------------------------------------------------------------- module Data.Geometry.Box( module Data.Geometry.Box.Internal - , topSide, leftSide, bottomSide, rightSide - , sides, sides' + , module Data.Geometry.Box.Corners + , module Data.Geometry.Box.Sides ) where import Control.DeepSeq +import Data.Geometry.Box.Corners import Data.Geometry.Box.Internal -import Data.Geometry.LineSegment +import Data.Geometry.Box.Sides import Data.Geometry.Vector -------------------------------------------------------------------------------- deriving instance (NFData p, NFData r, Arity d) => NFData (Box d p r) - - -topSide :: Num r => Rectangle p r -> LineSegment 2 p r -topSide = (\(l,r,_,_) -> ClosedLineSegment l r) . corners - --- | Oriented from *left to right* -bottomSide :: Num r => Rectangle p r -> LineSegment 2 p r -bottomSide = (\(_,_,r,l) -> ClosedLineSegment l r) . corners - --- -leftSide :: Num r => Rectangle p r -> LineSegment 2 p r -leftSide = (\(t,_,_,b) -> ClosedLineSegment b t) . corners - --- | The right side, oriented from *bottom* to top -rightSide :: Num r => Rectangle p r -> LineSegment 2 p r -rightSide = (\(_,t,b,_) -> ClosedLineSegment b t) . corners - - --- | The sides of the rectangle, in order (Top, Right, Bottom, Left). The sides --- themselves are also oriented in clockwise order. If, you want them in the --- same order as the functions `topSide`, `bottomSide`, `leftSide`, and --- `rightSide`, use `sides'` instead. -sides :: Num r => Rectangle p r -> ( LineSegment 2 p r - , LineSegment 2 p r - , LineSegment 2 p r - , LineSegment 2 p r - ) -sides = (\(t,r,b,l) -> (t,flipSegment r,flipSegment b,l)) . sides' - - --- | The sides of the rectangle. The order of the segments is (Top, Right, --- Bottom, Left). Note that the segments themselves, are oriented as described --- by the functions topSide, bottomSide, leftSide, rightSide (basically: from --- left to right, and from bottom to top). If you want the segments oriented --- along the boundary of the rectangle, use the `sides` function instead. -sides' :: Num r => Rectangle p r -> ( LineSegment 2 p r - , LineSegment 2 p r - , LineSegment 2 p r - , LineSegment 2 p r - ) -sides' r = (topSide r, rightSide r, bottomSide r, leftSide r) diff --git a/src/Data/Geometry/Box/Corners.hs b/src/Data/Geometry/Box/Corners.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Box/Corners.hs @@ -0,0 +1,70 @@ +{-# LANGUAGE TemplateHaskell #-} +module Data.Geometry.Box.Corners( Corners(Corners), northWest, northEast, southEast, southWest + , corners, cornersInDirection + ) where + +import Control.Lens (makeLenses,Ixed(..),Index, IxValue,(%~),(&),(^?!)) +import Data.Ext +import Data.Functor.Apply +import Data.Geometry.Box.Internal +import Data.Geometry.Directions +import Data.Geometry.Point +import Data.Semigroup.Foldable.Class +import Data.Semigroup.Traversable.Class +import Data.Util +import GHC.Generics (Generic) + +-------------------------------------------------------------------------------- + +-- | A Quadrant data type +data Corners a = Corners { _northWest :: !a + , _northEast :: !a + , _southEast :: !a + , _southWest :: !a + } deriving (Show,Eq,Ord,Generic,Functor,Foldable,Traversable) +makeLenses ''Corners + + +type instance Index (Corners a) = InterCardinalDirection +type instance IxValue (Corners a) = a + +instance Ixed (Corners a) where + ix = \case + NorthWest -> northWest + NorthEast -> northEast + SouthEast -> southEast + SouthWest -> southWest + +instance Foldable1 Corners +instance Traversable1 Corners where + traverse1 f (Corners a b c d) = Corners <$> f a <.> f b <.> f c <.> f d + +instance Applicative Corners where + pure x = Corners x x x x + (Corners f g h i) <*> (Corners a b c d) = Corners (f a) (g b) (h c) (i d) + +instance Semigroup a => Semigroup (Corners a) where + s <> s' = (<>) <$> s <*> s' +instance Monoid a => Monoid (Corners a) where + mempty = pure mempty + + +-------------------------------------------------------------------------------- + +-- | Get the corners of a rectangle, the order is: +-- (TopLeft, TopRight, BottomRight, BottomLeft). +-- The extra values in the Top points are taken from the Top point, +-- the extra values in the Bottom points are taken from the Bottom point +corners :: Num r => Rectangle p r -> Corners (Point 2 r :+ p) +corners r = let w = width r + p = (_maxP r)&core %~ _cwMax + q = (_minP r)&core %~ _cwMin + in Corners (p&core.xCoord %~ (subtract w)) p + (q&core.xCoord %~ (+ w)) q + + +-------------------------------------------------------------------------------- + +-- | Gets the corners in a particular direction +cornersInDirection :: CardinalDirection -> Corners p -> Two p +cornersInDirection d c = (\icd -> c^?!ix icd) <$> interCardinalsOf d 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 @@ -15,8 +15,11 @@ import Control.DeepSeq import Control.Lens +import Data.Bifoldable import Data.Bifunctor +import Data.Bitraversable import Data.Ext +import qualified Data.Foldable as F import Data.Geometry.Point import Data.Geometry.Properties import Data.Geometry.Transformation @@ -25,12 +28,11 @@ import qualified Data.List.NonEmpty as NE import qualified Data.Range as R import qualified Data.Semigroup.Foldable as F -import qualified Data.Foldable as F import qualified Data.Vector.Fixed as FV import Data.Vinyl.CoRec (asA) import GHC.Generics (Generic) import GHC.TypeLits -import Test.QuickCheck(Arbitrary(..)) +import Test.QuickCheck (Arbitrary(..)) -------------------------------------------------------------------------------- @@ -60,6 +62,10 @@ } deriving Generic makeLenses ''Box + + + + -- | Given the point with the lowest coordinates and the point with highest -- coordinates, create a box. box :: Point d r :+ p -> Point d r :+ p -> Box d p r @@ -108,12 +114,14 @@ r `intersect'` s = asA @(R.Range r) $ r `intersect` s instance Arity d => Bifunctor (Box d) where - bimap :: forall p q r s. (p -> q) -> (r -> s) -> Box d p r -> Box d q s - bimap f g (Box mi ma) = Box (bimap g' f mi) (bimap g' f ma) + bimap = bimapDefault +instance Arity d => Bifoldable (Box d) where + bifoldMap = bifoldMapDefault +instance Arity d => Bitraversable (Box d) where + bitraverse f g (Box mi ma) = Box <$> bitraverse (tr g) f mi <*> bitraverse (tr g) f ma where - g' :: Functor g => g (Point d r) -> g (Point d s) - g' = fmap (fmap g) - + tr :: (Traversable t, Applicative f) => (r -> f s) -> t (Point d r) -> f (t (Point d s)) + tr g' = traverse $ traverse g' -- -- In principle this should also just work for Boxes in higher dimensions. It is just -- -- that we need a better way to compute their corners @@ -240,23 +248,7 @@ height = widthIn (C :: C 2) --- | Get the corners of a rectangle, the order is: --- (TopLeft, TopRight, BottomRight, BottomLeft). --- The extra values in the Top points are taken from the Top point, --- the extra values in the Bottom points are taken from the Bottom point -corners :: Num r => Rectangle p r -> ( Point 2 r :+ p - , Point 2 r :+ p - , Point 2 r :+ p - , Point 2 r :+ p - ) -corners r = let w = width r - p = (_maxP r)&core %~ _cwMax - q = (_minP r)&core %~ _cwMin - in ( p&core.xCoord %~ (subtract w) - , p - , q&core.xCoord %~ (+ w) - , q - ) +-------------------------------------------------------------------------------- -------------------------------------------------------------------------------- -- * Constructing bounding boxes @@ -282,3 +274,6 @@ instance IsBoxable (Box d p r) where boundingBox (Box m m') = Box (m&extra .~ ()) (m'&extra .~ ()) + +instance IsBoxable c => IsBoxable (c :+ e) where + boundingBox = boundingBox . view core diff --git a/src/Data/Geometry/Box/Sides.hs b/src/Data/Geometry/Box/Sides.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Box/Sides.hs @@ -0,0 +1,91 @@ +{-# LANGUAGE TemplateHaskell #-} +module Data.Geometry.Box.Sides( Sides(Sides), north, east, south, west + , topSide, bottomSide, leftSide, rightSide + , sides, sides' + + , sideDirections + ) where + +import Data.Geometry.Directions +import Data.Geometry.Box.Internal +import Data.Geometry.Box.Corners +import Data.Geometry.LineSegment +import Data.Functor.Apply +import Data.Semigroup.Foldable.Class +import Data.Semigroup.Traversable.Class +import GHC.Generics (Generic) +import Control.Lens(makeLenses, Ixed(..), Index, IxValue) + +-------------------------------------------------------------------------------- + +-- | The four sides of a rectangle +data Sides a = Sides { _north :: !a + , _east :: !a + , _south :: !a + , _west :: !a + } deriving (Show,Read,Eq,Generic,Ord,Foldable,Functor,Traversable) +makeLenses ''Sides + +instance Applicative Sides where + pure x = Sides x x x x + (Sides f g h i) <*> (Sides a b c d) = Sides (f a) (g b) (h c) (i d) + +instance Foldable1 Sides +instance Traversable1 Sides where + traverse1 f (Sides a b c d) = Sides <$> f a <.> f b <.> f c <.> f d + +instance Semigroup a => Semigroup (Sides a) where + s <> s' = (<>) <$> s <*> s' +instance Monoid a => Monoid (Sides a) where + mempty = pure mempty + + +type instance Index (Sides a) = CardinalDirection +type instance IxValue (Sides a) = a + +instance Ixed (Sides a) where + ix = \case + North -> north + East -> east + South -> south + West -> west + +-- | Constructs a Sides value that indicates the appropriate +-- direction. +sideDirections :: Sides CardinalDirection +sideDirections = Sides North East South West + +-------------------------------------------------------------------------------- + +topSide :: Num r => Rectangle p r -> LineSegment 2 p r +topSide = (\(Corners l r _ _) -> ClosedLineSegment l r) . corners + +-- | Oriented from *left to right* +bottomSide :: Num r => Rectangle p r -> LineSegment 2 p r +bottomSide = (\(Corners _ _ r l) -> ClosedLineSegment l r) . corners + +-- +leftSide :: Num r => Rectangle p r -> LineSegment 2 p r +leftSide = (\(Corners t _ _ b) -> ClosedLineSegment b t) . corners + +-- | The right side, oriented from *bottom* to top +rightSide :: Num r => Rectangle p r -> LineSegment 2 p r +rightSide = (\(Corners _ t b _) -> ClosedLineSegment b t) . corners + + +-- | The sides of the rectangle, in order (Top, Right, Bottom, Left). The sides +-- themselves are also oriented in clockwise order. If, you want them in the +-- same order as the functions `topSide`, `bottomSide`, `leftSide`, and +-- `rightSide`, use `sides'` instead. +sides :: Num r => Rectangle p r -> Sides (LineSegment 2 p r) +sides r = let Corners nw ne se sw = corners r + in Sides (ClosedLineSegment nw ne) (ClosedLineSegment ne se) + (ClosedLineSegment se sw) (ClosedLineSegment sw nw) + +-- | The sides of the rectangle. The order of the segments is (Top, Right, +-- Bottom, Left). Note that the segments themselves, are oriented as described +-- by the functions topSide, bottomSide, leftSide, rightSide (basically: from +-- left to right, and from bottom to top). If you want the segments oriented +-- along the boundary of the rectangle, use the `sides` function instead. +sides' :: Num r => Rectangle p r -> Sides (LineSegment 2 p r) +sides' r = Sides (topSide r) (rightSide r) (bottomSide r) (leftSide r) diff --git a/src/Data/Geometry/Directions.hs b/src/Data/Geometry/Directions.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Directions.hs @@ -0,0 +1,49 @@ +{-# LANGUAGE TemplateHaskell #-} +module Data.Geometry.Directions( CardinalDirection(..) + , _North, _East, _South, _West + , oppositeDirection + + , InterCardinalDirection(..) + , _NorthWest, _NorthEast, _SouthEast, _SouthWest + + , interCardinalsOf + ) where + +import Control.Lens (makePrisms) +import Data.Util +import GHC.Generics (Generic) + +-------------------------------------------------------------------------------- + +data CardinalDirection = North | East | South | West deriving (Show,Read,Eq,Ord,Enum,Bounded) +makePrisms ''CardinalDirection + +-------------------------------------------------------------------------------- +-- * Functions on Cardinal Directions + +-- | Computes the direction opposite to the given one. +oppositeDirection :: CardinalDirection -> CardinalDirection +oppositeDirection = \case + North -> South + East -> West + South -> North + West -> East + +-------------------------------------------------------------------------------- + +-- | Intercardinal directions +data InterCardinalDirection = NorthWest | NorthEast | SouthEast | SouthWest + deriving (Show,Read,Eq,Ord,Enum,Generic) +makePrisms ''InterCardinalDirection + +-------------------------------------------------------------------------------- +-- * Functions on InterCardinal Directions + +-- | Get the two intercardinal directions, in increasing order, +-- corresponding to the cardinal direction. +interCardinalsOf :: CardinalDirection -> Two InterCardinalDirection +interCardinalsOf = \case + North -> Two NorthWest NorthEast + East -> Two NorthEast SouthEast + South -> Two SouthEast SouthWest + West -> Two SouthWest NorthWest diff --git a/src/Data/Geometry/Ellipse.hs b/src/Data/Geometry/Ellipse.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Ellipse.hs @@ -0,0 +1,56 @@ +{-# LANGUAGE TemplateHaskell #-} +module Data.Geometry.Ellipse( + Ellipse(Ellipse) + , affineTransformation + , ellipseMatrix + , unitEllipse + , circleToEllipse, ellipseToCircle, _EllipseCircle + ) where + +import Control.Lens +import Data.Ext +import Data.Geometry.Ball +import Data.Geometry.Matrix +import Data.Geometry.Transformation +import Data.Geometry.Point +import Data.Geometry.Properties +import Data.Geometry.Vector + +-------------------------------------------------------------------------------- + +-- | A typre representing planar ellipses +newtype Ellipse r = Ellipse { _affineTransformation :: Transformation 2 r } + deriving (Show,Eq,Functor,Foldable,Traversable) +makeLenses ''Ellipse + +type instance Dimension (Ellipse r) = 2 +type instance NumType (Ellipse r) = r + +instance Num r => IsTransformable (Ellipse r) where + transformBy t (Ellipse t') = Ellipse $ t |.| t' + + +ellipseMatrix :: Iso (Ellipse r) (Ellipse s) (Matrix 3 3 r) (Matrix 3 3 s) +ellipseMatrix = affineTransformation.transformationMatrix + +-- | Ellipse representing the unit circle +unitEllipse :: Num r => Ellipse r +unitEllipse = Ellipse $ Transformation identityMatrix + +-------------------------------------------------------------------------------- +-- | Converting between ellipses and circles + +_EllipseCircle :: (Floating r, Eq r) => Prism' (Ellipse r) (Circle () r) +_EllipseCircle = prism' circleToEllipse ellipseToCircle + +ellipseToCircle :: (Num r, Eq r) => Ellipse r -> Maybe (Circle () r) +ellipseToCircle e = case e^.ellipseMatrix of + Matrix (Vector3 (Vector3 sx 0 x) + (Vector3 0 sy y) + (Vector3 0 0 1) + ) + | sx == sy -> Just $ Circle (ext $ Point2 x y) (sx*sx) + _ -> Nothing + +circleToEllipse :: Floating r => Circle p r -> Ellipse r +circleToEllipse (Circle (Point v :+ _) rr) = Ellipse $ translation v |.| uniformScaling (sqrt rr) 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 @@ -31,7 +31,6 @@ makeLenses ''HalfLine deriving instance (Show r, Arity d) => Show (HalfLine d r) -deriving instance (Eq r, Arity d) => Eq (HalfLine d r) deriving instance (NFData r, Arity d) => NFData (HalfLine d r) deriving instance Arity d => Functor (HalfLine d) @@ -40,6 +39,11 @@ type instance Dimension (HalfLine d r) = d type instance NumType (HalfLine d r) = r + + +instance (Eq r, Fractional r, Arity d) => Eq (HalfLine d r) where + (HalfLine p u) == (HalfLine q v) = let lam = scalarMultiple u v + in p == q && (signum <$> lam) == Just 1 instance HasStart (HalfLine d r) where type StartCore (HalfLine d r) = Point d r 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 @@ -38,13 +38,14 @@ -------------------------------------------------------------------------------- --- | A Halfspace in \(d\) dimensions. +-- | A Halfspace in \(d\) dimensions. Note that the intended side of +-- the halfspace is already indicated by the normal vector of the +-- bounding plane. newtype HalfSpace d r = HalfSpace { _boundingPlane :: HyperPlane d r } deriving Generic makeLenses ''HalfSpace deriving instance (Arity d, Show r) => Show (HalfSpace d r) -deriving instance (Arity d, Eq r) => Eq (HalfSpace d r) -- deriving instance (NFData r, Arity d) => NFData (HalfSpace d r) deriving instance Arity d => Functor (HalfSpace d) deriving instance Arity d => Foldable (HalfSpace d) @@ -54,6 +55,12 @@ type instance Dimension (HalfSpace d r) = d deriving instance (Arity d, Arity (d + 1), Fractional r) => IsTransformable (HalfSpace d r) + +instance (Arity d, Eq r, Fractional r) => Eq (HalfSpace d r) where + (HalfSpace h) == (HalfSpace h') = let u = h^.normalVec + v = h'^.normalVec + d = quadrance (u ^+^ v) - (quadrance u) + in h == h' && signum d == 1 -------------------------------------------------------------------------------- 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 @@ -25,12 +25,15 @@ type instance NumType (HyperPlane d r) = r deriving instance (Arity d, Show r) => Show (HyperPlane d r) -deriving instance (Arity d, Eq r) => Eq (HyperPlane d r) deriving instance (NFData r, Arity d) => NFData (HyperPlane d r) deriving instance Arity d => Functor (HyperPlane d) deriving instance Arity d => Foldable (HyperPlane d) deriving instance Arity d => Traversable (HyperPlane d) +instance (Arity d, Eq r, Fractional r) => Eq (HyperPlane d r) where + (HyperPlane p u) == h@(HyperPlane _ v) = p `intersects` h && u `isScalarMultipleOf` v + + instance (Arity d, Arity (d + 1), Fractional r) => IsTransformable (HyperPlane d r) where transformBy t (HyperPlane p v) = HyperPlane (transformBy t p) (transformBy t v) @@ -68,12 +71,29 @@ pattern Plane :: Point 3 r -> Vector 3 r -> Plane r pattern Plane p n = HyperPlane p n +{-# COMPLETE Plane #-} +-- | Produces a plane. If r lies counter clockwise of q w.r.t. p then +-- the normal vector of the resulting plane is pointing "upwards". +-- +-- >>> from3Points origin (Point3 1 0 0) (Point3 0 1 0) +-- HyperPlane {_inPlane = Point3 [0,0,0], _normalVec = Vector3 [0,0,1]} from3Points :: Num r => Point 3 r -> Point 3 r -> Point 3 r -> HyperPlane 3 r from3Points p q r = let u = q .-. p v = r .-. p in HyperPlane p (u `cross` v) +instance OnSideUpDownTest (Plane r) where + -- >>> (Point3 5 5 5) `onSideUpDown` from3Points origin (Point3 1 0 0) (Point3 0 1 0) + -- Above + -- >>> (Point3 5 5 (-5)) `onSideUpDown` from3Points origin (Point3 1 0 0) (Point3 0 1 0) + -- Below + -- >>> (Point3 5 5 0) `onSideUpDown` from3Points origin (Point3 1 0 0) (Point3 0 1 0) + -- On + q `onSideUpDown` (Plane p n) = let v = q .-. p in case (n `dot` v) `compare` 0 of + LT -> Below + EQ -> On + GT -> Above type instance IntersectionOf (Line 3 r) (Plane r) = [NoIntersection, Point 3 r, Line 3 r] @@ -107,3 +127,23 @@ instance HasSupportingPlane (HyperPlane d r) where supportingPlane = id + + +-- | Given +-- * a plane, +-- * a unit vector in the plane that will represent the y-axis (i.e. the "view up" vector), and +-- * a point in the plane, +-- +-- computes the plane coordinates of the given point, using the +-- inPlane point as the origin, the normal vector of the plane as the +-- unit vector in the "z-direction" and the view up vector as the +-- y-axis. +-- +-- >>> planeCoordinatesWith (Plane origin (Vector3 0 0 1)) (Vector3 0 1 0) (Point3 10 10 0) +-- Point2 [10.0,10.0] +planeCoordinatesWith :: Fractional r => Plane r -> Vector 3 r -> Point 3 r -> Point 2 r +planeCoordinatesWith h vup = projectPoint . transformBy (planeCoordinatesTransform h vup) + +planeCoordinatesTransform :: Num r => Plane r -> Vector 3 r -> Transformation 3 r +planeCoordinatesTransform (HyperPlane o n) v = rotateTo (Vector3 (v `cross` n) v n) + |.| translation ((-1) *^ toVec o) 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 @@ -1,7 +1,10 @@ {-# LANGUAGE TemplateHaskell #-} module Data.Geometry.Interval( -- * 1 dimensional Intervals - Interval(..) + Interval + , fromRange, toRange + , _Range + , pattern OpenInterval , pattern ClosedInterval , pattern Interval @@ -12,12 +15,13 @@ , inInterval , shiftLeft' + , asProperInterval, flipInterval + , module Data.Range - ) - where + ) where import Control.DeepSeq -import Control.Lens (makeLenses, (^.),(%~),(&), Lens') +import Control.Lens (lens, (^.),(%~),(&), Lens') import Data.Bifunctor import Data.Bitraversable import Data.Ext @@ -34,20 +38,31 @@ -------------------------------------------------------------------------------- -- | An Interval is essentially a 'Data.Range' but with possible payload -newtype Interval a r = GInterval { _unInterval :: Range (r :+ a) } +-- +-- We can think of an interval being defined as: +-- +-- >>> data Interval a r = Interval (EndPoint (r :+ a)) (EndPoint (r :+ a)) +newtype Interval a r = GInterval { toRange :: Range (r :+ a) } deriving (Eq,Generic,Arbitrary) -makeLenses ''Interval +_Range :: Lens' (Interval a r) (Range (r :+ a)) +_Range = lens toRange (const GInterval) +{-# INLINE _Range #-} + +-- | Constrct an interval from a Range +fromRange :: Range (r :+ a) -> Interval a r +fromRange = GInterval + deriving instance (NFData a, NFData r) => NFData (Interval a r) instance (Show a, Show r) => Show (Interval a r) where show ~(Interval l u) = concat [ "Interval (", show l, ") (", show u,")"] instance Functor (Interval a) where - fmap = T.fmapDefault + fmap f (GInterval r) = GInterval $ fmap (first f) r instance F.Foldable (Interval a) where - foldMap = T.foldMapDefault + foldMap f (GInterval r) = foldMap (f . (^.core)) r instance T.Traversable (Interval a) where traverse f (GInterval r) = GInterval <$> T.traverse f' r @@ -63,7 +78,7 @@ -- inInterval and inRange is that the extra value is *not* used in the -- comparison with inInterval, whereas it is in inRange. inInterval :: Ord r => r -> Interval a r -> Bool -x `inInterval` r = x `inRange` (fmap (^.core) $ r^.unInterval ) +x `inInterval` r = x `inRange` (fmap (^.core) $ r^._Range ) pattern OpenInterval :: (r :+ a) -> (r :+ a) -> Interval a r @@ -87,7 +102,8 @@ instance HasStart (Interval a r) where type StartCore (Interval a r) = r type StartExtra (Interval a r) = a - start = unInterval.lower.unEndPoint + start = _Range.lower.unEndPoint + {-# INLINE start #-} class HasEnd t where type EndCore t @@ -97,7 +113,8 @@ instance HasEnd (Interval a r) where type EndCore (Interval a r) = r type EndExtra (Interval a r) = a - end = unInterval.upper.unEndPoint + end = _Range.upper.unEndPoint + {-# INLINE end #-} type instance Dimension (Interval a r) = 1 type instance NumType (Interval a r) = r @@ -121,6 +138,17 @@ g (Arg _ x) = x +-- | Shifts the interval to the left by delta +shiftLeft' :: Num r => r -> Interval a r -> Interval a r +shiftLeft' delta = fmap (subtract delta) -shiftLeft' :: Num r => r -> Interval a r -> Interval a r -shiftLeft' x = fmap (subtract x) + +-- | Makes sure the start and endpoint are oriented such that the +-- starting value is smaller than the ending value. +asProperInterval :: Ord r => Interval p r -> Interval p r +asProperInterval i | (i^.start.core) > (i^.end.core) = flipInterval i + | otherwise = i + +-- | Flips the start and endpoint of the interval. +flipInterval :: Interval a r -> Interval a r +flipInterval = _Range %~ \(Range s t) -> Range t s diff --git a/src/Data/Geometry/IntervalTree.hs b/src/Data/Geometry/IntervalTree.hs --- a/src/Data/Geometry/IntervalTree.hs +++ b/src/Data/Geometry/IntervalTree.hs @@ -55,7 +55,7 @@ => [i] -> IntervalTree i r fromIntervals is = foldr insert (createTree pts) is where - endPoints (toRange -> Range' a b) = [a,b] + endPoints (asRange -> Range' a b) = [a,b] pts = List.sort . concatMap endPoints $ is -- | Lists the intervals. We don't guarantee anything about the order @@ -100,7 +100,7 @@ => i -> IntervalTree i r -> IntervalTree i r insert i (IntervalTree t) = IntervalTree $ insert' t where - ri@(Range a b) = toRange i + ri@(Range a b) = asRange i insert' Nil = Nil insert' (Internal l nd@(_splitPoint -> m) r) @@ -119,7 +119,7 @@ => i -> IntervalTree i r -> IntervalTree i r delete i (IntervalTree t) = IntervalTree $ delete' t where - ri@(Range a b) = toRange i + ri@(Range a b) = asRange i delete' Nil = Nil delete' (Internal l nd@(_splitPoint -> m) r) @@ -137,15 +137,13 @@ -- | Anything that looks like an interval class IntervalLike i where - toRange :: i -> Range (NumType i) + asRange :: i -> Range (NumType i) instance IntervalLike (Range r) where - toRange = id + asRange = id instance IntervalLike (Interval p r) where - toRange = fmap (^.core) . _unInterval - - + asRange = fmap (^.core) . toRange -------------------------------------------------------------------------------- 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 @@ -80,7 +80,7 @@ _ -> coRec NoIntersection _ -> error "intersect; line x boundary rect; absurd" where - (t,r,b,l) = sides' rect + Sides t r b l = sides' rect ints = map (\s -> sl `intersect` toSL s) [t,r,b,l] nub' = map L.head . L.group . L.sort 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 @@ -53,6 +53,8 @@ instance (Arity d, Eq r, Fractional r) => Eq (Line d r) where l@(Line p _) == m = l `isParallelTo` m && p `onLine` m + + instance (Arbitrary r, Arity d, Num r, Eq r) => Arbitrary (Line d r) where arbitrary = do p <- arbitrary q <- suchThat arbitrary (/= p) @@ -215,29 +217,34 @@ :& (H $ \_ -> Nothing) -- l is a vertical line (through x=0) :& RNil + -- | Result of a side test data SideTestUpDown = Below | On | Above deriving (Show,Read,Eq,Ord) --- | Given a point q and a line l, compute to which side of l q lies. For --- vertical lines the left side of the line is interpeted as below. --- --- >>> Point2 10 10 `onSideUpDown` (lineThrough origin $ Point2 10 5) --- Above --- >>> Point2 10 10 `onSideUpDown` (lineThrough origin $ Point2 (-10) 5) --- Above --- >>> Point2 5 5 `onSideUpDown` (verticalLine 10) --- Below --- >>> Point2 5 5 `onSideUpDown` (lineThrough origin $ Point2 (-3) (-3)) --- On -onSideUpDown :: (Ord r, Num r) => Point 2 r -> Line 2 r -> SideTestUpDown -q `onSideUpDown` (Line p v) = let r = p .+^ v - f z = (z^.xCoord, -z^.yCoord) - minBy g a b = F.minimumBy (comparing g) [a,b] - maxBy g a b = F.maximumBy (comparing g) [a,b] - in case ccw (minBy f p r) (maxBy f p r) q of - CCW -> Above - CW -> Below - CoLinear -> On +class OnSideUpDownTest t where + onSideUpDown :: (d ~ Dimension t, r ~ NumType t, Ord r, Num r) + => Point d r -> t -> SideTestUpDown + +instance OnSideUpDownTest (Line 2 r) where + -- | Given a point q and a line l, compute to which side of l q lies. For + -- vertical lines the left side of the line is interpeted as below. + -- + -- >>> Point2 10 10 `onSideUpDown` (lineThrough origin $ Point2 10 5) + -- Above + -- >>> Point2 10 10 `onSideUpDown` (lineThrough origin $ Point2 (-10) 5) + -- Above + -- >>> Point2 5 5 `onSideUpDown` (verticalLine 10) + -- Below + -- >>> Point2 5 5 `onSideUpDown` (lineThrough origin $ Point2 (-3) (-3)) + -- On + q `onSideUpDown` (Line p v) = let r = p .+^ v + f z = (z^.xCoord, -z^.yCoord) + minBy g a b = F.minimumBy (comparing g) [a,b] + maxBy g a b = F.maximumBy (comparing g) [a,b] + in case ccw (minBy f p r) (maxBy f p r) q of + CCW -> Above + CW -> Below + CoLinear -> On -- | Result of a side test data SideTest = LeftSide | OnLine | RightSide deriving (Show,Read,Eq,Ord) diff --git a/src/Data/Geometry/LineSegment.hs b/src/Data/Geometry/LineSegment.hs --- a/src/Data/Geometry/LineSegment.hs +++ b/src/Data/Geometry/LineSegment.hs @@ -14,6 +14,7 @@ , pattern LineSegment , pattern LineSegment' , pattern ClosedLineSegment + , pattern OpenLineSegment , endPoints , _SubLine @@ -26,6 +27,8 @@ , segmentLength , sqDistanceToSeg, sqDistanceToSegArg , flipSegment + + , interpolate ) where import Control.Arrow ((&&&)) @@ -45,7 +48,7 @@ import Data.Vinyl import Data.Vinyl.CoRec import GHC.TypeLits -import Test.QuickCheck +import Test.QuickCheck(Arbitrary(..)) -------------------------------------------------------------------------------- -- * d-dimensional LineSegments @@ -77,12 +80,16 @@ pattern LineSegment' s t <- ((^.start) &&& (^.end) -> (s,t)) {-# COMPLETE LineSegment' #-} -pattern ClosedLineSegment :: Point d r :+ p - -> Point d r :+ p - -> LineSegment d p r +pattern ClosedLineSegment :: Point d r :+ p -> Point d r :+ p -> LineSegment d p r pattern ClosedLineSegment s t = GLineSegment (ClosedInterval s t) {-# COMPLETE ClosedLineSegment #-} +pattern OpenLineSegment :: Point d r :+ p -> Point d r :+ p -> LineSegment d p r +pattern OpenLineSegment s t = GLineSegment (OpenInterval s t) +{-# COMPLETE OpenLineSegment #-} + + + type instance Dimension (LineSegment d p r) = d type instance NumType (LineSegment d p r) = r @@ -187,7 +194,6 @@ :& (H $ coRec . subLineToSegment) :& RNil - instance (Ord r, Fractional r) => (LineSegment 2 p r) `IsIntersectableWith` (Line 2 r) where nonEmptyIntersection = defaultNonEmptyIntersection @@ -292,3 +298,18 @@ -- (q&unEndPoint %~ ff) -- ss'' = ss'^._SubLine + +-- | Linearly interpolate the two endpoints with a value in the range [0,1] +-- +-- >>> interpolate 0.5 $ ClosedLineSegment (ext $ origin) (ext $ Point2 10.0 10.0) +-- Point2 [5.0,5.0] +-- >>> interpolate 0.1 $ ClosedLineSegment (ext $ origin) (ext $ Point2 10.0 10.0) +-- Point2 [1.0,1.0] +-- >>> interpolate 0 $ ClosedLineSegment (ext $ origin) (ext $ Point2 10.0 10.0) +-- Point2 [0.0,0.0] +-- >>> interpolate 1 $ ClosedLineSegment (ext $ origin) (ext $ Point2 10.0 10.0) +-- Point2 [10.0,10.0] +interpolate :: (Fractional r, Arity d) => r -> LineSegment d p r -> Point d r +interpolate t (LineSegment' p q) = Point $ (asV p ^* (1-t)) ^+^ (asV q ^* t) + where + asV = (^.core.vector) diff --git a/src/Data/Geometry/Matrix.hs b/src/Data/Geometry/Matrix.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Matrix.hs @@ -0,0 +1,79 @@ +-------------------------------------------------------------------------------- +-- | +-- Module : Data.Geometry.Matrix +-- Copyright : (C) Frank Staals +-- License : see the LICENSE file +-- Maintainer : Frank Staals +-- +-- type-indexed matrices. +-- +-------------------------------------------------------------------------------- +module Data.Geometry.Matrix( + Matrix(Matrix) + , identityMatrix + + , multM + , mult + + , Invertible(..) + , HasDeterminant(..) + ) where + +import Control.Lens (imap) +import Data.Geometry.Matrix.Internal (mkRow) +import Data.Geometry.Vector +import Linear.Matrix ((!*),(!*!)) +import qualified Linear.Matrix as Lin +import Unsafe.Coerce (unsafeCoerce) + +-------------------------------------------------------------------------------- +-- * Matrices + +-- | a matrix of n rows, each of m columns, storing values of type r +newtype Matrix n m r = Matrix (Vector n (Vector m r)) + +deriving instance (Show r, Arity n, Arity m) => Show (Matrix n m r) +deriving instance (Eq r, Arity n, Arity m) => Eq (Matrix n m r) +deriving instance (Ord r, Arity n, Arity m) => Ord (Matrix n m r) +deriving instance (Arity n, Arity m) => Functor (Matrix n m) +deriving instance (Arity n, Arity m) => Foldable (Matrix n m) +deriving instance (Arity n, Arity m) => Traversable (Matrix n m) + +multM :: (Arity r, Arity c, Arity c', Num a) => Matrix r c a -> Matrix c c' a -> Matrix r c' a +(Matrix a) `multM` (Matrix b) = Matrix $ a !*! b + +mult :: (Arity m, Arity n, Num r) => Matrix n m r -> Vector m r -> Vector n r +(Matrix m) `mult` v = m !* v + +-- | Produces the Identity Matrix +identityMatrix :: (Arity d, Num r) => Matrix d d r +identityMatrix = Matrix $ imap mkRow (pure 1) + +class Invertible n r where + inverse' :: Matrix n n r -> Matrix n n r + +instance Fractional r => Invertible 2 r where + -- >>> inverse' $ Matrix $ Vector2 (Vector2 1 2) (Vector2 3 4.0) + -- Matrix Vector2 [Vector2 [-2.0,1.0],Vector2 [1.5,-0.5]] + inverse' (Matrix m) = Matrix . unsafeCoerce . Lin.inv22 . unsafeCoerce $ m + +instance Fractional r => Invertible 3 r where + -- >>> inverse' $ Matrix $ Vector3 (Vector3 1 2 4) (Vector3 4 2 2) (Vector3 1 1 1.0) + -- Matrix Vector3 [Vector3 [0.0,0.5,-1.0],Vector3 [-0.5,-0.75,3.5],Vector3 [0.5,0.25,-1.5]] + inverse' (Matrix m) = Matrix . unsafeCoerce . Lin.inv33 . unsafeCoerce $ m + +instance Fractional r => Invertible 4 r where + inverse' (Matrix m) = Matrix . unsafeCoerce . Lin.inv44 . unsafeCoerce $ m + + +class Arity d => HasDeterminant d where + det :: Num r => Matrix d d r -> r + +instance HasDeterminant 1 where + det (Matrix (Vector1 (Vector1 x))) = x +instance HasDeterminant 2 where + det = Lin.det22 . unsafeCoerce +instance HasDeterminant 3 where + det = Lin.det33 . unsafeCoerce +instance HasDeterminant 4 where + det = Lin.det44 . unsafeCoerce diff --git a/src/Data/Geometry/Matrix/Internal.hs b/src/Data/Geometry/Matrix/Internal.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Matrix/Internal.hs @@ -0,0 +1,14 @@ +{-# LANGUAGE Unsafe #-} +module Data.Geometry.Matrix.Internal where + +import Control.Lens (set) +import Data.Geometry.Vector +import qualified Data.Vector.Fixed as FV + +-------------------------------------------------------------------------------- +-- * Helper functions to easily create matrices + +-- | Creates a row with zeroes everywhere, except at position i, where the +-- value is the supplied value. +mkRow :: forall d r. (Arity d, Num r) => Int -> r -> Vector d r +mkRow i x = set (FV.element i) x zero 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 @@ -24,7 +24,6 @@ import Data.Geometry.PlanarSubdivision.Basic import Data.Geometry.PlanarSubdivision.Merge import Data.Geometry.Polygon -import qualified Data.PlaneGraph as PG import Data.Proxy 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 @@ -423,13 +423,23 @@ d'' = PG.nextIncidentEdge d' g in g^.dataOf d'' --- | All incoming edges incident to vertex v, in counterclockwise order around v. +-- | All edges incident to vertex v in incoming direction +-- (i.e. pointing into v) in counterclockwise order around v. +-- +-- running time: \(O(k)\), where \(k) is the total number of incident edges of v incomingEdges :: VertexId' s -> PlanarSubdivision s v e f r -> V.Vector (Dart s) -incomingEdges v ps = V.filter (not . isPositive) $ incidentEdges v ps +incomingEdges v ps = orient <$> incidentEdges v ps + where + orient d = if headOf d ps == v then d else twin d --- | All outgoing edges incident to vertex v, in counterclockwise order around v. +-- | All edges incident to vertex v in outgoing direction +-- (i.e. pointing away from v) in counterclockwise order around v. +-- +-- running time: \(O(k)\), where \(k) is the total number of incident edges of v outgoingEdges :: VertexId' s -> PlanarSubdivision s v e f r -> V.Vector (Dart s) -outgoingEdges v ps = V.filter isPositive $ incidentEdges v ps +outgoingEdges v ps = orient <$> incidentEdges v ps + where + orient d = if tailOf d ps == v then d else twin d -- | Gets the neighbours of a particular vertex, in counterclockwise order @@ -437,10 +447,7 @@ -- -- running time: \(O(k)\), where \(k\) is the output size neighboursOf :: VertexId' s -> PlanarSubdivision s v e f r -> V.Vector (VertexId' s) -neighboursOf v ps = otherVtx <$> incidentEdges v ps - where - otherVtx d = let u = tailOf d ps in if u == v then headOf d ps else u - +neighboursOf v ps = flip tailOf ps <$> incomingEdges v ps -- | The face to the left of the dart -- diff --git a/src/Data/Geometry/Point.hs b/src/Data/Geometry/Point.hs --- a/src/Data/Geometry/Point.hs +++ b/src/Data/Geometry/Point.hs @@ -13,18 +13,17 @@ module Data.Geometry.Point( Point(..) , origin, vector , pointFromList - - , coord , unsafeCoord - , projectPoint + , pattern Point1 , pattern Point2 , pattern Point3 , xCoord, yCoord, zCoord , PointFunctor(..) - , CCW(..), ccw, ccw' + , CCW, ccw, ccw' + , pattern CCW, pattern CW, pattern CoLinear , ccwCmpAround, cwCmpAround, ccwCmpAroundWith, cwCmpAroundWith , sortAround, insertIntoCyclicOrder @@ -34,394 +33,12 @@ , cmpByDistanceTo , squaredEuclideanDist, euclideanDist - ) where -import Control.DeepSeq -import Control.Lens -import Data.Aeson -import qualified Data.CircularList as C -import qualified Data.CircularList.Util as CU -import Data.Ext -import qualified Data.Foldable as F -import Data.Geometry.Properties -import Data.Geometry.Vector -import qualified Data.Geometry.Vector as Vec -import qualified Data.List as L -import Data.Ord (comparing) -import Data.Proxy -import GHC.Generics (Generic) -import GHC.TypeLits -import Test.QuickCheck (Arbitrary) -import Text.ParserCombinators.ReadP (ReadP, string,pfail) -import Text.ParserCombinators.ReadPrec (lift) -import Text.Read (Read(..),readListPrecDefault, readPrec_to_P,minPrec) - --------------------------------------------------------------------------------- --- $setup --- >>> :{ --- let myVector :: Vector 3 Int --- myVector = Vector3 1 2 3 --- myPoint = Point myVector --- :} - - --------------------------------------------------------------------------------- --- * A d-dimensional Point - --- | A d-dimensional point. -newtype Point d r = Point { toVec :: Vector d r } deriving (Generic) - -instance (Show r, Arity d) => Show (Point d r) where - show (Point v) = mconcat [ "Point", show $ F.length v , " " - , show $ F.toList v - ] -instance (Read r, Arity d) => Read (Point d r) where - readPrec = lift readPt - readListPrec = readListPrecDefault - -readPt :: forall d r. (Arity d, Read r) => ReadP (Point d r) -readPt = do let d = natVal (Proxy :: Proxy d) - _ <- string $ "Point" <> show d <> " " - rs <- readPrec_to_P readPrec minPrec - case pointFromList rs of - Just p -> pure p - _ -> pfail - -deriving instance (Eq r, Arity d) => Eq (Point d r) -deriving instance (Ord r, Arity d) => Ord (Point d r) -deriving instance Arity d => Functor (Point d) -deriving instance Arity d => Foldable (Point d) -deriving instance Arity d => Traversable (Point d) -deriving instance (Arity d, NFData r) => NFData (Point d r) -deriving instance (Arity d, Arbitrary r) => Arbitrary (Point d r) - -type instance NumType (Point d r) = r -type instance Dimension (Point d r) = d - -instance Arity d => Affine (Point d) where - type Diff (Point d) = Vector d - - p .-. q = toVec p ^-^ toVec q - p .+^ v = Point $ toVec p ^+^ v - -instance (FromJSON r, Arity d, KnownNat d) => FromJSON (Point d r) where - parseJSON = fmap Point . parseJSON - -instance (ToJSON r, Arity d) => ToJSON (Point d r) where - toJSON = toJSON . toVec - toEncoding = toEncoding . toVec - --- | Point representing the origin in d dimensions --- --- >>> origin :: Point 4 Int --- Point4 [0,0,0,0] -origin :: (Arity d, Num r) => Point d r -origin = Point $ pure 0 - - --- ** Accessing points - --- | Lens to access the vector corresponding to this point. --- --- >>> (Point3 1 2 3) ^. vector --- Vector3 [1,2,3] --- >>> origin & vector .~ Vector3 1 2 3 --- Point3 [1,2,3] -vector :: Lens' (Point d r) (Vector d r) -vector = lens toVec (const Point) - - --- | Get the coordinate in a given dimension. This operation is unsafe in the --- sense that no bounds are checked. Consider using `coord` instead. --- --- --- >>> Point3 1 2 3 ^. unsafeCoord 2 --- 2 -unsafeCoord :: Arity d => Int -> Lens' (Point d r) r -unsafeCoord i = vector . singular (ix (i-1)) - -- Points are 1 indexed, vectors are 0 indexed - --- | Get the coordinate in a given dimension --- --- >>> Point3 1 2 3 ^. coord (C :: C 2) --- 2 --- >>> Point3 1 2 3 & coord (C :: C 1) .~ 10 --- Point3 [10,2,3] --- >>> Point3 1 2 3 & coord (C :: C 3) %~ (+1) --- Point3 [1,2,4] -coord :: forall proxy i d r. (1 <= i, i <= d, ((i - 1) + 1) ~ i - , Arity (i - 1), Arity d - ) => proxy i -> Lens' (Point d r) r -coord _ = vector . Vec.element (Proxy :: Proxy (i-1)) -{-# INLINABLE coord #-} - - --- somehow these rules don't fire --- {-# SPECIALIZE coord :: C 1 -> Lens' (Point 2 r) r#-} --- {-# SPECIALIZE coord :: C 2 -> Lens' (Point 2 r) r#-} - - --- | Constructs a point from a list of coordinates --- --- >>> pointFromList [1,2,3] :: Maybe (Point 3 Int) --- Just Point3 [1,2,3] -pointFromList :: Arity d => [r] -> Maybe (Point d r) -pointFromList = fmap Point . Vec.vectorFromList - - --- | Project a point down into a lower dimension. -projectPoint :: (Arity i, Arity d, i <= d) => Point d r -> Point i r -projectPoint = Point . prefix . toVec - --------------------------------------------------------------------------------- --- * Convenience functions to construct 2 and 3 dimensional points - - --- | We provide pattern synonyms Point2 and Point3 for 2 and 3 dimensional points. i.e. --- we can write: --- --- >>> :{ --- let --- f :: Point 2 r -> r --- f (Point2 x y) = x --- in f (Point2 1 2) --- :} --- 1 --- --- if we want. -pattern Point2 :: r -> r -> Point 2 r -pattern Point2 x y = Point (Vector2 x y) -{-# COMPLETE Point2 #-} - --- | Similarly, we can write: --- --- >>> :{ --- let --- g :: Point 3 r -> r --- g (Point3 x y z) = z --- in g myPoint --- :} --- 3 -pattern Point3 :: r -> r -> r -> Point 3 r -pattern Point3 x y z = (Point (Vector3 x y z)) -{-# COMPLETE Point3 #-} - --- | Shorthand to access the first coordinate C 1 --- --- >>> Point3 1 2 3 ^. xCoord --- 1 --- >>> Point2 1 2 & xCoord .~ 10 --- Point2 [10,2] -xCoord :: (1 <= d, Arity d) => Lens' (Point d r) r -xCoord = coord (C :: C 1) -{-# INLINABLE xCoord #-} - --- | Shorthand to access the second coordinate C 2 --- --- >>> Point2 1 2 ^. yCoord --- 2 --- >>> Point3 1 2 3 & yCoord %~ (+1) --- Point3 [1,3,3] -yCoord :: (2 <= d, Arity d) => Lens' (Point d r) r -yCoord = coord (C :: C 2) -{-# INLINABLE yCoord #-} - --- | Shorthand to access the third coordinate C 3 --- --- >>> Point3 1 2 3 ^. zCoord --- 3 --- >>> Point3 1 2 3 & zCoord %~ (+1) --- Point3 [1,2,4] -zCoord :: (3 <= d, Arity d) => Lens' (Point d r) r -zCoord = coord (C :: C 3) -{-# INLINABLE zCoord #-} - - --------------------------------------------------------------------------------- --- * Point Functors - --- | Types that we can transform by mapping a function on each point in the structure -class PointFunctor g where - pmap :: (Point (Dimension (g r)) r -> Point (Dimension (g s)) s) -> g r -> g s - - -- pemap :: (d ~ Dimension (g r)) => (Point d r :+ p -> Point d s :+ p) -> g r -> g s - -- pemap = - -instance PointFunctor (Point d) where - pmap f = f - - --------------------------------------------------------------------------------- --- * Functions specific to Two Dimensional points - -data CCW = CCW | CoLinear | CW - deriving (Show,Eq) - --- | Given three points p q and r determine the orientation when going from p to r via q. -ccw :: (Ord r, Num r) => Point 2 r -> Point 2 r -> Point 2 r -> CCW -ccw p q r = case z `compare` 0 of - LT -> CW - GT -> CCW - EQ -> CoLinear - where - Vector2 ux uy = q .-. p - Vector2 vx vy = r .-. p - z = ux * vy - uy * vx - --- | Given three points p q and r determine the orientation when going from p to r via q. -ccw' :: (Ord r, Num r) => Point 2 r :+ a -> Point 2 r :+ b -> Point 2 r :+ c -> CCW -ccw' p q r = ccw (p^.core) (q^.core) (r^.core) - --- | Sort the points arround the given point p in counter clockwise order with --- respect to the rightward horizontal ray starting from p. If two points q --- and r are colinear with p, the closest one to p is reported first. --- running time: O(n log n) -sortAround :: (Ord r, Num r) - => Point 2 r :+ q -> [Point 2 r :+ p] -> [Point 2 r :+ p] -sortAround c = L.sortBy (ccwCmpAround c <> cmpByDistanceTo c) - - --- | Quadrants of two dimensional points. in CCW order -data Quadrant = TopRight | TopLeft | BottomLeft | BottomRight - deriving (Show,Read,Eq,Ord,Enum,Bounded) - --- | Quadrants around point c; quadrants are closed on their "previous" --- boundary (i..e the boundary with the previous quadrant in the CCW order), --- open on next boundary. The origin itself is assigned the topRight quadrant -quadrantWith :: (Ord r, 1 <= d, 2 <= d, Arity d) - => Point d r :+ q -> Point d r :+ p -> Quadrant -quadrantWith (c :+ _) (p :+ _) = case ( (c^.xCoord) `compare` (p^.xCoord) - , (c^.yCoord) `compare` (p^.yCoord) ) of - (EQ, EQ) -> TopRight - (LT, EQ) -> TopRight - (LT, LT) -> TopRight - (EQ, LT) -> TopLeft - (GT, LT) -> TopLeft - (GT, EQ) -> BottomLeft - (GT, GT) -> BottomLeft - (EQ, GT) -> BottomRight - (LT, GT) -> BottomRight - --- | Quadrants with respect to the origin -quadrant :: (Ord r, Num r, 1 <= d, 2 <= d, Arity d) => Point d r :+ p -> Quadrant -quadrant = quadrantWith (ext origin) - --- | Given a center point c, and a set of points, partition the points into --- quadrants around c (based on their x and y coordinates). The quadrants are --- reported in the order topLeft, topRight, bottomLeft, bottomRight. The points --- are in the same order as they were in the original input lists. --- Points with the same x-or y coordinate as p, are "rounded" to above. -partitionIntoQuadrants :: (Ord r, 1 <= d, 2 <= d, Arity d) - => Point d r :+ q - -> [Point d r :+ p] - -> ( [Point d r :+ p], [Point d r :+ p] - , [Point d r :+ p], [Point d r :+ p] - ) -partitionIntoQuadrants c pts = (topL, topR, bottomL, bottomR) - where - (below',above') = L.partition (on yCoord) pts - (bottomL,bottomR) = L.partition (on xCoord) below' - (topL,topR) = L.partition (on xCoord) above' - - on l q = q^.core.l < c^.core.l - - - --- | Given a zero vector z, a center c, and two points p and q, --- compute the ccw ordering of p and q around c with this vector as zero --- direction. --- --- pre: the points p,q /= c -ccwCmpAroundWith :: (Ord r, Num r) - => Vector 2 r - -> Point 2 r :+ c - -> Point 2 r :+ a -> Point 2 r :+ b - -> Ordering -ccwCmpAroundWith z@(Vector2 zx zy) (c :+ _) (q :+ _) (r :+ _) = - case (ccw c a q, ccw c a r) of - (CCW,CCW) -> cmp - (CCW,CW) -> LT - (CCW,CoLinear) | onZero r -> GT - | otherwise -> LT - - (CW, CCW) -> GT - (CW, CW) -> cmp - (CW, CoLinear) -> GT - - (CoLinear, CCW) | onZero q -> LT - | otherwise -> GT - - (CoLinear, CW) -> LT - (CoLinear,CoLinear) -> case (onZero q, onZero r) of - (True, True) -> EQ - (False, False) -> EQ - (True, False) -> LT - (False, True) -> GT - where - a = c .+^ z - b = c .+^ Vector2 (-zy) zx - -- b is on a perpendicular vector to z - - -- test if the point lies on the ray defined by z, starting in c - onZero d = case ccw c b d of - CCW -> False - CW -> True - CoLinear -> True -- this shouldh appen only when you ask for c itself - - cmp = case ccw c q r of - CCW -> LT - CW -> GT - CoLinear -> EQ - --- | Given a zero vector z, a center c, and two points p and q, --- compute the cw ordering of p and q around c with this vector as zero --- direction. --- --- pre: the points p,q /= c -cwCmpAroundWith :: (Ord r, Num r) - => Vector 2 r - -> Point 2 r :+ a - -> Point 2 r :+ b -> Point 2 r :+ c - -> Ordering -cwCmpAroundWith z c = flip (ccwCmpAroundWith z c) - - - --- | Compare by distance to the first argument -cmpByDistanceTo :: (Ord r, Num r, Arity d) - => Point d r :+ c -> Point d r :+ p -> Point d r :+ q -> Ordering -cmpByDistanceTo (c :+ _) p q = comparing (squaredEuclideanDist c) (p^.core) (q^.core) - - --- | Counter clockwise ordering of the points around c. Points are ordered with --- respect to the positive x-axis. -ccwCmpAround :: (Num r, Ord r) - => Point 2 r :+ qc -> Point 2 r :+ p -> Point 2 r :+ q -> Ordering -ccwCmpAround = ccwCmpAroundWith (Vector2 1 0) - --- | Clockwise ordering of the points around c. Points are ordered with --- respect to the positive x-axis. -cwCmpAround :: (Num r, Ord r) - => Point 2 r :+ qc -> Point 2 r :+ p -> Point 2 r :+ q -> Ordering -cwCmpAround = cwCmpAroundWith (Vector2 1 0) - - --- | Given a center c, a new point p, and a list of points ps, sorted in --- counter clockwise order around c. Insert p into the cyclic order. The focus --- of the returned cyclic list is the new point p. --- --- running time: O(n) -insertIntoCyclicOrder :: (Ord r, Num r) - => Point 2 r :+ q -> Point 2 r :+ p - -> C.CList (Point 2 r :+ p) -> C.CList (Point 2 r :+ p) -insertIntoCyclicOrder c = CU.insertOrdBy (ccwCmpAround c <> cmpByDistanceTo c) - - --- | Squared Euclidean distance between two points -squaredEuclideanDist :: (Num r, Arity d) => Point d r -> Point d r -> r -squaredEuclideanDist = qdA + , AsAPoint(..), coord, unsafeCoord, vector' + ) where --- | Euclidean distance between two points -euclideanDist :: (Floating r, Arity d) => Point d r -> Point d r -> r -euclideanDist = distanceA +import Data.Geometry.Point.Class +import Data.Geometry.Point.Internal hiding (coord, unsafeCoord) +import Data.Geometry.Point.Orientation.Degenerate +import Data.Geometry.Point.Quadrants diff --git a/src/Data/Geometry/Point/Class.hs b/src/Data/Geometry/Point/Class.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Point/Class.hs @@ -0,0 +1,66 @@ +module Data.Geometry.Point.Class where + +import Control.Lens +import Data.Geometry.Point.Internal (Point) +import qualified Data.Geometry.Point.Internal as Internal +import Data.Geometry.Vector +import GHC.TypeNats + +-------------------------------------------------------------------------------- + +-- $setup +-- >>> import Data.Geometry.Point.Internal (pattern Point2, pattern Point3) + +class ToAPoint point d r where + toPoint :: Prism' (point d r) (Point d r) + +class AsAPoint p where + asAPoint :: Lens (p d r) (p d' r') (Point d r) (Point d' r') + +vector' :: AsAPoint p => Lens (p d r) (p d r') (Vector d r) (Vector d r') +vector' = asAPoint . lens Internal.toVec (const Internal.Point) + +coord :: (1 <= i, i <= d, KnownNat i, Arity d, AsAPoint p) => proxy i -> Lens' (p d r) r +coord i = asAPoint.Internal.coord i + +unsafeCoord :: (Arity d, AsAPoint p) => Int -> Lens' (p d r) r +unsafeCoord i = asAPoint.Internal.unsafeCoord i + +instance ToAPoint Point d r where + toPoint = prism' id Just + +instance AsAPoint Point where + asAPoint = id + + + + +-- | Shorthand to access the first coordinate C 1 +-- +-- >>> Point3 1 2 3 ^. xCoord +-- 1 +-- >>> Point2 1 2 & xCoord .~ 10 +-- Point2 [10,2] +xCoord :: (1 <= d, Arity d, AsAPoint point) => Lens' (point d r) r +xCoord = coord (C :: C 1) +{-# INLINABLE xCoord #-} + +-- | Shorthand to access the second coordinate C 2 +-- +-- >>> Point2 1 2 ^. yCoord +-- 2 +-- >>> Point3 1 2 3 & yCoord %~ (+1) +-- Point3 [1,3,3] +yCoord :: (2 <= d, Arity d, AsAPoint point) => Lens' (point d r) r +yCoord = coord (C :: C 2) +{-# INLINABLE yCoord #-} + +-- | Shorthand to access the third coordinate C 3 +-- +-- >>> Point3 1 2 3 ^. zCoord +-- 3 +-- >>> Point3 1 2 3 & zCoord %~ (+1) +-- Point3 [1,2,4] +zCoord :: (3 <= d, Arity d,AsAPoint point) => Lens' (point d r) r +zCoord = coord (C :: C 3) +{-# INLINABLE zCoord #-} diff --git a/src/Data/Geometry/Point/Internal.hs b/src/Data/Geometry/Point/Internal.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Point/Internal.hs @@ -0,0 +1,252 @@ +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE UndecidableInstances #-} +-------------------------------------------------------------------------------- +-- | +-- Module : Data.Geometry.Point +-- Copyright : (C) Frank Staals +-- License : see the LICENSE file +-- Maintainer : Frank Staals +-- +-- \(d\)-dimensional points. +-- +-------------------------------------------------------------------------------- +module Data.Geometry.Point.Internal + ( Point(..) + , origin, vector + , pointFromList + + , coord , unsafeCoord + + , projectPoint + + , pattern Point1 + , pattern Point2 + , pattern Point3 + , PointFunctor(..) + + , cmpByDistanceTo + , squaredEuclideanDist, euclideanDist + ) where + +import Control.DeepSeq +import Control.Lens +import Data.Aeson +import Data.Ext +import qualified Data.Foldable as F +import Data.Geometry.Properties +import Data.Geometry.Vector +import qualified Data.Geometry.Vector as Vec +import Data.Hashable +import Data.Ord (comparing) +import Data.Proxy +import GHC.Generics (Generic) +import GHC.TypeLits +import System.Random (Random(..)) +import Test.QuickCheck (Arbitrary) +import Text.ParserCombinators.ReadP (ReadP, string,pfail) +import Text.ParserCombinators.ReadPrec (lift) +import Text.Read (Read(..),readListPrecDefault, readPrec_to_P,minPrec) + + +-------------------------------------------------------------------------------- +-- $setup +-- >>> :{ +-- let myVector :: Vector 3 Int +-- myVector = Vector3 1 2 3 +-- myPoint = Point myVector +-- :} + + +-------------------------------------------------------------------------------- +-- * A d-dimensional Point + +-- | A d-dimensional point. +newtype Point d r = Point { toVec :: Vector d r } deriving (Generic) + +instance (Show r, Arity d) => Show (Point d r) where + show (Point v) = mconcat [ "Point", show $ F.length v , " " + , show $ F.toList v + ] +instance (Read r, Arity d) => Read (Point d r) where + readPrec = lift readPt + readListPrec = readListPrecDefault + +readPt :: forall d r. (Arity d, Read r) => ReadP (Point d r) +readPt = do let d = natVal (Proxy :: Proxy d) + _ <- string $ "Point" <> show d <> " " + rs <- readPrec_to_P readPrec minPrec + case pointFromList rs of + Just p -> pure p + _ -> pfail + +deriving instance (Eq r, Arity d) => Eq (Point d r) +deriving instance (Ord r, Arity d) => Ord (Point d r) +deriving instance Arity d => Functor (Point d) +deriving instance Arity d => Foldable (Point d) +deriving instance Arity d => Traversable (Point d) +deriving instance (Arity d, NFData r) => NFData (Point d r) +deriving instance (Arity d, Arbitrary r) => Arbitrary (Point d r) +deriving instance (Arity d, Hashable r) => Hashable (Point d r) +deriving instance (Arity d, Random r) => Random (Point d r) + + +type instance NumType (Point d r) = r +type instance Dimension (Point d r) = d + +instance Arity d => Affine (Point d) where + type Diff (Point d) = Vector d + + p .-. q = toVec p ^-^ toVec q + p .+^ v = Point $ toVec p ^+^ v + +instance (FromJSON r, Arity d, KnownNat d) => FromJSON (Point d r) where + parseJSON = fmap Point . parseJSON + +instance (ToJSON r, Arity d) => ToJSON (Point d r) where + toJSON = toJSON . toVec + toEncoding = toEncoding . toVec + +-- | Point representing the origin in d dimensions +-- +-- >>> origin :: Point 4 Int +-- Point4 [0,0,0,0] +origin :: (Arity d, Num r) => Point d r +origin = Point $ pure 0 + + +-- ** Accessing points + +-- | Lens to access the vector corresponding to this point. +-- +-- >>> (Point3 1 2 3) ^. vector +-- Vector3 [1,2,3] +-- >>> origin & vector .~ Vector3 1 2 3 +-- Point3 [1,2,3] +vector :: Lens' (Point d r) (Vector d r) +vector = lens toVec (const Point) +{-# INLINABLE vector #-} + +-- | Get the coordinate in a given dimension. This operation is unsafe in the +-- sense that no bounds are checked. Consider using `coord` instead. +-- +-- +-- >>> Point3 1 2 3 ^. unsafeCoord 2 +-- 2 +unsafeCoord :: Arity d => Int -> Lens' (Point d r) r +unsafeCoord i = vector . singular (ix (i-1)) + -- Points are 1 indexed, vectors are 0 indexed +{-# INLINABLE unsafeCoord #-} + +-- | Get the coordinate in a given dimension +-- +-- >>> Point3 1 2 3 ^. coord (C :: C 2) +-- 2 +-- >>> Point3 1 2 3 & coord (C :: C 1) .~ 10 +-- Point3 [10,2,3] +-- >>> Point3 1 2 3 & coord (C :: C 3) %~ (+1) +-- Point3 [1,2,4] +coord :: forall proxy i d r. (1 <= i, i <= d, Arity d, KnownNat i) + => proxy i -> Lens' (Point d r) r +coord _ = unsafeCoord $ fromIntegral (natVal $ C @i) +{-# INLINABLE coord #-} + + -- somehow these rules don't fire +-- {-# SPECIALIZE coord :: C 1 -> Lens' (Point 2 r) r#-} +-- {-# SPECIALIZE coord :: C 2 -> Lens' (Point 2 r) r#-} +-- {-# SPECIALIZE coord :: C 3 -> Lens' (Point 3 r) r#-} + + +-- | Constructs a point from a list of coordinates. The length of the +-- list has to match the dimension exactly. +-- +-- >>> pointFromList [1,2,3] :: Maybe (Point 3 Int) +-- Just Point3 [1,2,3] +-- >>> pointFromList [1] :: Maybe (Point 3 Int) +-- Nothing +-- >>> pointFromList [1,2,3,4] :: Maybe (Point 3 Int) +-- Nothing +pointFromList :: Arity d => [r] -> Maybe (Point d r) +pointFromList = fmap Point . Vec.vectorFromList + + +-- | Project a point down into a lower dimension. +projectPoint :: (Arity i, Arity d, i <= d) => Point d r -> Point i r +projectPoint = Point . prefix . toVec + +-------------------------------------------------------------------------------- +-- * Convenience functions to construct 1, 2 and 3 dimensional points + +-- | We provide pattern synonyms for 1, 2 and 3 dimensional points. i.e. +-- we can write: +-- +-- +-- >>> :{ +-- let +-- f :: Num r => Point 1 r -> r +-- f (Point1 x) = x + 1 +-- in f (Point1 1) +-- :} +-- 2 +pattern Point1 :: r -> Point 1 r +pattern Point1 x = Point (Vector1 x) +{-# COMPLETE Point1 #-} + + +-- | Pattern synonym for 2 dimensional points +-- +-- >>> :{ +-- let +-- f :: Point 2 r -> r +-- f (Point2 x y) = x +-- in f (Point2 1 2) +-- :} +-- 1 +pattern Point2 :: r -> r -> Point 2 r +pattern Point2 x y = Point (Vector2 x y) +{-# COMPLETE Point2 #-} + +-- | Similarly, we can write: +-- +-- >>> :{ +-- let +-- g :: Point 3 r -> r +-- g (Point3 x y z) = z +-- in g myPoint +-- :} +-- 3 +pattern Point3 :: r -> r -> r -> Point 3 r +pattern Point3 x y z = (Point (Vector3 x y z)) +{-# COMPLETE Point3 #-} + +-------------------------------------------------------------------------------- +-- * Point Functors + +-- | Types that we can transform by mapping a function on each point in the structure +class PointFunctor g where + pmap :: (Point (Dimension (g r)) r -> Point (Dimension (g s)) s) -> g r -> g s + + -- pemap :: (d ~ Dimension (g r)) => (Point d r :+ p -> Point d s :+ p) -> g r -> g s + -- pemap = + +instance PointFunctor (Point d) where + pmap f = f + + +-------------------------------------------------------------------------------- +-- * Functions specific to Two Dimensional points + +-- | Compare by distance to the first argument +cmpByDistanceTo :: (Ord r, Num r, Arity d) + => Point d r :+ c -> Point d r :+ p -> Point d r :+ q -> Ordering +cmpByDistanceTo (c :+ _) p q = comparing (squaredEuclideanDist c) (p^.core) (q^.core) + + + + +-- | Squared Euclidean distance between two points +squaredEuclideanDist :: (Num r, Arity d) => Point d r -> Point d r -> r +squaredEuclideanDist = qdA + +-- | Euclidean distance between two points +euclideanDist :: (Floating r, Arity d) => Point d r -> Point d r -> r +euclideanDist = distanceA diff --git a/src/Data/Geometry/Point/Orientation.hs b/src/Data/Geometry/Point/Orientation.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Point/Orientation.hs @@ -0,0 +1,31 @@ +module Data.Geometry.Point.Orientation where + +import Algorithms.Geometry.SoS.Orientation +import Algorithms.Geometry.SoS.Sign +import Data.Ext +import Data.Geometry.Point.Internal +import Data.Geometry.Vector + +-------------------------------------------------------------------------------- + +-------------------------------------------------------------------------------- + +newtype StrictCCW = SCCW Sign deriving Eq + +pattern CCW :: StrictCCW +pattern CCW = SCCW Negative + +pattern CW :: StrictCCW +pattern CW = SCCW Positive +{-# COMPLETE CCW, CW #-} + +instance Show StrictCCW where + show = \case + CCW -> "CCW" + CW -> "CW" + + +-- | Given three points p q and r determine the orientation when going from p to r via q. +ccw :: (Ord r, Num r, Ord i) + => Point 2 r :+ i -> Point 2 r :+ i -> Point 2 r :+ i -> StrictCCW +ccw p q r = SCCW $ sideTest r (Vector2 p q) diff --git a/src/Data/Geometry/Point/Orientation/Degenerate.hs b/src/Data/Geometry/Point/Orientation/Degenerate.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Point/Orientation/Degenerate.hs @@ -0,0 +1,150 @@ +module Data.Geometry.Point.Orientation.Degenerate( + CCW(..) + , pattern CCW, pattern CW, pattern CoLinear + + , ccw, ccw' + + , sortAround + + , ccwCmpAroundWith, cwCmpAroundWith + , ccwCmpAround, cwCmpAround + + , insertIntoCyclicOrder + ) where + +import Control.Lens +import qualified Data.CircularList as C +import qualified Data.CircularList.Util as CU +import Data.Ext +import Data.Geometry.Point.Internal +import Data.Geometry.Vector +import qualified Data.List as L + +-------------------------------------------------------------------------------- + +-- | Data type for expressing the orientation of three points, with +-- the option of allowing Colinearities. +newtype CCW = CCWWrap Ordering deriving Eq + +pattern CCW :: CCW +pattern CCW = CCWWrap GT + +pattern CW :: CCW +pattern CW = CCWWrap LT + +pattern CoLinear :: CCW +pattern CoLinear = CCWWrap EQ +{-# COMPLETE CCW, CW, CoLinear #-} + +instance Show CCW where + show = \case + CCW -> "CCW" + CW -> "CW" + CoLinear -> "CoLinear" + + +-- | Given three points p q and r determine the orientation when going from p to r via q. +ccw :: (Ord r, Num r) => Point 2 r -> Point 2 r -> Point 2 r -> CCW +ccw p q r = CCWWrap $ z `compare` 0 + -- case z `compare` 0 of + -- LT -> CW + -- GT -> CCW + -- EQ -> CoLinear + where + Vector2 ux uy = q .-. p + Vector2 vx vy = r .-. p + z = ux * vy - uy * vx + +-- | Given three points p q and r determine the orientation when going from p to r via q. +ccw' :: (Ord r, Num r) => Point 2 r :+ a -> Point 2 r :+ b -> Point 2 r :+ c -> CCW +ccw' p q r = ccw (p^.core) (q^.core) (r^.core) + +-- | Sort the points arround the given point p in counter clockwise order with +-- respect to the rightward horizontal ray starting from p. If two points q +-- and r are colinear with p, the closest one to p is reported first. +-- running time: O(n log n) +sortAround :: (Ord r, Num r) + => Point 2 r :+ q -> [Point 2 r :+ p] -> [Point 2 r :+ p] +sortAround c = L.sortBy (ccwCmpAround c <> cmpByDistanceTo c) + + +-- | Given a zero vector z, a center c, and two points p and q, +-- compute the ccw ordering of p and q around c with this vector as zero +-- direction. +-- +-- pre: the points p,q /= c +ccwCmpAroundWith :: (Ord r, Num r) + => Vector 2 r + -> Point 2 r :+ c + -> Point 2 r :+ a -> Point 2 r :+ b + -> Ordering +ccwCmpAroundWith z@(Vector2 zx zy) (c :+ _) (q :+ _) (r :+ _) = + case (ccw c a q, ccw c a r) of + (CCW,CCW) -> cmp + (CCW,CW) -> LT + (CCW,CoLinear) | onZero r -> GT + | otherwise -> LT + + (CW, CCW) -> GT + (CW, CW) -> cmp + (CW, CoLinear) -> GT + + (CoLinear, CCW) | onZero q -> LT + | otherwise -> GT + + (CoLinear, CW) -> LT + (CoLinear,CoLinear) -> case (onZero q, onZero r) of + (True, True) -> EQ + (False, False) -> EQ + (True, False) -> LT + (False, True) -> GT + where + a = c .+^ z + b = c .+^ Vector2 (-zy) zx + -- b is on a perpendicular vector to z + + -- test if the point lies on the ray defined by z, starting in c + onZero d = case ccw c b d of + CCW -> False + CW -> True + CoLinear -> True -- this shouldh appen only when you ask for c itself + + cmp = case ccw c q r of + CCW -> LT + CW -> GT + CoLinear -> EQ + +-- | Given a zero vector z, a center c, and two points p and q, +-- compute the cw ordering of p and q around c with this vector as zero +-- direction. +-- +-- pre: the points p,q /= c +cwCmpAroundWith :: (Ord r, Num r) + => Vector 2 r + -> Point 2 r :+ a + -> Point 2 r :+ b -> Point 2 r :+ c + -> Ordering +cwCmpAroundWith z c = flip (ccwCmpAroundWith z c) + +-- | Counter clockwise ordering of the points around c. Points are ordered with +-- respect to the positive x-axis. +ccwCmpAround :: (Num r, Ord r) + => Point 2 r :+ qc -> Point 2 r :+ p -> Point 2 r :+ q -> Ordering +ccwCmpAround = ccwCmpAroundWith (Vector2 1 0) + +-- | Clockwise ordering of the points around c. Points are ordered with +-- respect to the positive x-axis. +cwCmpAround :: (Num r, Ord r) + => Point 2 r :+ qc -> Point 2 r :+ p -> Point 2 r :+ q -> Ordering +cwCmpAround = cwCmpAroundWith (Vector2 1 0) + + +-- | Given a center c, a new point p, and a list of points ps, sorted in +-- counter clockwise order around c. Insert p into the cyclic order. The focus +-- of the returned cyclic list is the new point p. +-- +-- running time: O(n) +insertIntoCyclicOrder :: (Ord r, Num r) + => Point 2 r :+ q -> Point 2 r :+ p + -> C.CList (Point 2 r :+ p) -> C.CList (Point 2 r :+ p) +insertIntoCyclicOrder c = CU.insertOrdBy (ccwCmpAround c <> cmpByDistanceTo c) diff --git a/src/Data/Geometry/Point/Quadrants.hs b/src/Data/Geometry/Point/Quadrants.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Point/Quadrants.hs @@ -0,0 +1,69 @@ +module Data.Geometry.Point.Quadrants where + +import Control.DeepSeq +import Control.Lens +import Data.Aeson +import Data.Ext +import qualified Data.Foldable as F +import Data.Geometry.Point.Class +import Data.Geometry.Point.Internal +import Data.Geometry.Properties +import Data.Geometry.Vector +import qualified Data.Geometry.Vector as Vec +import Data.Hashable +import qualified Data.List as L +import Data.Ord (comparing) +import Data.Proxy +import GHC.Generics (Generic) +import GHC.TypeLits +import System.Random (Random(..)) +import Test.QuickCheck (Arbitrary) +import Text.ParserCombinators.ReadP (ReadP, string,pfail) +import Text.ParserCombinators.ReadPrec (lift) +import Text.Read (Read(..),readListPrecDefault, readPrec_to_P,minPrec) + +-------------------------------------------------------------------------------- + +-- | Quadrants of two dimensional points. in CCW order +data Quadrant = TopRight | TopLeft | BottomLeft | BottomRight + deriving (Show,Read,Eq,Ord,Enum,Bounded) + +-- | Quadrants around point c; quadrants are closed on their "previous" +-- boundary (i..e the boundary with the previous quadrant in the CCW order), +-- open on next boundary. The origin itself is assigned the topRight quadrant +quadrantWith :: (Ord r, 1 <= d, 2 <= d, Arity d) + => Point d r :+ q -> Point d r :+ p -> Quadrant +quadrantWith (c :+ _) (p :+ _) = case ( (c^.xCoord) `compare` (p^.xCoord) + , (c^.yCoord) `compare` (p^.yCoord) ) of + (EQ, EQ) -> TopRight + (LT, EQ) -> TopRight + (LT, LT) -> TopRight + (EQ, LT) -> TopLeft + (GT, LT) -> TopLeft + (GT, EQ) -> BottomLeft + (GT, GT) -> BottomLeft + (EQ, GT) -> BottomRight + (LT, GT) -> BottomRight + +-- | Quadrants with respect to the origin +quadrant :: (Ord r, Num r, 1 <= d, 2 <= d, Arity d) => Point d r :+ p -> Quadrant +quadrant = quadrantWith (ext origin) + +-- | Given a center point c, and a set of points, partition the points into +-- quadrants around c (based on their x and y coordinates). The quadrants are +-- reported in the order topLeft, topRight, bottomLeft, bottomRight. The points +-- are in the same order as they were in the original input lists. +-- Points with the same x-or y coordinate as p, are "rounded" to above. +partitionIntoQuadrants :: (Ord r, 1 <= d, 2 <= d, Arity d) + => Point d r :+ q + -> [Point d r :+ p] + -> ( [Point d r :+ p], [Point d r :+ p] + , [Point d r :+ p], [Point d r :+ p] + ) +partitionIntoQuadrants c pts = (topL, topR, bottomL, bottomR) + where + (below',above') = L.partition (on yCoord) pts + (bottomL,bottomR) = L.partition (on xCoord) below' + (topL,topR) = L.partition (on xCoord) above' + + on l q = q^.core.l < c^.core.l diff --git a/src/Data/Geometry/PolyLine.hs b/src/Data/Geometry/PolyLine.hs --- a/src/Data/Geometry/PolyLine.hs +++ b/src/Data/Geometry/PolyLine.hs @@ -5,7 +5,9 @@ import Control.Lens import Data.Aeson +import Data.Bifoldable import Data.Bifunctor +import Data.Bitraversable import Data.Ext import qualified Data.Foldable as F import Data.Geometry.Box @@ -17,10 +19,18 @@ import Data.LSeq (LSeq, pattern (:<|)) import qualified Data.LSeq as LSeq import qualified Data.List.NonEmpty as NE -import GHC.Generics(Generic) +import GHC.Generics (Generic) import GHC.TypeLits -------------------------------------------------------------------------------- + +-------------------------------------------------------------------------------- +-- $setup +-- >>> :{ +-- let myPolyLine = fromPointsUnsafe $ map ext [origin, Point2 10.0 10.0, Point2 10.0 20.0] +-- :} + +-------------------------------------------------------------------------------- -- * d-dimensional Polygonal Lines (PolyLines) -- | A Poly line in R^d has at least 2 vertices @@ -50,25 +60,33 @@ pmap f = over points (fmap (first f)) instance Arity d => Bifunctor (PolyLine d) where - bimap f g (PolyLine pts) = PolyLine $ fmap (bimap (fmap g) f) pts + bimap = bimapDefault +instance Arity d => Bifoldable (PolyLine d) where + bifoldMap = bifoldMapDefault +instance Arity d => Bitraversable (PolyLine d) where + bitraverse f g (PolyLine pts) = PolyLine <$> traverse (bitraverse (traverse g) f) pts instance (ToJSON p, ToJSON r, Arity d) => ToJSON (PolyLine d p r) where toEncoding = genericToEncoding defaultOptions instance (FromJSON p, FromJSON r, Arity d, KnownNat d) => FromJSON (PolyLine d p r) +-- | Builds a Polyline from a list of points, if there are sufficiently many points +fromPoints :: [Point d r :+ p] -> Maybe (PolyLine d p r) +fromPoints = fmap PolyLine . LSeq.eval (C @ 2) . LSeq.fromList + -- | pre: The input list contains at least two points -fromPoints :: [Point d r :+ p] -> PolyLine d p r -fromPoints = PolyLine . LSeq.forceLSeq (C :: C 2) . LSeq.fromList +fromPointsUnsafe :: [Point d r :+ p] -> PolyLine d p r +fromPointsUnsafe = PolyLine . LSeq.forceLSeq (C @ 2) . LSeq.fromList -- | pre: The input list contains at least two points. All extra vields are -- initialized with mempty. -fromPoints' :: (Monoid p) => [Point d r] -> PolyLine d p r -fromPoints' = fromPoints . map (\p -> p :+ mempty) +fromPointsUnsafe' :: (Monoid p) => [Point d r] -> PolyLine d p r +fromPointsUnsafe' = fromPointsUnsafe . map (\p -> p :+ mempty) -- | We consider the line-segment as closed. fromLineSegment :: LineSegment d p r -> PolyLine d p r -fromLineSegment ~(LineSegment' p q) = fromPoints [p,q] +fromLineSegment ~(LineSegment' p q) = fromPointsUnsafe [p,q] -- | Convert to a closed line segment by taking the first two points. asLineSegment :: PolyLine d p r -> LineSegment d p r @@ -80,3 +98,23 @@ asLineSegment' (PolyLine pts) = case F.toList pts of [p,q] -> Just $ ClosedLineSegment p q _ -> Nothing + +-- | Computes the edges, as linesegments, of an LSeq +edgeSegments :: Arity d => PolyLine d p r -> LSeq 1 (LineSegment d p r) +edgeSegments pl = let vs = pl^.points + in LSeq.zipWith ClosedLineSegment (LSeq.init vs) (LSeq.tail vs) + + +-- | Linearly interpolate the polyline with a value in the range +-- \([0,n-1]\), where \(n\) is the number of vertices of the polyline. +-- +-- running time: \(O(\log n)\) +-- +-- >>> interpolatePoly 0.5 myPolyLine +-- Point2 [5.0,5.0] +-- >>> interpolatePoly 1.5 myPolyLine +-- Point2 [10.0,15.0] +interpolatePoly :: (RealFrac r, Arity d) => r -> PolyLine d p r -> Point d r +interpolatePoly t pl = let i = floor t in case edgeSegments pl^?ix i of + Nothing -> pl^.points.to LSeq.last.core + Just e -> interpolate (t-fromIntegral i) e diff --git a/src/Data/Geometry/Polygon/Core.hs b/src/Data/Geometry/Polygon/Core.hs --- a/src/Data/Geometry/Polygon/Core.hs +++ b/src/Data/Geometry/Polygon/Core.hs @@ -73,12 +73,15 @@ import Data.Util import Data.Vinyl.CoRec (asA) +-- import Data.RealNumber.Rational + -------------------------------------------------------------------------------- {- $setup +>>> import Data.RealNumber.Rational >>> :{ -- import qualified Data.CircularSeq as C -let simplePoly :: SimplePolygon () Rational +let simplePoly :: SimplePolygon () (RealNumber 10) simplePoly = SimplePolygon . C.fromList . map ext $ [ Point2 0 0 , Point2 10 0 , Point2 10 10 @@ -268,17 +271,17 @@ -- -- -- >>> mapM_ print . polygonVertices $ withIncidentEdges simplePoly --- Point2 [0 % 1,0 % 1] :+ SP LineSegment (Closed (Point2 [1 % 1,11 % 1] :+ ())) (Closed (Point2 [0 % 1,0 % 1] :+ ())) LineSegment (Closed (Point2 [0 % 1,0 % 1] :+ ())) (Closed (Point2 [10 % 1,0 % 1] :+ ())) --- Point2 [10 % 1,0 % 1] :+ SP LineSegment (Closed (Point2 [0 % 1,0 % 1] :+ ())) (Closed (Point2 [10 % 1,0 % 1] :+ ())) LineSegment (Closed (Point2 [10 % 1,0 % 1] :+ ())) (Closed (Point2 [10 % 1,10 % 1] :+ ())) --- Point2 [10 % 1,10 % 1] :+ SP LineSegment (Closed (Point2 [10 % 1,0 % 1] :+ ())) (Closed (Point2 [10 % 1,10 % 1] :+ ())) LineSegment (Closed (Point2 [10 % 1,10 % 1] :+ ())) (Closed (Point2 [5 % 1,15 % 1] :+ ())) --- Point2 [5 % 1,15 % 1] :+ SP LineSegment (Closed (Point2 [10 % 1,10 % 1] :+ ())) (Closed (Point2 [5 % 1,15 % 1] :+ ())) LineSegment (Closed (Point2 [5 % 1,15 % 1] :+ ())) (Closed (Point2 [1 % 1,11 % 1] :+ ())) --- Point2 [1 % 1,11 % 1] :+ SP LineSegment (Closed (Point2 [5 % 1,15 % 1] :+ ())) (Closed (Point2 [1 % 1,11 % 1] :+ ())) LineSegment (Closed (Point2 [1 % 1,11 % 1] :+ ())) (Closed (Point2 [0 % 1,0 % 1] :+ ())) +-- Point2 [0,0] :+ V2 LineSegment (Closed (Point2 [1,11] :+ ())) (Closed (Point2 [0,0] :+ ())) LineSegment (Closed (Point2 [0,0] :+ ())) (Closed (Point2 [10,0] :+ ())) +-- Point2 [10,0] :+ V2 LineSegment (Closed (Point2 [0,0] :+ ())) (Closed (Point2 [10,0] :+ ())) LineSegment (Closed (Point2 [10,0] :+ ())) (Closed (Point2 [10,10] :+ ())) +-- Point2 [10,10] :+ V2 LineSegment (Closed (Point2 [10,0] :+ ())) (Closed (Point2 [10,10] :+ ())) LineSegment (Closed (Point2 [10,10] :+ ())) (Closed (Point2 [5,15] :+ ())) +-- Point2 [5,15] :+ V2 LineSegment (Closed (Point2 [10,10] :+ ())) (Closed (Point2 [5,15] :+ ())) LineSegment (Closed (Point2 [5,15] :+ ())) (Closed (Point2 [1,11] :+ ())) +-- Point2 [1,11] :+ V2 LineSegment (Closed (Point2 [5,15] :+ ())) (Closed (Point2 [1,11] :+ ())) LineSegment (Closed (Point2 [1,11] :+ ())) (Closed (Point2 [0,0] :+ ())) withIncidentEdges :: Polygon t p r -> Polygon t (Two (LineSegment 2 p r)) r withIncidentEdges (SimplePolygon vs) = SimplePolygon $ C.zip3LWith f (C.rotateL vs) vs (C.rotateR vs) where - f p c n = c&extra .~ SP (ClosedLineSegment p c) (ClosedLineSegment c n) + f p c n = c&extra .~ Two (ClosedLineSegment p c) (ClosedLineSegment c n) withIncidentEdges (MultiPolygon vs hs) = MultiPolygon vs' hs' where (SimplePolygon vs') = withIncidentEdges $ SimplePolygon vs @@ -565,7 +568,7 @@ -- will be numbered last, in the same order. -- -- >>> numberVertices simplePoly --- SimplePolygon (CSeq [Point2 [0 % 1,0 % 1] :+ SP 0 (),Point2 [10 % 1,0 % 1] :+ SP 1 (),Point2 [10 % 1,10 % 1] :+ SP 2 (),Point2 [5 % 1,15 % 1] :+ SP 3 (),Point2 [1 % 1,11 % 1] :+ SP 4 ()]) +-- SimplePolygon (CSeq [Point2 [0,0] :+ SP 0 (),Point2 [10,0] :+ SP 1 (),Point2 [10,10] :+ SP 2 (),Point2 [5,15] :+ SP 3 (),Point2 [1,11] :+ SP 4 ()]) numberVertices :: Polygon t p r -> Polygon t (SP Int p) r numberVertices = snd . bimapAccumL (\a p -> (a+1,SP a p)) (\a r -> (a,r)) 0 -- TODO: Make sure that this does not have the same issues as foldl vs foldl' diff --git a/src/Data/Geometry/PrioritySearchTree.hs b/src/Data/Geometry/PrioritySearchTree.hs --- a/src/Data/Geometry/PrioritySearchTree.hs +++ b/src/Data/Geometry/PrioritySearchTree.hs @@ -26,6 +26,8 @@ import Data.Geometry.Point import Data.List.NonEmpty (NonEmpty(..)) import qualified Data.List.NonEmpty as NonEmpty +import Data.Measured.Class () +import Data.Measured.Size import Data.Ord (comparing, Down(..)) import Data.Range import qualified Data.Set as Set diff --git a/src/Data/Geometry/QuadTree.hs b/src/Data/Geometry/QuadTree.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/QuadTree.hs @@ -0,0 +1,203 @@ +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE TypeApplications #-} +module Data.Geometry.QuadTree-- ( module Data.Geometry.QuadTree.Cell + -- , module Data.Geometry.QuadTree.Quadrants + -- , module Data.Geometry.QuadTree.Split + -- , QuadTree(..) + -- , leaves + -- , withCells + -- ) + where + + +import Control.Lens (makeLenses, (^.), (.~), (&), (^?!), ix, view) +import Data.Ext +import qualified Data.Foldable as F +import Data.Geometry.Box +import Data.Geometry.Point +import Data.Geometry.QuadTree.Cell +import Data.Geometry.QuadTree.Quadrants +import Data.Geometry.QuadTree.Split +import Data.Geometry.QuadTree.Tree (Tree(..)) +import qualified Data.Geometry.QuadTree.Tree as Tree +import Data.Geometry.Vector +import Data.Intersection +import Data.List.NonEmpty (NonEmpty(..)) +import Data.Tree.Util (TreeNode(..), levels) +import GHC.Generics (Generic) +-------------------------------------------------------------------------------- + +-- | QuadTree on the starting cell +data QuadTree v p r = QuadTree { _startingCell :: !(Cell r) + , _tree :: !(Tree v p) + } + deriving (Show,Eq,Generic,Functor,Foldable,Traversable) +makeLenses ''QuadTree + +-------------------------------------------------------------------------------- +-- * Functions operating on the QuadTree (in terms of the 'Tree' type) + +withCells :: (Fractional r, Ord r) => QuadTree v p r -> QuadTree (v :+ Cell r) (p :+ Cell r) r +withCells qt = qt&tree .~ withCellsTree qt + +withCellsTree :: (Fractional r, Ord r) + => QuadTree v p r -> Tree (v :+ Cell r) (p :+ Cell r) +withCellsTree (QuadTree c t) = Tree.withCells c t + +leaves :: QuadTree v p r -> NonEmpty p +leaves = Tree.leaves . view tree + +perLevel :: QuadTree v p r -> NonEmpty (NonEmpty (TreeNode v p)) +perLevel = levels . Tree.toRoseTree . view tree + + +-------------------------------------------------------------------------------- + +-- | Given a starting cell, a Tree builder, and some input required by +-- the builder, constructs a quadTree. +buildOn :: Cell r -> (Cell r -> i -> Tree v p) -> i -> QuadTree v p r +buildOn c0 builder = QuadTree c0 . builder c0 + +-- | The Equivalent of Tree.build for constructing a QuadTree +build :: (Fractional r, Ord r) => (Cell r -> i -> Split i v p) -> Cell r -> i -> QuadTree v p r +build f c = buildOn c (Tree.build f) + +-- | Build a QuadtTree from a set of points. +-- +-- pre: the points lie inside the initial given cell. +-- +-- running time: \(O(nh)\), where \(n\) is the number of points and +-- \(h\) is the height of the resulting quadTree. +fromPointsBox :: (Fractional r, Ord r) + => Cell r -> [Point 2 r :+ p] -> QuadTree () (Maybe (Point 2 r :+ p)) r +fromPointsBox c = buildOn c Tree.fromPoints + +fromPoints :: (RealFrac r, Ord r) + => NonEmpty (Point 2 r :+ p) -> QuadTree () (Maybe (Point 2 r :+ p)) r +fromPoints pts = buildOn c Tree.fromPoints (F.toList pts) + where + c = fitsRectangle $ boundingBoxList (view core <$> pts) + +-- | Locates the cell containing the given point, if it exists. +-- +-- running time: \(O(h)\), where \(h\) is the height of the quadTree +findLeaf :: (Fractional r, Ord r) + => Point 2 r -> QuadTree v p r -> Maybe (p :+ Cell r) +findLeaf q (QuadTree c0 t) | q `intersects` c0 = Just $ findLeaf' c0 t + | otherwise = Nothing + where + -- | + -- pre: p intersects c + findLeaf' c = \case + Leaf p -> p :+ c + Node _ qs -> let quad = quadrantOf q c + in findLeaf' ((splitCell c)^?!ix quad) (qs^?!ix quad) + +-------------------------------------------------------------------------------- + + +fromZeros :: (Fractional r, Ord r, Num a, Eq a, v ~ Quadrants Sign) + => Cell r -> (Point 2 r -> a) -> QuadTree v (Either v Sign) r +fromZeros = fromZerosWith (limitWidthTo (-1)) + + +fromZerosWith :: (Fractional r, Ord r, Eq a, Num a) + => Limiter r (Corners Sign) (Corners Sign) Sign + -> Cell r + -> (Point 2 r -> a) + -> QuadTree (Quadrants Sign) (Signs Sign) r +fromZerosWith limit c0 f = fromZerosWith' limit c0 (fromSignum f) + + +type Signs sign = Either (Corners sign) sign + + +fromZerosWith' :: (Eq sign, Fractional r, Ord r) + => Limiter r (Corners sign) (Corners sign) sign + -> Cell r + -> (Point 2 r -> sign) + -> QuadTree (Quadrants sign) (Signs sign) r +fromZerosWith' limit c0 f = build (limit $ shouldSplitZeros f) c0 (f <$> cellCorners c0) + + + +-- type Sign = Ordering + +-- pattern Negative :: Sign +-- pattern Negative = LT +-- pattern Zero :: Sign +-- pattern Zero = EQ +-- pattern Positive :: Sign +-- pattern Positive = GT +-- {-# COMPLETE Negative, Zero, Positive #-} + +-- fromOrdering :: Ordering -> Sign +-- fromOrdering = id + + +data Sign = Negative | Zero | Positive deriving (Show,Eq,Ord) + + + +-- | Interpret an ordering result as a Sign +fromOrdering :: Ordering -> Sign +fromOrdering = \case + LT -> Negative + EQ -> Zero + GT -> Positive + +fromSignum :: (Num a, Eq a) => (b -> a) -> b -> Sign +fromSignum f = \x -> case signum (f x) of + -1 -> Negative + 0 -> Zero + 1 -> Positive + _ -> error "absurd: fromSignum" + +-- | Splitter that determines if we should split a cell based on the +-- sign of the corners. +shouldSplitZeros :: forall r sign. (Fractional r, Eq sign) + => (Point 2 r -> sign) -- ^ The function we are evaluating + -> Splitter r + (Quadrants sign) -- the input are the signs of the corners + (Quadrants sign) -- at internal nodes we store signs of corners + sign +shouldSplitZeros f (Cell w' p) qs@(Quadrants nw ne se sw) | all sameSign qs = No ne + | otherwise = Yes qs qs' + where + m = fAt rr rr + n = fAt rr ww + e = fAt ww rr + s = fAt rr 0 + w = fAt 0 rr + + sameSign = (== ne) + + -- signs at the new corners + qs' = Quadrants (Quadrants nw n m w) + (Quadrants n ne e m) + (Quadrants m e se s) + (Quadrants w m s sw) + + r = w' - 1 + rr = pow r + ww = pow w' + + fAt x y = f $ p .+^ Vector2 x y + + +isZeroCell :: (Eq sign) => sign -- ^ the zero value + -> Either v sign -> Bool +isZeroCell z = \case + Left _ -> True -- if we kept splitting then we must have a sign transition + Right s -> s == z + +-------------------------------------------------------------------------------- + + + +-- | Constructs an empty/complete tree from the starting width +completeTree :: (Fractional r, Ord r) => Cell r -> QuadTree () () r +completeTree c0 = + build (\_ w -> if w == 0 then No () else Yes () (pure $ w - 1)) c0 (c0^.cellWidthIndex) + +-------------------------------------------------------------------------------- diff --git a/src/Data/Geometry/QuadTree/Cell.hs b/src/Data/Geometry/QuadTree/Cell.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/QuadTree/Cell.hs @@ -0,0 +1,141 @@ +{-# LANGUAGE TemplateHaskell #-} +module Data.Geometry.QuadTree.Cell where + +import Control.Lens (makeLenses, (^.),(&),(%~),ix, to) +import Data.Ext +import Data.Geometry.Box +import Data.Geometry.Directions +import Data.Geometry.LineSegment +import Data.Geometry.Point +import Data.Geometry.Properties +import Data.Geometry.QuadTree.Quadrants +import Data.Geometry.Vector + +-------------------------------------------------------------------------------- + +-- | side lengths will be 2^i for some integer i +type WidthIndex = Int + +-- | A Cell corresponding to a node in the QuadTree +data Cell r = Cell { _cellWidthIndex :: {-# UNPACK #-} !WidthIndex + , _lowerLeft :: !(Point 2 r) + } deriving (Show,Eq,Functor,Foldable,Traversable) +makeLenses ''Cell + +-- | Computes a cell that contains the given rectangle +fitsRectangle :: (RealFrac r, Ord r) => Rectangle p r -> Cell r +fitsRectangle b = Cell w ((b^.to minPoint.core) .-^ Vector2 1 1) + where + w = lg' . ceiling . (1+) . maximum . size $ b + + -- "approximate log" that over approximates by a factor of at most two. + lg' :: Integer -> WidthIndex + lg' n = go 1 + where + go i | floor (pow i) <= n = go (i+1) -- note that the floor does not really do anything + -- since i is integral and >= 1. + | otherwise = i + +type instance Dimension (Cell r) = 2 +type instance NumType (Cell r) = r + +type instance IntersectionOf (Point 2 r) (Cell r) = '[ NoIntersection, Point 2 r] + +instance (Ord r, Fractional r) => (Point 2 r) `IsIntersectableWith` (Cell r) where + nonEmptyIntersection = defaultNonEmptyIntersection + p `intersect` c = p `intersect` toBox c + +pow :: Fractional r => WidthIndex -> r +pow i = case i `compare` 0 of + LT -> 1 / (2 ^ (-1*i)) + EQ -> 1 + GT -> 2 ^ i + +cellWidth :: Fractional r => Cell r -> r +cellWidth (Cell w _) = pow w + +toBox :: Fractional r => Cell r -> Box 2 () r +toBox (Cell w p) = box (ext $ p) (ext $ p .+^ Vector2 (pow w) (pow w)) + +inCell :: (Fractional r, Ord r) => Point 2 r :+ p -> Cell r -> Bool +inCell (p :+ _) c = p `inBox` toBox c + +cellCorners :: Fractional r => Cell r -> Quadrants (Point 2 r) +cellCorners = fmap (^.core) . corners . toBox + +-- | Sides are open +cellSides :: Fractional r => Cell r -> Sides (LineSegment 2 () r) +cellSides = fmap (\(ClosedLineSegment p q) -> OpenLineSegment p q) . sides . toBox + +splitCell :: (Num r, Fractional r) => Cell r -> Quadrants (Cell r) +splitCell (Cell w p) = Quadrants (Cell r $ f 0 rr) + (Cell r $ f rr rr) + (Cell r $ f rr 0) + (Cell r p) + where + r = w - 1 + rr = pow r + f x y = p .+^ Vector2 x y + + +midPoint :: Fractional r => Cell r -> Point 2 r +midPoint (Cell w p) = let rr = pow (w - 1) in p .+^ Vector2 rr rr + + +-------------------------------------------------------------------------------- + +-- | Partitions the points into quadrants. See 'quadrantOf' for the +-- precise rules. +partitionPoints :: (Fractional r, Ord r) + => Cell r -> [Point 2 r :+ p] -> Quadrants [Point 2 r :+ p] +partitionPoints c = foldMap (\p -> let q = quadrantOf (p^.core) c in mempty&ix q %~ (p:)) + +-- | Computes the quadrant of the cell corresponding to the current +-- point. Note that we decide the quadrant solely based on the +-- midpoint. If the query point lies outside the cell, it is still +-- assigned a quadrant. +-- +-- - The northEast quadrants includes its bottom and left side +-- - The southEast quadrant includes its left side +-- - The northWest quadrant includes its bottom side +-- - The southWest quadrants does not include any of its sides. +-- +-- +-- >>> quadrantOf (Point2 9 9) (Cell 4 origin) +-- NorthEast +-- >>> quadrantOf (Point2 8 9) (Cell 4 origin) +-- NorthEast +-- >>> quadrantOf (Point2 8 8) (Cell 4 origin) +-- NorthEast +-- >>> quadrantOf (Point2 8 7) (Cell 4 origin) +-- SouthEast +-- >>> quadrantOf (Point2 4 7) (Cell 4 origin) +-- SouthWest +-- >>> quadrantOf (Point2 4 10) (Cell 4 origin) +-- NorthWest +-- >>> quadrantOf (Point2 4 40) (Cell 4 origin) +-- NorthEast +-- >>> quadrantOf (Point2 4 40) (Cell 4 origin) +-- NorthWest +quadrantOf :: forall r. (Fractional r, Ord r) + => Point 2 r -> Cell r -> InterCardinalDirection +quadrantOf q c = let m = midPoint c + in case (q^.xCoord < m^.xCoord, q^.yCoord < m^.yCoord) of + (False,False) -> NorthEast + (False,True) -> SouthEast + (True,False) -> NorthWest + (True,True) -> SouthWest + + + +-- | Given two cells c and me, compute on which side of `me` the cell +-- `c` is. +-- +-- pre: c and me are non-overlapping +relationTo :: (Fractional r, Ord r) + => (p :+ Cell r) -> Cell r -> Sides (Maybe (p :+ Cell r)) +c `relationTo` me = f <$> Sides b l t r <*> cellSides me + where + Sides t r b l = cellSides (c^.extra) + f e e' | e `intersects` e' = Just c + | otherwise = Nothing diff --git a/src/Data/Geometry/QuadTree/Quadrants.hs b/src/Data/Geometry/QuadTree/Quadrants.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/QuadTree/Quadrants.hs @@ -0,0 +1,16 @@ +module Data.Geometry.QuadTree.Quadrants( pattern Quadrants + , Quadrants + , module Data.Geometry.Box.Corners + ) where + +import Data.Geometry.Box.Corners + +-------------------------------------------------------------------------------- + +type Quadrants = Corners + +pattern Quadrants :: a -> a -> a -> a -> Corners a +pattern Quadrants a b c d = Corners a b c d +{-# COMPLETE Quadrants #-} + +-------------------------------------------------------------------------------- diff --git a/src/Data/Geometry/QuadTree/Split.hs b/src/Data/Geometry/QuadTree/Split.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/QuadTree/Split.hs @@ -0,0 +1,32 @@ +{-# LANGUAGE TemplateHaskell #-} +module Data.Geometry.QuadTree.Split where + +import Control.Lens (makePrisms,(^.)) +import Data.Geometry.QuadTree.Cell +import Data.Geometry.QuadTree.Quadrants + +-------------------------------------------------------------------------------- + +-- | Data Type to Decide if we should continue splitting the current cell +data Split i v p = No !p | Yes !v (Quadrants i) deriving (Show,Eq,Ord) +makePrisms ''Split + +-- | A splitter is a function that determines weather or not we should the given cell +-- corresponding to the given input (i). +type Splitter r i v p = Cell r -> i -> Split i v p + +-- | Transformer that limits the depth of a splitter +type Limiter r i v p = Splitter r i v p + -> Splitter r i v (Either i p) + +-- | Split only when the Cell-width is at least wMin +limitWidthTo :: WidthIndex -- ^ smallest allowed width of a cell (i.e. width of a leaf) + -> Limiter r i v p +limitWidthTo wMin f = \c pts -> case f c pts of + No p -> No (Right p) + Yes v qs | wMin < c^.cellWidthIndex -> Yes v qs + | otherwise -> No (Left pts) + -- note that it is important that we still evaluate the function so + -- that we can distinguish at the last level i.e. between a regular + -- " we are done splitting (No (Right p))" and a "we are no longer + -- allowed to split further (No (Left p))" diff --git a/src/Data/Geometry/QuadTree/Tree.hs b/src/Data/Geometry/QuadTree/Tree.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/QuadTree/Tree.hs @@ -0,0 +1,116 @@ +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE TypeApplications #-} +module Data.Geometry.QuadTree.Tree where + + +import Control.Lens (makePrisms) +import Data.Bifoldable +import Data.Bifunctor +import Data.Bitraversable +import Data.Ext +import qualified Data.Foldable as F +import Data.Functor.Apply +import Data.Geometry.Point +import Data.Geometry.QuadTree.Cell +import Data.Geometry.QuadTree.Quadrants +import Data.Geometry.QuadTree.Split +import Data.List.NonEmpty (NonEmpty(..)) +import qualified Data.List.NonEmpty as NonEmpty +import Data.Semigroup.Foldable.Class +import Data.Semigroup.Traversable.Class +import qualified Data.Tree as RoseTree +import Data.Tree.Util (TreeNode(..)) + +-------------------------------------------------------------------------------- + +-- | Our cells use Rational numbers as their numeric type +-- type CellR = Cell (RealNumber 10) + +-- | The Actual Tree type representing a quadTree +data Tree v p = Leaf !p + | Node !v (Quadrants (Tree v p)) -- quadrants are stored lazily on purpose + deriving (Show,Eq) +makePrisms ''Tree + +instance Bifunctor Tree where + bimap = bimapDefault + +instance Bifoldable Tree where + bifoldMap = bifoldMapDefault + +instance Bitraversable Tree where + bitraverse f g = \case + Leaf p -> Leaf <$> g p + Node v qs -> Node <$> f v <*> traverse (bitraverse f g) qs + +instance Bifoldable1 Tree +instance Bitraversable1 Tree where + bitraverse1 f g = \case + Leaf p -> Leaf <$> g p + Node v qs -> Node <$> f v <.> traverse1 (bitraverse1 f g) qs + +-- | Fold on the Tree type. +foldTree :: (p -> b) -> (v -> Quadrants b -> b) -> Tree v p -> b +foldTree f g = go + where + go = \case + Leaf p -> f p + Node v qs -> g v (go <$> qs) + +-- | Produce a list of all leaves of a quad tree +leaves :: Tree v p -> NonEmpty p +leaves = NonEmpty.fromList . bifoldMap (const []) (:[]) + +-- | Converts into a RoseTree +toRoseTree :: Tree v p -> RoseTree.Tree (TreeNode v p) +toRoseTree = foldTree (\p -> RoseTree.Node (LeafNode p) []) + (\v qs -> RoseTree.Node (InternalNode v) (F.toList qs)) + +-- | Computes the height of the quadtree +height :: Tree v p -> Integer +height = foldTree (const 1) (\_ -> (1 +) . maximum) + + +-------------------------------------------------------------------------------- + +-------------------------------------------------------------------------------- +-- * Functions operating on the QuadTree (in temrs of the 'Tree' type) + +-- | Builds a QuadTree +build :: Fractional r + => Splitter r pts v p -> Cell r -> pts -> Tree v p +build shouldSplit = build' + where + build' cc pts = case shouldSplit cc pts of + No p -> Leaf p + Yes v qs -> Node v $ build' <$> splitCell cc <*> qs + +-- | Annotate the tree with its corresponing cells +withCells :: Fractional r => Cell r -> Tree v p -> Tree (v :+ Cell r) (p :+ Cell r) +withCells c0 = \case + Leaf p -> Leaf (p :+ c0) + Node v qs -> Node (v :+ c0) (withCells <$> splitCell c0 <*> qs) + + +-------------------------------------------------------------------------------- + + +-- | Build a QuadtTree from a set of points. +-- +-- pre: the points lie inside the initial given cell. +-- +-- running time: \(O(nh)\), where \(n\) is the number of points and +-- \(h\) is the height of the resulting quadTree. +fromPoints :: (Fractional r, Ord r) + => Cell r -> [Point 2 r :+ p] + -> Tree () (Maybe (Point 2 r :+ p)) +fromPoints = build fromPointsF + +-- | The function that can be used to build a quadTree 'fromPoints' +fromPointsF :: (Fractional r, Ord r) + => Splitter r [Point 2 r :+ p] () (Maybe (Point 2 r :+ p)) +fromPointsF c = \case + [] -> No Nothing + [p] -> No (Just p) + pts -> Yes () $ partitionPoints c pts + -- (\cell -> filter (`inCell` cell) pts) <$> splitCell c diff --git a/src/Data/Geometry/RangeTree.hs b/src/Data/Geometry/RangeTree.hs --- a/src/Data/Geometry/RangeTree.hs +++ b/src/Data/Geometry/RangeTree.hs @@ -2,7 +2,6 @@ module Data.Geometry.RangeTree where import Control.Lens hiding (element) -import Data.BinaryTree (Measured(..)) import Data.Ext import qualified Data.Foldable as F import Data.Geometry.Point @@ -11,6 +10,7 @@ import Data.Geometry.Vector import Data.List.NonEmpty (NonEmpty(..)) import qualified Data.List.NonEmpty as NonEmpty +import Data.Measured.Class import Data.Proxy import Data.Range import Data.Semigroup.Foldable diff --git a/src/Data/Geometry/RangeTree/Generic.hs b/src/Data/Geometry/RangeTree/Generic.hs --- a/src/Data/Geometry/RangeTree/Generic.hs +++ b/src/Data/Geometry/RangeTree/Generic.hs @@ -9,6 +9,8 @@ import Data.List.NonEmpty (NonEmpty(..)) import qualified Data.List.NonEmpty as NonEmpty import Data.Range +import Data.Measured.Class +import Data.Measured.Size import Data.Semigroup import Data.Semigroup.Foldable import qualified Data.Set as Set diff --git a/src/Data/Geometry/RangeTree/Measure.hs b/src/Data/Geometry/RangeTree/Measure.hs --- a/src/Data/Geometry/RangeTree/Measure.hs +++ b/src/Data/Geometry/RangeTree/Measure.hs @@ -1,6 +1,6 @@ module Data.Geometry.RangeTree.Measure where -import Data.BinaryTree(Measured(..)) +import Data.Measured.Class import Data.Functor.Product(Product(..)) import Data.Functor.Classes diff --git a/src/Data/Geometry/SegmentTree/Generic.hs b/src/Data/Geometry/SegmentTree/Generic.hs --- a/src/Data/Geometry/SegmentTree/Generic.hs +++ b/src/Data/Geometry/SegmentTree/Generic.hs @@ -34,6 +34,8 @@ import qualified Data.List as List import Data.List.NonEmpty (NonEmpty) import qualified Data.List.NonEmpty as NonEmpty +import Data.Measured.Class +import Data.Measured.Size import GHC.Generics (Generic) -------------------------------------------------------------------------------- @@ -73,7 +75,7 @@ -- AtomicRange -> OpenRange MinInfinity MaxInfinity -data BuildLeaf a = LeafSingleton a | LeafRange a a deriving (Show,Eq) +data BuildLeaf a = LeafSingleton !a | LeafRange !a !a deriving (Show,Eq) -- | Given a sorted list of endpoints, without duplicates, construct a segment tree -- @@ -118,7 +120,7 @@ -> NonEmpty (Interval p r) -> SegmentTree v r fromIntervals f is = foldr (insert . f) (createTree pts mempty) is where - endPoints (toRange -> Range' a b) = [a,b] + endPoints (asRange -> Range' a b) = [a,b] pts = nub' . NonEmpty.sort . NonEmpty.fromList . concatMap endPoints $ is nub' = fmap NonEmpty.head . NonEmpty.group1 @@ -186,7 +188,7 @@ => i -> SegmentTree v r -> SegmentTree v r insert i (SegmentTree t) = SegmentTree $ insertRoot t where - ri@(Range a b) = toRange i + ri@(Range a b) = asRange i insertRoot t' = maybe t' (flip insert' t') $ getRange t' insert' inR lf@(Leaf nd@(LeafData rr _)) @@ -209,7 +211,7 @@ => i -> SegmentTree v r -> SegmentTree v r delete i (SegmentTree t) = SegmentTree $ delete' t where - (Range _ b) = toRange i + (Range _ b) = asRange i delete' (Leaf ld) = Leaf $ ld&leafAssoc %~ deleteAssoc i delete' (Node l nd@(_splitPoint -> m) r) @@ -254,7 +256,7 @@ instance IntervalLike a => IntervalLike (I a) where - toRange = toRange . _unI + asRange = asRange . _unI fromIntervals' :: (Eq p, Ord r) 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 @@ -136,7 +136,7 @@ :& RNil where s' = (fixEndPoints sm)^.subRange - s'' = bimap (^.extra) id + s'' = asProperInterval . first (^.extra) $ s'&start.core .~ toOffset' (s'^.start.extra.core) l &end.core .~ toOffset' (s'^.end.extra.core) l 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,69 +1,37 @@ +{-# LANGUAGE Unsafe #-} {-# LANGUAGE UndecidableInstances #-} module Data.Geometry.Transformation where -import Control.Lens (lens,Lens',set) -import Unsafe.Coerce(unsafeCoerce) +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 qualified Data.Vector.Fixed as FV import GHC.TypeLits -import Linear.Matrix ((!*),(!*!)) -import qualified Linear.Matrix as Lin --------------------------------------------------------------------------------- --- * Matrices --- | a matrix of n rows, each of m columns, storing values of type r -newtype Matrix n m r = Matrix (Vector n (Vector m r)) - -deriving instance (Show r, Arity n, Arity m) => Show (Matrix n m r) -deriving instance (Eq r, Arity n, Arity m) => Eq (Matrix n m r) -deriving instance (Ord r, Arity n, Arity m) => Ord (Matrix n m r) -deriving instance (Arity n, Arity m) => Functor (Matrix n m) - -multM :: (Arity r, Arity c, Arity c', Num a) => Matrix r c a -> Matrix c c' a -> Matrix r c' a -(Matrix a) `multM` (Matrix b) = Matrix $ a !*! b - -mult :: (Arity m, Arity n, Num r) => Matrix n m r -> Vector m r -> Vector n r -(Matrix m) `mult` v = m !* v - - -class Invertible n r where - inverse' :: Matrix n n r -> Matrix n n r - -instance Fractional r => Invertible 2 r where - -- >>> inverse' $ Matrix $ Vector2 (Vector2 1 2) (Vector2 3 4.0) - -- Matrix Vector2 [Vector2 [-2.0,1.0],Vector2 [1.5,-0.5]] - inverse' (Matrix m) = Matrix . unsafeCoerce . Lin.inv22 . unsafeCoerce $ m - -instance Fractional r => Invertible 3 r where - -- >>> inverse' $ Matrix $ Vector3 (Vector3 1 2 4) (Vector3 4 2 2) (Vector3 1 1 1.0) - -- Matrix Vector3 [Vector3 [0.0,0.5,-1.0],Vector3 [-0.5,-0.75,3.5],Vector3 [0.5,0.25,-1.5]] - inverse' (Matrix m) = Matrix . unsafeCoerce . Lin.inv33 . unsafeCoerce $ m - -instance Fractional r => Invertible 4 r where - inverse' (Matrix m) = Matrix . unsafeCoerce . Lin.inv44 . unsafeCoerce $ m - -------------------------------------------------------------------------------- -- * Transformations -- | A type representing a Transformation for d dimensional objects newtype Transformation d r = Transformation { _transformationMatrix :: Matrix (d + 1) (d + 1) r } -transformationMatrix :: Lens' (Transformation d r) (Matrix (d + 1) (d + 1) r) -transformationMatrix = lens _transformationMatrix (const Transformation) +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 @@ -112,16 +80,17 @@ translation :: (Num r, Arity d, Arity (d + 1)) => Vector d r -> Transformation d r -translation v = Transformation . Matrix $ V.imap transRow (snoc v 1) +translation v = Transformation . Matrix $ imap transRow (snoc v 1) scaling :: (Num r, Arity d, Arity (d + 1)) => Vector d r -> Transformation d r -scaling v = Transformation . Matrix $ V.imap mkRow (snoc v 1) +scaling v = Transformation . Matrix $ imap mkRow (snoc v 1) uniformScaling :: (Num r, Arity d, Arity (d + 1)) => r -> Transformation d r uniformScaling = scaling . pure + -------------------------------------------------------------------------------- -- * Functions that execute transformations @@ -142,14 +111,6 @@ scaleUniformlyBy = transformBy . uniformScaling --------------------------------------------------------------------------------- --- * Helper functions to easily create matrices - --- | Creates a row with zeroes everywhere, except at position i, where the --- value is the supplied value. -mkRow :: forall d r. (Arity d, Num r) => Int -> r -> Vector d r -mkRow i x = set (FV.element i) x zero - -- | 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 @@ -171,3 +132,13 @@ (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) 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 @@ -3,8 +3,11 @@ {-# LANGUAGE UndecidableInstances #-} module Data.Geometry.Triangle where +import Control.DeepSeq import Control.Lens +import Data.Bifoldable import Data.Bifunctor +import Data.Bitraversable import Data.Either (partitionEithers) import Data.Ext import Data.Geometry.Ball (Disk, disk) @@ -19,28 +22,48 @@ import qualified Data.Geometry.Vector as V import qualified Data.List as List import Data.Maybe (mapMaybe) +import Data.Util import Data.Vinyl import Data.Vinyl.CoRec +import GHC.Generics (Generic) import GHC.TypeLits - -------------------------------------------------------------------------------- -- | Triangles in \(d\)-dimensional space. -data Triangle d p r = Triangle (Point d r :+ p) - (Point d r :+ p) - (Point d r :+ p) +data Triangle d p r = Triangle !(Point d r :+ p) + !(Point d r :+ p) + !(Point d r :+ p) + deriving (Generic) -deriving instance (Arity d, Show r, Show p) => Show (Triangle d p r) -deriving instance (Arity d, Read r, Read p) => Read (Triangle d p r) -deriving instance (Arity d, Eq r, Eq p) => Eq (Triangle d p r) +deriving instance (Arity d, Show r, Show p) => Show (Triangle d p r) +deriving instance (Arity d, Read r, Read p) => Read (Triangle d p r) +deriving instance (Arity d, Eq r, Eq p) => Eq (Triangle d p r) -instance Arity d => Functor (Triangle d p) where - fmap f (Triangle p q r) = let f' = first (fmap f) in Triangle (f' p) (f' q) (f' r) +instance (Arity d, NFData r, NFData p) => NFData (Triangle d p r) +instance Arity d => Bifunctor (Triangle d) where bimap = bimapDefault +instance Arity d => Bifoldable (Triangle d) where bifoldMap = bifoldMapDefault +instance Arity d => Bitraversable (Triangle d) where + bitraverse f g (Triangle p q r) = let tr = bitraverse (traverse g) f in + Triangle <$> tr p <*> tr q <*> tr r + +-- instance Arity d => Functor (Triangle d p) where +-- fmap f (Triangle p q r) = let f' = first (fmap f) in Triangle (f' p) (f' q) (f' r) + +instance Field1 (Triangle d p r) (Triangle d p r) (Point d r :+ p) (Point d r :+ p) where + _1 = lens (\(Triangle p _ _) -> p) (\(Triangle _ q r) p -> Triangle p q r) +instance Field2 (Triangle d p r) (Triangle d p r) (Point d r :+ p) (Point d r :+ p) where + _2 = lens (\(Triangle _ q _) -> q) (\(Triangle p _ r) q -> Triangle p q r) +instance Field3 (Triangle d p r) (Triangle d p r) (Point d r :+ p) (Point d r :+ p) where + _3 = lens (\(Triangle _ _ r) -> r) (\(Triangle p q _) r -> Triangle p q r) + type instance NumType (Triangle d p r) = r type instance Dimension (Triangle d p r) = d + +_TriangleThreePoints :: Iso' (Triangle d p r) (Three (Point d r :+ p)) +_TriangleThreePoints = iso (\(Triangle p q r) -> Three p q r) (\(Three p q r) -> Triangle p q r) instance PointFunctor (Triangle d p) where pmap f (Triangle p q r) = Triangle (p&core %~ f) (q&core %~ f) (r&core %~ f) diff --git a/src/Data/Geometry/Vector.hs b/src/Data/Geometry/Vector.hs --- a/src/Data/Geometry/Vector.hs +++ b/src/Data/Geometry/Vector.hs @@ -14,29 +14,29 @@ , module LV , C(..) , Affine(..) - , qdA, distanceA + , quadrance, qdA, distanceA , dot, norm, signorm , isScalarMultipleOf , scalarMultiple -- reexports , FV.replicate - , FV.imap , xComponent, yComponent, zComponent ) where import Control.Applicative (liftA2) -import Control.Lens(Lens') +import Control.Lens (Lens') +import Control.Monad.State import qualified Data.Foldable as F import Data.Geometry.Properties import Data.Geometry.Vector.VectorFamily import Data.Geometry.Vector.VectorFixed (C(..)) -import Data.Maybe import qualified Data.Vector.Fixed as FV import GHC.TypeLits import Linear.Affine (Affine(..), qdA, distanceA) -import Linear.Metric (dot,norm,signorm) -import Linear.Vector as LV -import Test.QuickCheck +import Linear.Metric (dot,norm,signorm,quadrance) +import Linear.Vector as LV hiding (E(..)) +import System.Random (Random(..)) +import Test.QuickCheck (Arbitrary(..),infiniteList) -------------------------------------------------------------------------------- @@ -46,13 +46,21 @@ instance (Arbitrary r, Arity d) => Arbitrary (Vector d r) where arbitrary = vectorFromListUnsafe <$> infiniteList +instance (Random r, Arity d) => Random (Vector d r) where + randomR (lows,highs) g0 = flip runState g0 $ + FV.zipWithM (\l h -> state $ randomR (l,h)) lows highs + random g0 = flip runState g0 $ FV.replicateM (state random) -- | 'isScalarmultipleof u v' test if v is a scalar multiple of u. -- -- >>> Vector2 1 1 `isScalarMultipleOf` Vector2 10 10 -- True +-- >>> Vector3 1 1 2 `isScalarMultipleOf` Vector3 10 10 20 +-- True -- >>> Vector2 1 1 `isScalarMultipleOf` Vector2 10 1 -- False +-- >>> Vector2 1 1 `isScalarMultipleOf` Vector2 (-1) (-1) +-- True -- >>> Vector2 1 1 `isScalarMultipleOf` Vector2 11.1 11.1 -- True -- >>> Vector2 1 1 `isScalarMultipleOf` Vector2 11.1 11.2 @@ -63,11 +71,20 @@ -- True -- >>> Vector2 2 1 `isScalarMultipleOf` Vector2 4 0 -- False +-- >>> Vector3 2 1 0 `isScalarMultipleOf` Vector3 4 0 5 +-- False +-- >>> Vector3 0 0 0 `isScalarMultipleOf` Vector3 4 0 5 +-- True isScalarMultipleOf :: (Eq r, Fractional r, Arity d) => Vector d r -> Vector d r -> Bool -u `isScalarMultipleOf` v = isJust $ scalarMultiple u v +u `isScalarMultipleOf` v = let d = u `dot` v + num = quadrance u * quadrance v + in num == 0 || 1 == d*d / num +-- u `isScalarMultipleOf` v = isJust $ scalarMultiple u v {-# SPECIALIZE isScalarMultipleOf :: (Eq r, Fractional r) => Vector 2 r -> Vector 2 r -> Bool #-} +{-# SPECIALIZE + isScalarMultipleOf :: (Eq r, Fractional r) => Vector 3 r -> Vector 3 r -> Bool #-} -- | scalarMultiple u v computes the scalar labmda s.t. v = lambda * u (if it exists) scalarMultiple :: (Eq r, Fractional r, Arity d) diff --git a/src/Data/Geometry/Vector/VectorFamily.hs b/src/Data/Geometry/Vector/VectorFamily.hs --- a/src/Data/Geometry/Vector/VectorFamily.hs +++ b/src/Data/Geometry/Vector/VectorFamily.hs @@ -39,6 +39,7 @@ import Text.ParserCombinators.ReadPrec (lift) import Text.Read (Read(..),readListPrecDefault, readPrec_to_P,minPrec) import Data.Proxy +import Data.Hashable -------------------------------------------------------------------------------- -- * d dimensional Vectors @@ -59,8 +60,11 @@ unV = lens _unV (const MKVector) {-# INLINE unV #-} -type Arity d = (ImplicitArity (Peano d), KnownNat d) +-- type Arity d = (ImplicitArity (Peano d), KnownNat d) +class (ImplicitArity (Peano d), KnownNat d) => Arity d +instance (ImplicitArity (Peano d), KnownNat d) => Arity d + deriving instance (Eq r, Arity d) => Eq (Vector d r) deriving instance (Ord r, Arity d) => Ord (Vector d r) @@ -69,6 +73,15 @@ deriving instance Arity d => Traversable (Vector d) deriving instance Arity d => Applicative (Vector d) + + +instance Arity d => FunctorWithIndex Int (Vector d) where + imap = V.imap +instance Arity d => FoldableWithIndex Int (Vector d) +instance Arity d => TraversableWithIndex Int (Vector d) where + itraverse = V.imapM + + deriving instance Arity d => Additive (Vector d) deriving instance Arity d => Metric (Vector d) instance Arity d => Affine (Vector d) where @@ -76,6 +89,8 @@ u .-. v = u ^-^ v p .+^ v = p ^+^ v +deriving instance (Arity d, Hashable r) => Hashable (Vector d r) + instance Arity d => Ixed (Vector d r) where ix = element' @@ -167,6 +182,9 @@ -------------------------------------------------------------------------------- -- * Snoccing and consindg + +cons :: (Arity d, Arity (d+1)) => r -> Vector d r -> Vector (d + 1) r +cons x = vectorFromListUnsafe . (x:) . F.toList -- | Add an element at the back of the vector snoc :: (Arity (d + 1), Arity d) => Vector d r -> r -> Vector (d + 1) r diff --git a/src/Data/Geometry/Vector/VectorFamilyPeano.hs b/src/Data/Geometry/Vector/VectorFamilyPeano.hs --- a/src/Data/Geometry/Vector/VectorFamilyPeano.hs +++ b/src/Data/Geometry/Vector/VectorFamilyPeano.hs @@ -5,7 +5,7 @@ import Control.Applicative (liftA2) import Control.DeepSeq import Control.Lens hiding (element) -import Data.Aeson(FromJSON(..),ToJSON(..)) +import Data.Aeson (FromJSON(..),ToJSON(..)) -- import Data.Aeson (ToJSON(..),FromJSON(..)) import qualified Data.Foldable as F import qualified Data.Geometry.Vector.VectorFixed as FV @@ -19,6 +19,7 @@ import qualified Linear.V3 as L3 import qualified Linear.V4 as L4 import Linear.Vector +import Data.Hashable -------------------------------------------------------------------------------- -- * Natural number stuff @@ -77,7 +78,9 @@ unVF = lens _unVF (const VectorFamily) {-# INLINE unVF #-} -type ImplicitArity d = (ImplicitPeano d, V.Arity (FromPeano d)) +-- type ImplicitArity d = (ImplicitPeano d, V.Arity (FromPeano d)) +class (ImplicitPeano d, V.Arity (FromPeano d)) => ImplicitArity d +instance (ImplicitPeano d, V.Arity (FromPeano d)) => ImplicitArity d instance (Eq r, ImplicitArity d) => Eq (VectorFamily d r) where (VectorFamily u) == (VectorFamily v) = case (implicitPeano :: SingPeano d) of @@ -187,6 +190,17 @@ (SS (SS (SS (SS SZ)))) -> rnf v (SS (SS (SS (SS (SS _))))) -> rnf v {-# INLINE rnf #-} + + +instance (ImplicitPeano d, Hashable r) => Hashable (VectorFamily d r) where + hashWithSalt = case (implicitPeano :: SingPeano d) of + SZ -> hashWithSalt + (SS SZ) -> hashWithSalt + (SS (SS SZ)) -> hashWithSalt + (SS (SS (SS SZ))) -> hashWithSalt + (SS (SS (SS (SS SZ)))) -> hashWithSalt + (SS (SS (SS (SS (SS _))))) -> hashWithSalt + instance ImplicitArity d => Ixed (VectorFamily d r) where ix = element' 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 @@ -13,7 +13,8 @@ -- embedding. -- -------------------------------------------------------------------------------- -module Data.PlaneGraph.Core( PlaneGraph(PlaneGraph), graph +module Data.PlaneGraph.Core( -- $setup + PlaneGraph(PlaneGraph), graph , PlanarGraph , VertexData(VertexData), vData, location, vtxDataToExt , fromSimplePolygon, fromConnectedSegments @@ -87,7 +88,6 @@ -------------------------------------------------------------------------------- - -- $setup -- >>> import Data.Proxy -- >>> import Data.PlaneGraph.AdjRep(Gr(Gr),Face(Face),Vtx(Vtx)) @@ -390,17 +390,26 @@ incidentEdges :: VertexId' s -> PlaneGraph s v e f r -> V.Vector (Dart s) incidentEdges v = PG.incidentEdges v . _graph --- | All incoming edges incident to vertex v, in counterclockwise order around v. + +-- | All edges incident to vertex v in incoming direction +-- (i.e. pointing into v) in counterclockwise order around v. -- +-- running time: \(O(k)\), where \(k) is the total number of incident edges of v +-- -- >>> incomingEdges (VertexId 1) smallG --- [Dart (Arc 1) -1] +-- [Dart (Arc 1) +1,Dart (Arc 4) -1,Dart (Arc 3) -1] incomingEdges :: VertexId' s -> PlaneGraph s v e f r -> V.Vector (Dart s) incomingEdges v = PG.incomingEdges v . _graph --- | All outgoing edges incident to vertex v, in counterclockwise order around v. + + +-- | All edges incident to vertex v in incoming direction +-- (i.e. pointing into v) in counterclockwise order around v. -- +-- running time: \(O(k)\), where \(k) is the total number of incident edges of v +-- -- >>> outgoingEdges (VertexId 1) smallG --- [Dart (Arc 4) +1,Dart (Arc 3) +1] +-- [Dart (Arc 1) -1,Dart (Arc 4) +1,Dart (Arc 3) +1] outgoingEdges :: VertexId' s -> PlaneGraph s v e f r -> V.Vector (Dart s) outgoingEdges v = PG.outgoingEdges v . _graph diff --git a/src/Graphics/Camera.hs b/src/Graphics/Camera.hs --- a/src/Graphics/Camera.hs +++ b/src/Graphics/Camera.hs @@ -22,9 +22,10 @@ ) where import Control.Lens +import Data.Geometry.Matrix import Data.Geometry.Point -import Data.Geometry.Vector import Data.Geometry.Transformation +import Data.Geometry.Vector -------------------------------------------------------------------------------- diff --git a/test/Data/Geometry/arrangement.ipe.out.ipe b/test/Data/Geometry/arrangement.ipe.out.ipe --- a/test/Data/Geometry/arrangement.ipe.out.ipe +++ b/test/Data/Geometry/arrangement.ipe.out.ipe @@ -133,7 +133,7 @@ -172.658536585365 637.215447154471 m +172.658536585365 637.215447154471 m 173.658536585365 638.048780487804 l 173.658536585365 638.048780487804 m 256.914285714285 707.428571428571 l