diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) 2012, Frank Staals + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Frank Staals nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/hgeometry.cabal b/hgeometry.cabal new file mode 100644 --- /dev/null +++ b/hgeometry.cabal @@ -0,0 +1,65 @@ +-- Initial hgeometry.cabal generated by cabal init. For further +-- documentation, see http://haskell.org/cabal/users-guide/ + +-- The name of the package. +name: hgeometry + +-- 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 + +version: 0.1.0.0 + +-- A short (one-line) description of the package. +synopsis: Geometry types in Haskell + +-- A longer description of the package. +description: Several basic geometry geometry types and functions on these types. + +-- The license under which the package is released. +license: BSD3 + +-- The file containing the license text. +license-file: LICENSE + +-- The package author(s). +author: Frank Staals + +-- An email address to which users can send suggestions, bug reports, and +-- patches. +maintainer: f.staals@uu.nl + +-- A copyright notice. +-- copyright: + +category: Geometry + +build-type: Simple + +-- Constraint on the version of Cabal needed to build this package. +cabal-version: >=1.8 + +library + -- Specify where we can find the source code + Hs-source-dirs: src + + GHC-Options: -Wall + + -- 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 + + + -- Modules included in this library but not exported. + -- other-modules: + + -- Other library packages from which modules are imported. + build-depends: base >=4.3 && < 5 diff --git a/src/Data/Geometry/BoundingBox.hs b/src/Data/Geometry/BoundingBox.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/BoundingBox.hs @@ -0,0 +1,90 @@ +{-# 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/Circle.hs b/src/Data/Geometry/Circle.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Circle.hs @@ -0,0 +1,89 @@ +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 + +-- | whether or not p lies in OR on the circle c +inCircle :: (Ord a, Floating a) => Point2' a -> Circle2' a -> Bool +p `inCircle` c = distanceToCenter p c <= radius c + +-- | whether or not p lies strictly inside the circle c +insideCircle :: (Floating a, Ord a) => Point2' a -> Circle2' a -> Bool +p `insideCircle` c = distanceToCenter p c < radius c + +-- | whether or not p lies on the circle +onCircle :: (Eq a, Floating a) => Point2' a -> Circle2' a -> Bool +p `onCircle` c = distanceToCenter p c == radius c + + +-- | whether or not a point lies in a disc: this includes its border +inDisc :: (Floating 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. +p `insideDisc` (Disc2 c) = p `insideCircle` c diff --git a/src/Data/Geometry/Geometry.hs b/src/Data/Geometry/Geometry.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Geometry.hs @@ -0,0 +1,123 @@ +{-# 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*z + 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/Line.hs b/src/Data/Geometry/Line.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Line.hs @@ -0,0 +1,122 @@ +{-# 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 + +toSimpleLineOption :: Polyline2' a -> Maybe (LineSegment2' a) +toSimpleLineOption p = if isSimpleLine p then Just (toSimpleLine p) else Nothing + + +--------------------------------------------------------------------- +-- | 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 + +onLineSegment :: (Eq a, Floating a) => Point2' a -> LineSegment2' a -> Bool +p `onLineSegment` l@(LineSegment2 s t) = p == pointAt lambda l + where + a = p |-| s + b = t |-| s + lambda = (a |@| b) / length l + +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) diff --git a/src/Data/Geometry/Point.hs b/src/Data/Geometry/Point.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Point.hs @@ -0,0 +1,37 @@ +module Data.Geometry.Point where + + +data Point2' a = Point2 (a,a) + deriving (Show,Read,Eq,Ord) + + +instance Functor Point2' where + fmap f (Point2 (x,y)) = Point2 (f x,f y) + +(|+|) :: Num a => Point2' a -> Point2' a -> Point2' a +(Point2 (x,y)) |+| (Point2 (a,b)) = Point2 (x+a,y+b) + +(|-|) :: Num a => Point2' a -> Point2' a -> Point2' a +(Point2 (x,y)) |-| (Point2 (a,b)) = Point2 (x-a,y-b) + + +-- | scalar multiplication +(|*|) :: Num a => a -> Point2' a -> Point2' a +s |*| (Point2 (x,y)) = Point2 (s*x,s*y) + + +-- | dot product +(|@|) :: Num a => Point2' a -> Point2' a -> a +(Point2 (x,y)) |@| (Point2 (a,b)) = x*a + y*b + + +getX :: Point2' a -> a +getX (Point2 (x,_)) = x + +getY :: Point2' a -> a +getY (Point2 (_,y)) = y + +-- | euclidean distance between p and q +dist :: Floating a => Point2' a -> Point2' a -> a +dist p q = let a = getX p - getX q + b = getY p - getY q in sqrt (a*a + b*b) diff --git a/src/Data/Geometry/Polygon.hs b/src/Data/Geometry/Polygon.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Polygon.hs @@ -0,0 +1,60 @@ +module Data.Geometry.Polygon where + +import Data.Geometry.BoundingBox +import Data.Geometry.Point +import Data.Geometry.Geometry + +--------------------------------------------------------------------- +-- | Polygons + +-- | 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 + + -- | 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]) + + isSimple :: p a -> Bool + containsHoles :: p a -> Bool + +--------------------------------------------------------------------- +-- | 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) + +instance HasPoints SimplePolygon' where + points (SimplePolygon pts) = pts + +instance IsPolygon SimplePolygon' where + isSimple = const True + containsHoles = const False + +instance IsPoint2Functor SimplePolygon' where + p2fmap f (SimplePolygon pts) = SimplePolygon (map f pts) + +--------------------------------------------------------------------- +-- | A multipolygon consists of several simple polygons +data MultiPolygon' a = MultiPolygon [SimplePolygon' a] + deriving (Show,Eq) + +instance HasPoints MultiPolygon' where + points (MultiPolygon pls) = concatMap points pls + +instance IsPolygon MultiPolygon' where + isSimple = const False + containsHoles = const False + +instance IsPoint2Functor MultiPolygon' where + p2fmap f (MultiPolygon polys) = MultiPolygon (map (p2fmap f) polys) + +--------------------------------------------------------------------- +-- | Bounding boxes can be used as polygons +instance IsPolygon BoundingBox2' where + isSimple = const True + containsHoles = const False diff --git a/src/Data/Geometry/SetOperations.hs b/src/Data/Geometry/SetOperations.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/SetOperations.hs @@ -0,0 +1,66 @@ +{-# 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