diff --git a/changelog.org b/changelog.org --- a/changelog.org +++ b/changelog.org @@ -2,7 +2,24 @@ * Changelog -** 0.9 (unreleased) +** 0.10 + +- renamed the smallest enclosing ball to RIC +- improved tangency finding on convex hulls/chains +- changes to how we order points in ccwCmpAround and cwCmpAround; + these will report EQ if points appear at the same angle from the + center point. +- new functions ccwCmpAroundWith and cwCmpAroundWith that allow you to + specify the direction corresponding to "zero". +- bugfixes, in particular triangulating a polygon with holes now works properly. +- removed some unused dependencies +- we are no longer depending on ghc-plugins; as a result hgeometry + now also compiles with ghcjs +- more ToJSON/FromJSON instances. +- removed the 'point2' and 'point3' functions in favor of the pattern + synonyms Point2 and Point3. + +** 0.9 - Implemented 2D Linear Programming using randomized incremental construction (in \(O(n)\) expected time). This allows us to solve diff --git a/docs/Data/PlaneGraph/small.png b/docs/Data/PlaneGraph/small.png new file mode 100644 Binary files /dev/null and b/docs/Data/PlaneGraph/small.png differ diff --git a/hgeometry.cabal b/hgeometry.cabal --- a/hgeometry.cabal +++ b/hgeometry.cabal @@ -1,5 +1,5 @@ name: hgeometry -version: 0.9.0.0 +version: 0.10.0.0 synopsis: Geometric Algorithms, Data structures, and Data types. description: HGeometry provides some basic geometry types, and geometric algorithms and @@ -44,7 +44,7 @@ extra-source-files: README.md changelog.org -Extra-doc-files: +Extra-doc-files: docs/Data/PlaneGraph/small.png -- docs/**/*.png cabal-version: 2.0 @@ -116,11 +116,13 @@ -- * Geometric Algorithms Algorithms.Geometry.ConvexHull.GrahamScan Algorithms.Geometry.ConvexHull.DivideAndConquer + Algorithms.Geometry.ConvexHull.QuickHull + -- Algorithms.Geometry.ConvexHull.JarvisMarch Algorithms.Geometry.LowerEnvelope.DualCH Algorithms.Geometry.SmallestEnclosingBall.Types - Algorithms.Geometry.SmallestEnclosingBall.RandomizedIncrementalConstruction + Algorithms.Geometry.SmallestEnclosingBall.RIC Algorithms.Geometry.SmallestEnclosingBall.Naive Algorithms.Geometry.DelaunayTriangulation.Types @@ -134,7 +136,7 @@ Algorithms.Geometry.WellSeparatedPairDecomposition.WSPD Algorithms.Geometry.WellSeparatedPairDecomposition.Types - Algorithms.Geometry.Diameter + Algorithms.Geometry.Diameter.Naive -- Algorithms.Geometry.Sweep @@ -160,6 +162,7 @@ Algorithms.Geometry.FrechetDistance.Discrete + -- * Embedded Planar Graphs Data.PlaneGraph Data.PlaneGraph.Core @@ -171,12 +174,14 @@ Graphics.Render other-modules: - + -- * Implementation Internals of Polygons + Data.Geometry.Polygon.Core + Data.Geometry.Polygon.Extremes -- other-extensions: build-depends: base >= 4.11 && < 5 - , hgeometry-combinatorial >= 0.9.0.0 + , hgeometry-combinatorial >= 0.10.0.0 , bifunctors >= 4.1 , bytestring >= 0.10 @@ -185,7 +190,7 @@ , lens >= 4.2 , semigroupoids >= 5 , semigroups >= 0.18 - , singletons >= 2.0 + -- , singletons >= 2.0 , linear >= 1.10 , fixed-vector >= 1.0 , vector-builder >= 0.3.7 @@ -196,7 +201,12 @@ , QuickCheck >= 2.5 , quickcheck-instances >= 0.3 , reflection >= 2.1 + , primitive >= 0.6.3.0 + -- , singleton-typelits >= 0.1.0.0 + -- , ghc-typelits-natnormalise >= 0.6 + -- , ghc-typelits-knownnat >= 0.6 + , vector >= 0.11 , data-clist >= 0.1.2.3 , text >= 1.1.1.0 @@ -204,12 +214,14 @@ , aeson >= 1.0 , yaml >= 0.8 - , mtl + , mtl >= 2.2 , template-haskell + , hspec, QuickCheck, quickcheck-instances - hs-source-dirs: src + hs-source-dirs: src test + default-language: Haskell2010 default-extensions: TypeFamilies @@ -234,7 +246,6 @@ , DeriveFoldable , DeriveTraversable , DeriveGeneric - , AutoDeriveTypeable , FlexibleInstances diff --git a/src/Algorithms/Geometry/ClosestPair/DivideAndConquer.hs b/src/Algorithms/Geometry/ClosestPair/DivideAndConquer.hs --- a/src/Algorithms/Geometry/ClosestPair/DivideAndConquer.hs +++ b/src/Algorithms/Geometry/ClosestPair/DivideAndConquer.hs @@ -14,15 +14,12 @@ , CP , CCP(..) , mergePairs - , mergeSortedBy - , mergeSortedListsBy ) where +import Algorithms.DivideAndConquer import Control.Lens -import Data.BinaryTree import Data.Ext -import qualified Data.Foldable as F import Data.Geometry.Point import Data.LSeq (LSeq) import qualified Data.LSeq as LSeq @@ -30,7 +27,7 @@ import Data.List.NonEmpty (NonEmpty(..)) import qualified Data.List.NonEmpty as NonEmpty import Data.Ord (comparing) -import Data.Semigroup.Foldable(foldMap1, toNonEmpty) +import Data.Semigroup.Foldable (toNonEmpty) import Data.UnBounded import Data.Util @@ -40,11 +37,11 @@ -- \(n\) points. -- -- running time: \(O(n)\) -closestPair :: ( Ord r, Num r) => LSeq 2 (Point 2 r :+ p) -> Two (Point 2 r :+ p) -closestPair = f . foldMap1 mkCCP . asBalancedBinLeafTree . toNonEmpty +closestPair :: (Ord r, Num r) => LSeq 2 (Point 2 r :+ p) -> Two (Point 2 r :+ p) +closestPair = f . divideAndConquer1 mkCCP . toNonEmpty . LSeq.unstableSortBy (comparing (^.core)) where - mkCCP (Elem p) = CCP (p :| []) Top + mkCCP p = CCP (p :| []) Top f = \case CCP _ (ValT (SP cp _)) -> cp CCP _ Top -> error "closestPair: absurd." @@ -68,24 +65,6 @@ cmp p q = comparing (^.core.yCoord) p q <> comparing (^.core.xCoord) p q --- | Given an ordering and two nonempty sequences ordered according to that --- ordering, merge them -mergeSortedBy :: (a -> a -> Ordering) -> NonEmpty a -> NonEmpty a -> NonEmpty a -mergeSortedBy cmp ls rs = NonEmpty.fromList - $ mergeSortedListsBy cmp (F.toList ls) (F.toList rs) - - --- | Given an ordering and two nonempty sequences ordered according to that --- ordering, merge them -mergeSortedListsBy :: (a -> a -> Ordering) -> [a] -> [a] -> [a] -mergeSortedListsBy cmp = go - where - go [] ys = ys - go xs [] = xs - go xs@(x:xs') ys@(y:ys') = case x `cmp` y of - LT -> x : go xs' ys - EQ -> x : go xs' ys - GT -> y : go xs ys' -- | Function that does the actual merging work mergePairs :: forall p r. (Ord r, Num r) diff --git a/src/Algorithms/Geometry/ClosestPair/Naive.hs b/src/Algorithms/Geometry/ClosestPair/Naive.hs --- a/src/Algorithms/Geometry/ClosestPair/Naive.hs +++ b/src/Algorithms/Geometry/ClosestPair/Naive.hs @@ -9,12 +9,16 @@ -- \(n\) points in \(\mathbb{R}^d\). -- -------------------------------------------------------------------------------- -module Algorithms.Geometry.ClosestPair.Naive where +module Algorithms.Geometry.ClosestPair.Naive( closestPair + , closestPairWith + , DistanceFunction + ) where +import Control.Lens ((^.),_1) import Data.Ext import qualified Data.Foldable as F -import Data.Geometry (qdA) import Data.Geometry.Point +import Data.Geometry.Properties (NumType) import Data.Geometry.Vector (Arity) import Data.LSeq (LSeq) import qualified Data.List.NonEmpty as NonEmpty @@ -23,24 +27,31 @@ -------------------------------------------------------------------------------- --- | Naive algorithm to compute the closest pair in \(d\) dimensions. Runs in --- \(O(n^2)\) time (for any constant \(d\)). Note that we need at least two elements --- for there to be a closest pair. -closestPair :: ( Ord r, Arity d, Num r - ) => LSeq 2 (Point d r :+ p) -> Two (Point d r :+ p) -closestPair = getVal . getMin . sconcat . fmap (uncurry' mkPair) . pairs - where - uncurry' f (Two a b) = f a b - getVal (Arg _ x) = x +-- | Naive algorithm to compute the closest pair according to the +-- (squared) Euclidean distance in \(d\) dimensions. Note that we need +-- at least two elements for there to be a closest pair. +-- +-- running time: \(O(dn^2)\) time. +closestPair :: ( Ord r, Arity d, Num r) + => LSeq 2 (Point d r :+ p) -> Two (Point d r :+ p) +closestPair = (^._1) . closestPairWith (\p q -> squaredEuclideanDist (p^.core) (q^.core)) --- | A pair of points -type PP d p r = ArgMin r (Two (Point d r :+ p)) --- | Create a pair of points -mkPair :: (Arity d, Num r) - => Point d r :+ p -> Point d r :+ p -> PP d p r -mkPair pp@(p :+ _) qq@(q :+ _) = let dst = qdA p q - in Min (Arg dst (Two pp qq)) +type DistanceFunction g = g -> g -> NumType g + +-- | Naive algorithm to compute the closest pair of points (and the +-- distance realized by those points) given a distance function. Note +-- that we need at least two elements for there to be a closest pair. +-- +-- running time: \(O(T(d)n^2)\), where \(T(d)\) is the time required +-- to evaluate the distance between two points in \(\mathbb{R}^d\). +closestPairWith :: Ord r + => DistanceFunction (Point d r :+ p) + -> LSeq 2 (Point d r :+ p) -> SP (Two (Point d r :+ p)) r +closestPairWith d = getVal . getMin . sconcat . fmap mkPair . pairs + where + getVal (Arg dist x) = SP x dist + mkPair (Two p q) = Min (Arg (d p q) (Two p q)) -- | Produce all lists from a vec of elements. Since the Vec contains at least two -- elements, the resulting list is non-empty 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 @@ -9,30 +9,74 @@ -- of a set of \(n\) points in \(\mathbb{R}^2\). -- -------------------------------------------------------------------------------- -module Algorithms.Geometry.ConvexHull.DivideAndConquer(convexHull) where +module Algorithms.Geometry.ConvexHull.DivideAndConquer( convexHull + , upperHull + , lowerHull + ) where -import Control.Lens ((^.)) -import Data.BinaryTree +import Algorithms.DivideAndConquer +import Control.Arrow ((&&&)) +import Control.Lens ((^.), to) import Data.Ext -import Data.Function (on) import Data.Geometry.Point import Data.Geometry.Polygon import Data.Geometry.Polygon.Convex +import Data.List.NonEmpty (NonEmpty(..)) import qualified Data.List.NonEmpty as NonEmpty -import Data.Semigroup.Foldable - +import Data.Util -------------------------------------------------------------------------------- -- | \(O(n \log n)\) time ConvexHull using divide and conquer. The resulting polygon is -- given in clockwise order. -convexHull :: (Ord r, Num r) - => NonEmpty.NonEmpty (Point 2 r :+ p) -> ConvexPolygon p r -convexHull = unMerge - . foldMap1 (Merge . ConvexPolygon . fromPoints . (:[]) . _unElem) - . asBalancedBinLeafTree - . NonEmpty.sortBy (compare `on` (^.core)) +convexHull :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> ConvexPolygon p r +convexHull (p :| []) = ConvexPolygon . fromPoints $ [p] +convexHull pts = combine . (upperHull' &&& lowerHull') . NonEmpty.sortBy incXdecY $ pts + where + combine (l:|uh,_:|lh) = ConvexPolygon . fromPoints $ l : uh <> reverse (init lh) -newtype Merge r p = Merge { unMerge :: ConvexPolygon p r } +---------------------------------------- +-- * Computing a lower hull -instance (Num r, Ord r) => Semigroup (Merge r p) where - (Merge lp) <> (Merge rp) = let (ch,_,_) = merge lp rp in Merge ch +-- | \(O(n \log n)\) time LowerHull using divide and conquer. The resulting Hull is +-- given from left to right, i.e. in counter clockwise order. +lowerHull :: (Ord r, Num r) + => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p) +lowerHull = lowerHull' . NonEmpty.sortBy incXdecY + +lowerHull' :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p) +lowerHull' = unLH . divideAndConquer1 (LH . (:|[])) + +newtype LH r p = LH { unLH :: NonEmpty (Point 2 r :+ p) } deriving (Eq,Show) + +instance (Num r, Ord r) => Semigroup (LH r p) where + (LH lh) <> (LH rh) = LH $ hull lowerTangent' lh rh + +---------------------------------------- +-- * Computing an upper hull + +-- | \(O(n \log n)\) time UpperHull using divide and conquer. The resulting Hull is +-- given from left to right, i.e. in clockwise order. +upperHull :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p) +upperHull = upperHull' . NonEmpty.sortBy incXdecY + +upperHull' :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p) +upperHull' = unUH . divideAndConquer1 (UH . (:|[])) + +newtype UH r p = UH { unUH :: NonEmpty (Point 2 r :+ p) } + +instance (Num r, Ord r) => Semigroup (UH r p) where + (UH lh) <> (UH rh) = UH $ hull upperTangent' lh rh + +---------------------------------------- + +-- | The function that does the actual merging part +hull :: (NonEmpty p -> NonEmpty p -> Two (p :+ [p])) + -> NonEmpty p -> NonEmpty p -> NonEmpty p +hull tangent lh rh = let Two (l :+ lh') (r :+ rh') = tangent (NonEmpty.reverse lh) rh + in NonEmpty.fromList $ (reverse lh') <> [l,r] <> rh' + +-------------------------------------------------------------------------------- + +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 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 @@ -22,10 +22,11 @@ lh = NonEmpty.tail . hull' $ reverse ps' in ConvexPolygon . fromPoints . reverse $ lh ++ uh +-- | Computes the upper hull. 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 = hull id - +upperHull = NonEmpty.reverse . hull id +-- | Computes the upper hull. The upper hull is given from left to right lowerHull :: (Ord r, Num r) => NonEmpty (Point 2 r :+ p) -> NonEmpty (Point 2 r :+ p) lowerHull = hull reverse diff --git a/src/Algorithms/Geometry/ConvexHull/QuickHull.hs b/src/Algorithms/Geometry/ConvexHull/QuickHull.hs new file mode 100644 --- /dev/null +++ b/src/Algorithms/Geometry/ConvexHull/QuickHull.hs @@ -0,0 +1,109 @@ +module Algorithms.Geometry.ConvexHull.QuickHull( convexHull ) where + +import Control.Lens ((^.),(&),(.~)) +import Data.Ext +import qualified Data.Foldable as F +import Data.Geometry.Line +import Data.Geometry.Point +import Data.Geometry.Polygon +import Data.Geometry.Polygon.Convex (ConvexPolygon(..)) +import Data.Geometry.Triangle +import qualified Data.List as List +import Data.List.NonEmpty (NonEmpty(..)) +import Data.Ord (comparing) +import Data.Util + + +-- import Data.Ratio +-- import qualified Data.List.NonEmpty as NonEmpty +-- import qualified Algorithms.Geometry.ConvexHull.GrahamScan as GC +-- import Debug.Trace +-------------------------------------------------------------------------------- + +-- | ConvexHull using Quickhull. The resulting polygon is given in +-- clockwise order. +-- +-- running time: \(O(n^2)\) +convexHull :: (Ord r, Fractional r, Show r, Show p) + => NonEmpty (Point 2 r :+ p) -> ConvexPolygon p r +convexHull (p :| []) = ConvexPolygon . fromPoints $ [p] +convexHull ps = ConvexPolygon . fromPoints + $ [l] <> hull l r above <> [r] <> (reverse $ hull l r below) + where + STR l r mids = findExtremes ps + m = lineThrough (l^.core) (r^.core) + (above,below) = List.partition (\(p :+ _) -> p `liesAbove` m) mids + +-- | Finds the leftmost and rightmost point in the list +findExtremes :: Ord r + => NonEmpty (Point 2 r :+ q) + -> STR (Point 2 r :+ q) (Point 2 r :+ q) [Point 2 r :+ q] +findExtremes (p :| pts ) = foldr f (STR p p []) pts + where + f q (STR l r ms) = case (incXdecY q l, incXdecY q r) of + (LT,_) -> STR q r (addIfNot r l ms) + (EQ,_) -> STR l r ms -- ditch q; it is the same as l + (GT,GT) -> STR l q (addIfNot l r ms) + (GT,EQ) -> STR l r ms -- ditch q; it is the same as r + (GT,LT) -> STR l r (q:ms) + + addIfNot y x xs | (x^.core) /= (y^.core) = x:xs + | otherwise = xs + +-- findExtremesBy :: (a -> a -> Ordering) +-- -> NonEmpty a +-- -> STR a a [a] +-- findExtremesBy cmp pts = let l = F.minimumBy cmp pts +-- r = F.maximumBy cmp pts +-- a /=. b = a `cmp` b /= EQ +-- in STR l r [p | p <- F.toList pts, p /=. l, p /=. r] + + +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 + +-- | include neigher left or right +-- +hull :: (Fractional r, Ord r) + => Point 2 r :+ p -> Point 2 r :+ p -> [Point 2 r :+ p] -> [Point 2 r :+ p] +hull _ _ [] = [] +hull l r pts = hull l mid ls <> [mid] <> hull mid r rs + where + m = lineThrough (l^.core) (r^.core) + mid = F.maximumBy (comparing dist) pts + + dist (p :+ _) = p `sqDistanceTo` m + t = Triangle l mid r + -- line through l and mid, which splits the remaining points in a left half and a right half. + splitL = lineThrough (l^.core) (mid^.core) + rightSide = (r^.core) `onSide` splitL -- define the side containing r the right side + + (ls,rs) = List.partition (\(p :+ _) -> p `onSide` splitL /= rightSide) + . filter (\(p :+ _) -> not $ p `onTriangle` t) $ pts + + +-- mPoint2 [x,y] = Point2 x y + +-- -- testPoints = NonEmpty.fromList +-- -- [ mPoint2 [22536303956634 % 7570647828779,(-5816376064439) % 1228319866920] :+ 1 +-- -- , mPoint2 [(-3136920648983) % 824638230353,(-14583744643665) % 9604445576558] :+ 2 +-- -- , mPoint2 [(-11653462784667) % 6525086575987,(-598434515815) % 1364557986096] :+ 3 +-- -- , mPoint2 [(-7841595901661) % 3282967141364,(-207167076115) % 482378191549] :+ 4 +-- -- ] + + +-- testPoints :: NonEmpty (Point 2 Rational :+ Int) +-- testPoints = read "(Point2 [(-11199966464450) % 1365514034959,4065659138075 % 2296468530516] :+ 1) :| [Point2 [86686001553073 % 2736621704548,(-63774454571048) % 1880665081093] :+ 2,Point2 [(-77322324895231) % 8260610289790,(-41165682123514) % 2291705829063] :+ 3,Point2 [2292642905947 % 2659329735076,87045289726355 % 2752214350419] :+ 4]" + + +-- toDouble :: Point 2 Rational :+ a -> Point 2 Double :+ a +-- toDouble (p :+ x) = (realToFrac <$> p) :+ x + + +-- -- Falsified (after 36 tests and 4 shrinks): + + + +-- testPoints2 :: NonEmpty (Point 2 Rational :+ Int) +-- testPoints2 = read "(Point2 [(-9876468593885) % 9254762894818,(-34982972348716) % 7450362538495] :+ 1) :| [Point2 [(-11974177119403) % 7705693443554,(-37634868551543) % 9311528788922] :+ 2,Point2 [(-32383659855458) % 9565531378857,20253950785876 % 8268868939819] :+ 3,Point2 [42425655100996 % 8786996213535,(-7972873491283) % 1604043452399] :+ 4]" 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 @@ -162,6 +162,7 @@ True -> (,False) <$> rotateR' l r r1 (pred' r1) False -> pure (r1,True) + -- | The code that does the actual rotating rotateR' :: (Ord r, Fractional r) => VertexID -> VertexID -> Vertex -> Vertex -> Merge p r Vertex @@ -231,7 +232,8 @@ . IM.adjustWithKey (insert'' u) v where -- inserts b into the adjacency list of a - insert'' bi ai = CU.insertOrdBy (cwCmpAround (ptMap V.! ai) `on` (ptMap V.!)) bi + insert'' bi ai = CU.insertOrdBy (cwCmpAround' (ptMap V.! ai) `on` (ptMap V.!)) bi + cwCmpAround' c p q = cwCmpAround c p q <> cmpByDistanceTo c p q -- | Deletes an edge @@ -241,21 +243,21 @@ delete' x = CL.filterL (/= x) -- should we rotate left or right if it is the focus? - - -- | Lifted version of Convex.IsLeftOf isLeftOf :: (Ord r, Num r) => VertexID -> (VertexID, VertexID) -> Merge p r Bool p `isLeftOf` (l,r) = withPtMap . snd . fst <$> ask where - withPtMap ptMap = (ptMap V.! p) `Convex.isLeftOf` (ptMap V.! l, ptMap V.! r) + withPtMap ptMap = (ptMap V.! p) `isLeftOf'` (ptMap V.! l, ptMap V.! r) + a `isLeftOf'` (b,c) = ccw' b c a == CCW -- | Lifted version of Convex.IsRightOf isRightOf :: (Ord r, Num r) => VertexID -> (VertexID, VertexID) -> Merge p r Bool p `isRightOf` (l,r) = withPtMap . snd . fst <$> ask where - withPtMap ptMap = (ptMap V.! p) `Convex.isRightOf` (ptMap V.! l, ptMap V.! r) + withPtMap ptMap = (ptMap V.! p) `isRightOf'` (ptMap V.! l, ptMap V.! r) + a `isRightOf'` (b,c) = ccw' b c a == CW -------------------------------------------------------------------------------- -- * Some Helper functions diff --git a/src/Algorithms/Geometry/DelaunayTriangulation/Naive.hs b/src/Algorithms/Geometry/DelaunayTriangulation/Naive.hs --- a/src/Algorithms/Geometry/DelaunayTriangulation/Naive.hs +++ b/src/Algorithms/Geometry/DelaunayTriangulation/Naive.hs @@ -59,7 +59,7 @@ -- running time: O(m log m), where m=|vs| is the number of vertices to sort. sortAround' :: (Num r, Ord r) => Mapping p r -> VertexID -> [VertexID] -> [VertexID] -sortAround' (_,ptsV) u vs = reverse . map (^.extra) $ sortArround (f u) (map f vs) +sortAround' (_,ptsV) u vs = reverse . map (^.extra) $ sortAround (f u) (map f vs) where f v = (ptsV V.! v)&extra .~ v diff --git a/src/Algorithms/Geometry/Diameter.hs b/src/Algorithms/Geometry/Diameter.hs deleted file mode 100644 --- a/src/Algorithms/Geometry/Diameter.hs +++ /dev/null @@ -1,35 +0,0 @@ -module Algorithms.Geometry.Diameter where - -import Control.Lens -import Data.Ext -import Data.Geometry -import Data.List(maximumBy) - --------------------------------------------------------------------------------- - - -diameterNaive :: (Ord r, Floating r, Arity d) => [Point d r :+ p] -> r -diameterNaive = maybe 0 (\(p,q) -> euclideanDist (p^.core) (q^.core)) - . diametralPairNaive - --- | Computes the Euclidean diametral pair by naively trying all pairs. --- --- running time: \(O(n^2)\) -diametralPairNaive :: (Ord r, Num r, Arity d) - => [Point d r :+ p] -> Maybe (Point d r :+ p, Point d r :+ p) -diametralPairNaive = diametralPairWithNaive squaredEuclideanDist - - --- | Given a distance function and a list of points pts, computes the diametral --- pair by naively trying all pairs. --- --- running time: \(O(n^2)\) -diametralPairWithNaive :: Ord r - => (Point d r -> Point d r -> r) - -> [Point d r :+ p] - -> Maybe (Point d r :+ p, Point d r :+ p) -diametralPairWithNaive f pts@(_:_:_) = Just $ maximumBy cmp [ (p,q) | p <- pts, q <- pts ] - where - f' (p,q) = f (p^.core) (q^.core) - tp `cmp` tq = f' tp `compare` f' tq -diametralPairWithNaive _ _ = Nothing diff --git a/src/Algorithms/Geometry/Diameter/Naive.hs b/src/Algorithms/Geometry/Diameter/Naive.hs new file mode 100644 --- /dev/null +++ b/src/Algorithms/Geometry/Diameter/Naive.hs @@ -0,0 +1,32 @@ +module Algorithms.Geometry.Diameter.Naive where + +import Control.Lens +import Data.Ext +import Data.Geometry +import Data.List(maximumBy) + +-------------------------------------------------------------------------------- + +diameter :: (Ord r, Floating r, Arity d) => [Point d r :+ p] -> r +diameter = maybe 0 (\(p,q) -> euclideanDist (p^.core) (q^.core)) . diametralPair + +-- | Computes the Euclidean diametral pair by naively trying all pairs. +-- +-- running time: \(O(n^2)\) +diametralPair :: (Ord r, Num r, Arity d) + => [Point d r :+ p] -> Maybe (Point d r :+ p, Point d r :+ p) +diametralPair = diametralPairWith squaredEuclideanDist + +-- | Given a distance function and a list of points pts, computes the diametral +-- pair by naively trying all pairs. +-- +-- running time: \(O(n^2)\) +diametralPairWith :: Ord r + => (Point d r -> Point d r -> r) + -> [Point d r :+ p] + -> Maybe (Point d r :+ p, Point d r :+ p) +diametralPairWith f pts@(_:_:_) = Just $ maximumBy cmp [ (p,q) | p <- pts, q <- pts ] + where + f' (p,q) = f (p^.core) (q^.core) + tp `cmp` tq = f' tp `compare` f' tq +diametralPairWith _ _ = Nothing diff --git a/src/Algorithms/Geometry/LowerEnvelope/DualCH.hs b/src/Algorithms/Geometry/LowerEnvelope/DualCH.hs --- a/src/Algorithms/Geometry/LowerEnvelope/DualCH.hs +++ b/src/Algorithms/Geometry/LowerEnvelope/DualCH.hs @@ -16,11 +16,11 @@ type Envelope a r = NonEmpty (Line 2 r :+ a) -- | Given a list of non-vertical lines, computes the lower envelope using --- duality. +-- duality. The lines are given in left to right order. -- -- \(O(n\log n)\) lowerEnvelope :: (Ord r, Fractional r) => NonEmpty (Line 2 r :+ a) -> Envelope a r -lowerEnvelope = lowerEnvelopeWith upperHull +lowerEnvelope = NonEmpty.reverse . lowerEnvelopeWith upperHull type UpperHullAlgorithm a r = NonEmpty (Point 2 r :+ a) -> NonEmpty (Point 2 r :+ a) diff --git a/src/Algorithms/Geometry/PolygonTriangulation/MakeMonotone.hs b/src/Algorithms/Geometry/PolygonTriangulation/MakeMonotone.hs --- a/src/Algorithms/Geometry/PolygonTriangulation/MakeMonotone.hs +++ b/src/Algorithms/Geometry/PolygonTriangulation/MakeMonotone.hs @@ -1,7 +1,13 @@ {-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE ScopedTypeVariables #-} -module Algorithms.Geometry.PolygonTriangulation.MakeMonotone where +module Algorithms.Geometry.PolygonTriangulation.MakeMonotone( makeMonotone + , computeDiagonals + + , VertexType(..) + , classifyVertices + ) where + import Algorithms.Geometry.LineSegmentIntersection.BentleyOttmann ( xCoordAt , ordAt) import Algorithms.Geometry.PolygonTriangulation.Types @@ -29,16 +35,15 @@ import qualified Data.Vector.Mutable as MV +-- import Debug.Trace +-- import qualified Data.CircularSeq as CC ---------------------------------------------------------------------------------- data VertexType = Start | Merge | Split | End | Regular deriving (Show,Read,Eq) - --- How about the hole vertices? - -- | assigns a vertex type to each vertex -- --- pre: the polygon is given in CCW order +-- pre: Both the outer boundary and the inner boundary of the polygon are given in CCW order. -- -- running time: \(O(n)\). classifyVertices :: (Num r, Ord r) @@ -119,7 +124,7 @@ f = first (\i -> vertexInfo^.ix' i._2) pg :: Polygon t (SP Int (p :+ VertexType)) r - pg = numberVertices . classifyVertices . toCounterClockWiseOrder $ p' + pg = numberVertices . holesToCW . classifyVertices . toCCW $ p' vertexInfo :: V.Vector (STR (Point 2 r) p VertexType) vertexInfo = let vs = polygonVertices pg n = F.length vs @@ -137,9 +142,18 @@ sweep'' :: NonEmpty.NonEmpty (Event r) -> Sweep p r () sweep'' = mapM_ handle + -- make everything counterclockwise + toCCW p = (toCounterClockWiseOrder' p)&polygonHoles'.traverse %~ toCounterClockWiseOrder' + -- make the holes clockwise: + holesToCW p = p&polygonHoles'.traverse %~ toClockwiseOrder' + + + -- | Computes a set of diagionals that decompose the polygon into y-monotone -- pieces. -- +-- pre: the polygon boundary is given in counterClockwise order. +-- -- running time: \(O(n\log n)\) makeMonotone :: (Fractional r, Ord r) => proxy s -> Polygon t p r @@ -277,21 +291,21 @@ -- testPolygon :: SimplePolygon Int Rational --- testPolygon = fromPoints [ point2 20 20 :+ 1 --- , point2 18 19 :+ 2 --- , point2 16 25 :+ 3 --- , point2 13 23 :+ 4 --- , point2 10 24 :+ 5 --- , point2 6 22 :+ 6 --- , point2 8 21 :+ 7 --- , point2 7 18 :+ 8 --- , point2 2 19 :+ 9 --- , point2 1 10 :+ 10 --- , point2 3 5 :+ 11 --- , point2 11 7 :+ 12 --- , point2 15 1 :+ 13 --- , point2 12 15 :+ 14 --- , point2 15 12 :+ 15 +-- testPolygon = fromPoints [ Point2 20 20 :+ 1 +-- , Point2 18 19 :+ 2 +-- , Point2 16 25 :+ 3 +-- , Point2 13 23 :+ 4 +-- , Point2 10 24 :+ 5 +-- , Point2 6 22 :+ 6 +-- , Point2 8 21 :+ 7 +-- , Point2 7 18 :+ 8 +-- , Point2 2 19 :+ 9 +-- , Point2 1 10 :+ 10 +-- , Point2 3 5 :+ 11 +-- , Point2 11 7 :+ 12 +-- , Point2 15 1 :+ 13 +-- , Point2 12 15 :+ 14 +-- , Point2 15 12 :+ 15 -- ] -- vertexTypes = [Start,Merge,Start,Merge,Start,Regular,Regular,Merge,Start,Regular,End,Split,End,Split,End] @@ -307,3 +321,11 @@ -- ] -- outFile = "/Users/frank/tmp/out.ipe" -- writeIpeFile outFile . singlePageFromContent $ out + + +-- myPoly :: Polygon Multi () Rational +-- myPoly = MultiPolygon (CC.fromList $ read "[Point2 [16 % 1,80 % 1] :+ (),Point2 [16 % 1,16 % 1] :+ (),Point2 [144 % 1,16 % 1] :+ (),Point2 [144 % 1,80 % 1] :+ ()]" +-- ) +-- [ fromPoints $ read "[Point2 [88 % 1,48 % 1] :+ (),Point2 [112 % 1,40 % 1] :+ (),Point2 [112 % 1,48 % 1] :+ (),Point2 [80 % 1,56 % 1] :+ ()]" +-- , fromPoints $ read "[Point2 [32 % 1,64 % 1] :+ (),Point2 [32 % 1,32 % 1] :+ (),Point2 [64 % 1,32 % 1] :+ (),Point2 [64 % 1,64 % 1] :+ ()]" +-- ] diff --git a/src/Algorithms/Geometry/PolygonTriangulation/Triangulate.hs b/src/Algorithms/Geometry/PolygonTriangulation/Triangulate.hs --- a/src/Algorithms/Geometry/PolygonTriangulation/Triangulate.hs +++ b/src/Algorithms/Geometry/PolygonTriangulation/Triangulate.hs @@ -31,8 +31,8 @@ -- -- running time: \(O(n \log n)\) triangulate' :: (Ord r, Fractional r) - => proxy s -> Polygon t p r - -> PlaneGraph s p PolygonEdgeType PolygonFaceData r + => proxy s -> Polygon t p r + -> PlaneGraph s p PolygonEdgeType PolygonFaceData r triangulate' px pg' = constructGraph px e es diags where (pg, diags) = computeDiagonals' pg' @@ -56,16 +56,17 @@ computeDiagonals' pg' = (pg, monotoneDiags <> extraDiags) where pg = toCounterClockWiseOrder pg' - monotoneP = MM.makeMonotone (Identity pg') pg -- use some arbitrary proxy type + monotoneP = MM.makeMonotone (Identity pg) pg -- use some arbitrary proxy type -- outerFaceId' = outerFaceId monotoneP monotoneDiags = map (^._2.core) . filter (\e' -> e'^._2.extra == Diagonal) . F.toList . edgeSegments $ monotoneP extraDiags = concatMap (TM.computeDiagonals . toCounterClockWiseOrder') . lefts . map (^._2.core) + . filter (\mp -> mp^._2.extra == Inside) -- triangulate only the insides -- . filter (\f -> f^._1 /= outerFaceId') . F.toList . rawFacePolygons $ monotoneP - -- we alredy know we get the polgyons in *clockwise* order, so skip the - -- check if it is counter clockwise - toCounterClockWiseOrder' = reverseOuterBoundary + -- -- we alredy know we get the polgyons in *clockwise* order, so skip the + -- -- check if it is counter clockwise + -- toCounterClockWiseOrder'' = reverseOuterBoundary diff --git a/src/Algorithms/Geometry/PolygonTriangulation/TriangulateMonotone.hs b/src/Algorithms/Geometry/PolygonTriangulation/TriangulateMonotone.hs --- a/src/Algorithms/Geometry/PolygonTriangulation/TriangulateMonotone.hs +++ b/src/Algorithms/Geometry/PolygonTriangulation/TriangulateMonotone.hs @@ -145,11 +145,11 @@ -------------------------------------------------------------------------------- --- testPolygon = fromPoints . map ext $ [ point2 10 10 --- , point2 5 20 --- , point2 3 14 --- , point2 1 1 --- , point2 8 8 ] +-- testPolygon = fromPoints . map ext $ [ Point2 10 10 +-- , Point2 5 20 +-- , Point2 3 14 +-- , Point2 1 1 +-- , Point2 8 8 ] diff --git a/src/Algorithms/Geometry/PolygonTriangulation/Types.hs b/src/Algorithms/Geometry/PolygonTriangulation/Types.hs --- a/src/Algorithms/Geometry/PolygonTriangulation/Types.hs +++ b/src/Algorithms/Geometry/PolygonTriangulation/Types.hs @@ -12,6 +12,7 @@ import qualified Data.PlaneGraph as PG import qualified Data.Vector as V import qualified Data.Vector.Mutable as MV + -------------------------------------------------------------------------------- data PolygonEdgeType = Original | Diagonal @@ -83,7 +84,7 @@ constructGraph px e origs diags = subdiv & PG.vertexData.traverse %~ NonEmpty.head & PG.faceData .~ faceData' - & PG.rawDartData.traverse %~ snd + & PG.rawDartData.traverse %~ snd where subdiv :: PG.PlaneGraph s (NonEmpty p) (Bool,PolygonEdgeType) () r subdiv = PG.fromConnectedSegments px $ e' : origs' <> diags' @@ -101,6 +102,7 @@ -- the interior faces intFaces = flip PG.leftFace subdiv <$> queryDarts + faceData' :: V.Vector PolygonFaceData faceData' = V.create $ do v' <- MV.replicate (PG.numFaces subdiv) Outside 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 @@ -9,7 +9,9 @@ -- points in \(\mathbb{R}^2\) -- -------------------------------------------------------------------------------- -module Algorithms.Geometry.SmallestEnclosingBall.Naive where +module Algorithms.Geometry.SmallestEnclosingBall.Naive( smallestEnclosingDisk + , enclosesAll + ) where -- just for the types import Control.Lens @@ -55,7 +57,6 @@ => [Point 2 r :+ p] -> [DiskResult p r] -> DiskResult p r smallestEnclosingDisk' pts = minimumBy (compare `on` (^.enclosingDisk.squaredRadius)) . filter (flip enclosesAll pts) - -- | check if a disk encloses all points enclosesAll :: (Num r, Ord r) => DiskResult p r -> [Point 2 r :+ q] -> Bool diff --git a/src/Algorithms/Geometry/SmallestEnclosingBall/RIC.hs b/src/Algorithms/Geometry/SmallestEnclosingBall/RIC.hs new file mode 100644 --- /dev/null +++ b/src/Algorithms/Geometry/SmallestEnclosingBall/RIC.hs @@ -0,0 +1,164 @@ +{-# LANGUAGE DeriveFunctor #-} +{-# LANGUAGE TemplateHaskell #-} +-------------------------------------------------------------------------------- +-- | +-- Module : Algorithms.Geometry.SmallestEnclosingBall.RIC +-- Copyright : (C) Frank Staals +-- License : see the LICENSE file +-- Maintainer : Frank Staals +-- +-- An randomized algorithm to compute the smallest enclosing disk of a set of +-- \(n\) points in \(\mathbb{R}^2\). The expected running time is \(O(n)\). +-- +-------------------------------------------------------------------------------- +module Algorithms.Geometry.SmallestEnclosingBall.RIC( + smallestEnclosingDisk' + , smallestEnclosingDisk + , smallestEnclosingDiskWithPoint + , smallestEnclosingDiskWithPoints + ) where + +import Algorithms.Geometry.SmallestEnclosingBall.Types +import Control.Lens +import Control.Monad.Random.Class +import Data.Ext +import qualified Data.Foldable as F +import Data.Geometry +import Data.Geometry.Ball +import qualified Data.List as List +import Data.List.NonEmpty(NonEmpty(..)) +import Data.Maybe (fromMaybe, mapMaybe, catMaybes) +import Data.Ord (comparing) +import System.Random.Shuffle (shuffle) + +import Debug.Trace + +-------------------------------------------------------------------------------- + +-- | Compute the smallest enclosing disk of a set of points, +-- implemented using randomized incremental construction. +-- +-- pre: the input has at least two points. +-- +-- running time: expected \(O(n)\) time, where \(n\) is the number of input points. +smallestEnclosingDisk :: (Ord r, Fractional r, MonadRandom m + -- , Show r, Show p + ) + => [Point 2 r :+ p] + -> m (DiskResult p r) + +smallestEnclosingDisk pts@(_:_:_) = ((\(p:q:pts') -> smallestEnclosingDisk' p q pts') + . F.toList) <$> shuffle pts +smallestEnclosingDisk _ = error "smallestEnclosingDisk: Too few points" + +-- | Smallest enclosing disk. +smallestEnclosingDisk' :: (Ord r, Fractional r + -- , Show r, Show p + + ) + => Point 2 r :+ p -> Point 2 r :+ p -> [Point 2 r :+ p] + -> DiskResult p r +smallestEnclosingDisk' a b = foldr addPoint (initial a b) . List.tails + where + -- The empty case occurs only initially + addPoint [] br = br + addPoint (p:pts) br@(DiskResult d _) + | (p^.core) `inClosedBall` d = br + | otherwise = fromJust' $ smallestEnclosingDiskWithPoint p (a :| (b : pts)) + fromJust' = fromMaybe (error "smallestEncosingDisk' : fromJust, absurd") + +-- | Smallest enclosing disk, given that p should be on it. +smallestEnclosingDiskWithPoint :: (Ord r, Fractional r + -- , Show r, Show p + ) + => Point 2 r :+ p -> NonEmpty (Point 2 r :+ p) + -> Maybe (DiskResult p r) +smallestEnclosingDiskWithPoint p (a :| pts) = foldr addPoint (Just $ initial p a) $ List.tails pts + where + addPoint [] br = br + addPoint (q:pts') br@(Just (DiskResult d _)) + | (q^.core) `inClosedBall` d = br + | otherwise = smallestEnclosingDiskWithPoints p q (a:pts') + addPoint _ br = br + + +-- | Smallest enclosing disk, given that p and q should be on it +-- +-- running time: \(O(n)\) +smallestEnclosingDiskWithPoints :: (Ord r, Fractional r + -- , Show r, Show p + ) + => Point 2 r :+ p -> Point 2 r :+ p -> [Point 2 r :+ p] + -> Maybe (DiskResult p r) +smallestEnclosingDiskWithPoints p q ps = minimumOn (^.enclosingDisk.squaredRadius) + $ catMaybes [mkEnclosingDisk dl, mkEnclosingDisk dr, mdc] + where + centers = mapMaybe disk' ps + -- generate a disk with p q and r + disk' r = (r:+) <$> disk (p^.core) (q^.core) (r^.core) + + -- partition the points in to those on the left and those on the + -- right. Note that centers still contains only those points (and + -- disks) for which the three points are not colinear. So the + -- points are either on the left or on the right. + (leftCenters,rightCenters) = List.partition (\(r :+ _) -> ccw' p q r == CCW) centers + -- note that we consider 'leftmost' with respect to going from p + -- to q. This does not really have a global meaning. + + -- we need to find the leftmost and rightmost center on the + -- bisector. In case there are left-centers, this means that among + -- the left centers we want to find the point that is furthest way + -- from p (or q). If there are no left-centers, we with to find + -- the closest one among the right-centers. + leftDist z = let c = z^.extra.center + s = if ccw' p q c == CCW then 1 else -1 + in s * squaredEuclideanDist (p^.core) (c^.core) + + dl = maximumOn leftDist leftCenters -- disk that has the "leftmost" center + dr = minimumOn leftDist rightCenters -- disk that has the "rightmost" center + + -- diameteral disk + dd = fromDiameter (p^.core) (q^.core) + mdc | isEnclosingDisk dd ps = Just $ DiskResult dd (Two p q) + | otherwise = Nothing + + -- test if d is an enclosing disk. + mkEnclosingDisk md = md >>= mkEnclosingDisk' + mkEnclosingDisk' (r :+ d) | isEnclosingDisk d ps = Just (DiskResult d (Three p q r)) + | otherwise = Nothing + + +isEnclosingDisk :: (Foldable t, Ord r, Num r) + => Disk p r -> t (Point 2 r :+ extra) -> Bool +isEnclosingDisk d = all (\s -> (s^.core) `inClosedBall` d) + +-- | Constructs the initial 'DiskResult' from two points +initial :: Fractional r => Point 2 r :+ p -> Point 2 r :+ p -> DiskResult p r +initial p q = DiskResult (fromDiameter (p^.core) (q^.core)) (Two p q) + +maximumOn :: Ord b => (a -> b) -> [a] -> Maybe a +maximumOn f = \case + [] -> Nothing + xs -> Just $ List.maximumBy (comparing f) xs + +minimumOn :: Ord b => (a -> b) -> [a] -> Maybe a +minimumOn f = \case + [] -> Nothing + xs -> Just $ List.minimumBy (comparing f) xs + + +-------------------------------------------------------------------------------- + +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] + +disk'' r = (r:+) <$> disk (p^.core) (q^.core) (r^.core) + where + p = ext $ Point2 0 (-6) + q = ext $ Point2 0 6 diff --git a/src/Algorithms/Geometry/SmallestEnclosingBall/RandomizedIncrementalConstruction.hs b/src/Algorithms/Geometry/SmallestEnclosingBall/RandomizedIncrementalConstruction.hs deleted file mode 100644 --- a/src/Algorithms/Geometry/SmallestEnclosingBall/RandomizedIncrementalConstruction.hs +++ /dev/null @@ -1,84 +0,0 @@ -{-# LANGUAGE DeriveFunctor #-} -{-# LANGUAGE TemplateHaskell #-} --------------------------------------------------------------------------------- --- | --- Module : Algorithms.Geometry.SmallestEnclosingBall.RandomizedIncrementalConstruction --- Copyright : (C) Frank Staals --- License : see the LICENSE file --- Maintainer : Frank Staals --- --- An randomized algorithm to compute the smallest enclosing disk of a set of --- \(n\) points in \(\mathbb{R}^2\). The expected running time is \(O(n)\). --- --------------------------------------------------------------------------------- -module Algorithms.Geometry.SmallestEnclosingBall.RandomizedIncrementalConstruction where - -import Algorithms.Geometry.SmallestEnclosingBall.Types -import Control.Lens -import Control.Monad.Random.Class -import Data.Ext -import qualified Data.Foldable as F -import Data.Geometry -import Data.Geometry.Ball -import qualified Data.List as L -import Data.List.NonEmpty -import Data.Maybe (fromMaybe) -import System.Random.Shuffle (shuffle) - --------------------------------------------------------------------------------- - --- | O(n) expected time algorithm to compute the smallest enclosing disk of a --- set of points. we need at least two points. --- implemented using randomized incremental construction -smallestEnclosingDisk :: (Ord r, Fractional r, MonadRandom m) - => [Point 2 r :+ p] - -> m (DiskResult p r) - -smallestEnclosingDisk pts@(_:_:_) = ((\(p:q:pts') -> smallestEnclosingDisk' p q pts') - . F.toList) <$> shuffle pts -smallestEnclosingDisk _ = error "smallestEnclosingDisk: Too few points" - --- | Smallest enclosing disk. -smallestEnclosingDisk' :: (Ord r, Fractional r) - => Point 2 r :+ p -> Point 2 r :+ p -> [Point 2 r :+ p] - -> DiskResult p r -smallestEnclosingDisk' a b = foldr addPoint (initial a b) . L.tails - where - -- The epty case occurs only initially - addPoint [] br = br - addPoint (p:pts) br@(DiskResult d _) - | (p^.core) `inClosedBall` d = br - | otherwise = smallestEnclosingDiskWithPoint p (a :| (b : pts)) - - --- | Smallest enclosing disk, given that p should be on it. -smallestEnclosingDiskWithPoint :: (Ord r, Fractional r) - => Point 2 r :+ p -> NonEmpty (Point 2 r :+ p) - -> DiskResult p r -smallestEnclosingDiskWithPoint p (a :| pts) = foldr addPoint (initial p a) $ L.tails pts - where - addPoint [] br = br - addPoint (q:pts') br@(DiskResult d _) - | (q^.core) `inClosedBall` d = br - | otherwise = smallestEnclosingDiskWithPoints p q (a:pts') - - - --- | Smallest enclosing disk, given that p and q should be on it -smallestEnclosingDiskWithPoints :: (Ord r, Fractional r) - => Point 2 r :+ p -> Point 2 r :+ p -> [Point 2 r :+ p] - -> DiskResult p r -smallestEnclosingDiskWithPoints p q = foldr addPoint (initial p q) - where - addPoint r br@(DiskResult d _) - | (r^.core) `inClosedBall` d = br - | otherwise = DiskResult (circle' r) (Three p q r) - - circle' r = fromMaybe degen $ disk (p^.core) (q^.core) (r^.core) - degen = error "smallestEnclosingDisk: Unhandled degeneracy, three points on a line" - -- TODO: handle degenerate case - - --- | Constructs the initial 'DiskResult' from two points -initial :: Fractional r => Point 2 r :+ p -> Point 2 r :+ p -> DiskResult p r -initial p q = DiskResult (fromDiameter (p^.core) (q^.core)) (Two p q) diff --git a/src/Algorithms/Geometry/SmallestEnclosingBall/Types.hs b/src/Algorithms/Geometry/SmallestEnclosingBall/Types.hs --- a/src/Algorithms/Geometry/SmallestEnclosingBall/Types.hs +++ b/src/Algorithms/Geometry/SmallestEnclosingBall/Types.hs @@ -1,4 +1,3 @@ -{-# LANGUAGE DeriveFunctor #-} {-# LANGUAGE TemplateHaskell #-} -------------------------------------------------------------------------------- -- | @@ -41,5 +40,5 @@ -- and the points defining it data DiskResult p r = DiskResult { _enclosingDisk :: Disk () r , _definingPoints :: TwoOrThree (Point 2 r :+ p) - } + } deriving (Show,Eq) makeLenses ''DiskResult 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 @@ -209,7 +209,7 @@ -- | Assigns the points to a level. Returns the list of levels used. The first -- level in the list is the level assigned to the rest of the nodes. Their -- level is actually still set to Nothing in the underlying array. -assignLevels :: (Fractional r, Ord r, Arity d, KnownNat d +assignLevels :: (Fractional r, Ord r, Arity d , Show r, Show p ) => Int -- ^ Number of items we need to collect @@ -362,9 +362,6 @@ -------------------------------------------------------------------------------- -- * Finding Well Separated Pairs --- type AlwaysTrueWSPD d = ( Arity d, KnownNat d --- , AlwaysTruePFT d, AlwaysTrueTransformation d) - findPairs :: (Floating r, Ord r, Arity d, Arity (d + 1)) => r -> SplitTree d p r a -> SplitTree d p r a -> [WSP d p r a] @@ -436,7 +433,7 @@ -- | Computes the maximum width of a splitTree -maxWidth :: (Arity d, KnownNat d, Num r) +maxWidth :: (Arity d, Num r) => SplitTree d p r a -> r maxWidth (Leaf _) = 0 maxWidth (Node _ (NodeData i b _) _) = fromJust $ widthIn' i b 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 @@ -84,11 +84,11 @@ -- | Test if a point lies strictly inside a ball -- --- >>> (point2 0.5 0.0) `insideBall` unitBall +-- >>> (Point2 0.5 0.0) `insideBall` unitBall -- True --- >>> (point2 1 0) `insideBall` unitBall +-- >>> (Point2 1 0) `insideBall` unitBall -- False --- >>> (point2 2 0) `insideBall` unitBall +-- >>> (Point2 2 0) `insideBall` unitBall -- False insideBall :: (Arity d, Ord r, Num r) => Point d r -> Ball d p r -> Bool @@ -104,9 +104,9 @@ -- | Test if a point lies on the boundary of a ball. -- --- >>> (point2 1 0) `onBall` unitBall +-- >>> (Point2 1 0) `onBall` unitBall -- True --- >>> (point3 1 1 0) `onBall` unitBall +-- >>> (Point3 1 1 0) `onBall` unitBall -- False onBall :: (Arity d, Ord r, Num r) => Point d r -> Ball d p r -> Bool @@ -145,7 +145,7 @@ -- | Given three points, get the disk through the three points. If the three -- input points are colinear we return Nothing -- --- >>> disk (point2 0 10) (point2 10 0) (point2 (-10) 0) +-- >>> disk (Point2 0 10) (Point2 10 0) (Point2 (-10) 0) -- Just (Ball {_center = Point2 [0.0,0.0] :+ (), _squaredRadius = 100.0}) disk :: (Eq r, Fractional r) => Point 2 r -> Point 2 r -> Point 2 r -> Maybe (Disk () r) @@ -177,7 +177,7 @@ ynom = det33 $ V3 (fy px py) (fy qx qy) (fy sx sy) denom = (2 *) . det33 $ V3 (V3 px py 1) (V3 qx qy 1) (V3 sx sy 1) - c = point2 (xnom / denom) (ynom / denom) + c = Point2 (xnom / denom) (ynom / denom) 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 @@ -25,6 +25,7 @@ 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) @@ -171,9 +172,9 @@ -- | Check if a point lies a box -- --- >>> origin `inBox` (boundingBoxList' [point3 1 2 3, point3 10 20 30] :: Box 3 () Int) +-- >>> origin `inBox` (boundingBoxList' [Point3 1 2 3, Point3 10 20 30] :: Box 3 () Int) -- False --- >>> origin `inBox` (boundingBoxList' [point3 (-1) (-2) (-3), point3 10 20 30] :: Box 3 () Int) +-- >>> origin `inBox` (boundingBoxList' [Point3 (-1) (-2) (-3), Point3 10 20 30] :: Box 3 () Int) -- True inBox :: (Arity d, Ord r) => Point d r -> Box d p r -> Bool p `inBox` b = FV.and . FV.zipWith R.inRange (toVec p) . extent $ b @@ -182,7 +183,7 @@ -- resulting vector is 0 indexed whereas one would normally count dimensions -- starting at zero. -- --- >>> extent (boundingBoxList' [point3 1 2 3, point3 10 20 30] :: Box 3 () Int) +-- >>> extent (boundingBoxList' [Point3 1 2 3, Point3 10 20 30] :: Box 3 () Int) -- Vector3 [Range (Closed 1) (Closed 10),Range (Closed 2) (Closed 20),Range (Closed 3) (Closed 30)] extent :: Arity d => Box d p r -> Vector d (R.Range r) @@ -191,16 +192,16 @@ -- | Get the size of the box (in all dimensions). Note that the resulting vector is 0 indexed -- whereas one would normally count dimensions starting at zero. -- --- >>> size (boundingBoxList' [origin, point3 1 2 3] :: Box 3 () Int) +-- >>> size (boundingBoxList' [origin, Point3 1 2 3] :: Box 3 () Int) -- Vector3 [1,2,3] size :: (Arity d, Num r) => Box d p r -> Vector d r size = fmap R.width . extent -- | Given a dimension, get the width of the box in that dimension. Dimensions are 1 indexed. -- --- >>> widthIn (C :: C 1) (boundingBoxList' [origin, point3 1 2 3] :: Box 3 () Int) +-- >>> widthIn (C :: C 1) (boundingBoxList' [origin, Point3 1 2 3] :: Box 3 () Int) -- 1 --- >>> widthIn (C :: C 3) (boundingBoxList' [origin, point3 1 2 3] :: Box 3 () Int) +-- >>> widthIn (C :: C 3) (boundingBoxList' [origin, Point3 1 2 3] :: Box 3 () Int) -- 3 widthIn :: forall proxy p i d r. (Arity d, Arity (i - 1), Num r, ((i-1)+1) <= d) => proxy i -> Box d p r -> r @@ -209,13 +210,13 @@ -- | Same as 'widthIn' but with a runtime int instead of a static dimension. -- --- >>> widthIn' 1 (boundingBoxList' [origin, point3 1 2 3] :: Box 3 () Int) +-- >>> widthIn' 1 (boundingBoxList' [origin, Point3 1 2 3] :: Box 3 () Int) -- Just 1 --- >>> widthIn' 3 (boundingBoxList' [origin, point3 1 2 3] :: Box 3 () Int) +-- >>> widthIn' 3 (boundingBoxList' [origin, Point3 1 2 3] :: Box 3 () Int) -- Just 3 --- >>> widthIn' 10 (boundingBoxList' [origin, point3 1 2 3] :: Box 3 () Int) +-- >>> widthIn' 10 (boundingBoxList' [origin, Point3 1 2 3] :: Box 3 () Int) -- Nothing -widthIn' :: (Arity d, KnownNat d, Num r) => Int -> Box d p r -> Maybe r +widthIn' :: (Arity d, Num r) => Int -> Box d p r -> Maybe r widthIn' i = preview (V.element' (i-1)) . size @@ -224,14 +225,14 @@ type Rectangle = Box 2 --- >>> width (boundingBoxList' [origin, point2 1 2] :: Rectangle () Int) +-- >>> width (boundingBoxList' [origin, Point2 1 2] :: Rectangle () Int) -- 1 -- >>> width (boundingBoxList' [origin] :: Rectangle () Int) -- 0 width :: Num r => Rectangle p r -> r width = widthIn (C :: C 1) --- >>> height (boundingBoxList' [origin, point2 1 2] :: Rectangle () Int) +-- >>> height (boundingBoxList' [origin, Point2 1 2] :: Rectangle () Int) -- 2 -- >>> height (boundingBoxList' [origin] :: Rectangle () Int) -- 0 @@ -270,9 +271,9 @@ -- | Unsafe version of boundingBoxList, that does not check if the list is non-empty -boundingBoxList' :: (IsBoxable g, Ord (NumType g), Arity (Dimension g)) - => [g] -> Box (Dimension g) () (NumType g) -boundingBoxList' = boundingBoxList . NE.fromList +boundingBoxList' :: (IsBoxable g, Foldable c, Ord (NumType g), Arity (Dimension g)) + => c g -> Box (Dimension g) () (NumType g) +boundingBoxList' = boundingBoxList . NE.fromList . F.toList ---------------------------------------- diff --git a/src/Data/Geometry/Duality.hs b/src/Data/Geometry/Duality.hs --- a/src/Data/Geometry/Duality.hs +++ b/src/Data/Geometry/Duality.hs @@ -13,7 +13,7 @@ -- | Returns Nothing if the input line is vertical -- Maps a line l: y = ax + b to a point (a,-b) dualPoint :: (Fractional r, Eq r) => Line 2 r -> Maybe (Point 2 r) -dualPoint l = (\(a,b) -> point2 a (-b)) <$> toLinearFunction l +dualPoint l = (\(a,b) -> Point2 a (-b)) <$> toLinearFunction l -- | Pre: the input line is not vertical dualPoint' :: (Fractional r, Eq r) => Line 2 r -> Point 2 r 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 @@ -68,10 +68,10 @@ lineThrough p q = Line p (q .-. p) verticalLine :: Num r => r -> Line 2 r -verticalLine x = Line (point2 x 0) (Vector2 0 1) +verticalLine x = Line (Point2 x 0) (Vector2 0 1) horizontalLine :: Num r => r -> Line 2 r -horizontalLine y = Line (point2 0 y) (Vector2 1 0) +horizontalLine y = Line (Point2 0 y) (Vector2 1 0) -- | Given a line l with anchor point p and vector v, get the line -- perpendicular to l that also goes through p. The resulting line m is @@ -94,9 +94,9 @@ -- | Test if the two lines are parallel. -- --- >>> lineThrough origin (point2 1 0) `isParallelTo` lineThrough (point2 1 1) (point2 2 1) +-- >>> lineThrough origin (Point2 1 0) `isParallelTo` lineThrough (Point2 1 1) (Point2 2 1) -- True --- >>> lineThrough origin (point2 1 0) `isParallelTo` lineThrough (point2 1 1) (point2 2 2) +-- >>> lineThrough origin (Point2 1 0) `isParallelTo` lineThrough (Point2 1 1) (Point2 2 2) -- False isParallelTo :: (Eq r, Fractional r, Arity d) => Line d r -> Line d r -> Bool @@ -106,11 +106,11 @@ -- | Test if point p lies on line l -- --- >>> origin `onLine` lineThrough origin (point2 1 0) +-- >>> origin `onLine` lineThrough origin (Point2 1 0) -- True --- >>> point2 10 10 `onLine` lineThrough origin (point2 2 2) +-- >>> Point2 10 10 `onLine` lineThrough origin (Point2 2 2) -- True --- >>> point2 10 5 `onLine` lineThrough origin (point2 2 2) +-- >>> Point2 10 5 `onLine` lineThrough origin (Point2 2 2) -- False onLine :: (Eq r, Fractional r, Arity d) => Point d r -> Line d r -> Bool p `onLine` (Line q v) = p == q || (p .-. q) `isScalarMultipleOf` v @@ -203,7 +203,7 @@ -- | Create a line from the linear function ax + b fromLinearFunction :: Num r => r -> r -> Line 2 r -fromLinearFunction a b = Line (point2 0 b) (Vector2 1 a) +fromLinearFunction a b = Line (Point2 0 b) (Vector2 1 a) -- | get values a,b s.t. the input line is described by y = ax + b. -- returns Nothing if the line is vertical @@ -216,34 +216,56 @@ :& RNil -- | Result of a side test -data SideTest = Below | On | Above deriving (Show,Read,Eq,Ord) +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 `onSide` (lineThrough origin $ point2 10 5) +-- >>> Point2 10 10 `onSideUpDown` (lineThrough origin $ Point2 10 5) -- Above --- >>> point2 10 10 `onSide` (lineThrough origin $ point2 (-10) 5) +-- >>> Point2 10 10 `onSideUpDown` (lineThrough origin $ Point2 (-10) 5) -- Above --- >>> point2 5 5 `onSide` (verticalLine 10) +-- >>> Point2 5 5 `onSideUpDown` (verticalLine 10) -- Below --- >>> point2 5 5 `onSide` (lineThrough origin $ point2 (-3) (-3)) +-- >>> Point2 5 5 `onSideUpDown` (lineThrough origin $ Point2 (-3) (-3)) -- On -onSide :: (Ord r, Num r) => Point 2 r -> Line 2 r -> SideTest -q `onSide` (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 +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 +-- | Result of a side test +data SideTest = LeftSide | OnLine | RightSide 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 `onSide` (lineThrough origin $ Point2 10 5) +-- LeftSide +-- >>> Point2 10 10 `onSide` (lineThrough origin $ Point2 (-10) 5) +-- RightSide +-- >>> Point2 5 5 `onSide` (verticalLine 10) +-- LeftSide +-- >>> Point2 5 5 `onSide` (lineThrough origin $ Point2 (-3) (-3)) +-- OnLine +onSide :: (Ord r, Num r) => Point 2 r -> Line 2 r -> SideTest +q `onSide` (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 p r q of + CCW -> LeftSide + CW -> RightSide + CoLinear -> OnLine -- | Test if the query point q lies (strictly) above line l liesAbove :: (Ord r, Num r) => Point 2 r -> Line 2 r -> Bool -q `liesAbove` l = q `onSide` l == Above +q `liesAbove` l = q `onSideUpDown` l == Above -- | Get the bisector between two points 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 @@ -203,30 +203,30 @@ -- | Test if a point lies on a line segment. -- --- >>> (point2 1 0) `onSegment` (ClosedLineSegment (origin :+ ()) (point2 2 0 :+ ())) +-- >>> (Point2 1 0) `onSegment` (ClosedLineSegment (origin :+ ()) (Point2 2 0 :+ ())) -- True --- >>> (point2 1 1) `onSegment` (ClosedLineSegment (origin :+ ()) (point2 2 0 :+ ())) +-- >>> (Point2 1 1) `onSegment` (ClosedLineSegment (origin :+ ()) (Point2 2 0 :+ ())) -- False --- >>> (point2 5 0) `onSegment` (ClosedLineSegment (origin :+ ()) (point2 2 0 :+ ())) +-- >>> (Point2 5 0) `onSegment` (ClosedLineSegment (origin :+ ()) (Point2 2 0 :+ ())) -- False --- >>> (point2 (-1) 0) `onSegment` (ClosedLineSegment (origin :+ ()) (point2 2 0 :+ ())) +-- >>> (Point2 (-1) 0) `onSegment` (ClosedLineSegment (origin :+ ()) (Point2 2 0 :+ ())) -- False --- >>> (point2 1 1) `onSegment` (ClosedLineSegment (origin :+ ()) (point2 3 3 :+ ())) +-- >>> (Point2 1 1) `onSegment` (ClosedLineSegment (origin :+ ()) (Point2 3 3 :+ ())) -- True -- -- Note that the segments are assumed to be closed. So the end points lie on the segment. -- --- >>> (point2 2 0) `onSegment` (ClosedLineSegment (origin :+ ()) (point2 2 0 :+ ())) +-- >>> (Point2 2 0) `onSegment` (ClosedLineSegment (origin :+ ()) (Point2 2 0 :+ ())) -- True --- >>> origin `onSegment` (ClosedLineSegment (origin :+ ()) (point2 2 0 :+ ())) +-- >>> origin `onSegment` (ClosedLineSegment (origin :+ ()) (Point2 2 0 :+ ())) -- True -- -- -- This function works for arbitrary dimensons. -- --- >>> (point3 1 1 1) `onSegment` (ClosedLineSegment (origin :+ ()) (point3 3 3 3 :+ ())) +-- >>> (Point3 1 1 1) `onSegment` (ClosedLineSegment (origin :+ ()) (Point3 3 3 3 :+ ())) -- True --- >>> (point3 1 2 1) `onSegment` (ClosedLineSegment (origin :+ ()) (point3 3 3 3 :+ ())) +-- >>> (Point3 1 2 1) `onSegment` (ClosedLineSegment (origin :+ ()) (Point3 3 3 3 :+ ())) -- False onSegment :: (Ord r, Fractional r, Arity d) => Point d r -> LineSegment d p r -> Bool @@ -277,7 +277,7 @@ in (s&start .~ q)&end .~ p -- testSeg :: LineSegment 2 () Rational --- testSeg = LineSegment (Open $ ext origin) (Closed $ ext (point2 10 0)) +-- testSeg = LineSegment (Open $ ext origin) (Closed $ ext (Point2 10 0)) -- horL' :: Line 2 Rational -- horL' = horizontalLine 0 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 @@ -180,7 +180,7 @@ oF@(FaceId (VertexId of')) = PG.leftFace ofD g mkFaceIdx i | i == 0 = Nothing - | otherwise = Just (c,mkFaceId i) + | otherwise = Just (c,mkFaceId . flipID $ i) -- at index i we are storing the outerface mkFaceData :: Int -> f -> FaceData (Dart s) f diff --git a/src/Data/Geometry/PlanarSubdivision/Merge.hs b/src/Data/Geometry/PlanarSubdivision/Merge.hs --- a/src/Data/Geometry/PlanarSubdivision/Merge.hs +++ b/src/Data/Geometry/PlanarSubdivision/Merge.hs @@ -15,10 +15,13 @@ , embedAsHolesIn ) where +import Algorithms.DivideAndConquer import Control.Lens hiding (holes) -import Data.BinaryTree (asBalancedBinLeafTree, foldUp, Elem(..)) +import Data.Ext import Data.Geometry.PlanarSubdivision.Basic import Data.Geometry.PlanarSubdivision.Raw +import Data.Geometry.Point +import Data.Geometry.Polygon import Data.PlanarGraph.Dart import Data.PlaneGraph ( Dart, VertexId(..), FaceId(..) , VertexId', FaceId' @@ -26,14 +29,7 @@ import qualified Data.PlaneGraph as PG import Data.Semigroup.Foldable import qualified Data.Vector as V - --- import Data.Coerce -import Unsafe.Coerce(unsafeCoerce) - - -import Data.Ext -import Data.Geometry.Point -import Data.Geometry.Polygon +import Unsafe.Coerce (unsafeCoerce) -------------------------------------------------------------------------------- -- * Embedding one subdivision in another one @@ -118,7 +114,7 @@ => (f -> f -> f) -> t (PlanarSubdivision s v e f r) -> PlanarSubdivision s v e f r -mergeAllWith f = foldUp (\l _ r -> mergeWith f l r) _unElem . asBalancedBinLeafTree . toNonEmpty +mergeAllWith f = divideAndConquer1With (mergeWith f) id . toNonEmpty -- | Merge a pair of *disjoint* planar subdivisions, unifying their -- outer face. For the outerface data it simply takes the data of the @@ -153,7 +149,7 @@ -- we have to shift the number of the *Arcs*. Since every dart -- consists of two arcs, we have to shift by numDarts / 2 -- Furthermore, we take numFaces - 1 since we want the first - -- *internal* face of p2 (the one with FaceId 1) to correspond with the first free + -- /internal/ face of p2 (the one with FaceId 1) to correspond with the first free -- position (at index numFaces) cs = p1^.components <> p2'^.components 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 @@ -10,8 +10,32 @@ -- \(d\)-dimensional points. -- -------------------------------------------------------------------------------- -module Data.Geometry.Point where +module Data.Geometry.Point( Point(..) + , origin, vector + , pointFromList + , coord , unsafeCoord + + , projectPoint + + , pattern Point2 + , pattern Point3 + , xCoord, yCoord, zCoord + + , PointFunctor(..) + + , CCW(..), ccw, ccw' + + , ccwCmpAround, cwCmpAround, ccwCmpAroundWith, cwCmpAroundWith + , sortAround, insertIntoCyclicOrder + + , Quadrant(..), quadrantWith, quadrant, partitionIntoQuadrants + + , cmpByDistanceTo + + , squaredEuclideanDist, euclideanDist + ) where + import Control.DeepSeq import Control.Lens import Data.Aeson @@ -23,13 +47,14 @@ 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) -import Test.QuickCheck(Arbitrary) -------------------------------------------------------------------------------- @@ -99,7 +124,7 @@ -- | Lens to access the vector corresponding to this point. -- --- >>> (point3 1 2 3) ^. vector +-- >>> (Point3 1 2 3) ^. vector -- Vector3 [1,2,3] -- >>> origin & vector .~ Vector3 1 2 3 -- Point3 [1,2,3] @@ -111,7 +136,7 @@ -- sense that no bounds are checked. Consider using `coord` instead. -- -- --- >>> point3 1 2 3 ^. unsafeCoord 2 +-- >>> Point3 1 2 3 ^. unsafeCoord 2 -- 2 unsafeCoord :: Arity d => Int -> Lens' (Point d r) r unsafeCoord i = vector . singular (ix (i-1)) @@ -119,11 +144,11 @@ -- | Get the coordinate in a given dimension -- --- >>> point3 1 2 3 ^. coord (C :: C 2) +-- >>> Point3 1 2 3 ^. coord (C :: C 2) -- 2 --- >>> point3 1 2 3 & coord (C :: C 1) .~ 10 +-- >>> 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 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 @@ -160,15 +185,13 @@ -- let -- f :: Point 2 r -> r -- f (Point2 x y) = x --- in f (point2 1 2) +-- in f (Point2 1 2) -- :} -- 1 -- -- if we want. pattern Point2 :: r -> r -> Point 2 r -pattern Point2 x y <- (_point2 -> (x,y)) - where - Point2 x y = point2 x y +pattern Point2 x y = Point (Vector2 x y) {-# COMPLETE Point2 #-} -- | Similarly, we can write: @@ -181,47 +204,14 @@ -- :} -- 3 pattern Point3 :: r -> r -> r -> Point 3 r -pattern Point3 x y z <- (_point3 -> (x,y,z)) - where - Point3 x y z = point3 x y z +pattern Point3 x y z = (Point (Vector3 x y z)) {-# COMPLETE Point3 #-} --- | Construct a 2 dimensional point --- --- >>> point2 1 2 --- Point2 [1,2] -point2 :: r -> r -> Point 2 r -point2 x y = Point $ Vector2 x y - --- | Destruct a 2 dimensional point --- --- >>> _point2 $ point2 1 2 --- (1,2) -_point2 :: Point 2 r -> (r,r) -_point2 = (\(Vector2 x y) -> (x,y)) . toVec - - - --- | Construct a 3 dimensional point --- --- >>> point3 1 2 3 --- Point3 [1,2,3] -point3 :: r -> r -> r -> Point 3 r -point3 x y z = Point $ Vector3 x y z - --- | Destruct a 3 dimensional point --- --- >>> _point3 $ point3 1 2 3 --- (1,2,3) -_point3 :: Point 3 r -> (r,r,r) -_point3 = (\(Vector3 x y z) -> (x,y,z)) . toVec - - -- | Shorthand to access the first coordinate C 1 -- --- >>> point3 1 2 3 ^. xCoord +-- >>> Point3 1 2 3 ^. xCoord -- 1 --- >>> point2 1 2 & xCoord .~ 10 +-- >>> Point2 1 2 & xCoord .~ 10 -- Point2 [10,2] xCoord :: (1 <= d, Arity d) => Lens' (Point d r) r xCoord = coord (C :: C 1) @@ -229,9 +219,9 @@ -- | Shorthand to access the second coordinate C 2 -- --- >>> point2 1 2 ^. yCoord +-- >>> Point2 1 2 ^. yCoord -- 2 --- >>> point3 1 2 3 & yCoord %~ (+1) +-- >>> 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) @@ -239,9 +229,9 @@ -- | Shorthand to access the third coordinate C 3 -- --- >>> point3 1 2 3 ^. zCoord +-- >>> Point3 1 2 3 ^. zCoord -- 3 --- >>> point3 1 2 3 & zCoord %~ (+1) +-- >>> 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) @@ -287,9 +277,9 @@ -- 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) -sortArround :: (Ord r, Num r) - => Point 2 r :+ q -> [Point 2 r :+ p] -> [Point 2 r :+ p] -sortArround c = L.sortBy (ccwCmpAround c) +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 @@ -336,38 +326,85 @@ 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. --- Points nearer to the center come before --- points further away. -ccwCmpAround :: (Num r, Ord r) - => Point 2 r :+ qc -> Point 2 r :+ p -> Point 2 r :+ q -> Ordering -ccwCmpAround c q r = case (quadrantWith c q `compare` quadrantWith c r) of - EQ -> case ccw (c^.core) (q^.core) (r^.core) of - CCW -> LT - CW -> GT - CoLinear -> qdA (c^.core) (q^.core) - `compare` - qdA (c^.core) (r^.core) - x -> x -- if the quadrant differs, use the order - -- specified by the quadrant. +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. Points nearer to the center come before --- points further away. -cwCmpAround :: (Num r, Ord r) - => Point 2 r :+ qc -> Point 2 r :+ p -> Point 2 r :+ q -> Ordering -cwCmpAround c q r = case (quadrantWith c q `compare` quadrantWith c r) of - EQ -> case ccw (c^.core) (q^.core) (r^.core) of - CCW -> GT - CW -> LT - CoLinear -> qdA (c^.core) (q^.core) - `compare` - qdA (c^.core) (r^.core) - LT -> GT - GT -> LT -- if the quadrant differs, use the order - -- specified by the quadrant. - +-- 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 @@ -378,7 +415,7 @@ 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) +insertIntoCyclicOrder c = CU.insertOrdBy (ccwCmpAround c <> cmpByDistanceTo c) -- | Squared Euclidean distance between two points 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 @@ -4,6 +4,7 @@ module Data.Geometry.PolyLine where import Control.Lens +import Data.Aeson import Data.Bifunctor import Data.Ext import qualified Data.Foldable as F @@ -16,13 +17,14 @@ import Data.LSeq (LSeq, pattern (:<|)) import qualified Data.LSeq as LSeq import qualified Data.List.NonEmpty as NE +import GHC.Generics(Generic) import GHC.TypeLits -------------------------------------------------------------------------------- -- * d-dimensional Polygonal Lines (PolyLines) -- | A Poly line in R^d has at least 2 vertices -newtype PolyLine d p r = PolyLine { _points :: LSeq 2 (Point d r :+ p) } +newtype PolyLine d p r = PolyLine { _points :: LSeq 2 (Point d r :+ p) } deriving (Generic) makeLenses ''PolyLine deriving instance (Show r, Show p, Arity d) => Show (PolyLine d p r) @@ -50,6 +52,9 @@ instance Arity d => Bifunctor (PolyLine d) where bimap f g (PolyLine pts) = PolyLine $ fmap (bimap (fmap 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) -- | pre: The input list contains at least two points fromPoints :: [Point d r :+ p] -> PolyLine d p r diff --git a/src/Data/Geometry/Polygon.hs b/src/Data/Geometry/Polygon.hs --- a/src/Data/Geometry/Polygon.hs +++ b/src/Data/Geometry/Polygon.hs @@ -1,3 +1,4 @@ +{-# LANGUAGE TemplateHaskell #-} -------------------------------------------------------------------------------- -- | -- Module : Data.Geometry.Polygon @@ -8,514 +9,58 @@ -- A Polygon data type and some basic functions to interact with them. -- -------------------------------------------------------------------------------- -module Data.Geometry.Polygon where - -import Algorithms.Geometry.LinearProgramming.LP2DRIC -import Algorithms.Geometry.LinearProgramming.Types -import Control.DeepSeq -import Control.Lens hiding (Simple) -import Control.Monad.Random.Class -import Data.Bifoldable -import Data.Bifunctor -import Data.Bitraversable -import qualified Data.CircularSeq as C -import Data.Ext -import qualified Data.Foldable as F -import Data.Geometry.Boundary -import Data.Geometry.Box -import Data.Geometry.Line -import Data.Geometry.HalfSpace(rightOf) -import Data.Geometry.LineSegment -import Data.Geometry.Point -import Data.Geometry.Properties -import Data.Geometry.Transformation -import Data.Geometry.Triangle (Triangle(..), inTriangle) -import Data.Geometry.Vector -import qualified Data.List as List -import Data.List.NonEmpty (NonEmpty(..)) -import qualified Data.List.NonEmpty as NonEmpty -import Data.Maybe (mapMaybe, catMaybes) -import Data.Ord (comparing) -import Data.Semigroup (sconcat) -import Data.Semigroup.Foldable -import qualified Data.Sequence as Seq -import Data.Util -import Data.Vinyl.CoRec (asA) - - --------------------------------------------------------------------------------- --- * Polygons - -{- $setup ->>> :{ --- import qualified Data.CircularSeq as C -let simplePoly :: SimplePolygon () Rational - simplePoly = SimplePolygon . C.fromList . map ext $ [ point2 0 0 - , point2 10 0 - , point2 10 10 - , point2 5 15 - , point2 1 11 - ] -:} -} - --- | We distinguish between simple polygons (without holes) and Polygons with holes. -data PolygonType = Simple | Multi - - -data Polygon (t :: PolygonType) p r where - SimplePolygon :: C.CSeq (Point 2 r :+ p) -> Polygon Simple p r - MultiPolygon :: C.CSeq (Point 2 r :+ p) -> [Polygon Simple p r] -> Polygon Multi p r - -instance Bifunctor (Polygon t) where - bimap = bimapDefault - -instance Bifoldable (Polygon t) where - bifoldMap = bifoldMapDefault - -instance Bitraversable (Polygon t) where - bitraverse f g p = case p of - SimplePolygon vs -> SimplePolygon <$> bitraverseVertices f g vs - MultiPolygon vs hs -> MultiPolygon <$> bitraverseVertices f g vs - <*> traverse (bitraverse f g) hs - -instance (NFData p, NFData r) => NFData (Polygon t p r) where - rnf (SimplePolygon vs) = rnf vs - rnf (MultiPolygon vs hs) = rnf (vs,hs) - -bitraverseVertices :: (Applicative f, Traversable t) => (p -> f q) -> (r -> f s) - -> t (Point 2 r :+ p) -> f (t (Point 2 s :+ q)) -bitraverseVertices f g = traverse (bitraverse (traverse g) f) - -type SimplePolygon = Polygon Simple - -type MultiPolygon = Polygon Multi - --- | Either a simple or multipolygon -type SomePolygon p r = Either (Polygon Simple p r) (Polygon Multi p r) - -type instance Dimension (SomePolygon p r) = 2 -type instance NumType (SomePolygon p r) = r - --- | Polygons are per definition 2 dimensional -type instance Dimension (Polygon t p r) = 2 -type instance NumType (Polygon t p r) = r - -instance (Show p, Show r) => Show (Polygon t p r) where - show (SimplePolygon vs) = "SimplePolygon " <> show vs - show (MultiPolygon vs hs) = "MultiPolygon " <> show vs <> " " <> show hs - -instance (Eq p, Eq r) => Eq (Polygon t p r) where - (SimplePolygon vs) == (SimplePolygon vs') = vs == vs' - (MultiPolygon vs hs) == (MultiPolygon vs' hs') = vs == vs' && hs == hs' - -instance PointFunctor (Polygon t p) where - pmap f (SimplePolygon vs) = SimplePolygon (fmap (first f) vs) - pmap f (MultiPolygon vs hs) = MultiPolygon (fmap (first f) vs) (map (pmap f) hs) - -instance Fractional r => IsTransformable (Polygon t p r) where - transformBy = transformPointFunctor - -instance IsBoxable (Polygon t p r) where - boundingBox = boundingBoxList' . toListOf (outerBoundary.traverse.core) - -type instance IntersectionOf (Line 2 r) (Boundary (Polygon t p r)) = - '[Seq.Seq (Either (Point 2 r) (LineSegment 2 () r))] - -type instance IntersectionOf (Point 2 r) (Polygon t p r) = [NoIntersection, Point 2 r] - -instance (Fractional r, Ord r) => (Point 2 r) `IsIntersectableWith` (Polygon t p r) where - nonEmptyIntersection = defaultNonEmptyIntersection - q `intersects` pg = q `inPolygon` pg /= Outside - q `intersect` pg | q `intersects` pg = coRec q - | otherwise = coRec NoIntersection - --- instance IsIntersectableWith (Line 2 r) (Boundary (Polygon t p r)) where --- nonEmptyIntersection _ _ (CoRec xs) = null xs --- l `intersect` (Boundary (SimplePolygon vs)) = --- undefined - -- l `intersect` (Boundary (MultiPolygon vs hs)) = coRec . - -- Seq.sortBy f . Seq.fromList - -- . concatMap (unpack . (l `intersect`) . Boundary) - -- $ SimplePolygon vs : hs - -- where - -- unpack (CoRec x) = x - -- f = undefined - - - - --- * Functions on Polygons - -outerBoundary :: forall t p r. Lens' (Polygon t p r) (C.CSeq (Point 2 r :+ p)) -outerBoundary = lens g s - where - g :: Polygon t p r -> C.CSeq (Point 2 r :+ p) - g (SimplePolygon vs) = vs - g (MultiPolygon vs _) = vs - - s :: Polygon t p r -> C.CSeq (Point 2 r :+ p) - -> Polygon t p r - s (SimplePolygon _) vs = SimplePolygon vs - s (MultiPolygon _ hs) vs = MultiPolygon vs hs - -polygonHoles :: forall p r. Lens' (Polygon Multi p r) [Polygon Simple p r] -polygonHoles = lens g s - where - g :: Polygon Multi p r -> [Polygon Simple p r] - g (MultiPolygon _ hs) = hs - s :: Polygon Multi p r -> [Polygon Simple p r] - -> Polygon Multi p r - s (MultiPolygon vs _) = MultiPolygon vs - - --- | Access the i^th vertex on the outer boundary -outerVertex :: Int -> Lens' (Polygon t p r) (Point 2 r :+ p) -outerVertex i = outerBoundary.C.item i - --- running time: \(O(\log i)\) -outerBoundaryEdge :: Int -> Polygon t p r -> LineSegment 2 p r -outerBoundaryEdge i p = let u = p^.outerVertex i - v = p^.outerVertex (i+1) - in LineSegment (Closed u) (Open v) - - --- | Get all holes in a polygon -holeList :: Polygon t p r -> [Polygon Simple p r] -holeList (SimplePolygon _) = [] -holeList (MultiPolygon _ hs) = hs - - --- | The vertices in the polygon. No guarantees are given on the order in which --- they appear! -polygonVertices :: Polygon t p r - -> NonEmpty.NonEmpty (Point 2 r :+ p) -polygonVertices (SimplePolygon vs) = toNonEmpty vs -polygonVertices (MultiPolygon vs hs) = - sconcat $ toNonEmpty vs NonEmpty.:| map polygonVertices hs - - --- | Creates a simple polygon from the given list of vertices. --- --- pre: the input list constains no repeated vertices. -fromPoints :: [Point 2 r :+ p] -> SimplePolygon p r -fromPoints = SimplePolygon . C.fromList - - --- | The edges along the outer boundary of the polygon. The edges are half open. --- --- running time: \(O(n)\) -outerBoundaryEdges :: Polygon t p r -> C.CSeq (LineSegment 2 p r) -outerBoundaryEdges = toEdges . (^.outerBoundary) - --- | Lists all edges. The edges on the outer boundary are given before the ones --- on the holes. However, no other guarantees are given on the order. --- --- running time: \(O(n)\) -listEdges :: Polygon t p r -> [LineSegment 2 p r] -listEdges pg = let f = F.toList . outerBoundaryEdges - in f pg <> concatMap f (holeList pg) - --- | Pairs every vertex with its incident edges. The first one is its --- predecessor edge, the second one its successor edge. --- --- >>> 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] :+ ())) -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) -withIncidentEdges (MultiPolygon vs hs) = MultiPolygon vs' hs' - where - (SimplePolygon vs') = withIncidentEdges $ SimplePolygon vs - hs' = map withIncidentEdges hs - --- -- | Gets the i^th edge on the outer boundary of the polygon, that is the edge ----- with vertices i and i+1 with respect to the current focus. All indices --- -- modulo n. --- -- - --- | Given the vertices of the polygon. Produce a list of edges. The edges are --- half-open. -toEdges :: C.CSeq (Point 2 r :+ p) -> C.CSeq (LineSegment 2 p r) -toEdges vs = C.zipLWith (\p q -> LineSegment (Closed p) (Open q)) vs (C.rotateR vs) - -- let vs' = F.toList vs in - -- C.fromList $ zipWith (\p q -> LineSegment (Closed p) (Open q)) vs' (tail vs' ++ vs') - - --- | Test if q lies on the boundary of the polygon. Running time: O(n) --- --- >>> point2 1 1 `onBoundary` simplePoly --- False --- >>> point2 0 0 `onBoundary` simplePoly --- True --- >>> point2 10 0 `onBoundary` simplePoly --- True --- >>> point2 5 13 `onBoundary` simplePoly --- False --- >>> point2 5 10 `onBoundary` simplePoly --- False --- >>> point2 10 5 `onBoundary` simplePoly --- True --- >>> point2 20 5 `onBoundary` simplePoly --- False --- --- TODO: testcases multipolygon -onBoundary :: (Fractional r, Ord r) => Point 2 r -> Polygon t p r -> Bool -q `onBoundary` pg = any (q `onSegment`) es - where - out = SimplePolygon $ pg^.outerBoundary - es = concatMap (F.toList . outerBoundaryEdges) $ out : holeList pg - --- | Check if a point lies inside a polygon, on the boundary, or outside of the polygon. --- Running time: O(n). --- --- >>> point2 1 1 `inPolygon` simplePoly --- Inside --- >>> point2 0 0 `inPolygon` simplePoly --- OnBoundary --- >>> point2 10 0 `inPolygon` simplePoly --- OnBoundary --- >>> point2 5 13 `inPolygon` simplePoly --- Inside --- >>> point2 5 10 `inPolygon` simplePoly --- Inside --- >>> point2 10 5 `inPolygon` simplePoly --- OnBoundary --- >>> point2 20 5 `inPolygon` simplePoly --- Outside --- --- TODO: Add some testcases with multiPolygons --- TODO: Add some more onBoundary testcases -inPolygon :: forall t p r. (Fractional r, Ord r) - => Point 2 r -> Polygon t p r - -> PointLocationResult -q `inPolygon` pg - | q `onBoundary` pg = OnBoundary - | odd kl && odd kr && not (any (q `inHole`) hs) = Inside - | otherwise = Outside - where - l = horizontalLine $ q^.yCoord - - -- Given a line segment, compute the intersection point (if a point) with the - -- line l - intersectionPoint = asA @(Point 2 r) . (`intersect` l) - - -- Count the number of intersections that the horizontal line through q - -- maxes with the polygon, that are strictly to the left and strictly to - -- the right of q. If these numbers are both odd the point lies within the polygon. - -- - -- - -- note that: - by the asA (Point 2 r) we ignore horizontal segments (as desired) - -- - by the filtering, we effectively limit l to an open-half line, starting - -- at the (open) point q. - -- - by using half-open segments as edges we avoid double counting - -- intersections that coincide with vertices. - -- - If the point is outside, and on the same height as the - -- minimum or maximum coordinate of the polygon. The number of - -- intersections to the left or right may be one. Thus - -- incorrectly classifying the point as inside. To avoid this, - -- we count both the points to the left *and* to the right of - -- p. Only if both are odd the point is inside. so that if - -- the point is outside, and on the same y-coordinate as one - -- of the extermal vertices (one ofth) - -- - -- See http://geomalgorithms.com/a03-_inclusion.html for more information. - SP kl kr = count (\p -> (p^.xCoord) `compare` (q^.xCoord)) - . mapMaybe intersectionPoint . F.toList . outerBoundaryEdges $ pg - - -- For multi polygons we have to test if we do not lie in a hole . - inHole = insidePolygon - hs = holeList pg - - count :: (a -> Ordering) -> [a] -> SP Int Int - count f = foldr (\x (SP lts gts) -> case f x of - LT -> SP (lts + 1) gts - EQ -> SP lts gts - GT -> SP lts (gts + 1)) (SP 0 0) - - --- | Test if a point lies strictly inside the polgyon. -insidePolygon :: (Fractional r, Ord r) => Point 2 r -> Polygon t p r -> Bool -q `insidePolygon` pg = q `inPolygon` pg == Inside - - --- testQ = map (`inPolygon` testPoly) [ point2 1 1 -- Inside --- , point2 0 0 -- OnBoundary --- , point2 5 14 -- Inside --- , point2 5 10 -- Inside --- , point2 10 5 -- OnBoundary --- , point2 20 5 -- Outside --- ] - --- testPoly :: SimplePolygon () Rational --- testPoly = SimplePolygon . C.fromList . map ext $ [ point2 0 0 --- , point2 10 0 --- , point2 10 10 --- , point2 5 15 --- , point2 1 11 --- ] - --- | Compute the area of a polygon -area :: Fractional r => Polygon t p r -> r -area poly@(SimplePolygon _) = abs $ signedArea poly -area (MultiPolygon vs hs) = area (SimplePolygon vs) - sum [area h | h <- hs] - - --- | Compute the signed area of a simple polygon. The the vertices are in --- clockwise order, the signed area will be negative, if the verices are given --- in counter clockwise order, the area will be positive. -signedArea :: Fractional r => SimplePolygon p r -> r -signedArea poly = x / 2 - where - x = sum [ p^.core.xCoord * q^.core.yCoord - q^.core.xCoord * p^.core.yCoord - | LineSegment' p q <- F.toList $ outerBoundaryEdges poly ] - - --- | Compute the centroid of a simple polygon. -centroid :: Fractional r => SimplePolygon p r -> Point 2 r -centroid poly = Point $ sum' xs ^/ (6 * signedArea poly) - where - xs = [ (toVec p ^+^ toVec q) ^* (p^.xCoord * q^.yCoord - q^.xCoord * p^.yCoord) - | LineSegment' (p :+ _) (q :+ _) <- F.toList $ outerBoundaryEdges poly ] - - sum' = F.foldl' (^+^) zero - - --- | Pick a point that is inside the polygon. --- --- (note: if the polygon is degenerate; i.e. has <3 vertices, we report a --- vertex of the polygon instead.) --- --- pre: the polygon is given in CCW order --- --- running time: \(O(n)\) -pickPoint :: (Ord r, Fractional r) => Polygon p t r -> Point 2 r -pickPoint pg | isTriangle pg = centroid . SimplePolygon $ pg^.outerBoundary - | otherwise = let LineSegment' (p :+ _) (q :+ _) = findDiagonal pg - in p .+^ (0.5 *^ (q .-. p)) - --- | Test if the polygon is a triangle --- --- running time: \(O(1)\) -isTriangle :: Polygon p t r -> Bool -isTriangle = \case - SimplePolygon vs -> go vs - MultiPolygon vs [] -> go vs - MultiPolygon _ _ -> False - where - go vs = case toNonEmpty vs of - (_ :| [_,_]) -> True - _ -> False - --- | Find a diagonal of the polygon. --- --- pre: the polygon is given in CCW order --- --- running time: \(O(n)\) -findDiagonal :: (Ord r, Fractional r) => Polygon t p r -> LineSegment 2 p r -findDiagonal pg = List.head . catMaybes . F.toList $ diags - -- note that a diagonal is guaranteed to exist, so the usage of head is safe. - where - vs = pg^.outerBoundary - diags = C.zip3LWith f (C.rotateL vs) vs (C.rotateR vs) - f u v w = case ccw (u^.core) (v^.core) (w^.core) of - CCW -> Just $ findDiag u v w - -- v is a convex vertex, so find a diagonal - -- (either uw) or from v to a point inside the - -- triangle - CW -> Nothing -- v is a reflex vertex - CoLinear -> Nothing -- colinear vertex!? - - -- we test if uw is a diagonal by figuring out if there is a vertex - -- strictly inside the triangle t. If there is no such vertex then uw must - -- be a diagonal (i.e. uw intersects the polygon boundary iff there is a - -- vtx inside t). If there are vertices inside the triangle, we find the - -- one z furthest from the line(segment) uw. It then follows that vz is a - -- diagonal. Indeed this is pretty much the argument used to prove that any - -- polygon can be triangulated. See BKOS Chapter 3 for details. - findDiag u v w = let t = Triangle u v w - uw = ClosedLineSegment u w - in maybe uw (ClosedLineSegment v) - . safeMaximumOn (distTo $ supportingLine uw) - . filter (\(z :+ _) -> z `inTriangle` t == Inside) - . F.toList . polygonVertices - $ pg - - distTo l (z :+ _) = sqDistanceTo z l - +module Data.Geometry.Polygon( PolygonType(..) + , Polygon(..) + , _SimplePolygon, _MultiPolygon + , SimplePolygon, MultiPolygon, SomePolygon -safeMaximumOn :: Ord b => (a -> b) -> [a] -> Maybe a -safeMaximumOn f = \case - [] -> Nothing - xs -> Just $ List.maximumBy (comparing f) xs + , fromPoints + , polygonVertices, listEdges --- | Test if the outer boundary of the polygon is in clockwise or counter --- clockwise order. --- --- running time: \(O(n)\) --- -isCounterClockwise :: (Eq r, Fractional r) => Polygon t p r -> Bool -isCounterClockwise = (\x -> x == abs x) . signedArea - . fromPoints . F.toList . (^.outerBoundary) + , outerBoundary, outerBoundaryEdges + , outerVertex, outerBoundaryEdge + , polygonHoles, polygonHoles' + , holeList --- | Orient the outer boundary to clockwise order -toClockwiseOrder :: (Eq r, Fractional r) => Polygon t p r -> Polygon t p r -toClockwiseOrder p - | isCounterClockwise p = reverseOuterBoundary p - | otherwise = p + , inPolygon, insidePolygon, onBoundary --- | Orient the outer boundary to counter clockwise order -toCounterClockWiseOrder :: (Eq r, Fractional r) => Polygon t p r -> Polygon t p r -toCounterClockWiseOrder p - | not $ isCounterClockwise p = reverseOuterBoundary p - | otherwise = p + , area, signedArea -reverseOuterBoundary :: Polygon t p r -> Polygon t p r -reverseOuterBoundary p = p&outerBoundary %~ C.reverseDirection + , centroid + , pickPoint + , isTriangle, isStarShaped --- | Convert a Polygon to a simple polygon by forgetting about any holes. -asSimplePolygon :: Polygon t p r -> SimplePolygon p r -asSimplePolygon poly@(SimplePolygon _) = poly -asSimplePolygon (MultiPolygon vs _) = SimplePolygon vs + , isCounterClockwise + , toCounterClockWiseOrder, toCounterClockWiseOrder' + , toClockwiseOrder, toClockwiseOrder' + , reverseOuterBoundary + , findDiagonal --- | Comparison that compares which point is 'larger' in the direction given by --- the vector u. -cmpExtreme :: (Num r, Ord r) - => Vector 2 r -> Point 2 r :+ p -> Point 2 r :+ q -> Ordering -cmpExtreme u p q = u `dot` (p^.core .-. q^.core) `compare` 0 + , withIncidentEdges, numberVertices + , asSimplePolygon + , extremesLinear, cmpExtreme + ) where --- | Finds the extreme points, minimum and maximum, in a given direction --- --- running time: \(O(n)\) -extremesLinear :: (Ord r, Num r) => Vector 2 r -> Polygon t p r - -> (Point 2 r :+ p, Point 2 r :+ p) -extremesLinear u p = let vs = p^.outerBoundary - f = cmpExtreme u - in (F.minimumBy f vs, F.maximumBy f vs) +import Algorithms.Geometry.LinearProgramming.LP2DRIC +import Algorithms.Geometry.LinearProgramming.Types +import Control.Lens hiding (Simple) +import Control.Monad.Random.Class +import Data.Ext +import qualified Data.Foldable as F +import Data.Geometry.HalfSpace (rightOf) +import Data.Geometry.Line +import Data.Geometry.Point +import Data.Geometry.Polygon.Core +import Data.Geometry.Polygon.Extremes --- | assigns unique integer numbers to all vertices. Numbers start from 0, and --- are increasing along the outer boundary. The vertices of holes --- 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 ()] -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' - -------------------------------------------------------------------------------- +-- * Polygons -- | Test if a Simple polygon is star-shaped. Returns a point in the kernel -- (i.e. from which the entire polygon is visible), if it exists. diff --git a/src/Data/Geometry/Polygon/Convex.hs b/src/Data/Geometry/Polygon/Convex.hs --- a/src/Data/Geometry/Polygon/Convex.hs +++ b/src/Data/Geometry/Polygon/Convex.hs @@ -12,8 +12,8 @@ -------------------------------------------------------------------------------- module Data.Geometry.Polygon.Convex( ConvexPolygon(..), simplePolygon , merge - , lowerTangent, upperTangent - , isLeftOf, isRightOf + , lowerTangent, lowerTangent' + , upperTangent, upperTangent' , extremes , maxInDirection @@ -26,22 +26,27 @@ import Control.DeepSeq import Control.Lens hiding ((:<), (:>)) -import Data.CircularSeq (focus,CSeq) +import Data.CircularSeq (CSeq) import qualified Data.CircularSeq as C import Data.Ext import qualified Data.Foldable as F -import Data.Function (on, ) +import Data.Function (on) import Data.Geometry.Box (IsBoxable(..)) import Data.Geometry.LineSegment import Data.Geometry.Point -import Data.Geometry.Polygon (fromPoints, SimplePolygon, cmpExtreme, outerBoundary) +import Data.Geometry.Polygon.Core (fromPoints, SimplePolygon, outerBoundary) +import Data.Geometry.Polygon.Extremes(cmpExtreme) import Data.Geometry.Properties import Data.Geometry.Transformation import Data.Geometry.Vector +import Data.List.NonEmpty (NonEmpty(..)) +import qualified Data.List.NonEmpty as NonEmpty import Data.Maybe (fromJust) import Data.Ord (comparing) +import Data.Semigroup.Foldable (Foldable1(..)) import Data.Sequence (viewl,viewr, ViewL(..), ViewR(..)) import qualified Data.Sequence as S +import Data.Util -- import Data.Geometry.Ipe -- import Debug.Trace @@ -88,7 +93,7 @@ -- -- print $ toPlaneGraph (Proxy :: Proxy DT) dt -- -- writeIpeFile outFile . singlePageFromContent $ out -- -- mapM_ (print . extremesNaive (v2 1 0)) polies --- pure $ map (flip rightTangent (point2 80 528)) polies +-- pure $ map (flip rightTangent (Point2 80 528)) polies @@ -220,11 +225,12 @@ rotateTo' :: Eq a => (a :+ b) -> CSeq (a :+ b) -> CSeq (a :+ b) rotateTo' x = fromJust . C.findRotateTo (coreEq x) - coreEq :: Eq a => (a :+ b) -> (a :+ b) -> Bool coreEq = (==) `on` (^.core) +-------------------------------------------------------------------------------- +-- * Computing Tangents -- | Compute the lower tangent of the two polgyons -- @@ -234,33 +240,45 @@ -- - The vertices of the polygons are given in clockwise order -- -- Running time: O(n+m), where n and m are the sizes of the two polygons respectively -lowerTangent :: (Num r, Ord r) - => ConvexPolygon p r - -> ConvexPolygon p r - -> LineSegment 2 p r -lowerTangent (getVertices -> l) (getVertices -> r) = rotate xx yy zz zz'' +lowerTangent :: (Num r, Ord r) + => ConvexPolygon p r + -> ConvexPolygon p r + -> LineSegment 2 p r +lowerTangent lp rp = ClosedLineSegment l r where - xx = rightMost l - yy = leftMost r + mkH f = NonEmpty.fromList . F.toList . f . getVertices + lh = mkH (C.rightElements . rightMost) lp + rh = mkH (C.leftElements . leftMost) rp + (Two (l :+ _) (r :+ _)) = lowerTangent' lh rh - zz = pred' yy - zz'' = succ' xx +-- | Compute the lower tangent of the two convex chains lp and rp +-- +-- pre: - the chains lp and rp have at least 1 vertex +-- - lp and rp are disjoint, and there is a vertical line +-- having lp on the left and rp on the right. +-- - The vertices in the left-chain are given in clockwise order, (right to left) +-- - The vertices in the right chain are given in counterclockwise order (left-to-right) +-- +-- The result returned is the two endpoints l and r of the tangents, +-- and the remainders lc and rc of the chains (i.e.) such that the lower hull +-- of both chains is: (reverse lc) ++ [l,h] ++ rc +-- +-- Running time: \(O(n+m)\), where n and m are the sizes of the two chains +-- respectively +lowerTangent' :: (Ord r, Num r, Foldable1 f) + => f (Point 2 r :+ p) -> f (Point 2 r :+ p) + -> Two ((Point 2 r :+ p) :+ [Point 2 r :+ p]) +lowerTangent' l0 r0 = go (toNonEmpty l0) (toNonEmpty r0) + where + ne = NonEmpty.fromList + isRight' [] _ _ = False + isRight' (x:_) l r = ccw' l r x /= CCW + go lh@(l:|ls) rh@(r:|rs) | isRight' rs l r = go lh (ne rs) + | isRight' ls l r = go (ne ls) rh + | otherwise = Two (l :+ ls) (r :+ rs) - rotate x y z z'' - | focus z `isRightOf` (focus x, focus y) = rotate x z (pred' z) z'' - -- rotate the right polygon CCW - | focus z'' `isRightOf` (focus x, focus y) = rotate z'' y z (succ' z'') - -- rotate the left polygon CW - | otherwise = ClosedLineSegment (focus x) - (focus y) -succ' :: CSeq a -> CSeq a -succ' = C.rotateR - -pred' :: CSeq a -> CSeq a -pred' = C.rotateL - -- | Compute the upper tangent of the two polgyons -- -- pre: - polygons lp and rp have at least 1 vertex @@ -269,33 +287,43 @@ -- - The vertices of the polygons are given in clockwise order -- -- Running time: O(n+m), where n and m are the sizes of the two polygons respectively -upperTangent :: (Num r, Ord r) - => ConvexPolygon p r - -> ConvexPolygon p r - -> LineSegment 2 p r -upperTangent (getVertices -> l) (getVertices -> r) = rotate xx yy zz zz' +upperTangent :: (Num r, Ord r) + => ConvexPolygon p r + -> ConvexPolygon p r + -> LineSegment 2 p r +upperTangent lp rp = ClosedLineSegment l r where - xx = rightMost l - yy = leftMost r - - zz = succ' yy - zz' = pred' xx - - rotate x y z z' - | focus z `isLeftOf` (focus x, focus y) = rotate x z (succ' z) z' - -- rotate the right polygon CW - | focus z' `isLeftOf` (focus x, focus y) = rotate z' y z (pred' z') - -- rotate the left polygon CCW - | otherwise = ClosedLineSegment (focus x) - (focus y) + mkH f = NonEmpty.fromList . F.toList . f . getVertices + lh = mkH (C.leftElements . rightMost) lp + rh = mkH (C.rightElements . leftMost) rp + (Two (l :+ _) (r :+ _)) = upperTangent' lh rh -isRightOf :: (Num r, Ord r) - => Point 2 r :+ p -> (Point 2 r :+ p', Point 2 r :+ p'') -> Bool -a `isRightOf` (b,c) = ccw (b^.core) (c^.core) (a^.core) == CW +-- | Compute the upper tangent of the two convex chains lp and rp +-- +-- pre: - the chains lp and rp have at least 1 vertex +-- - lp and rp are disjoint, and there is a vertical line +-- having lp on the left and rp on the right. +-- - The vertices in the left-chain are given in clockwise order, (right to left) +-- - The vertices in the right chain are given in counterclockwise order (left-to-right) +-- +-- The result returned is the two endpoints l and r of the tangents, +-- and the remainders lc and rc of the chains (i.e.) such that the upper hull +-- of both chains is: (reverse lc) ++ [l,h] ++ rc +-- +-- Running time: \(O(n+m)\), where n and m are the sizes of the two chains +-- respectively +upperTangent' :: (Ord r, Num r, Foldable1 f) + => f (Point 2 r :+ p) -> f (Point 2 r :+ p) + -> Two ((Point 2 r :+ p) :+ [Point 2 r :+ p]) +upperTangent' l0 r0 = go (toNonEmpty l0) (toNonEmpty r0) + where + ne = NonEmpty.fromList + isLeft' [] _ _ = False + isLeft' (x:_) l r = ccw' l r x /= CW -isLeftOf :: (Num r, Ord r) - => Point 2 r :+ p -> (Point 2 r :+ p', Point 2 r :+ p'') -> Bool -a `isLeftOf` (b,c) = ccw (b^.core) (c^.core) (a^.core) == CCW + go lh@(l:|ls) rh@(r:|rs) | isLeft' rs l r = go lh (ne rs) + | isLeft' ls l r = go (ne ls) rh + | otherwise = Two (l :+ ls) (r :+ rs) -------------------------------------------------------------------------------- @@ -363,13 +391,13 @@ -- in if p cur nxt then go xs' else xs -- test1 :: Num r => ConvexPolygon () r --- test1 = ConvexPolygon . fromPoints . map ext . reverse $ [origin, point2 1 4, point2 5 6, point2 10 3] +-- test1 = ConvexPolygon . fromPoints . map ext . reverse $ [origin, Point2 1 4, Point2 5 6, Point2 10 3] -- test2 :: Num r => ConvexPolygon () r --- test2 = ConvexPolygon . fromPoints . map ext . reverse $ [point2 11 6, point2 10 10, point2 15 18, point2 12 5] +-- test2 = ConvexPolygon . fromPoints . map ext . reverse $ [Point2 11 6, Point2 10 10, Point2 15 18, Point2 12 5] -- testA :: Num r => ConvexPolygon () r --- testA = ConvexPolygon . fromPoints . map ext $ [origin, point2 5 1, point2 2 2] +-- testA = ConvexPolygon . fromPoints . map ext $ [origin, Point2 5 1, Point2 2 2] -- testB :: Num r => ConvexPolygon () r --- testB = ConvexPolygon . fromPoints . map ext $ [origin, point2 5 3, point2 (-2) 2, point2 (-2) 1] +-- testB = ConvexPolygon . fromPoints . map ext $ [origin, Point2 5 3, Point2 (-2) 2, Point2 (-2) 1] diff --git a/src/Data/Geometry/Polygon/Core.hs b/src/Data/Geometry/Polygon/Core.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Polygon/Core.hs @@ -0,0 +1,571 @@ +-------------------------------------------------------------------------------- +-- | +-- Module : Data.Geometry.Polygon.Core +-- Copyright : (C) Frank Staals +-- License : see the LICENSE file +-- Maintainer : Frank Staals +-- +-- A Polygon data type and some basic functions to interact with them. +-- +-------------------------------------------------------------------------------- +module Data.Geometry.Polygon.Core( PolygonType(..) + , Polygon(..) + , _SimplePolygon, _MultiPolygon + , SimplePolygon, MultiPolygon, SomePolygon + + + , fromPoints + + , polygonVertices, listEdges + + , outerBoundary, outerBoundaryEdges + , outerVertex, outerBoundaryEdge + + , polygonHoles, polygonHoles' + , holeList + + , inPolygon, insidePolygon, onBoundary + + , area, signedArea + + , centroid + , pickPoint + + , isTriangle + + , isCounterClockwise + , toCounterClockWiseOrder, toCounterClockWiseOrder' + , toClockwiseOrder, toClockwiseOrder' + , reverseOuterBoundary + + , findDiagonal + + , withIncidentEdges, numberVertices + + , asSimplePolygon + ) where + +import Control.DeepSeq +import Control.Lens hiding (Simple) +import Data.Bifoldable +import Data.Bifunctor +import Data.Bitraversable +import qualified Data.CircularSeq as C +import Data.Ext +import qualified Data.Foldable as F +import Data.Geometry.Boundary +import Data.Geometry.Box +import Data.Geometry.Line +import Data.Geometry.LineSegment +import Data.Geometry.Point +import Data.Geometry.Properties +import Data.Geometry.Transformation +import Data.Geometry.Triangle (Triangle(..), inTriangle) +import Data.Geometry.Vector +import qualified Data.List as List +import Data.List.NonEmpty (NonEmpty(..)) +import qualified Data.List.NonEmpty as NonEmpty +import Data.Maybe (mapMaybe, catMaybes) +import Data.Ord (comparing) +import Data.Semigroup (sconcat) +import Data.Semigroup.Foldable +import qualified Data.Sequence as Seq +import Data.Util +import Data.Vinyl.CoRec (asA) + +-------------------------------------------------------------------------------- + +{- $setup +>>> :{ +-- import qualified Data.CircularSeq as C +let simplePoly :: SimplePolygon () Rational + simplePoly = SimplePolygon . C.fromList . map ext $ [ Point2 0 0 + , Point2 10 0 + , Point2 10 10 + , Point2 5 15 + , Point2 1 11 + ] +:} -} + +-- | We distinguish between simple polygons (without holes) and Polygons with holes. +data PolygonType = Simple | Multi + + +data Polygon (t :: PolygonType) p r where + SimplePolygon :: C.CSeq (Point 2 r :+ p) -> Polygon Simple p r + MultiPolygon :: C.CSeq (Point 2 r :+ p) -> [Polygon Simple p r] -> Polygon Multi p r + +-- | Prism to 'test' if we are a simple polygon +_SimplePolygon :: Prism' (Polygon Simple p r) (C.CSeq (Point 2 r :+ p)) +_SimplePolygon = prism' SimplePolygon (\(SimplePolygon vs) -> Just vs) + +-- | Prism to 'test' if we are a Multi polygon +_MultiPolygon :: Prism' (Polygon Multi p r) (C.CSeq (Point 2 r :+ p), [Polygon Simple p r]) +_MultiPolygon = prism' (uncurry MultiPolygon) (\(MultiPolygon vs hs) -> Just (vs,hs)) + +instance Bifunctor (Polygon t) where + bimap = bimapDefault + +instance Bifoldable (Polygon t) where + bifoldMap = bifoldMapDefault + +instance Bitraversable (Polygon t) where + bitraverse f g p = case p of + SimplePolygon vs -> SimplePolygon <$> bitraverseVertices f g vs + MultiPolygon vs hs -> MultiPolygon <$> bitraverseVertices f g vs + <*> traverse (bitraverse f g) hs + +instance (NFData p, NFData r) => NFData (Polygon t p r) where + rnf (SimplePolygon vs) = rnf vs + rnf (MultiPolygon vs hs) = rnf (vs,hs) + +bitraverseVertices :: (Applicative f, Traversable t) => (p -> f q) -> (r -> f s) + -> t (Point 2 r :+ p) -> f (t (Point 2 s :+ q)) +bitraverseVertices f g = traverse (bitraverse (traverse g) f) + +type SimplePolygon = Polygon Simple + +type MultiPolygon = Polygon Multi + +-- | Either a simple or multipolygon +type SomePolygon p r = Either (Polygon Simple p r) (Polygon Multi p r) + +type instance Dimension (SomePolygon p r) = 2 +type instance NumType (SomePolygon p r) = r + +-- | Polygons are per definition 2 dimensional +type instance Dimension (Polygon t p r) = 2 +type instance NumType (Polygon t p r) = r + +instance (Show p, Show r) => Show (Polygon t p r) where + show (SimplePolygon vs) = "SimplePolygon (" <> show vs <> ")" + show (MultiPolygon vs hs) = "MultiPolygon (" <> show vs <> ") (" <> show hs <> ")" + +-- instance (Read p, Read r) => Show (Polygon t p r) where +-- show (SimplePolygon vs) = "SimplePolygon (" <> show vs <> ")" +-- show (MultiPolygon vs hs) = "MultiPolygon (" <> show vs <> ") (" <> show hs <> ")" + + + +instance (Eq p, Eq r) => Eq (Polygon t p r) where + (SimplePolygon vs) == (SimplePolygon vs') = vs == vs' + (MultiPolygon vs hs) == (MultiPolygon vs' hs') = vs == vs' && hs == hs' + +instance PointFunctor (Polygon t p) where + pmap f (SimplePolygon vs) = SimplePolygon (fmap (first f) vs) + pmap f (MultiPolygon vs hs) = MultiPolygon (fmap (first f) vs) (map (pmap f) hs) + +instance Fractional r => IsTransformable (Polygon t p r) where + transformBy = transformPointFunctor + +instance IsBoxable (Polygon t p r) where + boundingBox = boundingBoxList' . toListOf (outerBoundary.traverse.core) + +type instance IntersectionOf (Line 2 r) (Boundary (Polygon t p r)) = + '[Seq.Seq (Either (Point 2 r) (LineSegment 2 () r))] + +type instance IntersectionOf (Point 2 r) (Polygon t p r) = [NoIntersection, Point 2 r] + +instance (Fractional r, Ord r) => (Point 2 r) `IsIntersectableWith` (Polygon t p r) where + nonEmptyIntersection = defaultNonEmptyIntersection + q `intersects` pg = q `inPolygon` pg /= Outside + q `intersect` pg | q `intersects` pg = coRec q + | otherwise = coRec NoIntersection + +-- instance IsIntersectableWith (Line 2 r) (Boundary (Polygon t p r)) where +-- nonEmptyIntersection _ _ (CoRec xs) = null xs +-- l `intersect` (Boundary (SimplePolygon vs)) = +-- undefined + -- l `intersect` (Boundary (MultiPolygon vs hs)) = coRec . + -- Seq.sortBy f . Seq.fromList + -- . concatMap (unpack . (l `intersect`) . Boundary) + -- $ SimplePolygon vs : hs + -- where + -- unpack (CoRec x) = x + -- f = undefined + + + + +-- * Functions on Polygons + +outerBoundary :: forall t p r. Lens' (Polygon t p r) (C.CSeq (Point 2 r :+ p)) +outerBoundary = lens g s + where + g :: Polygon t p r -> C.CSeq (Point 2 r :+ p) + g (SimplePolygon vs) = vs + g (MultiPolygon vs _) = vs + + s :: Polygon t p r -> C.CSeq (Point 2 r :+ p) + -> Polygon t p r + s (SimplePolygon _) vs = SimplePolygon vs + s (MultiPolygon _ hs) vs = MultiPolygon vs hs + +polygonHoles :: forall p r. Lens' (Polygon Multi p r) [Polygon Simple p r] +polygonHoles = lens g s + where + g :: Polygon Multi p r -> [Polygon Simple p r] + g (MultiPolygon _ hs) = hs + s :: Polygon Multi p r -> [Polygon Simple p r] + -> Polygon Multi p r + s (MultiPolygon vs _) = MultiPolygon vs + +polygonHoles' :: Traversal' (Polygon t p r) [Polygon Simple p r] +polygonHoles' = \f -> \case + p@(SimplePolygon _) -> pure p + (MultiPolygon vs hs) -> MultiPolygon vs <$> f hs + +-- | Access the i^th vertex on the outer boundary +outerVertex :: Int -> Lens' (Polygon t p r) (Point 2 r :+ p) +outerVertex i = outerBoundary.C.item i + +-- running time: \(O(\log i)\) +outerBoundaryEdge :: Int -> Polygon t p r -> LineSegment 2 p r +outerBoundaryEdge i p = let u = p^.outerVertex i + v = p^.outerVertex (i+1) + in LineSegment (Closed u) (Open v) + + +-- | Get all holes in a polygon +holeList :: Polygon t p r -> [Polygon Simple p r] +holeList (SimplePolygon _) = [] +holeList (MultiPolygon _ hs) = hs + + +-- | The vertices in the polygon. No guarantees are given on the order in which +-- they appear! +polygonVertices :: Polygon t p r + -> NonEmpty.NonEmpty (Point 2 r :+ p) +polygonVertices (SimplePolygon vs) = toNonEmpty vs +polygonVertices (MultiPolygon vs hs) = + sconcat $ toNonEmpty vs NonEmpty.:| map polygonVertices hs + + +-- | Creates a simple polygon from the given list of vertices. +-- +-- pre: the input list constains no repeated vertices. +fromPoints :: [Point 2 r :+ p] -> SimplePolygon p r +fromPoints = SimplePolygon . C.fromList + + +-- | The edges along the outer boundary of the polygon. The edges are half open. +-- +-- running time: \(O(n)\) +outerBoundaryEdges :: Polygon t p r -> C.CSeq (LineSegment 2 p r) +outerBoundaryEdges = toEdges . (^.outerBoundary) + +-- | Lists all edges. The edges on the outer boundary are given before the ones +-- on the holes. However, no other guarantees are given on the order. +-- +-- running time: \(O(n)\) +listEdges :: Polygon t p r -> [LineSegment 2 p r] +listEdges pg = let f = F.toList . outerBoundaryEdges + in f pg <> concatMap f (holeList pg) + +-- | Pairs every vertex with its incident edges. The first one is its +-- predecessor edge, the second one its successor edge (in terms of +-- the ordering along the boundary). +-- +-- +-- >>> 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] :+ ())) +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) +withIncidentEdges (MultiPolygon vs hs) = MultiPolygon vs' hs' + where + (SimplePolygon vs') = withIncidentEdges $ SimplePolygon vs + hs' = map withIncidentEdges hs + +-- -- | Gets the i^th edge on the outer boundary of the polygon, that is the edge +---- with vertices i and i+1 with respect to the current focus. All indices +-- -- modulo n. +-- -- + +-- | Given the vertices of the polygon. Produce a list of edges. The edges are +-- half-open. +toEdges :: C.CSeq (Point 2 r :+ p) -> C.CSeq (LineSegment 2 p r) +toEdges vs = C.zipLWith (\p q -> LineSegment (Closed p) (Open q)) vs (C.rotateR vs) + -- let vs' = F.toList vs in + -- C.fromList $ zipWith (\p q -> LineSegment (Closed p) (Open q)) vs' (tail vs' ++ vs') + + +-- | Test if q lies on the boundary of the polygon. Running time: O(n) +-- +-- >>> Point2 1 1 `onBoundary` simplePoly +-- False +-- >>> Point2 0 0 `onBoundary` simplePoly +-- True +-- >>> Point2 10 0 `onBoundary` simplePoly +-- True +-- >>> Point2 5 13 `onBoundary` simplePoly +-- False +-- >>> Point2 5 10 `onBoundary` simplePoly +-- False +-- >>> Point2 10 5 `onBoundary` simplePoly +-- True +-- >>> Point2 20 5 `onBoundary` simplePoly +-- False +-- +-- TODO: testcases multipolygon +onBoundary :: (Fractional r, Ord r) => Point 2 r -> Polygon t p r -> Bool +q `onBoundary` pg = any (q `onSegment`) es + where + out = SimplePolygon $ pg^.outerBoundary + es = concatMap (F.toList . outerBoundaryEdges) $ out : holeList pg + +-- | Check if a point lies inside a polygon, on the boundary, or outside of the polygon. +-- Running time: O(n). +-- +-- >>> Point2 1 1 `inPolygon` simplePoly +-- Inside +-- >>> Point2 0 0 `inPolygon` simplePoly +-- OnBoundary +-- >>> Point2 10 0 `inPolygon` simplePoly +-- OnBoundary +-- >>> Point2 5 13 `inPolygon` simplePoly +-- Inside +-- >>> Point2 5 10 `inPolygon` simplePoly +-- Inside +-- >>> Point2 10 5 `inPolygon` simplePoly +-- OnBoundary +-- >>> Point2 20 5 `inPolygon` simplePoly +-- Outside +-- +-- TODO: Add some testcases with multiPolygons +-- TODO: Add some more onBoundary testcases +inPolygon :: forall t p r. (Fractional r, Ord r) + => Point 2 r -> Polygon t p r + -> PointLocationResult +q `inPolygon` pg + | q `onBoundary` pg = OnBoundary + | odd kl && odd kr && not (any (q `inHole`) hs) = Inside + | otherwise = Outside + where + l = horizontalLine $ q^.yCoord + + -- Given a line segment, compute the intersection point (if a point) with the + -- line l + intersectionPoint = asA @(Point 2 r) . (`intersect` l) + + -- Count the number of intersections that the horizontal line through q + -- maxes with the polygon, that are strictly to the left and strictly to + -- the right of q. If these numbers are both odd the point lies within the polygon. + -- + -- + -- note that: - by the asA (Point 2 r) we ignore horizontal segments (as desired) + -- - by the filtering, we effectively limit l to an open-half line, starting + -- at the (open) point q. + -- - by using half-open segments as edges we avoid double counting + -- intersections that coincide with vertices. + -- - If the point is outside, and on the same height as the + -- minimum or maximum coordinate of the polygon. The number of + -- intersections to the left or right may be one. Thus + -- incorrectly classifying the point as inside. To avoid this, + -- we count both the points to the left *and* to the right of + -- p. Only if both are odd the point is inside. so that if + -- the point is outside, and on the same y-coordinate as one + -- of the extermal vertices (one ofth) + -- + -- See http://geomalgorithms.com/a03-_inclusion.html for more information. + SP kl kr = count (\p -> (p^.xCoord) `compare` (q^.xCoord)) + . mapMaybe intersectionPoint . F.toList . outerBoundaryEdges $ pg + + -- For multi polygons we have to test if we do not lie in a hole . + inHole = insidePolygon + hs = holeList pg + + count :: (a -> Ordering) -> [a] -> SP Int Int + count f = foldr (\x (SP lts gts) -> case f x of + LT -> SP (lts + 1) gts + EQ -> SP lts gts + GT -> SP lts (gts + 1)) (SP 0 0) + + +-- | Test if a point lies strictly inside the polgyon. +insidePolygon :: (Fractional r, Ord r) => Point 2 r -> Polygon t p r -> Bool +q `insidePolygon` pg = q `inPolygon` pg == Inside + + +-- testQ = map (`inPolygon` testPoly) [ Point2 1 1 -- Inside +-- , Point2 0 0 -- OnBoundary +-- , Point2 5 14 -- Inside +-- , Point2 5 10 -- Inside +-- , Point2 10 5 -- OnBoundary +-- , Point2 20 5 -- Outside +-- ] + +-- testPoly :: SimplePolygon () Rational +-- testPoly = SimplePolygon . C.fromList . map ext $ [ Point2 0 0 +-- , Point2 10 0 +-- , Point2 10 10 +-- , Point2 5 15 +-- , Point2 1 11 +-- ] + +-- | Compute the area of a polygon +area :: Fractional r => Polygon t p r -> r +area poly@(SimplePolygon _) = abs $ signedArea poly +area (MultiPolygon vs hs) = area (SimplePolygon vs) - sum [area h | h <- hs] + + +-- | Compute the signed area of a simple polygon. The the vertices are in +-- clockwise order, the signed area will be negative, if the verices are given +-- in counter clockwise order, the area will be positive. +signedArea :: Fractional r => SimplePolygon p r -> r +signedArea poly = x / 2 + where + x = sum [ p^.core.xCoord * q^.core.yCoord - q^.core.xCoord * p^.core.yCoord + | LineSegment' p q <- F.toList $ outerBoundaryEdges poly ] + + +-- | Compute the centroid of a simple polygon. +centroid :: Fractional r => SimplePolygon p r -> Point 2 r +centroid poly = Point $ sum' xs ^/ (6 * signedArea poly) + where + xs = [ (toVec p ^+^ toVec q) ^* (p^.xCoord * q^.yCoord - q^.xCoord * p^.yCoord) + | LineSegment' (p :+ _) (q :+ _) <- F.toList $ outerBoundaryEdges poly ] + + sum' = F.foldl' (^+^) zero + + +-- | Pick a point that is inside the polygon. +-- +-- (note: if the polygon is degenerate; i.e. has <3 vertices, we report a +-- vertex of the polygon instead.) +-- +-- pre: the polygon is given in CCW order +-- +-- running time: \(O(n)\) +pickPoint :: (Ord r, Fractional r) => Polygon p t r -> Point 2 r +pickPoint pg | isTriangle pg = centroid . SimplePolygon $ pg^.outerBoundary + | otherwise = let LineSegment' (p :+ _) (q :+ _) = findDiagonal pg + in p .+^ (0.5 *^ (q .-. p)) + +-- | Test if the polygon is a triangle +-- +-- running time: \(O(1)\) +isTriangle :: Polygon p t r -> Bool +isTriangle = \case + SimplePolygon vs -> go vs + MultiPolygon vs [] -> go vs + MultiPolygon _ _ -> False + where + go vs = case toNonEmpty vs of + (_ :| [_,_]) -> True + _ -> False + +-- | Find a diagonal of the polygon. +-- +-- pre: the polygon is given in CCW order +-- +-- running time: \(O(n)\) +findDiagonal :: (Ord r, Fractional r) => Polygon t p r -> LineSegment 2 p r +findDiagonal pg = List.head . catMaybes . F.toList $ diags + -- note that a diagonal is guaranteed to exist, so the usage of head is safe. + where + vs = pg^.outerBoundary + diags = C.zip3LWith f (C.rotateL vs) vs (C.rotateR vs) + f u v w = case ccw (u^.core) (v^.core) (w^.core) of + CCW -> Just $ findDiag u v w + -- v is a convex vertex, so find a diagonal + -- (either uw) or from v to a point inside the + -- triangle + CW -> Nothing -- v is a reflex vertex + CoLinear -> Nothing -- colinear vertex!? + + -- we test if uw is a diagonal by figuring out if there is a vertex + -- strictly inside the triangle t. If there is no such vertex then uw must + -- be a diagonal (i.e. uw intersects the polygon boundary iff there is a + -- vtx inside t). If there are vertices inside the triangle, we find the + -- one z furthest from the line(segment) uw. It then follows that vz is a + -- diagonal. Indeed this is pretty much the argument used to prove that any + -- polygon can be triangulated. See BKOS Chapter 3 for details. + findDiag u v w = let t = Triangle u v w + uw = ClosedLineSegment u w + in maybe uw (ClosedLineSegment v) + . safeMaximumOn (distTo $ supportingLine uw) + . filter (\(z :+ _) -> z `inTriangle` t == Inside) + . F.toList . polygonVertices + $ pg + + distTo l (z :+ _) = sqDistanceTo z l + + +safeMaximumOn :: Ord b => (a -> b) -> [a] -> Maybe a +safeMaximumOn f = \case + [] -> Nothing + xs -> Just $ List.maximumBy (comparing f) xs + + +-- | Test if the outer boundary of the polygon is in clockwise or counter +-- clockwise order. +-- +-- running time: \(O(n)\) +-- +isCounterClockwise :: (Eq r, Fractional r) => Polygon t p r -> Bool +isCounterClockwise = (\x -> x == abs x) . signedArea + . fromPoints . F.toList . (^.outerBoundary) + + +-- | Make sure that every edge has the polygon's interior on its +-- right, by orienting the outer boundary into clockwise order, and +-- the inner borders (i.e. any holes, if they exist) into +-- counter-clockwise order. +-- +-- running time: \(O(n)\) +-- | Orient the outer boundary of the polygon to clockwise order +toClockwiseOrder :: (Eq r, Fractional r) => Polygon t p r -> Polygon t p r +toClockwiseOrder p = (toClockwiseOrder' p)&polygonHoles'.traverse %~ toCounterClockWiseOrder' + +-- | Orient the outer boundary into clockwise order. Leaves any holes +-- as they are. +-- +toClockwiseOrder' :: (Eq r, Fractional r) => Polygon t p r -> Polygon t p r +toClockwiseOrder' pg + | isCounterClockwise pg = reverseOuterBoundary pg + | otherwise = pg + +-- | Make sure that every edge has the polygon's interior on its left, +-- by orienting the outer boundary into counter-clockwise order, and +-- the inner borders (i.e. any holes, if they exist) into clockwise order. +-- +-- running time: \(O(n)\) +toCounterClockWiseOrder :: (Eq r, Fractional r) => Polygon t p r -> Polygon t p r +toCounterClockWiseOrder p = + (toCounterClockWiseOrder' p)&polygonHoles'.traverse %~ toClockwiseOrder' + +-- | Orient the outer boundary into counter-clockwise order. Leaves +-- any holes as they are. +-- +toCounterClockWiseOrder' :: (Eq r, Fractional r) => Polygon t p r -> Polygon t p r +toCounterClockWiseOrder' p + | not $ isCounterClockwise p = reverseOuterBoundary p + | otherwise = p + +reverseOuterBoundary :: Polygon t p r -> Polygon t p r +reverseOuterBoundary p = p&outerBoundary %~ C.reverseDirection + + +-- | Convert a Polygon to a simple polygon by forgetting about any holes. +asSimplePolygon :: Polygon t p r -> SimplePolygon p r +asSimplePolygon poly@(SimplePolygon _) = poly +asSimplePolygon (MultiPolygon vs _) = SimplePolygon vs + + +-- | assigns unique integer numbers to all vertices. Numbers start from 0, and +-- are increasing along the outer boundary. The vertices of holes +-- 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 ()]) +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/Polygon/Extremes.hs b/src/Data/Geometry/Polygon/Extremes.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Polygon/Extremes.hs @@ -0,0 +1,39 @@ +{-# LANGUAGE TemplateHaskell #-} +-------------------------------------------------------------------------------- +-- | +-- Module : Data.Geometry.Polygon.Extremes +-- Copyright : (C) Frank Staals +-- License : see the LICENSE file +-- Maintainer : Frank Staals +-- +-- Finding the Extremal vertex of a polygon in a given direction. +-- +-------------------------------------------------------------------------------- +module Data.Geometry.Polygon.Extremes( cmpExtreme + , extremesLinear + ) where + +import Control.Lens hiding (Simple) +import Data.Ext +import qualified Data.Foldable as F +import Data.Geometry.Point +import Data.Geometry.Polygon.Core +import Data.Geometry.Vector + +-------------------------------------------------------------------------------- + +-- | Comparison that compares which point is 'larger' in the direction given by +-- the vector u. +cmpExtreme :: (Num r, Ord r) + => Vector 2 r -> Point 2 r :+ p -> Point 2 r :+ q -> Ordering +cmpExtreme u p q = u `dot` (p^.core .-. q^.core) `compare` 0 + + +-- | Finds the extreme points, minimum and maximum, in a given direction +-- +-- running time: \(O(n)\) +extremesLinear :: (Ord r, Num r) => Vector 2 r -> Polygon t p r + -> (Point 2 r :+ p, Point 2 r :+ p) +extremesLinear u p = let vs = p^.outerBoundary + f = cmpExtreme u + in (F.minimumBy f vs, F.maximumBy f vs) 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 @@ -19,7 +19,7 @@ , queryRange ) where -import Algorithms.Geometry.ClosestPair.DivideAndConquer (mergeSortedListsBy) +import Algorithms.DivideAndConquer (mergeSortedListsBy) import Control.Lens import Data.BinaryTree import Data.Ext diff --git a/src/Data/Geometry/Slab.hs b/src/Data/Geometry/Slab.hs --- a/src/Data/Geometry/Slab.hs +++ b/src/Data/Geometry/Slab.hs @@ -69,8 +69,8 @@ (Slab h) `intersect` (Slab v) = coRec $ box low high where - low = point2 (v^.start.core) (h^.start.core) :+ (v^.start.extra, h^.start.extra) - high = point2 (v^.end.core) (h^.end.core) :+ (v^.end.extra, h^.end.extra) + low = Point2 (v^.start.core) (h^.start.core) :+ (v^.start.extra, h^.start.extra) + high = Point2 (v^.end.core) (h^.end.core) :+ (v^.end.extra, h^.end.extra) 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,10 +1,4 @@ -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE UndecidableInstances #-} -{-# LANGUAGE StandaloneDeriving #-} -{-# LANGUAGE GeneralizedNewtypeDeriving #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE DeriveFunctor #-} module Data.Geometry.Transformation where import Control.Lens (lens,Lens',set) @@ -130,8 +124,6 @@ -------------------------------------------------------------------------------- -- * Functions that execute transformations - --- type AlwaysTrueTransformation d = (Arity (1 + d), AlwaysTrueSnoc d, Arity d, Index' (1+d-1) (1+d)) translateBy :: ( IsTransformable g, Num (NumType g) , Arity (Dimension g), Arity (Dimension g + 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 @@ -26,6 +26,7 @@ -------------------------------------------------------------------------------- +-- | Triangles in \(d\)-dimensional space. data Triangle d p r = Triangle (Point d r :+ p) (Point d r :+ p) (Point d r :+ p) @@ -47,10 +48,11 @@ instance (Fractional r, Arity d, Arity (d + 1)) => IsTransformable (Triangle d p r) where transformBy = transformPointFunctor - -- | convenience function to construct a triangle without associated data. -triangle' :: Point d r -> Point d r -> Point d r -> Triangle d () r -triangle' p q r = Triangle (ext p) (ext q) (ext r) +pattern Triangle' :: Point d r -> Point d r -> Point d r -> Triangle d () r +pattern Triangle' p q r <- Triangle (p :+ ()) (q :+ ()) (r :+ ()) + where + Triangle' p q r = Triangle (ext p) (ext q) (ext r) sideSegments :: Triangle d p r -> [LineSegment d p r] 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 @@ -41,13 +41,13 @@ -------------------------------------------------------------------------------- type instance Dimension (Vector d r) = d -type instance NumType (Vector d r) = r +type instance NumType (Vector d r) = r instance (Arbitrary r, Arity d) => Arbitrary (Vector d r) where arbitrary = vectorFromListUnsafe <$> infiniteList --- | Test if v is a scalar multiple of u. +-- | 'isScalarmultipleof u v' test if v is a scalar multiple of u. -- -- >>> Vector2 1 1 `isScalarMultipleOf` Vector2 10 10 -- True @@ -69,7 +69,7 @@ {-# SPECIALIZE isScalarMultipleOf :: (Eq r, Fractional r) => Vector 2 r -> Vector 2 r -> Bool #-} --- | Get the scalar labmda s.t. v = lambda * u (if it exists) +-- | scalarMultiple u v computes the scalar labmda s.t. v = lambda * u (if it exists) scalarMultiple :: (Eq r, Fractional r, Arity d) => Vector d r -> Vector d r -> Maybe r scalarMultiple u v 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 @@ -282,6 +282,6 @@ destruct v = (head $ F.toList v, vectorFromListUnsafe . tail $ F.toList v) -- FIXME: this implementaion of tail is not particularly nice -snoc :: (ImplicitArity d, ImplicitArity (S d), (1 + FromPeano d) ~ (FromPeano d + 1)) - => VectorFamily d r -> r -> VectorFamily (S d) r -snoc = flip V.snoc +-- snoc :: (ImplicitArity d, ImplicitArity (S d)) +-- => VectorFamily d r -> r -> VectorFamily (S d) r +-- snoc = flip V.snoc 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 @@ -223,7 +223,7 @@ sing x = x NonEmpty.:| [] - vts = map (\(p,sp) -> (p,map (^.extra) . sortArround (ext p) <$> sp)) + vts = map (\(p,sp) -> (p,map (^.extra) . sortAround (ext p) <$> sp)) . M.assocs $ pts -- vertex Data vxData = V.fromList . map (\(p,sp) -> VertexData p (sp^._1)) $ vts @@ -270,6 +270,7 @@ vertices' = PG.vertices' . _graph -- | Enumerate all vertices, together with their vertex data +-- -- >>> mapM_ print $ vertices smallG -- (VertexId 0,VertexData {_location = Point2 [0,0], _vData = 0}) -- (VertexId 1,VertexData {_location = Point2 [2,2], _vData = 1}) diff --git a/src/Data/PlaneGraph/IO.hs b/src/Data/PlaneGraph/IO.hs --- a/src/Data/PlaneGraph/IO.hs +++ b/src/Data/PlaneGraph/IO.hs @@ -91,9 +91,8 @@ -------------------------------------------------------------------------------- --- | Transforms the planar graph into a format taht can be easily converted --- into JSON format. For every vertex, the adjacent vertices are given in --- counter clockwise order. +-- | Transforms the plane graph into adjacency lists. For every +-- vertex, the adjacent vertices are given in counter clockwise order. -- -- See 'toAdjacencyLists' for notes on how we handle self-loops. -- @@ -102,6 +101,11 @@ toAdjRep = first (\(PGA.Vtx v aj (VertexData p x)) -> Vtx v p aj x) . PGIO.toAdjRep . view graph +-- | Given the AdjacencyList representation of a plane graph, +-- construct the plane graph representing it. All the adjacencylists +-- should be in counter clockwise order. +-- +-- running time: \(O(n)\) fromAdjRep :: proxy s -> Gr (Vtx v e r) (Face f) -> PlaneGraph s v e f r fromAdjRep px = PlaneGraph . PGIO.fromAdjRep px . first (\(Vtx v p aj x) -> PGA.Vtx v aj $ VertexData p x) @@ -124,5 +128,8 @@ -- sort the adjacencies around every vertex v sort' (Vtx v p ajs x) = Vtx v p (List.sortBy (around p) ajs) x around p (a,_) (b,_) = ccwCmpAround (ext p) (location' V.! a) (location' V.! b) + -- note: since the graph is planar, there should not be + -- any pairs of points for which ccwCmpAround returns EQ + -- hence, no need to pick a secondary comparison -------------------------------------------------------------------------------- diff --git a/test/Data/PlaneGraph/myPlaneGraph.yaml b/test/Data/PlaneGraph/myPlaneGraph.yaml --- a/test/Data/PlaneGraph/myPlaneGraph.yaml +++ b/test/Data/PlaneGraph/myPlaneGraph.yaml @@ -1,21 +1,4 @@ -faces: -- incidentEdge: - - 0 - - 4 - fData: [] -- incidentEdge: - - 0 - - 2 - fData: [] -- incidentEdge: - - 0 - - 1 - fData: [] -- incidentEdge: - - 0 - - 3 - fData: [] -ajacencies: +adjacencies: - adj: - - 4 - [] @@ -88,3 +71,20 @@ vData: - [] - [] +faces: +- fData: [] + incidentEdge: + - 0 + - 4 +- fData: [] + incidentEdge: + - 0 + - 2 +- fData: [] + incidentEdge: + - 0 + - 1 +- fData: [] + incidentEdge: + - 0 + - 3 diff --git a/test/Data/PlaneGraph/small.yaml b/test/Data/PlaneGraph/small.yaml --- a/test/Data/PlaneGraph/small.yaml +++ b/test/Data/PlaneGraph/small.yaml @@ -1,16 +1,3 @@ -faces: -- incidentEdge: - - 0 - - 2 - fData: OuterFace -- incidentEdge: - - 0 - - 1 - fData: A -- incidentEdge: - - 0 - - 3 - fData: B ajacencies: - adj: - - 2 @@ -56,3 +43,16 @@ - -1 - 4 vData: 3 +faces: +- fData: OuterFace + incidentEdge: + - 0 + - 2 +- fData: A + incidentEdge: + - 0 + - 1 +- fData: B + incidentEdge: + - 0 + - 3