Add new files

src/Geometry/ConvexPoly.hs

@@ -0,0 +1,77 @@
+{-# LANGUAGE TemplateHaskell #-}
+{-# LANGUAGE BangPatterns #-}
+module Geometry.ConvexPoly
+    ( ConvexPoly (..)
+    , cpPoints
+    , cpCen
+    , cpRad
+    , centroid
+    , pointsToPoly
+    , convexPolysOverlap
+    ) where
+import Geometry.Data
+import Geometry.Vector
+import Geometry.LHS
+import Geometry.Intersect
+
+import Data.Maybe
+import Control.Lens
+import qualified Control.Foldl as L
+data ConvexPoly = ConvexPoly
+    { _cpPoints :: [Point2]
+    , _cpCen :: Point2
+    , _cpRad    :: Float
+    }
+pointsToPoly :: [Point2] -> ConvexPoly
+pointsToPoly xs = ConvexPoly
+    { _cpPoints = xs
+    , _cpCen = cen
+    , _cpRad    = minimum $ map (dist cen) xs
+    }
+  where
+    cen = centroid xs
+
+-- | Test whether two polygons intersect or if one is contained in the other.
+convexPolysOverlap :: ConvexPoly -> ConvexPoly -> Bool
+convexPolysOverlap cp1 cp2 = dist (_cpCen cp1) (_cpCen cp2) < _cpRad cp1 + _cpRad cp2
+    && polyPointsOverlap (_cpPoints cp1) (_cpPoints cp2)
+
+--pointInConvexPoly :: Point2 -> ConvexPoly -> Bool
+--pointInConvexPoly p cp = dist p (_cpCen cp) < _cpRad cp
+--    && pointInPolyPoints p (_cpPoints cp)
+
+-- | Test whether two polygons intersect or if one is contained in the other.
+polyPointsOverlap :: [Point2] -> [Point2] -> Bool
+polyPointsOverlap (p:ps) (q:qs) = pointInPolyPoints p (q:qs)
+    || pointInPolyPoints q (p:ps)
+    || polyPointsIntersect (p:ps) (q:qs) 
+polyPointsOverlap _ _ = False
+
+-- | Test whether a point is strictly inside a polygon.
+-- Supposes the points in the polygon are listed in anticlockwise order.
+pointInPolyPoints :: Point2 -> [Point2] -> Bool
+pointInPolyPoints !p (x:xs) = all (\l -> uncurry isLHS l p) $ zip (x:xs) (xs ++ [x])
+pointInPolyPoints _ [] = False
+
+polyPointsIntersect :: [Point2] -> [Point2] -> Bool
+polyPointsIntersect (a:b:xs) ps = go a (a:b:xs) ps
+  where
+    go x' (a':b':xs') ps' = pairPolyPointsIntersect a' b' ps' || go x' (b':xs') ps'
+    go b' (a':[]) ps' = pairPolyPointsIntersect a' b' ps'
+    go _ _ _ = False
+polyPointsIntersect _ _ = False
+
+pairPolyPointsIntersect :: Point2 -> Point2 -> [Point2] -> Bool
+pairPolyPointsIntersect a' b' (c':d':xs') = go c' a' b' (c':d':xs')
+  where
+    go x a b (c:d:xs) 
+        | isJust $ myIntersectSegSeg a b c d = True
+        | otherwise = go x a b (d:xs)
+    go d a b (c:[]) = isJust $ myIntersectSegSeg a b c d
+    go _ _ _ _ = False
+pairPolyPointsIntersect _ _ _ = False
+
+centroid :: Foldable t => t Point2 -> Point2
+centroid = L.fold $ (/) <$> L.Fold (+.+) (V2 0 0) id <*> L.genericLength
+
+makeLenses ''ConvexPoly

src/Geometry/Intersect.hs

@@ -87,6 +87,7 @@ myIntersectSegSeg
     -> Point2
     -> Point2
     -> Maybe Point2
+{-# INLINE myIntersectSegSeg #-}
 myIntersectSegSeg a@(V2 ax ay) b@(V2 bx by) c@(V2 cx cy) d@(V2 dx dy) = case ratIntersectLineLine a b c d of
     Nothing    -> Nothing
     Just (V2 x y) -> if inbetween x && inbetween' y 
@@ -99,6 +100,7 @@ myIntersectSegSeg a@(V2 ax ay) b@(V2 bx by) c@(V2 cx cy) d@(V2 dx dy) = case rat
         &&  ((cy <= y && y <= dy) || (dy <= y && y <= cy))
 -- | Polymorphic intersection of fractional line points.
 myIntersectLineLine :: (Eq a,Fractional a) => V2 a -> V2 a -> V2 a -> V2 a -> Maybe (V2 a)
+{-# INLINE myIntersectLineLine #-}
 myIntersectLineLine a@(V2 ax _) b c@(V2 cx _) d
     | isNothing (linGrad a b) = V2 ax <$> axisInt (c *-* V2 ax 0) (d *-* V2 ax 0)
     | isNothing (linGrad c d) = V2 cx <$> axisInt (a *-* V2 cx 0) (b *-* V2 cx 0)
@@ -115,6 +117,7 @@ myIntersectLineLine a@(V2 ax _) b c@(V2 cx _) d
     (*-*) (V2 ax' ay) (V2 bx by) = V2 (ax'-bx) (ay-by)
 -- | Transforms floating points to rationals then performs line intersection.
 ratIntersectLineLine :: Point2 -> Point2 -> Point2 -> Point2 -> Maybe Point2
+{-# INLINE ratIntersectLineLine #-}
 ratIntersectLineLine a b c d = toNumPoint2 
     <$> myIntersectLineLine (toRatPoint2 a) (toRatPoint2 b) (toRatPoint2 c) (toRatPoint2 d)
   where 
@@ -131,18 +134,19 @@ but is symmetric around 0:
 
 >>> roundPoint2 (0.5,-0.5)
 (0.0,0.0)
-
 -}
 roundPoint2 :: Point2 -> Point2
 roundPoint2 (V2 x y) = V2 (fromIntegral (round x :: Int)) (fromIntegral (round y :: Int))
 -- | Given two points, finds the linear gradient if it is non-infinite.
 linGrad :: (Eq a,Fractional a) => V2 a -> V2 a -> Maybe a
+{-# INLINE linGrad #-}
 linGrad (V2 x y) (V2 a b) 
     | x-a == 0  = Nothing
     | otherwise = Just $ (y-b)/(x-a)
 -- | Given two points, finds the intersection with the y axis if it exists.
 axisInt :: (Eq a,Fractional a) => V2 a -> V2 a -> Maybe a
-axisInt p (V2 a b) =  (\lg -> b - (a*lg)) <$> linGrad p (V2 a b)
+{-# INLINE axisInt #-}
+axisInt p (V2 a b) = (\lg -> b - (a*lg)) <$> linGrad p (V2 a b)
 -- | Placeholder, undefined.
 intersectSegsSeg :: [Point2] -> Point2 -> Point2 -> Maybe Point2
 intersectSegsSeg = undefined

src/Geometry/LHS.hs

@@ -0,0 +1,44 @@
+module Geometry.LHS
+    ( isLHS
+    , isRHS
+    ) where
+import Geometry.Data
+import Geometry.Vector
+-- | Test whether a point is on the LHS of a line.
+-- Returns False if the line is of zero length.
+isLHS 
+    :: Point2  -- ^ First line point. 
+    -> Point2  -- ^ Second line point.
+    -> Point2  -- ^ Point not on line.
+    -> Bool
+{-# INLINE isLHS #-}
+isLHS 
+    (V2 x y)      
+    (V2 x' y')   
+    (V2 x'' y'') 
+    | (x,y) == (x',y')  = False
+    | otherwise = a1 * b2 - a2 * b1 > 0
+  where 
+    a1 = x' - x
+    a2 = y' - y
+    b1 = x'' - x
+    b2 = y'' - y
+-- | Test whether a point is on the LHS of a line.
+-- Returns False if the line is of zero length.
+isRHS 
+    :: Point2 -- ^ First line point. 
+    -> Point2 -- ^ Second line point.
+    -> Point2 -- ^ Point not on line.
+    -> Bool
+{-# INLINE isRHS #-}
+isRHS 
+    (V2 x y)       
+    (V2 x' y')    
+    (V2 x'' y'')  
+    | (x,y) == (x',y')  = False
+    | otherwise = a1 * b2 - a2 * b1 < 0
+  where 
+    a1 = x' - x
+    a2 = y' - y
+    b1 = x'' - x
+    b2 = y'' - y

src/Geometry/Vector.hs

@@ -122,4 +122,7 @@ projV fromv onv
     | otherwise = (fromv `dotV` onv) / den *.* onv
   where
     den = onv `dotV` onv
-
+-- | Return distance between two points.
+dist :: Point2 -> Point2 -> Float
+{-# INLINE dist #-}
+dist !p1 !p2 = magV (p2 -.- p1)