diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright (c) 2012, Frank Staals +Copyright (c) 2014, Frank Staals All rights reserved. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,64 @@ +HGeometry +========= + +[![Build Status](https://travis-ci.org/noinia/hgeometry.svg?branch=master)](https://travis-ci.org/noinia/hgeometry) +[![Hackage](https://img.shields.io/hackage/v/hgeometry.svg)](https://hackage.haskell.org/package/hgeometry) + +HGeometry provides some basic geometry types, and geometric algorithms and data +structures for them. The main two focusses are: (1) Strong type safety, and (2) +implementations of geometric algorithms and data structures with good +asymptotic running time guarantees. Design choices showing these aspects are +for example: + +- we provide a data type `Point d r` parameterized by a +type-level natural number `d`, representing d-dimensional points (in all cases +our type parameter `r` represents the (numeric) type for the (real)-numbers): + +```haskell +newtype Point (d :: Nat) (r :: *) = Point { toVec :: Vector d r } +``` +- the vertices of a `PolyLine d p r` are stored in a `Data.Seq2` which enforces +that a polyline is a proper polyline, and thus has at least two vertices. + +Please note that aspect (2), implementing good algorithms, is much work in +progress. HGeometry currently has only very basic types, and implements only +two algorithms: an (optimal) $O(n \log n)$ time algorithm for convex hull, and +an $O(n)$ expected time algorithm for smallest enclosing disk (both in $R^2$). + +Current work is on implementing $O(n \log n + k)$ time red-blue line segment +intersection. This would also allow for efficient polygon intersection and map +overlay. + +A Note on the Ext (:+) data type +--------------------------------- + +In many applications we do not just have geometric data, e.g. `Point d r`s or +`Polygon r`s, but instead, these types have some additional properties, like a +color, size, thickness, elevation, or whatever. Hence, we would like that our +library provides functions that also allow us to work with `ColoredPolygon r`s +etc. The typical Haskell approach would be to construct type-classes such as +`PolygonLike` and define functions that work with any type that is +`PolygonLike`. However, geometric algorithms are often hard enough by +themselves, and thus we would like all the help that the type-system/compiler +can give us. Hence, we choose to work with concrete types. + +To still allow for some extensibility our types will use the Ext (:+) type. For +example, our `Polygon` data type, has an extra type parameter `p` that allows +the vertices of the polygon to cary some extra information of type `p` (for +example a color, a size, or whatever). + +```haskell +data Polygon (t :: PolygonType) p r where + SimplePolygon :: C.CList (Point 2 r :+ p) -> Polygon Simple p r + MultiPolygon :: C.CList (Point 2 r :+ p) -> [Polygon Simple p r] -> Polygon Multi p r + ``` + +In all places this extra data is accessable by the (:+) type in Data.Ext, which +is essentially just a pair. + +Reading and Writing Ipe files +----------------------------- + +Appart from geometric types, HGeometry provides some interface for reading and +writing Ipe (http://ipe7.sourceforge.net). However, this is all very work in +progress. Hence, the API is experimental and may change at any time! diff --git a/doctests.hs b/doctests.hs new file mode 100644 --- /dev/null +++ b/doctests.hs @@ -0,0 +1,35 @@ +import Test.DocTest + +import Data.Monoid + +main = doctest $ ["-isrc" ] ++ ghcExts ++ files + +ghcExts = map ("-X" ++) + [ "TypeFamilies" + , "GADTs" + , "KindSignatures" + , "DataKinds" + , "TypeOperators" + , "ConstraintKinds" + , "PolyKinds" + , "RankNTypes" + + , "PatternSynonyms" + , "ViewPatterns" + + , "StandaloneDeriving" + , "GeneralizedNewtypeDeriving" + , "FlexibleInstances" + , "FlexibleContexts" + ] + +files = geomModules + +geomModules = map ("src/Data/Geometry/" <>) [ "Point.hs" + , "Vector.hs" + , "Line.hs" + , "LineSegment.hs" + , "PolyLine.hs" + , "Ball.hs" + , "Box.hs" + ] diff --git a/examples/bapc_examples.hs b/examples/bapc_examples.hs new file mode 100644 --- /dev/null +++ b/examples/bapc_examples.hs @@ -0,0 +1,73 @@ +module Main where + + +import qualified BAPC2012.Gunslinger +import qualified BAPC2014.Armybase + + + +import Control.Applicative +import Control.Monad(unless) +import System.Exit + +type Algorithm = String -> String + +examplesPrefixPath = "examples/" + + +data BAPCTest = BAPC { name :: String + , directory :: FilePath + , algo :: Algorithm + , files :: [(FilePath,FilePath)] + } + + +bapcTests = [ BAPC "Armybase" "BAPC2014" BAPC2014.Armybase.armybase + [ ("sample.in", "sample.out") + , ("testdata.in", "testdata.out") + ] + -- , BAPC "Gunslinger" "BAPC2012" BAPC2012.Gunslinger.gunslinger + -- [ ("sampleG.in", "sampleG.out") + -- , ("G.in", "G.out") + -- ] + ] + +main :: IO () +main = mapM_ runBAPCTest bapcTests + +-------------------------------------------------------------------------------- + +runBAPCTest :: BAPCTest -> IO () +runBAPCTest (BAPC n p alg fs) = do + let dash = replicate 80 '-' + putStrLn dash + putStrLn $ "Running tests for " ++ n + putStrLn dash + b <- runTests p alg fs + unless b exitFailure + putStrLn $ "Tests for " ++ n ++ " PASSED." + putStrLn dash + + +runTests :: FilePath -> Algorithm -> [(FilePath, FilePath)] -> IO Bool +runTests p f = runTests' f . map (both (p' ++)) + where + p' = concat [examplesPrefixPath, p, "/"] + both g (a,b) = (g a, g b) + +-- | Given an algorithm and a list of pairs: (inputFile,solutionFile), run all tests +runTests' :: Algorithm -> [(FilePath, FilePath)] -> IO Bool +runTests' f = allM (uncurry $ runTest f) + +runTest :: Algorithm + -> FilePath -> FilePath + -> IO Bool +runTest f inFile solutionFile = (\input solution -> f input == solution) + <$> readFile inFile + <*> readFile solutionFile + +-------------------------------------------------------------------------------- + + +allM :: (Functor m , Monad m) => (a -> m Bool) -> [a] -> m Bool +allM f = fmap and . mapM f diff --git a/hgeometry.cabal b/hgeometry.cabal --- a/hgeometry.cabal +++ b/hgeometry.cabal @@ -1,68 +1,179 @@ -- Initial hgeometry.cabal generated by cabal init. For further -- documentation, see http://haskell.org/cabal/users-guide/ --- The name of the package. name: hgeometry +version: 0.4.0.0 +synopsis: Data types for geometric objects, geometric algorithms, and data structures. +description: + HGeometry provides some basic geometry types, and geometric algorithms and + data structures for them. The main two focusses are: (1) Strong type safety, + and (2) implementations of geometric algorithms and data structures with good + asymptotic running time guarantees. --- The package version. See the Haskell package versioning policy (PVP) --- for standards guiding when and how versions should be incremented. --- http://www.haskell.org/haskellwiki/Package_versioning_policy --- PVP summary: +-+------- breaking API changes --- | | +----- non-breaking API additions --- | | | +--- code changes with no API change +homepage: http://fstaals.net/software/hgeometry +license: BSD3 +license-file: LICENSE +author: Frank Staals +maintainer: f.staals@uu.nl +-- copyright: -version: 0.1.1.1 +category: Geometry +build-type: Simple +extra-source-files: README.md +cabal-version: >=1.10 +source-repository head + type: git + location: http://github.com/noinia/hgeometry --- A short (one-line) description of the package. -synopsis: Geometry types in Haskell +library + exposed-modules: Data.Geometry + Data.Geometry.Properties + Data.Geometry.Vector + -- Data.Geometry.Vector.Vinyl + Data.Geometry.Transformation + Data.Geometry.Vector.VectorFixed + Data.Geometry.Interval --- A longer description of the package. -description: Several basic geometry types and functions on these types. + Data.Geometry.Point + Data.Geometry.Line + Data.Geometry.Line.Internal + Data.Geometry.LineSegment + Data.Geometry.HalfLine + Data.Geometry.PolyLine --- The license under which the package is released. -license: BSD3 + Data.Geometry.Triangle --- The file containing the license text. -license-file: LICENSE --- The package author(s). -author: Frank Staals + -- Data.Geometry.Plane --- An email address to which users can send suggestions, bug reports, and --- patches. -maintainer: f.staals@uu.nl + Data.Geometry.Box + Data.Geometry.Ball + Data.Geometry.Polygon --- A copyright notice. --- copyright: -category: Geometry + Algorithms.Geometry.ConvexHull.GrahamScan + Algorithms.Geometry.SmallestEnclosingBall.RandomizedIncrementalConstruction -build-type: Simple + -- Data.TypeLevel.Common + -- Data.TypeLevel.Filter --- Constraint on the version of Cabal needed to build this package. -cabal-version: >=1.8 + -- Data.Type.Nat + -- Data.Type.List --- -- Specify where we can find the source code --- Hs-source-dirs: src + -- Data.Vinyl.Universe.Geometry + -- Data.Vinyl.Show + -- Data.Vinyl.Extra -library - -- Specify where we can find the source code - Hs-source-dirs: src +-- Data.Geometry.Vector - GHC-Options: -Wall + Data.Geometry.Ipe + Data.Geometry.Ipe.Attributes + Data.Geometry.Ipe.Types + Data.Geometry.Ipe.Writer + Data.Geometry.Ipe.Reader + Data.Geometry.Ipe.PathParser - -- Modules exported by the library. - exposed-modules: Data.Geometry.Geometry - , Data.Geometry.Point - , Data.Geometry.BoundingBox - , Data.Geometry.Line - , Data.Geometry.Polygon - , Data.Geometry.Circle - , Data.Geometry.SetOperations + Data.Ext - -- Modules included in this library but not exported. - -- other-modules: + Data.Seq2 + System.Random.Shuffle + Control.Monad.State.Persistent - -- Other library packages from which modules are imported. - build-depends: base >=4.3 && < 5 + + other-modules: + Data.Geometry.Ipe.ParserPrimitives + + + + + -- other-extensions: + build-depends: base >= 4.7 && < 5 + , containers >= 0.5.5 + , vinyl >= 0.5 && < 0.6 + , linear >= 1.10 + , lens >= 4.2 + , singletons >= 1.0 && < 1.1 + , bifunctors >= 4.0 + + , text >= 0.11 + , bytestring >= 0.10 + + , semigroups >= 0.15 + , bifunctors >= 4.1 + , validation >= 0.4 + + , containers >= 0.5 + , parsec >= 3 + -- , tranformers >= 0.3 + + , vector >= 0.10 + , fixed-vector >= 0.6.4.0 && < 0.7 + , data-clist >= 0.0.7.2 + -- , HList >= 0.3 && < 0.4 + + , hexpat >= 0.20.7 + , mtl + , random + + + hs-source-dirs: src + + default-language: Haskell2010 + + default-extensions: TypeFamilies + , GADTs + , KindSignatures + , DataKinds + , TypeOperators + , ConstraintKinds + , PolyKinds + , RankNTypes + + , PatternSynonyms + , ViewPatterns + + , StandaloneDeriving + , GeneralizedNewtypeDeriving + , FlexibleInstances + , FlexibleContexts + + + + + +test-suite doctests + type: exitcode-stdio-1.0 + ghc-options: -threaded + main-is: doctests.hs + build-depends: base, doctest >= 0.8 + + default-language: Haskell2010 + + -- Extensions are enabled in doctests.hs + -- default-extensions: + +test-suite bapc_examples + type: exitcode-stdio-1.0 + ghc-options: -O2 + main-is: bapc_examples.hs + hs-source-dirs: examples + build-depends: base + , doctest >= 0.8 + , array >= 0.5 + , hgeometry + , lens + , data-clist + , linear + + + default-language: Haskell2010 + + default-extensions: TypeFamilies + , GADTs + , DataKinds + , TypeOperators + , ConstraintKinds + , PolyKinds + , PatternSynonyms + , ViewPatterns diff --git a/src/Algorithms/Geometry/ConvexHull/GrahamScan.hs b/src/Algorithms/Geometry/ConvexHull/GrahamScan.hs new file mode 100644 --- /dev/null +++ b/src/Algorithms/Geometry/ConvexHull/GrahamScan.hs @@ -0,0 +1,74 @@ +module Algorithms.Geometry.ConvexHull.GrahamScan( ConvexHull(..) + , DegenerateCH + , convexHull + , upperHull + , lowerHull + ) where + +import Control.Lens((^.)) +import Data.Ext +import Data.Geometry.Point +import Data.Geometry.Polygon +import qualified Data.List as L +import Data.Monoid + + +-- | Two dimensional convex hulls +newtype ConvexHull p r = ConvexHull { _hull :: (SimplePolygon p r) } +-- deriving (Show,Eq) + +type DegenerateCH p r = Maybe (Point 2 r :+ p) + + +-- | O(n log n) time ConvexHull using Graham-Scan +convexHull :: (Ord r, Num r) + => [Point 2 r :+ p] -> Either (DegenerateCH p r) (ConvexHull p r) +convexHull [] = Left Nothing +convexHull [p] = Left (Just p) +convexHull ps = let ps' = L.sortBy incXdecY ps + uh = tail . hull' $ ps' + lh = tail . hull' $ reverse ps' + in Right . ConvexHull . fromPoints $ lh ++ uh + +upperHull :: (Ord r, Num r) + => [Point 2 r :+ p] -> Either (DegenerateCH p r) [Point 2 r :+ p] +upperHull = hull id + + +lowerHull :: (Ord r, Num r) + => [Point 2 r :+ p] -> Either (DegenerateCH p r) [Point 2 r :+ p] +lowerHull = hull reverse + + +-- | Helper function so that that can compute both the upper or the lower hull, depending +-- on the function f +hull :: (Ord r, Num r) + => ([Point 2 r :+ p] -> [Point 2 r :+ p]) + -> [Point 2 r :+ p] + -> Either (DegenerateCH p r) [Point 2 r :+ p] +hull _ [] = Left Nothing +hull _ [p] = Left (Just p) +hull f ps = Right . hull' . f . L.sortBy incXdecY $ ps + + + +incXdecY :: Ord a => Ext t (Point 2 a) -> Ext t1 (Point 2 a) -> Ordering +incXdecY (Point2 px py :+ _) (Point2 qx qy :+ _) = + compare px qx <> compare qy py + + +-- | Precondition: The list of input points is sorted +hull' :: (Ord r, Num r) => [Point 2 r :+ p] -> [Point 2 r :+ p] +hull' (a:b:ps) = hull'' [b,a] ps + where + hull'' h [] = h + hull'' h (p:ps) = hull'' (cleanMiddle (p:h)) ps + + cleanMiddle [b,a] = [b,a] + cleanMiddle h@(c:b:a:rest) + | rightTurn (a^.core) (b^.core) (c^.core) = h + | otherwise = cleanMiddle (c:a:rest) + + +rightTurn :: (Ord r, Num r) => Point 2 r -> Point 2 r -> Point 2 r -> Bool +rightTurn a b c = ccw a b c == CW diff --git a/src/Algorithms/Geometry/SmallestEnclosingBall/RandomizedIncrementalConstruction.hs b/src/Algorithms/Geometry/SmallestEnclosingBall/RandomizedIncrementalConstruction.hs new file mode 100644 --- /dev/null +++ b/src/Algorithms/Geometry/SmallestEnclosingBall/RandomizedIncrementalConstruction.hs @@ -0,0 +1,88 @@ +{-# LANGUAGE DeriveFunctor #-} +{-# LANGUAGE TemplateHaskell #-} +module Algorithms.Geometry.SmallestEnclosingBall.RandomizedIncrementalConstruction where + +import Control.Lens +import Data.Ext +import qualified Data.Foldable as F +import Data.Geometry.Ball +import Data.Geometry.Point +import qualified Data.List as L +import Data.List.NonEmpty +import Data.Maybe(fromMaybe) +import Data.Monoid + +import System.Random +import System.Random.Shuffle(shuffle) + + + +-- | List of two or three elements +data TwoOrThree a = Two !a !a | Three !a !a !a deriving (Show,Read,Eq,Ord,Functor) + +instance F.Foldable TwoOrThree where + foldMap f (Two a b) = f a <> f b + foldMap f (Three a b c) = f a <> f b <> f c + +-- | The result of a smallest enclosing disk computation: The smallest ball +-- and the points defining it +data DiskResult p r = DiskResult { _enclosingDisk :: Circle () r + , _definingPoints :: TwoOrThree (Point 2 r :+ p) + } +makeLenses ''DiskResult + +-- | 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, RandomGen g) + => g + -> Point 2 r :+ p -> Point 2 r :+ p -> [Point 2 r :+ p] + -> DiskResult p r +smallestEnclosingDisk g p q pts = let (p':q':pts') = shuffle g (p:q:pts) + in smallestEnclosingDisk' p' q' pts' + + +-- | 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 $ circle (p^.core) (q^.core) (r^.core) + degen = error "smallestEnclosingDisk: Unhandled degeneracy" + -- 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/Control/Monad/State/Persistent.hs b/src/Control/Monad/State/Persistent.hs new file mode 100644 --- /dev/null +++ b/src/Control/Monad/State/Persistent.hs @@ -0,0 +1,56 @@ +{-# LANGUAGE DeriveFunctor #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +module Control.Monad.State.Persistent( PersistentStateT + , PersistentState + , store + , runPersistentStateT + , runPersistentState + ) where + + +import Control.Applicative +import Control.Monad.State +import Control.Monad.Identity(Identity(..)) + +import Data.List.NonEmpty(NonEmpty(..),(<|),toList) + + +-------------------------------------------------------------------------------- + +-- | A State monad that can store earlier versions of the state. +newtype PersistentStateT s m a = + PersistentStateT { runPersistentStateT' :: StateT (NonEmpty s) m a } + deriving (Functor,Applicative,Monad) + -- We store all the versions in reverse order + +-- | Create a snapshot of the current state and add it to the list of states +-- that we store. +store :: Monad m => PersistentStateT s m () +store = PersistentStateT $ do + ss@(s :| _) <- get + put (s <| ss) + + +instance Monad m => MonadState s (PersistentStateT s m) where + state f = PersistentStateT $ do + (s :| os) <- get + let (x,s') = f s + put (s' :| os) + return x + +-- | run a persistentStateT, returns a triplet with the value, the last state +-- and a list of all states (including the last one) in chronological order +runPersistentStateT :: Functor m => PersistentStateT s m a -> s -> m (a,s,[s]) +runPersistentStateT (PersistentStateT act) initS = f <$> runStateT act (initS :| []) + where + f (x,ss@(s :| _)) = (x, s, reverse $ toList ss) + + +-------------------------------------------------------------------------------- + +type PersistentState s = PersistentStateT s Identity + +runPersistentState :: PersistentState s a -> s -> (a,s,[s]) +runPersistentState act = runIdentity . runPersistentStateT act diff --git a/src/Data/Ext.hs b/src/Data/Ext.hs new file mode 100644 --- /dev/null +++ b/src/Data/Ext.hs @@ -0,0 +1,40 @@ +module Data.Ext where + +import Control.Applicative +import Control.Lens +import Data.Semigroup + +-------------------------------------------------------------------------------- + +data Ext extra core = core :+ extra deriving (Show,Read,Eq,Ord) + +instance Functor (Ext e) where + fmap f (c :+ e) = f c :+ e + +instance Monoid e => Applicative (Ext e) where + pure x = x :+ mempty + -- | This implementation ignores any extra values f may have + (f :+ _) <*> (c :+ce) = f c :+ ce + +instance (Semigroup core, Semigroup extra) => Semigroup (Ext extra core) where + (c :+ e) <> (c' :+ e') = c <> c' :+ e <> e' + + +type core :+ extra = Ext extra core +infixr 1 :+ + +_core :: (core :+ extra) -> core +_core (c :+ _) = c + +_extra :: (core :+ extra) -> extra +_extra (_ :+ e) = e + +core :: Lens (core :+ extra) (core' :+ extra) core core' +core = lens _core (\(_ :+ e) c -> (c :+ e)) + +extra :: Lens (core :+ extra) (core :+ extra') extra extra' +extra = lens _extra (\(c :+ _) e -> (c :+ e)) + + +only :: a -> a :+ () +only x = x :+ () diff --git a/src/Data/Geometry.hs b/src/Data/Geometry.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry.hs @@ -0,0 +1,18 @@ +module Data.Geometry( module Prop + , module T + , module P + , module V + , module L + , module Linear.Affine + ) where + + +import Data.Geometry.Properties as Prop +import Data.Geometry.Transformation as T + +import Data.Geometry.Point as P +import Data.Geometry.Vector as V +import Data.Geometry.Line as L + +import Linear.Affine hiding (Point, origin) +import Linear.Vector as V diff --git a/src/Data/Geometry/Ball.hs b/src/Data/Geometry/Ball.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Ball.hs @@ -0,0 +1,174 @@ +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE DeriveFunctor #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE UndecidableInstances #-} +module Data.Geometry.Ball where + +import Data.Ext +import Control.Lens hiding (only) +import qualified Data.List as L +import Data.Geometry.Line +import Data.Geometry.LineSegment +import Data.Geometry.Point +import Data.Geometry.Properties +import Data.Geometry.Vector +import GHC.TypeLits +import Linear.Affine(qdA, (.-.), (.+^)) +import Linear.Vector((^/),(*^),(^+^)) + +-------------------------------------------------------------------------------- +-- * A d-dimensional ball + +data Ball d p r = Ball { _center :: Point d r :+ p + , _squaredRadius :: r + } +makeLenses ''Ball + +deriving instance (Show r, Show p, Arity d) => Show (Ball d p r) +deriving instance (Eq r, Eq p, Arity d) => Eq (Ball d p r) +deriving instance Arity d => Functor (Ball d p) + +type instance NumType (Ball d p r) = r +type instance Dimension (Ball d p r) = d + +-- * Constructing Balls + +-- | Given two points on the diameter of the ball, construct a ball. +fromDiameter :: (Arity d, Fractional r) => Point d r -> Point d r -> Ball d () r +fromDiameter p q = let c = p .+^ ((q .-. p) ^/ 2) in Ball (only c) (qdA c p) + +-- | Construct a ball given the center point and a point p on the boundary. +fromCenterAndPoint :: (Arity d, Num r) => Point d r :+ p -> Point d r :+ p -> Ball d p r +fromCenterAndPoint c p = Ball c $ qdA (c^.core) (p^.core) + +-- | A d dimensional unit ball centered at the origin. +unitBall :: (Arity d, Num r) => Ball d () r +unitBall = Ball (only origin) 1 + +-- * Querying if a point lies in a ball + +-- | Result of a inBall query +data PointBallQueryResult = Inside | On | Outside deriving (Show,Read,Eq) + +inBall :: (Arity d, Ord r, Num r) + => Point d r -> Ball d p r -> PointBallQueryResult +p `inBall` (Ball c sr) = case qdA p (c^.core) `compare` sr of + LT -> Inside + EQ -> On + GT -> Outside + +-- | Test if a point lies strictly inside a ball +-- +-- >>> (point2 0.5 0) `insideBall` unitBall +-- True +-- >>> (point2 1 0) `insideBall` unitBall +-- False +-- >>> (point2 2 0) `insideBall` unitBall +-- False +insideBall :: (Arity d, Ord r, Num r) + => Point d r -> Ball d p r -> Bool +p `insideBall` b = p `inBall` b == Inside + +-- | Test if a point lies in or on the ball +-- +inClosedBall :: (Arity d, Ord r, Num r) + => Point d r -> Ball d p r -> Bool +p `inClosedBall` b = p `inBall` b /= Outside + +-- TODO: Add test cases + +-- | Test if a point lies on the boundary of a ball. +-- +-- >>> (point2 1 0) `onBall` unitBall +-- True +-- >>> (point3 1 1 0) `onBall` unitBall +-- False +onBall :: (Arity d, Ord r, Num r) + => Point d r -> Ball d p r -> Bool +p `onBall` b = p `inBall` b == On + + +-------------------------------------------------------------------------------- +-- * Circles, aka 2-dimensional Balls + +type Circle = Ball 2 + +-- | Given three points, get the circle through the three points. If the three +-- input points are colinear we return Nothing +-- +-- >>> circle (point2 0 10) (point2 10 0) (point2 (-10) 0) +-- Just (Ball {_center = Point {toVec = Vector {_unV = fromList [0.0,0.0]}} :+ (), _squaredRadius = 100.0}) +circle :: (Eq r, Fractional r) + => Point 2 r -> Point 2 r -> Point 2 r -> Maybe (Circle () r) +circle p q r = case f p `intersect` f q of + LineLineIntersection c -> Just $ Ball (only c) (qdA c p) + _ -> Nothing -- The two lines f p and f q are + -- parallel, that means the three + -- input points where colinear. + where + -- Given a point p', get the line perpendicular, and through the midpoint + -- of the line segment p'r + f p' = let v = r .-. p' + midPoint = p' .+^ (v ^/ 2) + in perpendicularTo (Line midPoint v) + + +instance (Ord r, Floating r) => (Line 2 r) `IsIntersectableWith` (Circle p r) where + + data Intersection (Line 2 r) (Circle p r) = NoLineCircleIntersection + | LineTouchesCircle (Point 2 r) + | LineCircleIntersection (Point 2 r) (Point 2 r) + deriving (Show,Eq) + + nonEmptyIntersection NoLineCircleIntersection = False + nonEmptyIntersection _ = True + + (Line p' v) `intersect` (Ball (c :+ _) r) = case discr `compare` 0 of + LT -> NoLineCircleIntersection + EQ -> LineTouchesCircle $ q' (lambda (+)) + GT -> let [l1,l2] = L.sort [lambda (-), lambda (+)] + in LineCircleIntersection (q' l1) (q' l2) + where + (Vector2 vx vy) = v + -- (px, py) is the vector/point after translating the circle s.t. it is centered at the + -- origin + pv@(Vector2 px py) = p' .-. c + + -- q alpha is a point on the translated line + q alpha = Point $ pv ^+^ alpha *^ v + -- a point q alpha after translating it back in the situation where c is the center of the circle. + q' alpha = q alpha .+^ toVec c + + -- let q lambda be the intersection point. We solve the following equation + -- solving the equation (q_x)^2 + (q_y)^2 = r^2 then yields the equation + -- L^2(vx^2 + vy^2) + L2(px*vx + py*vy) + px^2 + py^2 = 0 + -- where L = \lambda + aa = vx^2 + vy^2 + bb = 2 * (px * vx + py * vy) + cc = px^2 + py^2 - r^2 + discr = bb^2 - 4*aa*cc + discr' = sqrt discr + -- This thus gives us the following value(s) for lambda + lambda (|+-|) = (-bb |+-| discr') / (2*aa) + + +instance (Ord r, Floating r) => (LineSegment 2 p r) `IsIntersectableWith` (Circle q r) where + + data Intersection (LineSegment 2 p r) (Circle q r) = NoLineSegmentCircleIntersection + | LineSegmentTouchesCircle (Point 2 r) + | LineSegmentIntersectsCircle (Point 2 r) + | LineSegmentCrossesCircle (Point 2 r) (Point 2 r) + deriving (Show,Eq) + + nonEmptyIntersection NoLineSegmentCircleIntersection = False + nonEmptyIntersection _ = True + + s `intersect` c = case supportingLine s `intersect` c of + NoLineCircleIntersection -> NoLineSegmentCircleIntersection + LineTouchesCircle p -> if p `onSegment` s then LineSegmentTouchesCircle p + else NoLineSegmentCircleIntersection + LineCircleIntersection p q -> case (p `onSegment` s, q `onSegment` s) of + (False,False) -> NoLineSegmentCircleIntersection + (False,True) -> LineSegmentIntersectsCircle q + (True, False) -> LineSegmentIntersectsCircle p + (True, True) -> LineSegmentCrossesCircle p q diff --git a/src/Data/Geometry/BoundingBox.hs b/src/Data/Geometry/BoundingBox.hs deleted file mode 100644 --- a/src/Data/Geometry/BoundingBox.hs +++ /dev/null @@ -1,90 +0,0 @@ -{-# Language - FlexibleInstances, - UndecidableInstances - #-} -module Data.Geometry.BoundingBox( - -- ** BoundingBoxes - BoundingBox2'(..) - , IsBoxable(..) - , mergeBoxes - , bbFromPoints - , bbLeft - , bbRight - , bbTop - , bbBottom - , width - , height - ) where - - -import Data.Geometry.Point -import Data.Geometry.Geometry - ---------------------------------------------------------------------- --- | Bounding boxes - --- | Note that a bounding box is always axis parallel, so rotating may have not --- | the expected effect - - -data BoundingBox2' a = BoundingBox2 { lowerLeft :: Point2' a - , upperRight :: Point2' a - } - deriving (Show,Eq,Read) - - -instance IsPoint2Functor BoundingBox2' where - p2fmap f (BoundingBox2 p q) = BoundingBox2 (f p) (f q) - -instance HasPoints BoundingBox2' where - points (BoundingBox2 p@(Point2 (x,y)) q@(Point2 (a,b))) = - [p,Point2 (x,b), q, Point2 (a,y)] - ---------------------------------------------------------------------- --- | - --- | A class of objects for which we can compute a boundingbox -class IsBoxable g where - boundingBox :: Ord a => g a -> BoundingBox2' a - - -- TODO: it would be nice if we can use some similar trick as used in - -- show, to get an instance of IsBoxable for things of type [g a] - bbFromList :: Ord a => [g a] -> BoundingBox2' a - bbFromList = mergeBoxes . map boundingBox - - -bbFromPoints :: Ord a => [Point2' a] -> BoundingBox2' a -bbFromPoints pts = BoundingBox2 (Point2 (llx,lly)) (Point2 (urx,ury)) - where - xs = map getX pts - ys = map getY pts - llx = minimum xs - lly = minimum ys - urx = maximum xs - ury = maximum ys - --- | get the bounding box of a list of things -mergeBoxes :: Ord a => [BoundingBox2' a] -> BoundingBox2' a -mergeBoxes = bbFromPoints . concatMap points - -instance HasPoints p => IsBoxable p where - boundingBox = bbFromPoints . points - - -width :: Num a => BoundingBox2' a -> a -width b = bbRight b - bbLeft b - -height :: Num a => BoundingBox2' a -> a -height b = bbTop b - bbBottom b - -bbLeft :: BoundingBox2' a -> a -bbLeft = getX . lowerLeft - -bbRight :: BoundingBox2' a -> a -bbRight = getX . upperRight - -bbTop :: BoundingBox2' a -> a -bbTop = getY . upperRight - -bbBottom :: BoundingBox2' a -> a -bbBottom = getY . lowerLeft diff --git a/src/Data/Geometry/Box.hs b/src/Data/Geometry/Box.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Box.hs @@ -0,0 +1,184 @@ +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE DeriveFunctor #-} +module Data.Geometry.Box( + -- * d-dimensional boxes + Box(..) + + -- * Constructing bounding boxes + , IsBoxable(..) + , IsAlwaysTrueBoundingBox + , boundingBoxList + + -- * Functions on d-dimensonal boxes + , minPoint, maxPoint + , extent, size, widthIn + , inBox + + -- * Rectangles, aka 2-dimensional boxes + , Rectangle + , width , height + ) where + +import Control.Applicative +import qualified Data.Foldable as F +import Data.Maybe(catMaybes, maybe) +import Data.Semigroup + +import Control.Lens(Getter,to,(^.),view) +import Data.Ext +import Data.Geometry.Point +import qualified Data.Geometry.Interval as I +import Data.Geometry.Properties +import Data.Geometry.Transformation +import qualified Data.Geometry.Vector as V +import Data.Geometry.Vector(Vector, Arity, Index',C(..)) +import Linear.Affine((.-.)) + +import qualified Data.Vector.Fixed as FV + +import GHC.TypeLits +-------------------------------------------------------------------------------- + +data Box d p r = Empty + | Box { _minP :: Min (Point d r) :+ p + , _maxP :: Max (Point d r) :+ p + } + +deriving instance (Show r, Show p, Arity d) => Show (Box d p r) +deriving instance (Eq r, Eq p, Arity d) => Eq (Box d p r) +deriving instance (Ord r, Ord p, Arity d) => Ord (Box d p r) + +instance (Arity d, Ord r, Semigroup p) => Semigroup (Box d p r) where + Empty <> b = b + b <> Empty = b + (Box mi ma) <> (Box mi' ma') = Box (mi <> mi') (ma <> ma') + + +instance (Arity d, Ord r, Semigroup p) => Monoid (Box d p r) where + mempty = Empty + b `mappend` b' = b <> b' + + +instance (Arity d, Ord r) => (Box d p r) `IsIntersectableWith` (Box d p r) where + data Intersection (Box d p r) (Box d p r) = BoxBoxIntersection !(Box d () r) + + nonEmptyIntersection (BoxBoxIntersection Empty) = True + nonEmptyIntersection _ = False + + (Box a b) `intersect` (Box c d) = BoxBoxIntersection $ Box (mi :+ ()) (ma :+ ()) + where + mi = (a^.core) `max` (c^.core) + ma = (b^.core) `min` (d^.core) + +deriving instance (Show r, Show p, Arity d) => Show (Intersection (Box d p r) (Box d p r)) +deriving instance (Eq r, Eq p, Arity d) => Eq (Intersection (Box d p r) (Box d p r)) +deriving instance (Ord r, Ord p, Arity d) => Ord (Intersection (Box d p r) (Box d p r)) + + + + +-- Note that this does not guarantee the box is still a proper box +instance PointFunctor (Box d p) where + pmap f (Box mi ma) = Box (fmap f <$> mi) (fmap f <$> ma) + +instance (Num r, AlwaysTruePFT d) => IsTransformable (Box d p r) where + transformBy = transformPointFunctor + + +type instance Dimension (Box d p r) = d +type instance NumType (Box d p r) = r + +to' :: (m -> Point d r) -> (Box d p r -> m :+ p) -> + Getter (Box d p r) (Maybe (Point d r :+ p)) +to' f g = to $ \x -> case x of + Empty -> Nothing + b -> Just . fmap f . g $ b + +minPoint :: Getter (Box d p r) (Maybe (Point d r :+ p)) +minPoint = to' getMin _minP + +maxPoint :: Getter (Box d p r) (Maybe (Point d r :+ p)) +maxPoint = to' getMax _maxP + +-- | Check if a point lies a box +-- +-- >>> 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) +-- True +inBox :: (Arity d, Ord r) => Point d r -> Box d p r -> Bool +p `inBox` b = maybe False f $ extent b + where + f = FV.and . FV.zipWith I.inInterval (toVec p) + + + +-- | Get a vector with the extent of the box in each dimension. Note that the +-- 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) +-- Just (Vector {_unV = fromList [Interval {_start = 1 :+ (), _end = 10 :+ ()},Interval {_start = 2 :+ (), _end = 20 :+ ()},Interval {_start = 3 :+ (), _end = 30 :+ ()}]}) +extent :: (Arity d) + => Box d p r -> Maybe (Vector d (I.Interval p r)) +extent Empty = Nothing +extent (Box (Min a :+ p) (Max b :+ q)) = Just $ FV.zipWith f (toVec a) (toVec b) + where + f x y = I.Interval (x :+ p) (y :+ q) + +-- | 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) +-- Vector {_unV = fromList [1,2,3]} +size :: (Arity d, Num r) => Box d p r -> Vector d r +size = maybe (pure 0) (fmap I.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) +-- 1 +-- >>> widthIn (C :: C 3) (boundingBoxList [origin, point3 1 2 3] :: Box 3 () Int) +-- 3 +widthIn :: forall proxy p i d r. (Arity d, Num r, Index' (i-1) d) => proxy i -> Box d p r -> r +widthIn _ = view (V.element (C :: C (i - 1))) . size + +---------------------------------------- + +type Rectangle = Box 2 + +-- >>> 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) +-- 2 +-- >>> height (boundingBoxList [origin] :: Rectangle () Int) +-- 0 +height :: Num r => Rectangle p r -> r +height = widthIn (C :: C 2) + + +-------------------------------------------------------------------------------- + +class IsBoxable g where + boundingBox :: (Monoid p, Semigroup p, Ord (NumType g)) + => g -> Box (Dimension g) p (NumType g) + +type IsAlwaysTrueBoundingBox g p = (Semigroup p, Arity (Dimension g)) + + +boundingBoxList :: (IsBoxable g, Monoid p, F.Foldable c, Ord (NumType g) + , IsAlwaysTrueBoundingBox g p + ) => c g -> Box (Dimension g) p (NumType g) +boundingBoxList = F.foldMap boundingBox + +---------------------------------------- + +instance IsBoxable (Point d r) where + boundingBox p = Box (Min p :+ mempty) (Max p :+ mempty) diff --git a/src/Data/Geometry/Circle.hs b/src/Data/Geometry/Circle.hs deleted file mode 100644 --- a/src/Data/Geometry/Circle.hs +++ /dev/null @@ -1,95 +0,0 @@ -module Data.Geometry.Circle( Circle2'(..) - , Disc2'(..) - , IsCircleLike(..) - , inCircle - , insideCircle - , onCircle - , inDisc - , insideDisc - ) where - -import Data.Geometry.Point -import Data.Geometry.Geometry - ---------------------------------------------------------------------- --- | A circle in the plane - -data Circle2' a = Circle2 (Point2' a) a - deriving (Eq,Ord,Show,Read) - -instance HasPoints Circle2' where - points (Circle2 p _) = [p] - --- TODO: instance for transformable - ---------------------------------------------------------------------- --- | A disc in the plane (i.e. a circle inclusiding its contents) - -newtype Disc2' a = Disc2 { border :: Circle2' a } - deriving (Show,Eq,Ord,Read) - - -instance HasPoints Disc2' where - points = points . border - --- TODO: instance for transformable - --------------------------------------------------------------------------------- --- | functions on circles - --- | Class expressing functions that circlelike objects all have. Like a center --- and a radius. Minimal implementation is either getCircle or center and radius -class IsCircleLike t where - getCircle :: t a -> Circle2' a - getCircle x = Circle2 (center x) (radius x) - - center :: t a -> Point2' a - center = center . getCircle - - radius :: t a -> a - radius = radius . getCircle - - distance :: Floating a => Point2' a -> t a -> a - distance p = distance p . getCircle - - distanceToCenter :: Floating a => Point2' a -> t a -> a - distanceToCenter p = distanceToCenter p . getCircle - -instance IsCircleLike Circle2' where - center (Circle2 p _) = p - radius (Circle2 _ r) = r - distanceToCenter p (Circle2 q _) = dist p q - - distance p c@(Circle2 _ r) = distanceToCenter p c - r - -instance IsCircleLike Disc2' where - getCircle = border - --------------------------------------------------------------------------------- --- | Checking if points lie in or on a circle/disc - - --- | Squared distance to the center -l22ToCenter :: Num a => Point2' a -> Circle2' a -> a -l22ToCenter p (Circle2 q _) = l22dist p q - --- | whether or not p lies in OR on the circle c -inCircle :: (Ord a, Num a) => Point2' a -> Circle2' a -> Bool -p `inCircle` c = l22ToCenter p c <= (radius c)^2 - --- | whether or not p lies strictly inside the circle c -insideCircle :: (Num a, Ord a) => Point2' a -> Circle2' a -> Bool -p `insideCircle` c = l22ToCenter p c < (radius c)^2 - --- | whether or not p lies on the circle -onCircle :: (Eq a, Num a) => Point2' a -> Circle2' a -> Bool -p `onCircle` c = l22ToCenter p c == (radius c)^2 - - --- | whether or not a point lies in a disc: this includes its border -inDisc :: (Num a, Ord a) => Point2' a -> Disc2' a -> Bool -p `inDisc` (Disc2 c) = p `inCircle` c - --- | whether or not a point lies strictly inside a disc. -insideDisc :: (Num a, Ord a) => Point2' a -> Disc2' a -> Bool -p `insideDisc` (Disc2 c) = p `insideCircle` c diff --git a/src/Data/Geometry/Geometry.hs b/src/Data/Geometry/Geometry.hs deleted file mode 100644 --- a/src/Data/Geometry/Geometry.hs +++ /dev/null @@ -1,123 +0,0 @@ -{-# Language FlexibleInstances, - UndecidableInstances, - OverlappingInstances - #-} -module Data.Geometry.Geometry( - -- ** Point based geometries - IsPoint2Functor(..) - , IsTransformable(..) - , HasPoints(..) - , Vec3(..) - , Matrix3(..) - , identityMatrix3 - , matrix3FromLists - , matrix3FromList - , matrix3ToList - , matrix3ToLists - ) where - - -import Data.Geometry.Point - ---------------------------------------------------------------------- --- | Point based geometries - --- | A class that defines a point2 functor. This defines that every operation that --- we can do on a point we can also do on instances of this class. i.e. by --- applying the operation on the underlying points. - -class IsPoint2Functor g where - p2fmap :: (Point2' a -> Point2' b) -> g a -> g b - -class HasPoints g where - points :: g a -> [Point2' a] - -instance HasPoints Point2' where - points p = [p] - ---------------------------------------------------------------------- --- | Basic linear algebra to support affine transformations in 2D - --- | Type to represent a matrix, form is: --- [ [ a11, a12, a13 ] --- [ a21, a22, a23 ] --- [ a31, a32, a33 ] ] - -newtype Vec3 a = Vec3 (a,a,a) - deriving (Show,Eq) -newtype Matrix3 a = Matrix3 (Vec3 (Vec3 a)) - deriving (Show,Eq) - -instance Functor Vec3 where - fmap f (Vec3 (a,b,c)) = Vec3 (f a, f b, f c) - -instance Functor Matrix3 where - fmap f (Matrix3 (Vec3 (a,b,c))) = Matrix3 $ Vec3 (fmap f a, fmap f b, fmap f c) - -v3FromList :: [a] -> Vec3 a -v3FromList [a,b,c] = Vec3 (a,b,c) -v3FromList _ = error "v3FromList needs exactly 3 elements." - - --- | given a single list of 9 elements, construct a Matrix3 -matrix3FromList :: [a] -> Matrix3 a -matrix3FromList = matrix3FromLists . rtake 3 - where - rtake k ss = let (xs,ys) = splitAt k ss in - xs : rtake k ys - --- | Given a 3x3 matrix as a list of lists, convert it to a Matrix3 -matrix3FromLists :: [[a]] -> Matrix3 a -matrix3FromLists = Matrix3 . v3FromList . map v3FromList - -v3ToList :: Vec3 a -> [a] -v3ToList (Vec3 (a,b,c)) = [a,b,c] - - --- | Gather the elements of the matrix in one long list (in row by row order) -matrix3ToList :: Matrix3 a -> [a] -matrix3ToList = concat . matrix3ToLists - -matrix3ToLists :: Matrix3 a -> [[a]] -matrix3ToLists (Matrix3 (Vec3 (a,b,c))) = map v3ToList [a,b,c] - - -identityMatrix3 :: Num a => Matrix3 a -identityMatrix3 = matrix3FromLists [ [ 1, 0, 0 ] - , [ 0, 1, 0 ] - , [ 0, 0, 1 ] ] - -multiply3 :: Num a => Matrix3 a -> Vec3 a -> Vec3 a -multiply3 (Matrix3 (Vec3 ( Vec3 (a,b,c) - , Vec3 (d,e,f) - , Vec3 (g,h,i)))) (Vec3 (x,y,z)) = - Vec3 ( a*x + b*y + c*z, - d*x + e*y + f*z, - g*x + h*y + i*z) - - - - - --- | Class that indicates that something can be transformable using --- an affine transformation - -class IsTransformable g where - transformWith :: Num a => Matrix3 a -> g a -> g a - --- | Points are transformable -instance IsTransformable Point2' where - transformWith = transformPoint - -transformPoint :: Num a => Matrix3 a -> Point2' a -> Point2' a -transformPoint m (Point2 (x,y)) = Point2 (x',y') - where - v = Vec3 (x,y,1) - Vec3 (x',y',_) = multiply3 m v - --- | Everything that is built from points is transformable -tranformPointBased :: (Num a, IsPoint2Functor g) => Matrix3 a -> g a -> g a -tranformPointBased m = p2fmap (transformWith m) - -instance IsPoint2Functor g => IsTransformable g where - transformWith = tranformPointBased diff --git a/src/Data/Geometry/HalfLine.hs b/src/Data/Geometry/HalfLine.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/HalfLine.hs @@ -0,0 +1,131 @@ +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE DeriveFunctor #-} +{-# LANGUAGE UndecidableInstances #-} +module Data.Geometry.HalfLine where + +import Control.Applicative +import Control.Lens +import Data.Ext +import Data.Geometry.Interval +import Data.Geometry.Point +import Data.Geometry.Properties +import Data.Geometry.Transformation +import Data.Geometry.Vector +import Data.Geometry.Line +import Data.Geometry.LineSegment +import Linear.Vector((*^)) +import Linear.Affine(Affine(..),distanceA) + +-------------------------------------------------------------------------------- +-- * d-dimensional Half-Lines + +-- | d-dimensional Half-Lines +data HalfLine d r = HalfLine { _startPoint :: Point d r + , _halfLineDirection :: Vector d r + } +makeLenses ''HalfLine + +deriving instance (Show r, Arity d) => Show (HalfLine d r) +deriving instance (Eq r, Arity d) => Eq (HalfLine d r) +deriving instance Arity d => Functor (HalfLine d) + +type instance Dimension (HalfLine d r) = d +type instance NumType (HalfLine d r) = r + +instance HasStart (HalfLine d r) where + type StartCore (HalfLine d r) = Point d r + type StartExtra (HalfLine d r) = () + + start = lens ((:+ ()) . _startPoint) (\(HalfLine _ v) p -> HalfLine (p^.core) v) + +instance HasSupportingLine (HalfLine d r) where + supportingLine ~(HalfLine p v) = Line p v + +-- Half-Lines are transformable +instance (Num r, AlwaysTruePFT d) => IsTransformable (HalfLine d r) where + transformBy t = toHalfLine . transformPointFunctor t . toLineSegment' + where + toLineSegment' :: (Num r, Arity d) => HalfLine d r -> LineSegment d () r + toLineSegment' (HalfLine p v) = LineSegment (p :+ ()) ((p .+^ v) :+ ()) + +instance (Ord r, Fractional r) => (HalfLine 2 r) `IsIntersectableWith` (Line 2 r) where + data Intersection (HalfLine 2 r) (Line 2 r) = NoHalfLineLineIntersection + | HalfLineLineIntersection !(Point 2 r) + | HalfLineLineOverlap !(HalfLine 2 r) + deriving (Show,Eq) + + nonEmptyIntersection NoHalfLineLineIntersection = False + nonEmptyIntersection _ = True + + hl `intersect` l = case supportingLine hl `intersect` l of + SameLine _ -> HalfLineLineOverlap hl + LineLineIntersection p -> if p `onHalfLine` hl then HalfLineLineIntersection p + else NoHalfLineLineIntersection + ParallelLines -> NoHalfLineLineIntersection + + +instance (Ord r, Fractional r) => (HalfLine 2 r) `IsIntersectableWith` (HalfLine 2 r) where + data Intersection (HalfLine 2 r) (HalfLine 2 r) = NoHalfLineHalfLineIntersection + | HLHLIntersectInPoint !(Point 2 r) + | HLHLIntersectInSegment !(LineSegment 2 () r) + | HLHLIntersectInHalfLine !(HalfLine 2 r) + deriving (Show,Eq) + + nonEmptyIntersection NoHalfLineHalfLineIntersection = False + nonEmptyIntersection _ = True + + hl' `intersect` hl = case supportingLine hl' `intersect` supportingLine hl of + ParallelLines -> NoHalfLineHalfLineIntersection + LineLineIntersection p -> if p `onHalfLine` hl' && p `onHalfLine` hl then HLHLIntersectInPoint p + else NoHalfLineHalfLineIntersection + SameLine _ -> let p = _startPoint hl' + q = _startPoint hl + seg = LineSegment (p :+ ()) (q :+ ()) + in case (p `onHalfLine` hl, q `onHalfLine` hl') of + (False,False) -> NoHalfLineHalfLineIntersection + (False,True) -> HLHLIntersectInHalfLine hl + (True, False) -> HLHLIntersectInHalfLine hl' + (True, True) -> if hl == hl' then HLHLIntersectInHalfLine hl + else HLHLIntersectInSegment seg + + + +instance (Ord r, Fractional r) => (LineSegment 2 p r) `IsIntersectableWith` (HalfLine 2 r) where + data Intersection (LineSegment 2 p r) (HalfLine 2 r) = NoSegmentHalfLineIntersection + | SegmentHalfLineIntersection !(Point 2 r) + | SegmentOnHalfLine !(LineSegment 2 () r) + + nonEmptyIntersection NoSegmentHalfLineIntersection = False + nonEmptyIntersection _ = True + + s `intersect` hl = case supportingLine s `intersect` supportingLine hl of + ParallelLines -> NoSegmentHalfLineIntersection + LineLineIntersection p -> if p `onSegment` s && p `onHalfLine` hl then SegmentHalfLineIntersection p + else NoSegmentHalfLineIntersection + SameLine _ -> let p = s ^.start.core + q = s ^.end.core + r = hl ^.start.core + seg a b = LineSegment (a :+ ()) (b :+ ()) + in case (p `onHalfLine` hl, q `onHalfLine` hl) of + (False, False) -> NoSegmentHalfLineIntersection + (False, True) -> SegmentOnHalfLine $ seg r q + (True, False) -> SegmentOnHalfLine $ seg p r + (True, True) -> SegmentOnHalfLine $ seg p q + + + + +-- | Test if a point lies on a half-line +onHalfLine :: (Ord r, Fractional r, Arity d) => Point d r -> HalfLine d r -> Bool +p `onHalfLine` (HalfLine q v) = maybe False (>= 0) $ scalarMultiple (p .-. q) v + + + + + +-- | Transform a LineSegment into a half-line, by forgetting the second endpoint. +toHalfLine :: (Num r, Arity d) => LineSegment d p r -> HalfLine d r +toHalfLine (LineSegment p' q') = let p = p' ^.core + q = q' ^.core + in HalfLine p (q .-. p) diff --git a/src/Data/Geometry/Interval.hs b/src/Data/Geometry/Interval.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Interval.hs @@ -0,0 +1,123 @@ +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE DeriveFunctor #-} +module Data.Geometry.Interval( + -- * 1 dimensional Intervals + Interval(..) + , Intersection(..) + + -- * querying the start and end of intervals + , HasStart(..), HasEnd(..) + -- * Working with intervals + , width + , inInterval + ) where + +import Control.Lens +import Data.Ext +import Data.Geometry.Properties + +-------------------------------------------------------------------------------- + +data EndPointType = Open | Closed deriving (Show, Eq, Ord) + +newtype EndPoint (t :: EndPointType) a = EndPoint a deriving (Show,Read,Eq,Ord) + + +data Range l u a = Range { _lower :: EndPoint l a + , _upper :: EndPoint u a + } deriving (Show,Read,Eq,Ord) + + +-- newtype GLineSegment s t d p r = GLineSegment (GInterval s t p (Point d r)) + +-- newtype GInterval s t a r = GInterval (Range s t (r :+ a)) + +-- data CurveSegment c s t r = SubLine { _support :: c , _range :: Range s t r } + + + +-------------------------------------------------------------------------------- + +data Interval a r = Interval { _start :: r :+ a + , _end :: r :+ a + } + deriving (Show,Read,Eq,Functor) + + +class HasStart t where + type StartCore t + type StartExtra t + start :: Lens' t (StartCore t :+ StartExtra t) + +instance HasStart (Interval a r) where + type StartCore (Interval a r) = r + type StartExtra (Interval a r) = a + start = lens _start (\(Interval _ e) s -> Interval s e) + +class HasEnd t where + type EndCore t + type EndExtra t + end :: Lens' t (EndCore t :+ EndExtra t) + +instance HasEnd (Interval a r) where + type EndCore (Interval a r) = r + type EndExtra (Interval a r) = a + end = lens _end (\(Interval s _) e -> Interval s e) + + +-- | When comparing intervals, compare them lexicographically on +-- (start^.core,end^.core,start^.extra,end^.extra) +instance (Ord a, Ord r) => Ord (Interval a r) where + (Interval (s :+ a) (t :+ b)) `compare` (Interval (p :+ c) (q :+ d)) = + (s,t,a,b) `compare` (p,q,c,d) + +type instance Dimension (Interval a r) = 1 +type instance NumType (Interval a r) = r + + +instance Ord r => IsIntersectableWith (Interval a r) (Interval a r) where + + data Intersection (Interval a r) (Interval a r) = IntervalIntersection (Interval a r) + | NoOverlap + deriving (Show,Read,Eq,Ord) + + nonEmptyIntersection NoOverlap = False + nonEmptyIntersection _ = True + + (Interval a b) `intersect` (Interval c d) + | s^.core <= t^.core = IntervalIntersection $ Interval s t + | otherwise = NoOverlap + where + + s = a `maxOnCore` c + t = b `minOnCore` d + + +instance Ord r => IsUnionableWith (Interval a r) (Interval a r) where + + data Union (Interval a r) (Interval a r) = DisjointIntervals (Interval a r) (Interval a r) + | OneInterval (Interval a r) + deriving (Show,Read,Eq,Ord) + + i@(Interval a b) `union` j@(Interval c d) + | i `intersects` j = OneInterval $ Interval s t + | otherwise = DisjointIntervals i j + where + s = a `minOnCore` c + t = b `maxOnCore` d + + +maxOnCore :: Ord c => c :+ e -> c :+ e -> c :+ e +l@(lc :+ _) `maxOnCore` r@(rc :+ _) = if lc >= rc then l else r + +minOnCore :: Ord c => c :+ e -> c :+ e -> c :+ e +l@(lc :+ _) `minOnCore` r@(rc :+ _) = if lc <= rc then l else r + + +-- | Get the width of the interval +width :: Num r => Interval a r -> r +width i = i^.end.core - i^.start.core + +inInterval :: Ord r => r -> Interval a r -> Bool +x `inInterval` i = i^.start.core <= x && x <= i^.end.core diff --git a/src/Data/Geometry/Ipe.hs b/src/Data/Geometry/Ipe.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Ipe.hs @@ -0,0 +1,8 @@ +module Data.Geometry.Ipe( module I + ) where + +import qualified Data.Geometry.Ipe.Types as I +import qualified Data.Geometry.Ipe.Writer as I +import qualified Data.Geometry.Ipe.Reader as I +import qualified Data.Geometry.Ipe.Relations as I +import qualified Data.Geometry.Ipe.Attributes as I diff --git a/src/Data/Geometry/Ipe/Attributes.hs b/src/Data/Geometry/Ipe/Attributes.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Ipe/Attributes.hs @@ -0,0 +1,143 @@ +{-# LANGUAGE TemplateHaskell #-} +module Data.Geometry.Ipe.Attributes where + +import Control.Lens +import Data.Text(Text) +import Data.Geometry.Transformation(Matrix) +import Data.Singletons +import Data.Singletons.TH + + +-------------------------------------------------------------------------------- +-- | Common Attributes + +-- IpeObjects may have attributes. Essentially attributes are (key,value) +-- pairs. The key is some name. Which attributes an object can have depends on +-- the type of the object. However, all ipe objects support the following +-- 'common attributes': +data CommonAttributeUniverse = Layer | Matrix | Pin | Transformations + deriving (Show,Read,Eq) + + +-- | The CommonAttrElf family lists names of the the common attributes +-- (i.e. the keys in the key value pairs), and specifies the types that the +-- values should have. For example, it specifies that 'Matrix' is a (common) +-- attribute, and that the values of this attribute are of type 'Matrix 3 3 r'. +type family CommonAttrElf (f :: CommonAttributeUniverse) (r :: *) where + CommonAttrElf 'Layer r = Text + CommonAttrElf 'Matrix r = Matrix 3 3 r + CommonAttrElf Pin r = PinType + CommonAttrElf Transformations r = TransformationTypes + +-- | A wrapper type around common attributes. +newtype CommonAttribute r s = CommonAttribute (CommonAttrElf s r) + +-- | Possible values for Pin +data PinType = No | Yes | Horizontal | Vertical + deriving (Eq,Show,Read) + +-- | Possible values for Transformation +data TransformationTypes = Affine | Rigid | Translations deriving (Show,Read,Eq) + +-------------------------------------------------------------------------------- +-- Text Attributes + +-- these Attributes are speicifc to IpeObjects representing TextLabels and +-- MiniPages. The same structure as for the `CommonAttributes' applies here. + +-- | TODO +newtype TextLabelAttribute s r = TextLabelAttribute (CommonAttribute s r) +newtype MiniPageAttribute s r = MiniPageAttribute (CommonAttribute s r) + + +-------------------------------------------------------------------------------- +-- | Symbol Attributes + +-- | The optional Attributes for a symbol +data SymbolAttributeUniverse = SymbolStroke | SymbolFill | SymbolPen | Size + deriving (Show,Eq) + +-- | And the corresponding types +type family SymbolAttrElf (s :: SymbolAttributeUniverse) (r :: *) :: * where + SymbolAttrElf SymbolStroke r = IpeColor + SymbolAttrElf SymbolPen r = IpePen r + SymbolAttrElf SymbolFill r = IpeColor + SymbolAttrElf Size r = IpeSize r + + +-- | Wrapper around the possible SymbolAttributes +newtype SymbolAttribute r s = SymbolAttribute (SymbolAttrElf s r) + + +-- | Many types either consist of a symbolc value, or a value of type v +data IpeValue v = Named Text | Valued v deriving (Show,Eq,Ord) + +type Colour = Text -- TODO: Make this a Colour.Colour + +newtype IpeSize r = IpeSize (IpeValue r) deriving (Show,Eq,Ord) +newtype IpePen r = IpePen (IpeValue r) deriving (Show,Eq,Ord) +newtype IpeColor = IpeColor (IpeValue Colour) deriving (Show,Eq,Ord) + +------------------------------------------------------------------------------- +-- | Path Attributes + +-- | Possible attributes for a path +data PathAttributeUniverse = Stroke | Fill | Dash | Pen | LineCap | LineJoin + | FillRule | Arrow | RArrow | Opacity | Tiling | Gradient + deriving (Show,Eq) + +-- | and their types +type family PathAttrElf (s :: PathAttributeUniverse) (r :: *) :: * where + PathAttrElf Stroke r = IpeColor + PathAttrElf Fill r = IpeColor + PathAttrElf Dash r = IpeDash r + PathAttrElf Pen r = IpePen r + PathAttrElf LineCap r = Int + PathAttrElf LineJoin r = Int + PathAttrElf FillRule r = FillType + PathAttrElf Arrow r = IpeArrow r + PathAttrElf RArrow r = IpeArrow r + PathAttrElf Opacity r = IpeOpacity + PathAttrElf Tiling r = IpeTiling + PathAttrElf Gradient r = IpeGradient + +-- | Wrapper type around possible PathAttributes +newtype PathAttribute r s = PathAttribute (PathAttrElf s r) + +-- | Possible values for Dash +data IpeDash r = DashNamed Text + | DashPattern [r] r + +-- | Allowed Fill types +data FillType = Wind | EOFill deriving (Show,Read,Eq) + +-- | IpeOpacity, IpeTyling, and IpeGradient are all symbolic values +type IpeOpacity = Text +type IpeTiling = Text +type IpeGradient = Text + +-- | Possible values for an ipe arrow +data IpeArrow r = IpeArrow { _arrowName :: Text + , _arrowSize :: IpeSize r + } deriving (Show,Eq) +makeLenses ''IpeArrow + +-------------------------------------------------------------------------------- +-- | Group Attributes + + +-- | The only group attribute is a Clip +data GroupAttributeUniverse = Clip deriving (Show,Read,Eq,Ord) + +-- A clipping path is a Path. Which is defined in Data.Geometry.Ipe.Types. To +-- avoid circular imports, we define GroupAttrElf and GroupAttribute there. + + +-------------------------------------------------------------------------------- + +-- | Generate singletons for all the Universes +genSingletons [ ''CommonAttributeUniverse + , ''SymbolAttributeUniverse + , ''PathAttributeUniverse + , ''GroupAttributeUniverse + ] diff --git a/src/Data/Geometry/Ipe/ParserPrimitives.hs b/src/Data/Geometry/Ipe/ParserPrimitives.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Ipe/ParserPrimitives.hs @@ -0,0 +1,144 @@ +{-# Language FlexibleContexts #-} +{-# Language OverloadedStrings #-} +module Data.Geometry.Ipe.ParserPrimitives( runP, runP' + , pMany, pMany1, pChoice + , pChar, pSpace, pWhiteSpace, pInteger, pNatural + , (<*><>) , (<*><) + , (<***>) , (<***) , (***>) + , pMaybe , pCount , pSepBy + , Parser(..) , ParseError + , pNotFollowedBy + ) where + + +-- import Data.Functor.Identity + +import Control.Applicative hiding (many,(<|>)) + +import Text.Parsec(try) +import Text.Parsec(ParsecT(..),Parsec, Stream(..)) +import Text.Parsec.Text +import Text.ParserCombinators.Parsec hiding (Parser,try) + +import qualified Data.Text as T + + +runP' :: Parser a -> T.Text -> (a, T.Text) +runP' p s = case runP p s of + Left e -> error $ show e + Right x -> x + +runP :: Parser a -> T.Text -> Either ParseError (a,T.Text) +runP p = parse ((,) <$> p <*> getInput) "" + + +---------------------------------------------------------------------------- +-- | reexporting some standard combinators + +pMany :: Parser a -> Parser [a] +pMany = many + +pMany1 :: Parser a -> Parser [a] +pMany1 = many1 + + +pChoice :: [Parser a] -> Parser a +pChoice = choice . map try + + +pNatural :: Parser Integer +pNatural = read <$> pMany1 digit + +pInteger :: Parser Integer +pInteger = pNatural + <|> + negate <$> (pChar '-' *> pNatural) + +pChar :: Char -> Parser Char +pChar = char + +pSpace :: Parser Char +pSpace = pChar ' ' + +pWhiteSpace :: Parser Char +pWhiteSpace = space + +pMaybe :: Parser a -> Parser (Maybe a) +pMaybe = optionMaybe + +pCount :: Int -> Parser a -> Parser [a] +pCount = count + +pSepBy :: Parser a -> Parser b -> Parser [a] +pSepBy = sepBy + + +-- | infix variant of notfollowed by +p `pNotFollowedBy` q = do { x <- p ; notFollowedBy' q ; return x } + where + -- | copy of the original notFollowedBy but replaced the error message + -- to get rid of the Show dependency + notFollowedBy' p = try (do{ c <- try p; unexpected "not followed by" } + <|> return () + ) + +---------------------------------------------------------------------------- +-- | Running parsers in reverse + +infix 1 <*><>, <*>< + +-- | Runs parser q ``in reverse'' on the end of the input stream +(<*><>) :: (Reversable s, Stream s m t) + => ParsecT s u m (a -> b) -> ParsecT s u m a -> ParsecT s u m b +p <*><> q = do + rev + x <- q + rev + f <- p + return $ f x + + +p <*>< q = (\a _ -> a) <$> p <*><> q + + +rev :: (Reversable s, Stream s m t) => ParsecT s u m () +rev = getInput >>= (setInput . reverseS) + +-- as :: Parser String +-- as = many (char 'a') + +-- foo :: Parser String +-- foo = reverse <$> (string . reverse $ "foo") + +-- prs :: Parser (String,String) +-- prs = (,) <$> as <*>< foo' + +-- foo' :: Parser String +-- foo' = spaces ***> foo + +-- (<***>) :: Parser (a -> b) + +infixr 2 <***>, ***>, <*** + +-- | run the parsers in reverse order, first q, then p +(<***>) :: Monad m => m (t -> b) -> m t -> m b +p <***> q = do + x <- q + f <- p + return $ f x + +-- | the variants with missing brackets +(***>) :: (Functor m, Monad m) => m a -> m b -> m b +p ***> q = (\_ s -> s) <$> p <***> q + +(<***) :: (Functor m, Monad m) => m b -> m t -> m b +p <*** q = (\s _ -> s) <$> p <***> q + +class Reversable s where + reverseS :: s -> s + +instance Reversable [c] where + reverseS = reverse + +instance Reversable T.Text where + reverseS = T.reverse diff --git a/src/Data/Geometry/Ipe/PathParser.hs b/src/Data/Geometry/Ipe/PathParser.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Ipe/PathParser.hs @@ -0,0 +1,150 @@ +{-# Language FlexibleInstances #-} +{-# Language OverloadedStrings #-} +module Data.Geometry.Ipe.PathParser where + +import Numeric + +import Control.Applicative +import Control.Monad + +import Data.Bifunctor +import Data.Monoid(mconcat) +import Data.Semigroup +import Data.Validation(AccValidation(..)) + +import Data.Char(isSpace) +import Data.Ratio + +import Text.Parsec.Error(messageString, errorMessages) +import Data.Geometry.Point +import Data.Geometry.Vector +import Data.Geometry.Transformation +import Data.Geometry.Ipe.ParserPrimitives +import Data.Geometry.Ipe.Types(Operation(..)) +import Data.Text(Text) +import qualified Data.Text as T + +-- type Matrix d m r = () + + +----------------------------------------------------------------------- +-- | Represent stuff that can be used as a coordinate in ipe. (similar to show/read) + +class Num r => Coordinate r where + fromSeq :: Integer -> Maybe Integer -> r + +defaultFromSeq :: (Ord r, Fractional r) => Integer -> Maybe Integer -> r +defaultFromSeq x Nothing = fromInteger x +defaultFromSeq x (Just y) = let x' = fromInteger x + y' = fromInteger y + asDecimal = head . dropWhile (>= 1) . iterate (* 0.1) + in signum x' * (abs x' + asDecimal y') + +instance Coordinate Double where + fromSeq = defaultFromSeq + +instance Coordinate (Ratio Integer) where + fromSeq x Nothing = fromInteger x + fromSeq x (Just y) = fst . head $ readSigned readFloat (show x ++ "." ++ show y) + +----------------------------------------------------------------------- +-- | Running the parsers + +readPoint :: Coordinate r => Text -> Either Text (Point 2 r) +readPoint = bimap errorText fst . runP pPoint + +-- Collect errors +data Either' l r = Left' l | Right' r deriving (Show,Eq) +instance (Semigroup l, Semigroup r, Monoid r) => Monoid (Either' l r) where + mempty = Right' mempty + (Left' l) `mappend` (Left' l') = Left' $ l <> l' + (Left' l) `mappend` _ = Left' l + _ `mappend` (Left' l') = Left' l' + (Right' r) `mappend` (Right' r') = Right' $ r <> r' + +either' :: (l -> a) -> (r -> a) -> Either' l r -> a +either' lf _ (Left' l) = lf l +either' _ rf (Right' r) = rf r +-- TODO: Use Validation instead of this home-brew one + +readPathOperations :: Coordinate r => Text -> Either Text [Operation r] +readPathOperations = unWrap . mconcat . map (wrap . runP pOperation) + . clean . splitKeepDelims "mlcqeasuh" + where + -- Unwrap the Either'. If it is a Left containing all our errors, + -- combine them into one error. Otherwise just ReWrap it in an proper Either + unWrap = either' (Left . combineErrors) Right + -- for the lefts: wrap the error in a list, for the rights: we only care + -- about the result, so wrap that in a list as well. Collecting the + -- results is done using the Semigroup instance of Either' + wrap = either (Left' . (:[])) (Right' . (:[]) . fst) + -- Split the input string in pieces, each piece represents one operation + trim = T.dropWhile isSpace + clean = filter (not . T.null) . map trim + -- TODO: Do the splitting on the Text rather than unpacking and packing + -- the thing + +errorText :: ParseError -> Text +errorText = T.pack . unlines . map messageString . errorMessages + +combineErrors :: [ParseError] -> Text +combineErrors = T.unlines . map errorText + + +splitKeepDelims :: [Char] -> Text -> [Text] +splitKeepDelims delims t = maybe mPref continue $ T.uncons rest + where + mPref = if T.null pref then [] else [pref] + (pref,rest) = T.break (`elem` delims) t + continue (c,t') = pref `T.snoc` c : splitKeepDelims delims t' + + +readMatrix :: Coordinate r => Text -> Either Text (Matrix 3 3 r) +readMatrix = bimap errorText fst . runP pMatrix + +----------------------------------------------------------------------- +-- | The parsers themselves + + +pOperation :: Coordinate r => Parser (Operation r) +pOperation = pChoice [ MoveTo <$> pPoint *>> 'm' + , LineTo <$> pPoint *>> 'l' + , CurveTo <$> pPoint <*> pPoint' <*> pPoint' *>> 'c' + , QCurveTo <$> pPoint <*> pPoint' *>> 'q' + , Ellipse <$> pMatrix *>> 'e' + , ArcTo <$> pMatrix <*> pPoint' *>> 'a' + , Spline <$> pPoint `pSepBy` pWhiteSpace *>> 's' + , ClosedSpline <$> pPoint `pSepBy` pWhiteSpace *>> 'u' + , pChar 'h' *> pure ClosePath + ] + where + pPoint' = pWhiteSpace *> pPoint + p *>> c = p <*>< pWhiteSpace ***> pChar c + + +pPoint :: Coordinate r => Parser (Point 2 r) +pPoint = point2 <$> pCoordinate <* pWhiteSpace <*> pCoordinate + + +pCoordinate :: Coordinate r => Parser r +pCoordinate = fromSeq <$> pInteger <*> pDecimal + where + pDecimal = pMaybe (pChar '.' *> pInteger) + +pMatrix :: Coordinate r => Parser (Matrix 3 3 r) +pMatrix = (\a b -> mkMatrix (a:b)) <$> pCoordinate + <*> pCount 5 (pWhiteSpace *> pCoordinate) + + +-- | Generate a matrix from a list of 6 coordinates. +mkMatrix :: Coordinate r => [r] -> Matrix 3 3 r +mkMatrix [a,b,c,d,e,f] = Matrix $ v3 (v3 a c e) + (v3 b d f) + (v3 0 0 1) + -- We need the matrix in the following order: + -- 012 + -- 345 + -- + -- But ipe uses the following order: + -- 024 + -- 135 diff --git a/src/Data/Geometry/Ipe/Reader.hs b/src/Data/Geometry/Ipe/Reader.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Ipe/Reader.hs @@ -0,0 +1,119 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE FlexibleInstances #-} +module Data.Geometry.Ipe.Reader where + +import Data.Either(rights) +import Control.Applicative +import Control.Lens + +import Data.Ext +import qualified Data.Foldable as F +import Data.Vinyl + +import Data.Validation + +import Data.Geometry.Point +import Data.Geometry.Line +import Data.Geometry.LineSegment +import Data.Geometry.PolyLine +import Data.Geometry.Ipe.Types +import Data.Geometry.Ipe.Attributes + +import qualified Data.ByteString as B +import Data.Monoid +import Data.Text(Text) + +import Data.Geometry.Ipe.PathParser + +import Text.XML.Expat.Tree + +import qualified Data.Text as T + +-- fromIpeFile :: (Coordinate r, IpeRead t) => FilePath -> IO [PolyLine 2 () r] +-- fromIpeFile + + +fromIpeXML :: (Coordinate r, IpeRead t) => B.ByteString -> Either ConversionError (t r) +fromIpeXML b = (bimap (T.pack . show) id $ parse' defaultParseOptions b) >>= ipeRead + +class IpeReadText t where + ipeReadText :: Coordinate r => Text -> Either ConversionError (t r) + +type ConversionError = Text + +-- TODO: We also want to do something with the attributes + +class IpeRead t where + ipeRead :: Coordinate r => Node Text Text -> Either ConversionError (t r) + +-- instance IpeRead IpeSymbol where +-- ipeRead (Element "use" ats _) = case extract ["pos","name"] ats of + +-- given a list of keys, and a list of attributes. Extracts the values matching +-- the keys. The result is a pair (vs,others), where others are the remaining +-- attributes, and vs are the values corresponding to the keys. Sorted on +-- increasing order of their keys. +extract :: [Text] -> [(Text,Text)] -> ([Text],[(Text,Text)]) +extract = undefined + +instance IpeReadText (PolyLine 2 ()) where + ipeReadText t = readPathOperations t >>= fromOps + where + fromOps (MoveTo p:LineTo q:ops) = (\ps -> fromPoints $ [p,q] ++ ps) + <$> validateAll "Expected LineTo p" _LineTo ops + fromOps _ = Left "Expected MoveTo p:LineTo q:... " + +validateAll :: ConversionError -> Prism' (Operation r) (Point 2 r) -> [Operation r] + -> Either ConversionError [Point 2 r] +validateAll err fld = bimap T.unlines id . validateAll' err fld + + +validateAll' :: err -> Prism' (Operation r) (Point 2 r) -> [Operation r] + -> Either [err] [Point 2 r] +validateAll' err field = toEither . foldr (\op res -> f op <> res) (Right' []) + where + f op = maybe (Left' [err]) (\p -> Right' [p]) $ op ^? field + toEither = either' Left Right + +-- This is a bit of a hack +instance IpeRead (PolyLine 2 ()) where + ipeRead (Element "path" ats ts) = ipeReadText . T.unlines . map unText $ ts + -- apparently hexpat already splits the text into lines + ipeRead _ = Left "iperead: no polyline." + +unText (Text t) = t + + +instance IpeRead PathSegment where + ipeRead = fmap PolyLineSegment . ipeRead + +testP :: B.ByteString +testP = "\n128 656 m\n224 768 l\n304 624 l\n432 752 l\n" + +testO :: Text +testO = "\n128 656 m\n224 768 l\n304 624 l\n432 752 l\n" + +testPoly :: Either Text (PolyLine 2 () Double) +testPoly = fromIpeXML testP + +-- ipeRead' :: [Element Text Text] +-- ipeRead' = map ipeRead + +-- instance IpeRead (IpePage gs) where +-- ipeRead (Element "page" ats chs) = Right . IpePage [] [] . fromList' . rights $ map ipeRead chs +-- where +-- fromList' = Group' . foldr (\x r -> (IpeObject x :& RNil) :& r) RNil +-- ipeRead _ = Left "ipeRead: Not a page" + +readPolyLines :: Coordinate r => Node Text Text -> [PolyLine 2 () r] +readPolyLines (Element "ipe" _ chs) = concatMap readPolyLines' chs + + +readPolyLines' :: Coordinate r => Node Text Text -> [PolyLine 2 () r] +readPolyLines' (Element "page" _ chs) = rights $ map ipeRead chs +readPolyLines' _ = [] + +polylinesFromIpeFile :: (Coordinate r) => FilePath -> IO [PolyLine 2 () r] +polylinesFromIpeFile = fmap readPolies . B.readFile + where + readPolies = either (const []) readPolyLines . parse' defaultParseOptions diff --git a/src/Data/Geometry/Ipe/Types.hs b/src/Data/Geometry/Ipe/Types.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Ipe/Types.hs @@ -0,0 +1,347 @@ +{-# LANGUAGE DeriveFunctor #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE TemplateHaskell #-} + +{-# LANGUAGE OverloadedStrings #-} +module Data.Geometry.Ipe.Types where + +import Control.Applicative +import Control.Lens +import Data.Proxy +import Data.Vinyl + +import Linear.Affine((.-.), qdA) + +import Data.Ext +import Data.Geometry.Ball +import Data.Geometry.Point +import Data.Geometry.Properties +import Data.Geometry.Transformation(Matrix) +import Data.Geometry.Box(Rectangle) +import Data.Geometry.Line +import Data.Geometry.PolyLine + +import Data.Geometry.Ipe.Attributes +import Data.Text(Text) +import Data.TypeLevel.Filter + +import GHC.Exts + +import GHC.TypeLits + +import qualified Data.Sequence as S +import qualified Data.Seq2 as S2 + +-------------------------------------------------------------------------------- + + +-------------------------------------------------------------------------------- +-- | Image Objects + + +data Image r = Image { _imageData :: () + , _rect :: Rectangle () r + } deriving (Show,Eq,Ord) +makeLenses ''Image + +-------------------------------------------------------------------------------- +-- | Text Objects + +data TextLabel r = Label Text (Point 2 r) + deriving (Show,Eq,Ord) + +data MiniPage r = MiniPage Text (Point 2 r) r + deriving (Show,Eq,Ord) + +width :: MiniPage t -> t +width (MiniPage _ _ w) = w + +-------------------------------------------------------------------------------- +-- | Ipe Symbols, i.e. Points + +-- | A symbol (point) in ipe +data IpeSymbol r = Symbol { _symbolPoint :: Point 2 r + , _symbolName :: Text + } + deriving (Show,Eq,Ord) +makeLenses ''IpeSymbol + +type instance NumType (IpeSymbol r) = r + + +-- | Example of an IpeSymbol. I.e. A symbol that expresses that the size is 'large' +sizeSymbol :: SymbolAttribute Int Size +sizeSymbol = SymbolAttribute . IpeSize $ Named "large" + +-------------------------------------------------------------------------------- +-- | Paths + +-- | Paths consist of Path Segments. PathSegments come in the following forms: +data PathSegment r = PolyLineSegment (PolyLine 2 () r) + -- TODO + | PolygonPath + | CubicBezierSegment -- (CubicBezier 2 r) + | QuadraticBezierSegment -- (QuadraticBezier 2 r) + | EllipseSegment (Matrix 3 3 r) + | ArcSegment + | SplineSegment -- (Spline 2 r) + | ClosedSplineSegment -- (ClosedSpline 2 r) + deriving (Show,Eq,Ord) +makePrisms ''PathSegment + +-- | A path is a non-empty sequence of PathSegments. +newtype Path r = Path { _pathSegments :: S2.ViewL1 (PathSegment r) } + deriving (Show,Eq,Ord) +makeLenses ''Path + + +type instance NumType (Path r) = r + + +-- | type that represents a path in ipe. +data Operation r = MoveTo (Point 2 r) + | LineTo (Point 2 r) + | CurveTo (Point 2 r) (Point 2 r) (Point 2 r) + | QCurveTo (Point 2 r) (Point 2 r) + | Ellipse (Matrix 3 3 r) + | ArcTo (Matrix 3 3 r) (Point 2 r) + | Spline [Point 2 r] + | ClosedSpline [Point 2 r] + | ClosePath + deriving (Eq, Show) +makePrisms ''Operation + + +-------------------------------------------------------------------------------- +-- | Group Attributes + +-- | Now that we know what a Path is we can define the Attributes of a Group. +type family GroupAttrElf (s :: GroupAttributeUniverse) (r :: *) :: * where + GroupAttrElf Clip r = Path r -- strictly we event want this to be a closed path I guess + +newtype GroupAttribute r s = GroupAttribute (GroupAttrElf s r) + +-------------------------------------------------------------------------------- +-- | Groups + +-- | To define groups, we need some poly kinded, type-level only, 2-tuples. +data (a :: ka) :.: (b :: kb) + +-- | An IpeGroup can store IpeObjects. We distinguish the following different +-- ipeObjects. The parameter t will cary additional information about the +-- particular object. In particular, we will use it to keep track of the +-- attributes each item has. +-- +-- Note: We will use this type on the type-level only! In particular, we will +-- use it as a Label in a Vinyl Rec. +data IpeObjectType t = IpeGroup t + | IpeImage t + | IpeTextLabel t + | IpeMiniPage t + | IpeUse t + | IpePath t + deriving (Show,Read,Eq) + +-- | A group is essentially a hetrogenious list of IpeObjects. We represent a +-- group by means of a Vinyl Rec. +type Group gt r = Rec (IpeObject r) gt + +type instance NumType (Group gt r) = r + +-- | This type family links each 'IpeObjectType' to the type (in Haskell) that +-- we use to represent such an IpeObject. In principle we represent each object +-- by means of an Ext (:+), in which the 'core' is one of the previously seen types +-- (i.e. an Image, TextLabel, IpeSymbol, Path, etc), and the extra is a Vinyl record +-- storing the attributes. +-- +-- The different ipe types use a different Universe to draw the labels form, to +-- make sure that i.e. a Path only has attributes applicable to a Path. +-- +-- We also see the parameter of the IpeObjectType here: in all but the group case it is +-- a type level list that tells which attributes the object has. In case of a group it is a +-- poly-kinded pair (i.e. a `gt :.: gs`) where the gt is a type level list that captures +-- *ALL* type information of the objects stored in this group, and the *gs* is a type level +-- list that specifies which attributes this group itself has. +type IpeObjectElF r (f :: IpeObjectType k) = IpeObjectValueElF r f :+ IpeObjectAttrElF r f + +type family IpeObjectValueElF r (f :: IpeObjectType k) :: * where + IpeObjectValueElF r (IpeGroup (gt :.: gs)) = Group gt r + IpeObjectValueElF r (IpeImage is) = Image r + IpeObjectValueElF r (IpeTextLabel ts) = TextLabel r + IpeObjectValueElF r (IpeMiniPage mps) = MiniPage r + IpeObjectValueElF r (IpeUse ss) = IpeSymbol r + IpeObjectValueElF r (IpePath ps) = Path r + + +type family RevIpeObjectValueElF (t :: *) :: (k -> IpeObjectType k) where + RevIpeObjectValueElF (Group gt r) = IpeGroup + RevIpeObjectValueElF (Image r) = IpeImage + RevIpeObjectValueElF (TextLabel r) = IpeTextLabel + RevIpeObjectValueElF (MiniPage r) = IpeMiniPage + RevIpeObjectValueElF (IpeSymbol r ) = IpeUse + RevIpeObjectValueElF (Path r) = IpePath + + + + +type family IpeObjectAttrElF r (f :: IpeObjectType k) :: * where + IpeObjectAttrElF r (IpeGroup (gt :.: gs)) = Rec (GroupAttribute r) gs + IpeObjectAttrElF r (IpeImage is) = Rec (CommonAttribute r) is + IpeObjectAttrElF r (IpeTextLabel ts) = Rec (TextLabelAttribute r) ts + IpeObjectAttrElF r (IpeMiniPage mps) = Rec (MiniPageAttribute r) mps + IpeObjectAttrElF r (IpeUse ss) = Rec (SymbolAttribute r) ss + IpeObjectAttrElF r (IpePath ps) = Rec (PathAttribute r) ps + + + +type family IpeObjectAttrFunctorElF (f :: IpeObjectType k) :: (* -> u -> *) where + IpeObjectAttrFunctorElF (IpeGroup (gt :.: gs)) = GroupAttribute + IpeObjectAttrFunctorElF (IpeImage is) = CommonAttribute + IpeObjectAttrFunctorElF (IpeTextLabel ts) = TextLabelAttribute + IpeObjectAttrFunctorElF (IpeMiniPage mps) = MiniPageAttribute + IpeObjectAttrFunctorElF (IpeUse ss) = SymbolAttribute + IpeObjectAttrFunctorElF (IpePath ps) = PathAttribute + + + +-- TODO: Maybe split this TF into two TFS' one that determines the core type, the other +-- that gives the Attribute wrapper type +-- It would be nice if we could tell taht IpeObjecTELF was injective ... + +-- | An ipe Object is then simply a thin wrapper around the IpeObjectELF type family. +newtype IpeObject r (fld :: IpeObjectType k) = + IpeObject { _ipeObject :: IpeObjectElF r fld } + +makeLenses ''IpeObject + +type instance NumType (IpeObject r t) = r + +-------------------------------------------------------------------------------- + + +symb'' :: IpeObjectElF Int (IpeUse '[Size]) +symb'' = Symbol origin "myLargesymbol" :+ ( sizeSymbol :& RNil ) + +symb :: IpeObjectElF Int (IpeUse ('[] :: [SymbolAttributeUniverse])) +symb = Symbol origin "foo" :+ RNil + +symb' :: IpeObject Int (IpeUse '[Size]) +symb' = IpeObject symb'' + +gr :: Group '[IpeUse '[Size]] Int +gr = symb' :& RNil + +grr :: IpeObjectElF Int (IpeGroup ('[IpeUse '[Size]] + :.: + ('[] :: [GroupAttributeUniverse]) + ) + ) +grr = gr :+ RNil + + +grrr :: IpeObject Int (IpeGroup ('[IpeUse '[Size]] :.: + ('[] :: [GroupAttributeUniverse]) + ) + ) +grrr = IpeObject grr + + +points' :: forall gt r. Group gt r -> [Point 2 r] +points' = fmap (^.ipeObject.core.symbolPoint) . filterRec' + +filterRec' :: forall gt r fld. (fld ~ IpeUse '[Size]) => + Rec (IpeObject r) gt -> [IpeObject r fld] +filterRec' = undefined +-- filterRec' = filterRec (Proxy :: Proxy fld) + + + + + +-------------------------------------------------------------------------------- + + + +type XmlTree = Text + + +newtype Layer = Layer {_layerName :: Text } deriving (Show,Read,Eq,Ord,IsString) + + +-- | The definition of a view +-- make active layer into an index ? +data View = View { _layerNames :: [Layer] + , _activeLayer :: Layer + } + deriving (Eq, Ord, Show) +makeLenses ''View + + +-- | for now we pretty much ignore these +data IpeStyle = IpeStyle { _styleName :: Maybe Text + , _styleData :: XmlTree + } + deriving (Eq,Show,Read,Ord) +makeLenses ''IpeStyle + +-- | The maybe string is the encoding +data IpePreamble = IpePreamble { _encoding :: Maybe Text + , _preambleData :: XmlTree + } + deriving (Eq,Read,Show,Ord) +makeLenses ''IpePreamble + +type IpeBitmap = XmlTree + + + +-------------------------------------------------------------------------------- +-- Ipe Pages + + + + + + + +-- | An IpePage is essentially a Group, together with a list of layers and a +-- list of views. +data IpePage gs r = IpePage { _layers :: [Layer] + , _views :: [View] + , _pages :: Group gs r + } + -- deriving (Eq, Show) +makeLenses ''IpePage + + + + + +-- pGr :: IpePage '[IpeUse '[Size]] Int +pGr = IpePage [] [] gr + + + + +newtype Page r gs = Page { _unP :: IpePage gs r } +makeLenses ''Page + + +ppGr = Page pGr + +type IpePages gss r = Rec (Page r) gss + + + +-- | A complete ipe file +data IpeFile gs r = IpeFile { _preamble :: Maybe IpePreamble + , _styles :: [IpeStyle] + , _ipePages :: IpePages gs r + } + -- deriving (Eq,Show) + +-- ifP :: IpeFile '[ '[IpeUse '[Size]]] Int +ifP = IpeFile Nothing [] (ppGr :& RNil) + +makeLenses ''IpeFile diff --git a/src/Data/Geometry/Ipe/Writer.hs b/src/Data/Geometry/Ipe/Writer.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Ipe/Writer.hs @@ -0,0 +1,407 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TupleSections #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE UndecidableInstances #-} +module Data.Geometry.Ipe.Writer where + +import Control.Applicative hiding (Const(..)) +import Control.Lens((^.),(^..),(.~),(&), Prism', (#)) +import Data.Ext +import qualified Data.Foldable as F +import Data.Geometry.Ipe.Types +import qualified Data.Geometry.Ipe.Types as IT +import Data.Geometry.LineSegment +import Data.Geometry.PolyLine +import qualified Data.Geometry.Transformation as GT +import Data.Geometry.Point +import Data.Geometry.Vector +import Data.Maybe(catMaybes, mapMaybe, fromMaybe) +import Data.Monoid +import Data.Proxy +import qualified Data.Traversable as Tr +import Data.Vinyl +import Data.Vinyl.Functor +import Data.Vinyl.TypeLevel + +import Data.Geometry.Ipe.Attributes +import GHC.Exts + +import qualified Data.ByteString as B +import qualified Data.ByteString.Char8 as C +import Data.List(nub) +import qualified Data.Seq2 as S2 +import Data.Text(Text) + +import Text.XML.Expat.Tree +import Text.XML.Expat.Format(format') + +import System.IO(hPutStrLn,stderr) + +import qualified Data.Text as T + +-------------------------------------------------------------------------------- + +-- | Given a prism to convert something of type g into an ipe file, a file path, +-- and a g. Convert the geometry and write it to file. +writeIpe :: ( RecAll (Page r) gs IpeWrite + , IpeWriteText r + ) => Prism' (IpeFile gs r) g -> FilePath -> g -> IO () +writeIpe p fp g = writeIpeFile (p # g) fp + +-- | Write an IpeFiele to file. +writeIpeFile :: ( RecAll (Page r) gs IpeWrite + , IpeWriteText r + ) => IpeFile gs r -> FilePath -> IO () +writeIpeFile = writeIpeFile' + +-- | Convert the input to ipeXml, and prints it to standard out in such a way +-- that the copied text can be pasted into ipe as a geometry object. +printAsIpeSelection :: IpeWrite t => t -> IO () +printAsIpeSelection = C.putStrLn . fromMaybe "" . toIpeSelectionXML + +-- | Convert input into an ipe selection. +toIpeSelectionXML :: IpeWrite t => t -> Maybe B.ByteString +toIpeSelectionXML = fmap (format' . ipeSelection) . ipeWrite + where + ipeSelection x = Element "ipeselection" [] [x] + + +-- | Convert to Ipe xml +toIpeXML :: IpeWrite t => t -> Maybe B.ByteString +toIpeXML = fmap format' . ipeWrite + + +-- | Convert to ipe XML and write the output to a file. +writeIpeFile' :: IpeWrite t => t -> FilePath -> IO () +writeIpeFile' i fp = maybe err (B.writeFile fp) . toIpeXML $ i + where + err = hPutStrLn stderr $ + "writeIpeFile: error converting to xml. File '" <> fp <> "'not written" + +-------------------------------------------------------------------------------- + +-- | For types that can produce a text value +class IpeWriteText t where + ipeWriteText :: t -> Maybe Text + +-- | Types that correspond to an XML Element. All instances should produce an +-- Element. If the type should produce a Node with the Text constructor, use +-- the `IpeWriteText` typeclass instead. +class IpeWrite t where + ipeWrite :: t -> Maybe (Node Text Text) + + +-- | Ipe write for Exts +ipeWriteExt :: ( IpeWrite t + , RecAll a as IpeWriteText, AllSatisfy IpeAttrName as + ) => t :+ Rec a as -> Maybe (Node Text Text) +ipeWriteExt (x :+ arec) = ipeWrite x `mAddAtts` ipeWriteAttrs arec + + +-- | For the types representing attribute values we can get the name/key to use +-- when serializing to ipe. +class IpeAttrName a where + attrName :: Proxy a -> Text + + +type Attr = (Text,Text) + +-- | Functon to write all attributes in a Rec +ipeWriteAttrs :: ( RecAll f rs IpeWriteText + , AllSatisfy IpeAttrName rs + ) => Rec f rs -> [Attr] +ipeWriteAttrs rs = mapMaybe (\(x,my) -> (x,) <$> my) $ zip (writeAttrNames rs) + (writeAttrValues rs) + + +-- | Writing the attribute values +writeAttrValues :: RecAll f rs IpeWriteText => Rec f rs -> [Maybe Text] +writeAttrValues = recordToList + . rmap (\(Compose (Dict x)) -> Const $ ipeWriteText x) + . reifyConstraint (Proxy :: Proxy IpeWriteText) + +-- | Writing Attribute names +writeAttrNames :: AllSatisfy IpeAttrName rs => Rec f rs -> [Text] +writeAttrNames RNil = [] +writeAttrNames (x :& xs) = write'' x : writeAttrNames xs + where + write'' :: forall f s. IpeAttrName s => f s -> Text + write'' _ = attrName (Proxy :: Proxy s) + +-- | Function that states that all elements in xs satisfy a given constraint c +type family AllSatisfy (c :: k -> Constraint) (xs :: [k]) :: Constraint where + AllSatisfy c '[] = () + AllSatisfy c (x ': xs) = (c x, AllSatisfy c xs) + + +instance IpeWriteText Text where + ipeWriteText = Just + +-- | Add attributes to a node +addAtts :: Node Text Text -> [(Text,Text)] -> Node Text Text +n `addAtts` ats = n { eAttributes = ats ++ eAttributes n } + +-- | Same as `addAtts` but then for a Maybe node +mAddAtts :: Maybe (Node Text Text) -> [(Text, Text)] -> Maybe (Node Text Text) +mn `mAddAtts` ats = fmap (`addAtts` ats) mn + + +-------------------------------------------------------------------------------- + +instance IpeWriteText Double where + ipeWriteText = writeByShow + +instance IpeWriteText Int where + ipeWriteText = writeByShow + + +writeByShow :: Show t => t -> Maybe Text +writeByShow = ipeWriteText . T.pack . show + + + +unwords' :: [Maybe Text] -> Maybe Text +unwords' = fmap T.unwords . sequence + +unlines' :: [Maybe Text] -> Maybe Text +unlines' = fmap T.unlines . sequence + + +instance IpeWriteText r => IpeWriteText (Point 2 r) where + ipeWriteText (Point2 x y) = unwords' [ipeWriteText x, ipeWriteText y] + + +-------------------------------------------------------------------------------- + +instance IpeWriteText v => IpeWriteText (IpeValue v) where + ipeWriteText (Named t) = ipeWriteText t + ipeWriteText (Valued v) = ipeWriteText v + +deriving instance IpeWriteText r => IpeWriteText (IpeSize r) +deriving instance IpeWriteText r => IpeWriteText (IpePen r) +deriving instance IpeWriteText IpeColor + +-------------------------------------------------------------------------------- +instance IpeWriteText r => IpeWrite (IpeSymbol r) where + ipeWrite (Symbol p n) = f <$> ipeWriteText p + where + f ps = Element "use" [ ("pos", ps) + , ("name", n) + ] [] + +instance IpeWriteText (SymbolAttrElf rs r) => IpeWriteText (SymbolAttribute r rs) where + ipeWriteText (SymbolAttribute x) = ipeWriteText x + + +-- CommonAttributeUnivers +instance IpeAttrName Layer where attrName _ = "layer" +instance IpeAttrName Matrix where attrName _ = "matrix" +instance IpeAttrName Pin where attrName _ = "pin" +instance IpeAttrName Transformations where attrName _ = "transformations" + +-- IpeSymbolAttributeUniversre +instance IpeAttrName SymbolStroke where attrName _ = "stroke" +instance IpeAttrName SymbolFill where attrName _ = "fill" +instance IpeAttrName SymbolPen where attrName _ = "pen" +instance IpeAttrName Size where attrName _ = "size" + +-- PathAttributeUniverse +instance IpeAttrName Stroke where attrName _ = "stroke" +instance IpeAttrName Fill where attrName _ = "fill" +instance IpeAttrName Dash where attrName _ = "dash" +instance IpeAttrName Pen where attrName _ = "pen" +instance IpeAttrName LineCap where attrName _ = "cap" +instance IpeAttrName LineJoin where attrName _ = "join" +instance IpeAttrName FillRule where attrName _ = "fillrule" +instance IpeAttrName Arrow where attrName _ = "arrow" +instance IpeAttrName RArrow where attrName _ = "rarrow" +instance IpeAttrName Opacity where attrName _ = "opacity" +instance IpeAttrName Tiling where attrName _ = "tiling" +instance IpeAttrName Gradient where attrName _ = "gradient" + +-- GroupAttributeUniverse +instance IpeAttrName Clip where attrName _ = "clip" + +-------------------------------------------------------------------------------- + +instance IpeWriteText r => IpeWriteText (GT.Matrix 3 3 r) where + ipeWriteText (GT.Matrix m) = unwords' [a,b,c,d,e,f] + where + (Vector3 r1 r2 _) = m + + (Vector3 a c e) = ipeWriteText <$> r1 + (Vector3 b d f) = ipeWriteText <$> r2 + -- TODO: The third row should be (0,0,1) I guess. + + +instance IpeWriteText r => IpeWriteText (Operation r) where + ipeWriteText (MoveTo p) = unwords' [ ipeWriteText p, Just "m"] + ipeWriteText (LineTo p) = unwords' [ ipeWriteText p, Just "l"] + ipeWriteText (CurveTo p q r) = unwords' [ ipeWriteText p + , ipeWriteText q + , ipeWriteText r, Just "m"] + ipeWriteText (Ellipse m) = unwords' [ ipeWriteText m, Just "e"] + -- TODO: The rest + ipeWriteText ClosePath = Just "h" + + +instance IpeWriteText r => IpeWriteText (PolyLine 2 () r) where + ipeWriteText pl = case pl^..points.Tr.traverse.core of + (p : rest) -> unlines' . map ipeWriteText $ MoveTo p : map LineTo rest + -- the polyline type guarantees that there is at least one point + + +instance IpeWriteText r => IpeWriteText (PathSegment r) where + ipeWriteText (PolyLineSegment p) = ipeWriteText p + ipeWriteText (EllipseSegment m) = ipeWriteText $ Ellipse m + +instance IpeWriteText (PathAttrElf rs r) => IpeWriteText (PathAttribute r rs) where + ipeWriteText (PathAttribute x) = ipeWriteText x + +instance IpeWriteText r => IpeWrite (Path r) where + ipeWrite (Path segs) = (\t -> Element "path" [] [Text t]) <$> mt + where + concat' = F.foldr1 (\t t' -> t <> "\n" <> t') + mt = fmap concat' . Tr.sequence . fmap ipeWriteText $ segs + +-------------------------------------------------------------------------------- + + +instance ( IpeObjectElF r fld ~ (g :+ Rec f ats) + , IpeWrite g + , RecAll f ats IpeWriteText, AllSatisfy IpeAttrName ats + ) => IpeWrite (IpeObject r fld) where + ipeWrite (IpeObject (g :+ ats)) = ipeWrite g `mAddAtts` ipeWriteAttrs ats + + +ipeWriteRec :: RecAll f rs IpeWrite => Rec f rs -> [Node Text Text] +ipeWriteRec = catMaybes . recordToList + . rmap (\(Compose (Dict x)) -> Const $ ipeWrite x) + . reifyConstraint (Proxy :: Proxy IpeWrite) + +instance RecAll (IpeObject r) gt IpeWrite => IpeWrite (Group gt r) where + -- basically the same implementation as ipeWriteAttrs: convert the rec to a Rec Const + -- then turn that into a list. If the list is non-empty we construct a new group element. + ipeWrite = fmap (Element "group" []) + . wrap . ipeWriteRec + where + wrap [] = Nothing + wrap xs = Just xs + +instance IpeWriteText (GroupAttrElf rs r) => IpeWriteText (GroupAttribute r rs) where + ipeWriteText (GroupAttribute x) = ipeWriteText x + + +-------------------------------------------------------------------------------- + +instance IpeWrite IT.Layer where + ipeWrite (IT.Layer l) = Just $ Element "layer" [("name", l)] [] + +instance IpeWrite View where + ipeWrite (View lrs act) = Just $ Element "view" [ ("layers", ls) + , ("active", _layerName act) + ] [] + where + ls = T.unwords . map _layerName $ lrs + +instance IpeWrite (Group gs r) => IpeWrite (IpePage gs r) where + ipeWrite (IpePage lrs vs pgs) = Just . Element "page" [] . catMaybes . concat $ + [ map ipeWrite lrs + , map ipeWrite vs + , [ipeWrite pgs] + ] + +instance RecAll (Page r) gs IpeWrite => IpeWrite (IpeFile gs r) where + ipeWrite (IpeFile p s pgs) = Just $ Element "ipe" ipeAtts chs + where + ipeAtts = [("version","70005"),("creator", "HGeometry")] + -- TODO: Add preamble and styles + chs = ipeWriteRec pgs + + +-------------------------------------------------------------------------------- + +type Atts = [(Text,Text)] + +ipeWritePolyLines :: IpeWriteText r + => [(PolyLine 2 () r, Atts)] -> Node Text Text +ipeWritePolyLines pls = Element "ipe" ipeAtts [Element "page" [] chs] + where + chs = layers pls ++ mapMaybe f pls + ipeAtts = [("version","70005"),("creator", "HGeometry 0.4.0.0")] + + f (pl,ats) = ipeWrite (mkPath pl) `mAddAtts` ats + mkPath = Path . S2.l1Singleton . PolyLineSegment + layers = map mkLayer . nub . mapMaybe (lookup "layer" . snd) + mkLayer n = Element "layer" [("name",n)] [] + + +writePolyLineFile :: IpeWriteText r => FilePath -> [(PolyLine 2 () r, Atts)] -> IO () +writePolyLineFile fp = B.writeFile fp . format' . ipeWritePolyLines + + +instance (IpeWriteText r, IpeWrite p) => IpeWrite (PolyLine 2 p r) where + ipeWrite p = ipeWrite path + where + path = fromPolyLine $ p & points.Tr.traverse.extra .~ () + -- TODO: Do something with the p's + +fromPolyLine = Path . S2.l1Singleton . PolyLineSegment + + +instance (IpeWriteText r) => IpeWrite (LineSegment 2 p r) where + ipeWrite (LineSegment p q) = ipeWrite . fromPolyLine . fromPoints . map (^.core) $ [p,q] + + +instance IpeWrite () where + ipeWrite = const Nothing + +-- -- | slightly clever instance that produces a group if there is more than one +-- -- element and just an element if there is only one value produced +-- instance IpeWrite a => IpeWrite [a] where +-- ipeWrite = combine . mapMaybe ipeWrite + + +combine :: [Node Text Text] -> Maybe (Node Text Text) +combine [] = Nothing +combine [n] = Just n +combine ns = Just $ Element "group" [] ns + +-- instance (IpeWrite a, IpeWrite b) => IpeWrite (a,b) where +-- ipeWrite (a,b) = combine . catMaybes $ [ipeWrite a, ipeWrite b] + + + +-- -- | The default symbol for a point +-- ipeWritePoint :: IpeWriteText r => Point 2 r -> Maybe (Node Text Text) +-- ipeWritePoint = ipeWrite . flip Symbol "mark/disk(sx)" + + +-- instance (IpeWriteText r, Floating r) => IpeWrite (Circle r) where +-- ipeWrite = ipeWrite . Path . S2.l1Singleton . fromCircle + + + +-------------------------------------------------------------------------------- + + + +testPoly :: PolyLine 2 () Double +testPoly = fromPoints [origin, point2 0 10, point2 10 10, point2 100 100] + + + + +testWriteUse :: Maybe (Node Text Text) +testWriteUse = ipeWriteExt sym + where + sym :: IpeSymbol Double :+ (Rec (SymbolAttribute Double) [Size, SymbolStroke]) + sym = Symbol origin "mark" :+ ( SymbolAttribute (IpeSize $ Named "normal") + :& SymbolAttribute (IpeColor $ Named "green") + :& RNil + ) + + + +foo = ipeWrite grrr diff --git a/src/Data/Geometry/Line.hs b/src/Data/Geometry/Line.hs --- a/src/Data/Geometry/Line.hs +++ b/src/Data/Geometry/Line.hs @@ -1,133 +1,15 @@ -{-# Language - TypeFamilies - #-} -module Data.Geometry.Line where - -import Prelude hiding(length) - -import Data.Geometry.Point -import Data.Geometry.Geometry - - -import qualified Data.List as L - ---------------------------------------------------------------------- --- | A simple line segment in 2D consisint of a start and an end-point -data LineSegment2' a = LineSegment2 { startPoint :: Point2' a - , endPoint :: Point2' a - } - deriving (Eq, Ord, Show, Read) - -instance Functor LineSegment2' where - fmap f (LineSegment2 p q) = LineSegment2 (fmap f p) (fmap f q) - -instance IsPoint2Functor LineSegment2' where - p2fmap f (LineSegment2 p q) = LineSegment2 (f p) (f q) - -instance HasPoints LineSegment2' where - points (LineSegment2 p q) = [p,q] - ---------------------------------------------------------------------- --- | An infinite line - -newtype Line2' a = Line2 (LineSegment2' a) - deriving (Eq,Ord,Show,Read) - - -instance Functor Line2' where - fmap f (Line2 l) = Line2 $ fmap f l - -instance IsPoint2Functor Line2' where - p2fmap f (Line2 l) = Line2 $ p2fmap f l - -instance HasPoints Line2' where - points (Line2 l) = points l - ---------------------------------------------------------------------- --- | Polylines -newtype Polyline2' a = Polyline2 [LineSegment2' a] - deriving (Eq, Show, Read) - -instance IsPoint2Functor Polyline2' where - p2fmap f (Polyline2 ls) = Polyline2 (map (p2fmap f) ls) - -instance HasPoints Polyline2' where - points (Polyline2 ls) = case ls of - [] -> [] - (l:ls') -> points l ++ map endPoint ls' - ---------------------------------------------------------------------- --- | Constructing polylines - -polyLine :: [Point2' a] -> Polyline2' a -polyLine = Polyline2 . makeLines - where - makeLines :: [Point2' a] -> [LineSegment2' a] - makeLines [] = - error "Polyline consists of at least two points. No points given." - makeLines [_] = - error "Polyline consists of at least two points. Only one point given." - makeLines pts = zipWith LineSegment2 pts (tail pts) - ---------------------------------------------------------------------- --- | functions on Linesegments and Polylines - -isSimpleLine :: Polyline2' a -> Bool -isSimpleLine (Polyline2 []) = error "polyline without line segments" -isSimpleLine (Polyline2 [_]) = True -isSimpleLine _ = False - -toSimpleLine :: Polyline2' a -> LineSegment2' a -toSimpleLine (Polyline2 ls) = head ls +module Data.Geometry.Line(module I + ) where -toSimpleLineOption :: Polyline2' a -> Maybe (LineSegment2' a) -toSimpleLineOption p = if isSimpleLine p then Just (toSimpleLine p) else Nothing +import Data.Geometry.Line.Internal as I +import Data.Geometry.LineSegment +import Data.Geometry.Transformation +import Data.Geometry.Vector ---------------------------------------------------------------------- --- | Linear interpolation / points on line segments etc. - --- | simple linear interpolation, assuming t in [0,1] -linear :: Num a => a -> a -> a -> a -linear t x y = (1-t)*x + t*y - -inRange :: Ord a => a -> (a,a) -> Bool -x `inRange` (a,b) = a <= x && x <= b - -onLineSegment :: (Ord a, Fractional a) => Point2' a -> LineSegment2' a -> Bool -p `onLineSegment` l@(LineSegment2 s t) = - if t == s then p == s else (lambda `inRange` (0,1) && p == pointAt lambda l) +-- | Lines are transformable, via line segments +instance (Num r, AlwaysTruePFT d) => IsTransformable (Line d r) where + transformBy t = supportingLine . transformPointFunctor t . toLineSegment' where - a = p |-| s -- the vector from s to p - b = t |-| s -- the vector from s to t - lambda = (a |@| b) / (len b) - -- we translate such that s corresponds with the origin. In this coord system - -- b represents the input line segment. - -- We orthoganally project a onto b. Let c be this point (on the vector b) - -- then : d = a |@| b / length b denotes the distance between (0,0) and c - -- We can now get the lambda such that : c = linear (0,0) b by dividing - -- d / length b. Hence in total we divide through (length b)^2. This means - -- we can avoid computing the square root. - len (Point2 (x,y)) = x^2 + y^2 - -class HasLength c where - type PM c -- the precision model - -- | The length of the line-like segment - length :: c -> PM c - -instance Floating a => HasLength (LineSegment2' a) where - type PM (LineSegment2' a) = a - length (LineSegment2 s t) = dist s t - -instance Floating a => HasLength (Polyline2' a) where - type PM (Polyline2' a) = a - length (Polyline2 ls) = sum . map length $ ls - - -class LineLike c where - -- | get the point at `time' t (t in [0,1]) - pointAt :: Num a => a -> c a -> Point2' a - -instance LineLike LineSegment2' where - pointAt t (LineSegment2 (Point2 (px,py)) (Point2 (qx,qy))) = - Point2 (linear t px qx, linear t py qy) + toLineSegment' :: (Num r, Arity d) => Line d r -> LineSegment d () r + toLineSegment' = toLineSegment diff --git a/src/Data/Geometry/Line/Internal.hs b/src/Data/Geometry/Line/Internal.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Line/Internal.hs @@ -0,0 +1,131 @@ +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE DeriveFunctor #-} +module Data.Geometry.Line.Internal where + +import Control.Applicative +import Control.Lens +import Data.Ext +import qualified Data.Foldable as F +import Data.Geometry.Box +import Data.Geometry.Interval +import Data.Geometry.Point +import Data.Geometry.Properties +import Data.Geometry.Vector +import Linear.Affine(Affine(..)) +import Linear.Vector((*^)) + + +-------------------------------------------------------------------------------- +-- * d-dimensional Lines + +-- | A line is given by an anchor point and a vector indicating the +-- direction. +data Line d r = Line { _anchorPoint :: Point d r + , _direction :: Vector d r + } +makeLenses ''Line + +deriving instance (Show r, Arity d) => Show (Line d r) +deriving instance Arity d => Functor (Line d) + +type instance Dimension (Line d r) = d +type instance NumType (Line d r) = r + +-- ** Functions on lines + +-- | A line may be constructed from two points. +lineThrough :: (Num r, Arity d) => Point d r -> Point d r -> Line d r +lineThrough p q = Line p (q .-. p) + +verticalLine :: Num r => r -> Line 2 r +verticalLine x = Line (point2 x 0) (v2 0 1) + +horizontalLine :: Num r => r -> Line 2 r +horizontalLine y = Line (point2 0 y) (v2 1 0) + +perpendicularTo :: Num r => Line 2 r -> Line 2 r +perpendicularTo (Line p (Vector2 vx vy)) = Line p (v2 (-vy) vx) + + + + + +-- | Test if two lines are identical, meaning; if they have exactly the same +-- anchor point and directional vector. +isIdenticalTo :: (Eq r, Arity d) => Line d r -> Line d r -> Bool +(Line p u) `isIdenticalTo` (Line q v) = (p,u) == (q,v) + + +-- | Test if the two lines are parallel. +-- +-- >>> 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) +-- False +isParallelTo :: (Eq r, Fractional r, Arity d) + => Line d r -> Line d r -> Bool +(Line _ u) `isParallelTo` (Line _ v) = u `isScalarMultipleOf` v + -- TODO: Maybe use a specialize pragma for 2D (see intersect instance for two lines.) + + +-- | Test if point p lies on line l +-- +-- >>> origin `onLine` lineThrough origin (point2 1 0) +-- True +-- >>> point2 10 10 `onLine` lineThrough origin (point2 2 2) +-- True +-- >>> 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 + -- TODO: Maybe use a specialize pragma for 2D with an implementation using ccw + + + +instance (Eq r, Fractional r) => (Line 2 r) `IsIntersectableWith` (Line 2 r) where + + data Intersection (Line 2 r) (Line 2 r) = SameLine !(Line 2 r) + | LineLineIntersection !(Point 2 r) + | ParallelLines -- ^ No intersection + deriving (Show) + + nonEmptyIntersection ParallelLines = False + nonEmptyIntersection _ = True + + l@(Line p (Vector2 ux uy)) `intersect` m@(Line q v@(Vector2 vx vy)) + | areParallel = if q `onLine` l then SameLine l else ParallelLines + | otherwise = LineLineIntersection r + where + r = q .+^ alpha *^ v + + denom = vy * ux - vx * uy + areParallel = denom == 0 + -- Instead of using areParallel, we can also use the generic 'isParallelTo' function + -- for lines of arbitrary dimension, but this is a bit more efficient. + + alpha = (ux * (py - qy) + uy * (qx - px)) / denom + + Point2 px py = p + Point2 qx qy = q + +instance (Eq r, Fractional r) => Eq (Intersection (Line 2 r) (Line 2 r)) where + (SameLine l) == (SameLine m) = case (l `intersect` m) of + SameLine _ -> True + _ -> False + (LineLineIntersection p) == (LineLineIntersection q) = p == q + ParallelLines == ParallelLines = True + _ == _ = False + + + +-------------------------------------------------------------------------------- +-- * Supporting Lines + +-- | Types for which we can compute a supporting line, i.e. a line that contains the thing of type t. +class HasSupportingLine t where + supportingLine :: t -> Line (Dimension t) (NumType t) + +instance HasSupportingLine (Line d r) where + supportingLine = id diff --git a/src/Data/Geometry/LineSegment.hs b/src/Data/Geometry/LineSegment.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/LineSegment.hs @@ -0,0 +1,196 @@ +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE DeriveFunctor #-} +module Data.Geometry.LineSegment where + +import Control.Applicative +import Control.Lens +import Data.Ext +import Data.Geometry.Box +import Data.Geometry.Interval +import Data.Geometry.Point +import Data.Geometry.Line.Internal +import Data.Geometry.Properties +import Data.Geometry.Transformation +import Data.Geometry.Vector +import Linear.Vector((*^)) +import Linear.Affine(Affine(..),distanceA) +import qualified Data.List as L +import Data.Maybe(maybe) +import Data.Ord(comparing) +import Data.Semigroup + +-------------------------------------------------------------------------------- +-- * d-dimensional LineSegments + +-- | Line segments. LineSegments have a start and end point, both of which may +-- contain additional data of type p. +newtype LineSegment d p r = LineSeg { _unLineSeg :: Interval p (Point d r) } +pattern LineSegment s t = LineSeg (Interval s t) + +instance HasStart (LineSegment d p r) where + type StartCore (LineSegment d p r) = Point d r + type StartExtra (LineSegment d p r) = p + start = lens (_start . _unLineSeg) (\(LineSegment _ t) s -> LineSegment s t) + +instance HasEnd (LineSegment d p r) where + type EndCore (LineSegment d p r) = Point d r + type EndExtra (LineSegment d p r) = p + end = lens (_end . _unLineSeg) (\(LineSegment s _) t -> LineSegment s t) + +instance (Num r, Arity d) => HasSupportingLine (LineSegment d p r) where + supportingLine (LineSegment (p :+ _) (q :+ _)) = lineThrough p q + +deriving instance (Show r, Show p, Arity d) => Show (LineSegment d p r) +deriving instance (Eq r, Eq p, Arity d) => Eq (LineSegment d p r) +deriving instance (Ord r, Ord p, Arity d) => Ord (LineSegment d p r) +deriving instance Arity d => Functor (LineSegment d p) +type instance Dimension (LineSegment d p r) = d +type instance NumType (LineSegment d p r) = r + +instance PointFunctor (LineSegment d p) where + pmap f (LineSegment s e) = LineSegment (f <$> s) (f <$> e) + +instance Arity d => IsBoxable (LineSegment d p r) where + boundingBox l = boundingBoxList [l^.start.core, l^.end.core] + +instance (Num r, AlwaysTruePFT d) => IsTransformable (LineSegment d p r) where + transformBy = transformPointFunctor + + + + + + + + +-- ** Converting between Lines and LineSegments + +toLineSegment :: (Monoid p, Num r, Arity d) => Line d r -> LineSegment d p r +toLineSegment (Line p v) = LineSegment (p :+ mempty) + (p .+^ v :+ mempty) + +-- *** Intersecting LineSegments + +instance (Ord r, Fractional r) => + (LineSegment 2 p r) `IsIntersectableWith` (LineSegment 2 p r) where + + data Intersection (LineSegment 2 p r) (LineSegment 2 p r) = + OverlappingSegment !(LineSegment 2 p r) + | LineSegLineSegIntersection !(Point 2 r) + | NoIntersection + deriving (Show,Eq) + + nonEmptyIntersection NoIntersection = False + nonEmptyIntersection _ = True + + a@(LineSegment p q) `intersect` b@(LineSegment s t) = case la `intersect` lb of + SameLine _ -> + maybe NoIntersection OverlappingSegment $ overlap a b + LineLineIntersection r | onBothSegments r -> LineSegLineSegIntersection r + _ -> NoIntersection + where + la = supportingLine a + lb = supportingLine b + onBothSegments r = onSegment r a && onSegment r b + +instance (Ord r, Fractional r) => + (LineSegment 2 p r) `IsIntersectableWith` (Line 2 r) where + data Intersection (LineSegment 2 p r) (Line 2 r) = + LineContainsSegment !(LineSegment 2 p r) + | LineLineSegmentIntersection !(Point 2 r) + | NoLineLineSegmentIntersection + deriving (Show,Eq) + + nonEmptyIntersection NoLineLineSegmentIntersection = False + nonEmptyIntersection _ = True + + s `intersect` l = case (supportingLine s) `intersect` l of + SameLine _ -> LineContainsSegment s + LineLineIntersection p | p `onSegment` s -> LineLineSegmentIntersection p + _ -> NoLineLineSegmentIntersection + + +-- * Functions on LineSegments + +-- | Test if a point lies on a line segment. +-- +-- >>> (point2 1 0) `onSegment` (LineSegment (origin :+ ()) (point2 2 0 :+ ())) +-- True +-- >>> (point2 1 1) `onSegment` (LineSegment (origin :+ ()) (point2 2 0 :+ ())) +-- False +-- >>> (point2 5 0) `onSegment` (LineSegment (origin :+ ()) (point2 2 0 :+ ())) +-- False +-- >>> (point2 (-1) 0) `onSegment` (LineSegment (origin :+ ()) (point2 2 0 :+ ())) +-- False +-- >>> (point2 1 1) `onSegment` (LineSegment (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` (LineSegment (origin :+ ()) (point2 2 0 :+ ())) +-- True +-- >>> origin `onSegment` (LineSegment (origin :+ ()) (point2 2 0 :+ ())) +-- True +-- +-- +-- This function works for arbitrary dimensons. +-- +-- >>> (point3 1 1 1) `onSegment` (LineSegment (origin :+ ()) (point3 3 3 3 :+ ())) +-- True +-- >>> (point3 1 2 1) `onSegment` (LineSegment (origin :+ ()) (point3 3 3 3 :+ ())) +-- False +onSegment :: (Ord r, Fractional r, Arity d) + => Point d r -> LineSegment d p r -> Bool +p `onSegment` l = let s = l^.start.core + t = l^.end.core + inRange x = 0 <= x && x <= 1 + in maybe False inRange $ scalarMultiple (p .-. s) (t .-. s) + + +-- | Compute the overlap between the two segments (if they overlap/intersect) +overlap :: (Ord r, Fractional r, Arity d) + => LineSegment d p r -> LineSegment d p r -> Maybe (LineSegment d p r) +overlap l m = mim >>= \im -> case il `intersect` im of + IntervalIntersection (Interval s e) -> Just $ LineSegment (s^.extra) (e^.extra) + NoOverlap -> Nothing + where + p = l^.start + q = l^.end + r = m^.start + s = m^.end + + u = q^.core .-. p^.core + + -- lineseg l corresp to an interval from 0 1 + il = Interval (0 :+ p) (1 :+ q) + + + -- let lambda x denote the scalar s.t. x = p + (lambda x) *^ u + -- + -- lambda' computes lambda' and pairs it with the associated point. + -- lambda :: (Point d r :+ extra) -> Maybe (r :+ (Point d r :+ extra) + lambda' x = (:+ x) <$> scalarMultiple u (x^.core .-. p^.core) + + -- lineseg m corresponds to an interval + -- [min (lambda r, lambda s), max (lambda r, lambda s)] + -- + -- mim denotes this interval, assuming it exists (i.e. that is, + -- assuming r and s are indeed colinear with pq. + mim = mapM lambda' [r, s] >>= (f . L.sortBy (comparing (^.core))) + + -- Make sure we have two elems, s.t. both r and s are colinear with pq + f [a,b] = Just $ Interval a b + f _ = Nothing + + +-- | The left and right end point (or left below right if they have equal x-coords) +orderedEndPoints :: Ord r => LineSegment 2 p r -> (Point 2 r :+ p, Point 2 r :+ p) +orderedEndPoints s = if pc <= qc then (p, q) else (q,p) + where + p@(pc :+ _) = s^.start + q@(qc :+ _) = s^.end + + +-- | Length of the line segment +segmentLength :: (Arity d, Floating r) => LineSegment d p r -> r +segmentLength (LineSegment p q) = distanceA (p^.core) (q^.core) 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 @@ -1,40 +1,230 @@ +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE AutoDeriveTypeable #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE DeriveFunctor #-} module Data.Geometry.Point where +import Data.Proxy -data Point2' a = Point2 (a,a) - deriving (Show,Read,Eq,Ord) +import Control.Applicative +import Control.Lens +-- import Data.TypeLevel.Common +import Data.Typeable +import Data.Vinyl hiding (Nat) +import qualified Data.Vinyl as V +import Data.Vinyl.Functor(Const(..)) +import qualified Data.Vinyl.TypeLevel as TV -instance Functor Point2' where - fmap f (Point2 (x,y)) = Point2 (f x,f y) +import Data.Vinyl.TypeLevel hiding (Nat) -(|+|) :: Num a => Point2' a -> Point2' a -> Point2' a -(Point2 (x,y)) |+| (Point2 (a,b)) = Point2 (x+a,y+b) +import qualified Data.Vector.Fixed as FV +-- import qualified Data.Vector.Fixed.Cont as C -(|-|) :: Num a => Point2' a -> Point2' a -> Point2' a -(Point2 (x,y)) |-| (Point2 (a,b)) = Point2 (x-a,y-b) +import Data.Geometry.Properties +import Data.Geometry.Vector +import GHC.TypeLits --- | scalar multiplication -(|*|) :: Num a => a -> Point2' a -> Point2' a -s |*| (Point2 (x,y)) = Point2 (s*x,s*y) +import qualified Data.Geometry.Vector as Vec +import Linear.Vector(Additive(..)) +import Linear.Affine hiding (Point(..), origin) --- | dot product -(|@|) :: Num a => Point2' a -> Point2' a -> a -(Point2 (x,y)) |@| (Point2 (a,b)) = x*a + y*b +-------------------------------------------------------------------------------- +-------------------------------------------------------------------------------- +-- $setup +-- >>> :{ +-- let myVector :: Vector 3 Int +-- myVector = v3 1 2 3 +-- myPoint = Point myVector +-- :} -getX :: Point2' a -> a -getX (Point2 (x,_)) = x -getY :: Point2' a -> a -getY (Point2 (_,y)) = y +-------------------------------------------------------------------------------- +-- * A d-dimensional Point --- | euclidean distance between p and q -dist :: Floating a => Point2' a -> Point2' a -> a -dist p q = sqrt $ l22dist p q +-- | A d-dimensional point. +newtype Point d r = Point { toVec :: Vector d r } --- | Squared euclidean distance between p and q -l22dist :: Num a => Point2' a -> Point2' a -> a -l22dist p q = let (Point2 (a,b)) = q |-| p in a*a + b*b +deriving instance (Show r, Arity d) => Show (Point d r) +deriving instance (Eq r, Arity d) => Eq (Point d r) +deriving instance (Ord r, Arity d) => Ord (Point d r) +deriving instance Arity d => Functor (Point d) + + +type instance NumType (Point d r) = r +type instance Dimension (Point d r) = d + + + +instance Arity d => Affine (Point d) where + type Diff (Point d) = Vector d + + p .-. q = toVec p ^-^ toVec q + p .+^ v = Point $ toVec p ^+^ v + + +-- | Point representing the origin in d dimensions +-- +-- >>> origin :: Point 4 Int +-- Point {toVec = Vector {_unV = fromList [0,0,0,0]}} +origin :: (Arity d, Num r) => Point d r +origin = Point $ pure 0 + + +-- ** Accessing points + +-- | Lens to access the vector corresponding to this point. +-- +-- >>> (point3 1 2 3) ^. vector +-- Vector {_unV = fromList [1,2,3]} +-- >>> origin & vector .~ v3 1 2 3 +-- Point {toVec = Vector {_unV = fromList [1,2,3]}} +vector :: Lens' (Point d r) (Vector d r) +vector = lens toVec (const Point) + + +-- | Get the coordinate in a given dimension. This operation is unsafe in the +-- sense that no bounds are checked. Consider using `coord` instead. +-- +-- +-- >>> point3 1 2 3 ^. unsafeCoord 2 +-- 2 +unsafeCoord :: Arity d => Int -> Lens' (Point d r) r +unsafeCoord i = vector . FV.element (i-1) + -- Points are 1 indexed, vectors are 0 indexed + +-- | Get the coordinate in a given dimension +-- +-- >>> point3 1 2 3 ^. coord (C :: C 2) +-- 2 +-- >>> point3 1 2 3 & coord (C :: C 1) .~ 10 +-- Point {toVec = Vector {_unV = fromList [10,2,3]}} +-- >>> point3 1 2 3 & coord (C :: C 3) %~ (+1) +-- Point {toVec = Vector {_unV = fromList [1,2,4]}} +coord :: forall proxy i d r. (Index' (i-1) d, Arity d) => proxy i -> Lens' (Point d r) r +coord _ = vector . Vec.element (Proxy :: Proxy (i-1)) + + + + +-------------------------------------------------------------------------------- +-- * Convenience functions to construct 2 and 3 dimensional points + + +-- | We provide pattern synonyms Point2 and Point3 for 2 and 3 dimensional points. i.e. +-- we can write: +-- +-- >>> :{ +-- let +-- f :: Point 2 r -> r +-- f (Point2 x y) = x +-- in f (point2 1 2) +-- :} +-- 1 +-- +-- if we want. +pattern Point2 x y <- (_point2 -> (x,y)) + +-- | Similarly, we can write: +-- +-- >>> :{ +-- let +-- g :: Point 3 r -> r +-- g (Point3 x y z) = z +-- in g myPoint +-- :} +-- 3 +pattern Point3 x y z <- (_point3 -> (x,y,z)) + +-- | Construct a 2 dimensional point +-- +-- >>> point2 1 2 +-- Point {toVec = Vector {_unV = fromList [1,2]}} +point2 :: r -> r -> Point 2 r +point2 x y = Point $ v2 x y + +-- | Destruct a 2 dimensional point +-- +-- >>> _point2 $ point2 1 2 +-- (1,2) +_point2 :: Point 2 r -> (r,r) +_point2 p = (p^.xCoord, p^.yCoord) + + +-- | Construct a 3 dimensional point +-- +-- >>> point3 1 2 3 +-- Point {toVec = Vector {_unV = fromList [1,2,3]}} +point3 :: r -> r -> r -> Point 3 r +point3 x y z = Point $ v3 x y z + +-- | Destruct a 3 dimensional point +-- +-- >>> _point3 $ point3 1 2 3 +-- (1,2,3) +_point3 :: Point 3 r -> (r,r,r) +_point3 p = (p^.xCoord, p^.yCoord, p^.zCoord) + + +type i <=. d = (Index' (i-1) d, Arity d) + +-- | Shorthand to access the first coordinate C 1 +-- +-- >>> point3 1 2 3 ^. xCoord +-- 1 +-- >>> point2 1 2 & xCoord .~ 10 +-- Point {toVec = Vector {_unV = fromList [10,2]}} +xCoord :: (1 <=. d) => Lens' (Point d r) r +xCoord = coord (C :: C 1) + +-- | Shorthand to access the second coordinate C 2 +-- +-- >>> point2 1 2 ^. yCoord +-- 2 +-- >>> point3 1 2 3 & yCoord %~ (+1) +-- Point {toVec = Vector {_unV = fromList [1,3,3]}} +yCoord :: (2 <=. d) => Lens' (Point d r) r +yCoord = coord (C :: C 2) + +-- | Shorthand to access the third coordinate C 3 +-- +-- >>> point3 1 2 3 ^. zCoord +-- 3 +-- >>> point3 1 2 3 & zCoord %~ (+1) +-- Point {toVec = Vector {_unV = fromList [1,2,4]}} +zCoord :: (3 <=. d) => Lens' (Point d r) r +zCoord = coord (C :: C 3) + +-------------------------------------------------------------------------------- +-- * Point Functors + +-- | Types that we can transform by mapping a function on each point in the structure +class PointFunctor g where + pmap :: (Point (Dimension (g r)) r -> Point (Dimension (g s)) s) -> g r -> g s + + -- pemap :: (d ~ Dimension (g r)) => (Point d r :+ p -> Point d s :+ p) -> g r -> g s + -- pemap = + +instance PointFunctor (Point d) where + pmap f = f + + +-------------------------------------------------------------------------------- +-- * Functions specific to Two Dimensional points + +data CCW = CCW | CoLinear | CW + deriving (Show,Eq) + +-- | Given three points p q and r determine the orientation when going from p to r via q. +ccw :: (Ord r, Num r) => Point 2 r -> Point 2 r -> Point 2 r -> CCW +ccw p q r = case z `compare` 0 of + LT -> CW + GT -> CCW + EQ -> CoLinear + where + Vector2 ux uy = q .-. p + Vector2 vx vy = r .-. p + z = ux * vy - uy * vx diff --git a/src/Data/Geometry/PolyLine.hs b/src/Data/Geometry/PolyLine.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/PolyLine.hs @@ -0,0 +1,47 @@ +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE DeriveFunctor #-} +module Data.Geometry.PolyLine where + +import Control.Applicative +import Control.Lens +import Data.Ext +import qualified Data.Foldable as F +import Data.Geometry.Box +import Data.Geometry.Interval +import Data.Geometry.Point +import Data.Geometry.Properties +import Data.Geometry.Transformation +import Data.Geometry.Vector +import qualified Data.Seq2 as S2 +import Data.Semigroup + +-------------------------------------------------------------------------------- +-- * d-dimensional Polygonal Lines (PolyLines) + +-- | A Poly line in R^d +newtype PolyLine d p r = PolyLine { _points :: S2.Seq2 (Point d r :+ p) } +makeLenses ''PolyLine + +deriving instance (Show r, Show p, Arity d) => Show (PolyLine d p r) +deriving instance (Eq r, Eq p, Arity d) => Eq (PolyLine d p r) +deriving instance (Ord r, Ord p, Arity d) => Ord (PolyLine d p r) +deriving instance Arity d => Functor (PolyLine d p) +type instance Dimension (PolyLine d p r) = d +type instance NumType (PolyLine d p r) = r + +instance Semigroup (PolyLine d p r) where + (PolyLine pts) <> (PolyLine pts') = PolyLine $ pts <> pts' + +instance Arity d => IsBoxable (PolyLine d p r) where + boundingBox = boundingBoxList . toListOf (points.traverse.core) + +instance (Num r, AlwaysTruePFT d) => IsTransformable (PolyLine d p r) where + transformBy = transformPointFunctor + +instance PointFunctor (PolyLine d p) where + pmap f = over points (fmap (fmap f)) + + +-- | pre: The input list contains at least two points +fromPoints :: (Monoid p) => [Point d r] -> PolyLine d p r +fromPoints = PolyLine . S2.fromList . map (\p -> p :+ mempty) 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,60 +1,68 @@ +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE ScopedTypeVariables #-} module Data.Geometry.Polygon where -import Data.Geometry.BoundingBox -import Data.Geometry.Point -import Data.Geometry.Geometry ---------------------------------------------------------------------- --- | Polygons +import Control.Applicative +import Control.Lens hiding (Simple) +import qualified Data.Foldable as F --- | Class that defines what a polygon is. Note that it is assumed that the --- first and the last point are *NOT* the same point. --- -class (HasPoints p, IsTransformable p) => IsPolygon p where - vertices :: p a -> [Point2' a] - vertices = points +import Data.Ext +import Data.Geometry.Point +import Data.Geometry.Properties +import Data.Geometry.Transformation +import Data.Geometry.Box - -- | default implementation assumes points are in order - edges :: p a -> [(Point2' a, Point2' a)] - edges p = let pts = points p in - zip pts (tail pts ++ [head pts]) +import qualified Data.CircularList as C - isSimple :: p a -> Bool - containsHoles :: p a -> Bool +-------------------------------------------------------------------------------- +-- * Polygons ---------------------------------------------------------------------- --- | Simple polygons, i.e. polygons consisting of a sequence of points (vertices) --- | such that the edges do not intersect. Simple polygons do not contain holes -data SimplePolygon' a = SimplePolygon [Point2' a] - deriving (Show,Eq) +-- | We distinguish between simple polygons (without holes) and Polygons with holes. +data PolygonType = Simple | Multi -instance HasPoints SimplePolygon' where - points (SimplePolygon pts) = pts -instance IsPolygon SimplePolygon' where - isSimple = const True - containsHoles = const False +data Polygon (t :: PolygonType) p r where + SimplePolygon :: C.CList (Point 2 r :+ p) -> Polygon Simple p r + MultiPolygon :: C.CList (Point 2 r :+ p) -> [Polygon Simple p r] -> Polygon Multi p r -instance IsPoint2Functor SimplePolygon' where - p2fmap f (SimplePolygon pts) = SimplePolygon (map f pts) +type SimplePolygon = Polygon Simple ---------------------------------------------------------------------- --- | A multipolygon consists of several simple polygons -data MultiPolygon' a = MultiPolygon [SimplePolygon' a] - deriving (Show,Eq) +type MultiPolygon = Polygon Multi -instance HasPoints MultiPolygon' where - points (MultiPolygon pls) = concatMap points pls +-- | Polygons are per definition 2 dimensional +type instance Dimension (Polygon t p r) = 2 +type instance NumType (Polygon t p r) = r -instance IsPolygon MultiPolygon' where - isSimple = const False - containsHoles = const False +-- * Functions on Polygons -instance IsPoint2Functor MultiPolygon' where - p2fmap f (MultiPolygon polys) = MultiPolygon (map (p2fmap f) polys) +outerBoundary :: forall t p r. Lens' (Polygon t p r) (C.CList (Point 2 r :+ p)) +outerBoundary = lens get set + where + get :: Polygon t p r -> C.CList (Point 2 r :+ p) + get (SimplePolygon vs) = vs + get (MultiPolygon vs _) = vs ---------------------------------------------------------------------- --- | Bounding boxes can be used as polygons -instance IsPolygon BoundingBox2' where - isSimple = const True - containsHoles = const False + set :: Polygon t p r -> C.CList (Point 2 r :+ p) -> Polygon t p r + set (SimplePolygon _) vs = SimplePolygon vs + set (MultiPolygon _ hs) vs = MultiPolygon vs hs + +holes :: forall p r. Lens' (Polygon Multi p r) [Polygon Simple p r] +holes = lens get set + where + get :: Polygon Multi p r -> [Polygon Simple p r] + get (MultiPolygon _ hs) = hs + set :: Polygon Multi p r -> [Polygon Simple p r] -> Polygon Multi p r + set (MultiPolygon vs _) hs = MultiPolygon vs hs + + +-- | The vertices in the polygon. No guarantees are given on the order in which +-- they appear! +vertices :: Polygon t p r -> [Point 2 r :+ p] +vertices (SimplePolygon vs) = C.toList vs +vertices (MultiPolygon vs hs) = C.toList vs ++ concatMap vertices hs + + + +fromPoints :: [Point 2 r :+ p] -> SimplePolygon p r +fromPoints = SimplePolygon . C.fromList diff --git a/src/Data/Geometry/Properties.hs b/src/Data/Geometry/Properties.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Properties.hs @@ -0,0 +1,32 @@ +{-# LANGUAGE MultiParamTypeClasses #-} +module Data.Geometry.Properties where + +-------------------------------------------------------------------------------- +import GHC.TypeLits + +-- | A type family for types that are associated with a dimension. +type family Dimension t :: Nat + +-- | A type family for types that have an associated numeric type. +type family NumType t :: * + + +class IsIntersectableWith g h where + data Intersection g h + intersect :: g -> h -> Intersection g h + + -- | g `intersects` h <=> The intersection of g and h is non-empty. + -- + -- The default implementation computes the intersection of g and h, + -- and uses nonEmptyIntersection to determine if the intersection is + -- non-empty. + intersects :: g -> h -> Bool + g `intersects` h = nonEmptyIntersection $ g `intersect` h + + -- | Helper to implement `intersects`. + nonEmptyIntersection :: Intersection g h -> Bool + {-# MINIMAL intersect , nonEmptyIntersection #-} + +class IsUnionableWith g h where + data Union g h + union :: g -> h -> Union g h diff --git a/src/Data/Geometry/SetOperations.hs b/src/Data/Geometry/SetOperations.hs deleted file mode 100644 --- a/src/Data/Geometry/SetOperations.hs +++ /dev/null @@ -1,66 +0,0 @@ -{-# Language - MultiParamTypeClasses - , FlexibleInstances - , UndecidableInstances - #-} -module Data.Geometry.SetOperations( AreIntersectable(..) - ) where - -import Data.List -import Data.Geometry.Point -import Data.Geometry.Line -import Data.Geometry.Circle - -import Debug.Trace - --------------------------------------------------------------------------------- - --- | A class to represent that a pair of geometry objects (both parameterized --- over a) can be intersected. -class AreIntersectable g h a where - intersectionPoints :: g a -> h a -> [Point2' a] - --- | Intersection is symetrical --- instance AreIntersectable g h a => AreIntersectable h g a where --- intersectionPoints h g = intersectionPoints g h - - --- instance AreIntersectable LineSegment2' LineSegment2' a where --- intersectionPoints _ _ = [] - -instance (Ord a, Floating a) => AreIntersectable Circle2' LineSegment2' a where - -- | The intersection points, ordered along the line segment - intersectionPoints c l = filter (`onLineSegment` l) $ circleAndLine c l - -instance (Ord a, Floating a) => AreIntersectable Circle2' Polyline2' a where - -- | The intersection points, ordered along the polyline - intersectionPoints c (Polyline2 ls) = concatMap (intersectionPoints c) ls - -instance (Ord a, Floating a) => AreIntersectable Circle2' Line2' a where - intersectionPoints c (Line2 l) = circleAndLine c l - - --- | represents the intersection of a circle and an infinite line (as LineSegment ) -circleAndLine :: (Ord a, Floating a) => Circle2' a -> LineSegment2' a -> [Point2' a] -circleAndLine (Circle2 p r) (LineSegment2 s t) = map (|+| p) $ pts discr - where - s'@(Point2 (sx,sy)) = s |-| p -- translate so the circle is centred at (0,0) - t' = t |-| p - d@(Point2 (dx,dy)) = t' |-| s' -- vector from s' to t' - ---- - -- any point on the line segment is given as: q = s' + \lambda * d - q lamb = s' |+| (lamb |*| d) - -- solving the equation (q_x)^2 + (q_y)^2 = r^2 then yields the equation - -- L^2(dx^2 + dy^2) + L2(sx*dx + sy*dy) + sx^2 + sy^2 = 0 - -- where L = \lambda - a = dx^2 + dy^2 - b = 2*(sx*dx + sy*dy) - c = sx^2 + sy^2 - r^2 - discr = b^2 - 4*a*c - discr' = sqrt discr - lambda (|+-|) = (-b |+-| discr') / (2*a) - pts dscr - | dscr < 0 = [] -- no intersections - | dscr == 0 = [q $ lambda (+)] -- line tangent - | otherwise = let lambdas = sort [lambda (-), lambda (+)] in - map q lambdas -- intersection diff --git a/src/Data/Geometry/Transformation.hs b/src/Data/Geometry/Transformation.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Transformation.hs @@ -0,0 +1,140 @@ +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE DeriveFunctor #-} +module Data.Geometry.Transformation where + + +import Control.Applicative +import Control.Lens(lens,Lens',set) +import Data.Geometry.Point +import Data.Geometry.Properties +import Data.Geometry.Vector +import Data.Proxy + +import GHC.TypeLits +import Linear.Matrix((!*),(!*!)) +import Linear.Vector(zero) + +import qualified Data.Vector.Fixed as FV +import qualified Data.Geometry.Vector as V + +import Data.Vinyl.TypeLevel hiding (Nat) + +-------------------------------------------------------------------------------- +-- * Matrices + +-- | a matrix of n rows, each of m columns, storing values of type r +newtype Matrix n m r = Matrix (Vector n (Vector m r)) + +deriving instance (Show r, Arity n, Arity m) => Show (Matrix n m r) +deriving instance (Eq r, Arity n, Arity m) => Eq (Matrix n m r) +deriving instance (Ord r, Arity n, Arity m) => Ord (Matrix n m r) +deriving instance (Arity n, Arity m) => Functor (Matrix n m) + +multM :: (Arity r, Arity c, Arity c', Num a) => Matrix r c a -> Matrix c c' a -> Matrix r c' a +(Matrix a) `multM` (Matrix b) = Matrix $ a !*! b + +mult :: (Arity m, Arity n, Num r) => Matrix n m r -> Vector m r -> Vector n r +(Matrix m) `mult` v = m !* v + +-------------------------------------------------------------------------------- +-- * Transformations + +-- | A type representing a Transformation for d dimensional objects +newtype Transformation d r = Transformation { _transformationMatrix :: Matrix (1 + d) (1 + d) r } + +transformationMatrix :: Lens' (Transformation d r) (Matrix (1 + d) (1 + d) r) +transformationMatrix = lens _transformationMatrix (const Transformation) + +deriving instance (Show r, Arity (1 + d)) => Show (Transformation d r) +deriving instance (Eq r, Arity (1 + d)) => Eq (Transformation d r) +deriving instance (Ord r, Arity (1 + d)) => Ord (Transformation d r) +deriving instance Arity (1 + d) => Functor (Transformation d) + +type instance NumType (Transformation d r) = r + + +-- | Compose transformations (right to left) +(|.|) :: (Num r, Arity (1 + d)) => Transformation d r -> Transformation d r -> Transformation d r +(Transformation f) |.| (Transformation g) = Transformation $ f `multM` g + +-------------------------------------------------------------------------------- +-- * Transformable geometry objects + +-- | A class representing types that can be transformed using a transformation +class IsTransformable g where + transformBy :: Transformation (Dimension g) (NumType g) -> g -> g + +transformAllBy :: (Functor c, IsTransformable g) + => Transformation (Dimension g) (NumType g) -> c g -> c g +transformAllBy t = fmap (transformBy t) + + +type AlwaysTruePFT d = AlwaysTrueDestruct d (1 + d) + + +transformPointFunctor :: ( PointFunctor g, Num r, d ~ Dimension (g r) + , AlwaysTruePFT d + ) => Transformation d r -> g r -> g r +transformPointFunctor t = pmap (transformBy t) + +instance ( Num r + , Arity d, AlwaysTrueDestruct d (1 + d) + ) => IsTransformable (Point d r) where + transformBy (Transformation m) (Point v) = Point . V.init $ m `mult` v' + where + v' = snoc v 0 + + +-------------------------------------------------------------------------------- +-- * Common transformations + +translation :: ( Num r, Arity (1 + d) + , AlwaysTrueSnoc d, Arity d, Index' (1+d-1) (1+d)) + => Vector d r -> Transformation d r +translation v = Transformation . Matrix $ V.imap transRow (snoc v 1) + + +scaling :: (Num r, Arity (1 + d), AlwaysTrueSnoc d, Arity d) => Vector d r -> Transformation d r +scaling v = Transformation . Matrix $ V.imap mkRow (snoc v 1) + +uniformScaling :: (Num r, Arity (1 + d), AlwaysTrueSnoc d, Arity d) => r -> Transformation d r +uniformScaling = scaling . pure + +-------------------------------------------------------------------------------- +-- * 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) + , AlwaysTrueTransformation (Dimension g) + ) => Vector (Dimension g) (NumType g) -> g -> g +translateBy = transformBy . translation + +scaleBy :: ( IsTransformable g, Num (NumType g) + , AlwaysTrueTransformation (Dimension g) + ) => Vector (Dimension g) (NumType g) -> g -> g +scaleBy = transformBy . scaling + + +scaleUniformlyBy :: ( IsTransformable g, Num (NumType g) + , AlwaysTrueTransformation (Dimension g) + ) => NumType g -> g -> g +scaleUniformlyBy = transformBy . uniformScaling + + +-------------------------------------------------------------------------------- +-- * Helper functions to easily create matrices + +-- | Creates a row with zeroes everywhere, except at position i, where the +-- value is the supplied value. +mkRow :: forall d r. (Arity d, Num r) => Int -> r -> Vector d r +mkRow i x = set (FV.element i) x zero + +-- | Row in a translation matrix +transRow :: forall n r. (Arity n, Index' (n-1) n, Num r) => Int -> r -> Vector n r +transRow i x = set (V.element (Proxy :: Proxy (n-1))) x $ mkRow i 1 diff --git a/src/Data/Geometry/Triangle.hs b/src/Data/Geometry/Triangle.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Triangle.hs @@ -0,0 +1,45 @@ +{-# LANGUAGE DeriveFunctor #-} +module Data.Geometry.Triangle where + +import Control.Lens +import Data.Ext +import Data.Geometry.Point +import Data.Geometry.Ball +import Data.Geometry.Properties +import Data.Geometry.Transformation + +data Triangle p r = Triangle (Point 2 r :+ p) + (Point 2 r :+ p) + (Point 2 r :+ p) + deriving (Show,Eq,Functor) + +type instance NumType (Triangle p r) = r +type instance Dimension (Triangle p r) = 2 + +instance PointFunctor (Triangle p) where + pmap f (Triangle p q r) = Triangle (p&core %~ f) (q&core %~ f) (r&core %~ f) + +instance Num r => IsTransformable (Triangle d r) where + transformBy = transformPointFunctor + + +-- | Compute the area of a triangle +area :: Fractional r => Triangle p r -> r +area t = doubleArea t / 2 + +-- | 2*the area of a triangle. +doubleArea :: Num r => Triangle p r -> r +doubleArea (Triangle a b c) = abs $ ax*by - ax*cy + + bx*cy - bx*ay + + cx*ay - cx*by + -- Based on determinant of a 3x3 matrix (shoelace formula) + where + Point2 ax ay = a^.core + Point2 bx by = b^.core + Point2 cx cy = c^.core + + +-- | get the inscribed circle. Returns Nothing if the triangle is degenerate, +-- i.e. if the points are colinear. +inscribedCircle :: (Eq r, Fractional r) => Triangle p r -> Maybe (Circle () r) +inscribedCircle (Triangle p q r) = circle (p^.core) (q^.core) (r^.core) diff --git a/src/Data/Geometry/Vector.hs b/src/Data/Geometry/Vector.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Vector.hs @@ -0,0 +1,48 @@ +module Data.Geometry.Vector( module GV + , module FV + , isScalarMultipleOf + , scalarMultiple + ) where + +import qualified Data.Vector.Fixed as FV +import qualified Data.Foldable as F +import Data.Geometry.Vector.VectorFixed as GV + + +-- | Test if v is a scalar multiple of u. +-- +-- >>> v2 1 1 `isScalarMultipleOf` v2 10 10 +-- True +-- >>> v2 1 1 `isScalarMultipleOf` v2 10 1 +-- False +-- >>> v2 1 1 `isScalarMultipleOf` v2 11.1 11.1 +-- True +-- >>> v2 1 1 `isScalarMultipleOf` v2 11.1 11.2 +-- False +-- >>> v2 2 1 `isScalarMultipleOf` v2 11.1 11.2 +-- False +-- >>> v2 2 1 `isScalarMultipleOf` v2 4 2 +-- True +isScalarMultipleOf :: (Eq r, Fractional r, GV.Arity d) + => Vector d r -> Vector d r -> Bool +u `isScalarMultipleOf` v = fst $ scalarMultiple' u v + +-- | Get the scalar labmda s.t. v = lambda * u (if it exists) +scalarMultiple :: (Eq r, Fractional r, GV.Arity d) + => Vector d r -> Vector d r -> Maybe r +scalarMultiple u v = case scalarMultiple' u v of + (False,_) -> Nothing + (_, mm@(Just lambda)) -> mm + (_, _) -> Just . fromIntegral $ 0 + +-- | Helper function for computing the scalar multiple. The result is a pair +-- (b,mm), where b indicates if v is a scalar multiple of u, and mm is a Maybe +-- scalar multiple. If the result is Nothing, the scalar multiple is zero. +scalarMultiple' :: (Eq r, Fractional r, GV.Arity d) + => Vector d r -> Vector d r -> (Bool,Maybe r) +scalarMultiple' u v = F.foldr allLambda (True,Nothing) $ FV.zipWith f u v + where + f ui vi = (ui == 0 && vi == 0, ui / vi) + allLambda (True,_) x = x + allLambda (_, myLambda) (b,Nothing) = (b,Just myLambda) -- no lambda yet + allLambda (_, myLambda) (b,Just lambda) = (myLambda == lambda && b, Just lambda) diff --git a/src/Data/Geometry/Vector/VectorFixed.hs b/src/Data/Geometry/Vector/VectorFixed.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Vector/VectorFixed.hs @@ -0,0 +1,166 @@ +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE UndecidableInstances #-} +module Data.Geometry.Vector.VectorFixed where + +import Control.Applicative + +import Control.Lens + +import Data.Proxy + +import Data.Foldable +import Data.Traversable + +import Data.Vector.Fixed.Boxed +import Data.Vector.Fixed.Cont(Z(..),S(..),ToPeano(..),ToNat(..)) + +import GHC.TypeLits + +import Linear.Affine +import Linear.Metric +import Linear.Vector + +import qualified Data.Vector.Fixed as V + +import qualified Linear.V3 as L3 + + +-------------------------------------------------------------------------------- + +-- | A proxy which can be used for the coordinates. +data C (n :: Nat) = C deriving (Show,Read,Eq,Ord) + +-------------------------------------------------------------------------------- +-- * d dimensional Vectors + +-- | Datatype representing d dimensional vectors. Our implementation wraps the +-- implementation provided by fixed-vector. +newtype Vector (d :: Nat) (r :: *) = Vector { _unV :: Vec (ToPeano d) r } + + +unV :: Lens' (Vector d r) (Vec (ToPeano d) r) +unV = lens _unV (const Vector) + +---------------------------------------- +type Arity (n :: Nat) = V.Arity (ToPeano n) + +type Index' i d = V.Index (ToPeano i) (ToPeano d) + +element :: forall proxy i d r. (Arity d, Index' i d) => proxy i -> Lens' (Vector d r) r +element _ = V.elementTy (undefined :: (ToPeano i)) + + + +deriving instance (Show r, Arity d) => Show (Vector d r) +deriving instance (Eq r, Arity d) => Eq (Vector d r) +deriving instance (Ord r, Arity d) => Ord (Vector d r) +deriving instance Arity d => Functor (Vector d) + +deriving instance Arity d => Foldable (Vector d) +deriving instance Arity d => Applicative (Vector d) + +instance Arity d => Traversable (Vector d) where + traverse f (Vector v) = Vector <$> traverse f v + + +instance Arity d => Additive (Vector d) where + zero = pure 0 + (Vector u) ^+^ (Vector v) = Vector $ V.zipWith (+) u v + +instance Arity d => Affine (Vector d) where + type Diff (Vector d) = Vector d + + u .-. v = u ^-^ v + p .+^ v = p ^+^ v + + +instance Arity d => Metric (Vector d) + +type instance V.Dim (Vector d) = ToPeano d + +instance Arity d => V.Vector (Vector d) r where + construct = Vector <$> V.construct + inspect v = V.inspect (_unV v) + basicIndex v = V.basicIndex (_unV v) + +-- ---------------------------------------- + +type AlwaysTrueDestruct pd d = (Arity pd, ToPeano d ~ S (ToPeano pd)) + + +-- | Get the head and tail of a vector +destruct :: AlwaysTrueDestruct predD d + => Vector d r -> (r, Vector predD r) +destruct (Vector v) = (V.head v, Vector $ V.tail v) + + +-- | Cross product of two three-dimensional vectors +cross :: Num r => Vector 3 r -> Vector 3 r -> Vector 3 r +u `cross` v = fromV3 $ (toV3 u) `L3.cross` (toV3 v) + + +-------------------------------------------------------------------------------- + +-- | Conversion to a Linear.V3 +toV3 :: Vector 3 a -> L3.V3 a +toV3 (Vector3 a b c) = L3.V3 a b c + +-- | Conversion from a Linear.V3 +fromV3 :: L3.V3 a -> Vector 3 a +fromV3 (L3.V3 a b c) = v3 a b c + +---------------------------------------------------------------------------------- + + +type AlwaysTrueSnoc d = ToPeano (1 + d) ~ S (ToPeano d) + +-- | Add an element at the back of the vector +snoc :: (AlwaysTrueSnoc d, Arity d) => Vector d r -> r -> Vector (1 + d) r +snoc = flip V.snoc + +-- | Get a vector of the first d - 1 elements. +init :: AlwaysTrueDestruct predD d => Vector d r -> Vector predD r +init = Vector . V.reverse . V.tail . V.reverse . _unV + +-- | Get a prefix of i elements of a vector +prefix :: (Prefix (ToPeano i) (ToPeano d)) => Vector d r -> Vector i r +prefix (Vector v) = Vector $ prefix' v + +class Prefix i d where + prefix' :: Vec d r -> Vec i r + +instance Prefix Z d where + prefix' _ = V.vector V.empty + +instance (V.Arity i, V.Arity d, V.Index i d, Prefix i d) => Prefix (S i) (S d) where + prefix' v = V.vector $ V.head v `V.cons` (prefix' $ V.tail v) + + +-- | Map with indices +imap :: Arity d => (Int -> r -> s ) -> Vector d r -> Vector d s +imap = V.imap + +-------------------------------------------------------------------------------- +-- * Functions specific to two and three dimensional vectors. + +-- | Construct a 2 dimensional vector +v2 :: r -> r -> Vector 2 r +v2 a b = Vector $ V.mk2 a b + +-- | Construct a 3 dimensional vector +v3 :: r -> r -> r -> Vector 3 r +v3 a b c = Vector $ V.mk3 a b c + + +-- | Destruct a 2 dim vector into a pair +_unV2 :: Vector 2 r -> (r,r) +_unV2 v = let [x,y] = V.toList v in (x,y) + +_unV3 :: Vector 3 r -> (r,r,r) +_unV3 v = let [x,y,z] = V.toList v in (x,y,z) + + +-- | Pattern synonym for two and three dim vectors +pattern Vector2 x y <- (_unV2 -> (x,y)) +pattern Vector3 x y z <- (_unV3 -> (x,y,z)) diff --git a/src/Data/Seq2.hs b/src/Data/Seq2.hs new file mode 100644 --- /dev/null +++ b/src/Data/Seq2.hs @@ -0,0 +1,150 @@ +module Data.Seq2 where + +import Prelude hiding (foldr,foldl,head,tail,last) + +import Control.Applicative +import Data.List.NonEmpty +import Data.Semigroup + + +import qualified Data.Traversable as T +import qualified Data.Foldable as F +import qualified Data.Sequence as S + +-- | Basically Data.Sequence but with the guarantee that the list contains at +-- least two elements. +data Seq2 a = Seq2 a (S.Seq a) a + deriving (Eq,Ord,Show,Read) + + +instance T.Traversable Seq2 where + -- Applicative f => (a -> f b) -> t a -> f (t b) + traverse f ~(Seq2 l s r) = Seq2 <$> f l <*> T.traverse f s <*> f r + +instance Functor Seq2 where + fmap = T.fmapDefault + +instance F.Foldable Seq2 where + foldMap = T.foldMapDefault + +instance Semigroup (Seq2 a) where + l <> r = l >< r + +duo :: a -> a -> Seq2 a +duo a b = Seq2 a S.empty b + +length :: Seq2 a -> Int +length ~(Seq2 _ s _) = 2 + S.length s + + +-- | get the element with index i, counting from the left and starting at 0. +-- O(log(min(i,n-i))) +index :: Seq2 a -> Int -> a +index ~(Seq2 l s r) i + | i == 0 = l + | i < 1 + sz = S.index s (i+1) + | i == sz + 1 = r + | otherwise = error "index: index out of bounds." + where + sz = S.length s + + +(<|) :: a -> Seq2 a -> Seq2 a +x <| ~(Seq2 l s r) = Seq2 x (l S.<| s) r + + +(|>) :: Seq2 a -> a -> Seq2 a +~(Seq2 l s r) |> x = Seq2 l (s S.|> r) x + + +-- | Concatenate two sequences. O(log(min(n1,n2))) +~(Seq2 ll ls lr) >< ~(Seq2 rl rs rr) = Seq2 ll ((ls S.|> lr) S.>< (rl S.<| rs)) rr + + +-- | pre: the list contains at least two elements +fromList :: [a] -> Seq2 a +fromList (a:b:xs) = F.foldl' (\s x -> s |> x) (duo a b) xs +fromList _ = error "Seq2.fromList: Not enough values" + +-------------------------------------------------------------------------------- +-- | Left views + +data ViewL2 a = a :<< ViewR1 a deriving (Show,Read,Eq,Ord) + +-- | At least two elements +instance T.Traversable ViewL2 where + traverse f ~(a :<< s) = (:<<) <$> f a <*> T.traverse f s + +instance Functor ViewL2 where + fmap = T.fmapDefault + +instance F.Foldable ViewL2 where + foldMap = T.foldMapDefault + + +-- | At least one element +data ViewL1 a = a :< S.Seq a deriving (Show,Read,Eq,Ord) + +instance T.Traversable ViewL1 where + traverse f ~(a :< s) = (:<) <$> f a <*> T.traverse f s + +instance Functor ViewL1 where + fmap = T.fmapDefault + +instance F.Foldable ViewL1 where + foldMap = T.foldMapDefault + +-- | We throw away information here; namely that the combined list contains two elements. +instance Semigroup (ViewL1 a) where + ~(a :< s) <> ~(b :< t) = a :< (s <> S.singleton b <> t) + + +toNonEmpty :: ViewL1 a -> NonEmpty a +toNonEmpty ~(a :< s) = (a :| F.toList s) + + +-- | O(1) get a left view +viewl :: Seq2 a -> ViewL2 a +viewl ~(Seq2 l s r) = l :<< (s :> r) + + +l1Singleton :: a -> ViewL1 a +l1Singleton = (:< S.empty) + +-------------------------------------------------------------------------------- +-- | Right views + +-- | A view of the right end of the seq, with the guarantee that it +-- has at least two elements +data ViewR2 a = ViewL1 a :>> a deriving (Show,Read,Eq,Ord) + +instance T.Traversable ViewR2 where + traverse f ~(s :>> a) = (:>>) <$> T.traverse f s <*> f a + +instance Functor ViewR2 where + fmap = T.fmapDefault + +instance F.Foldable ViewR2 where + foldMap = T.foldMapDefault + + +-- | A view of the right end of the sequence, with the guarantee that it has at +-- least one element. +data ViewR1 a = S.Seq a :> a deriving (Show,Read,Eq,Ord) + +instance T.Traversable ViewR1 where + traverse f ~(s :> a) = (:>) <$> T.traverse f s <*> f a + +instance Functor ViewR1 where + fmap = T.fmapDefault + +instance F.Foldable ViewR1 where + foldMap = T.foldMapDefault + + +-- | O(1) get a right view +viewr :: Seq2 a -> ViewR2 a +viewr ~(Seq2 l s r) = (l :< s) :>> r + +r1Singleton :: a -> ViewR1 a +r1Singleton = (S.empty :>) diff --git a/src/System/Random/Shuffle.hs b/src/System/Random/Shuffle.hs new file mode 100644 --- /dev/null +++ b/src/System/Random/Shuffle.hs @@ -0,0 +1,28 @@ +module System.Random.Shuffle(shuffle) where + +import Control.Monad +import qualified Data.Foldable as F +import System.Random + +import qualified Data.Vector.Mutable as MV +import qualified Data.Vector as V + + + +-- | Fisher–Yates shuffle, which shuffles in O(n) time. +shuffle :: (F.Foldable f, RandomGen g) => g -> f a -> [a] +shuffle gen = V.toList . V.modify (\v -> + do + let n = MV.length v + forM_ (rands gen $ n - 1) $ \(SP i j) -> MV.swap v i j + ) . V.fromList . F.toList + +-- | Strict pair +data SP a b = SP !a !a + +-- | Generate a list of indices in decreasing order, coupled with a random +-- value in the range [0,i]. +rands :: RandomGen g => g -> Int -> [SP Int Int] +rands g i + | i <= 0 = [] + | otherwise = let (j,g') = randomR (0,i) g in SP i j : rands g' (i-1)